Моделирование жидкости методом SPH

Материал из Department of Theoretical and Applied Mechanics
Перейти к: навигация, поиск
Курсовые работы по МДС: 2016-2017 > Моделирование жидкости методом SPH

Курсовой проект по Механике дискретных сред

Исполнитель: Шварёв Николай

Группа: 09 (43604/1)

Семестр: осень 2016

Описание метода

SPH (Smoothed Particle Hydrodynamics, Гидродинамика сглаженных частиц) - вычислительный лагранжев метод для симуляции жидкостей и газа, который используется в биологии, астрофизике, баллистике, вулканографии, океанографии.

Сутью метода SPH является разбиение жидкости на дискретные элементы-частицы. Любая физическая величина любой частицы может быть получена путём суммирования соответствующих величин всех частиц которые находятся в пределах двух сглаженных с помощью функции ядра длин.

Значение любой величины A на любом расстоянии r задается формулой:

[math]A(\pmb{r})= \sum_{j}m_{j}\frac{A_{j}}{ρ_{j}}W(|\pmb r -\pmb r_{j}|,h) \qquad \qquad \qquad \qquad \qquad j = 1...N,[/math]

где [math]m_{j}[/math] - масса частицы j, [math]A_{j}[/math] - значение величины A для частицы j, [math]ρ_{j}[/math] - плотность частицы j, W - функция ядра, h - радиус обрезания.

Таким образом, плотность любой частицы высчитывается по формуле:

[math]ρ(\pmb{r})= \sum_{j}m_{j}W(|\pmb r -\pmb r_{j}|,h) \qquad \qquad \qquad \qquad \qquad i,j = 1...N,[/math]

В нашем случае роль функции ядра исполняет функция Люси:

[math] W(r \lt h) = \frac{5}{π{b}^{2}}(1+3\frac{r}{h}){(1-\frac{r}{h})}^{3}, [/math]

где [math]r = |\pmb r - \pmb r_{j}|[/math]

Используя уравнение Эйлера, уравнения изменения плотности и положения

[math] \begin{cases} \frac{d\pmb v}{dt} = -\frac{1}{ρ}∇P + \pmb g \\ \frac{dρ}{dt} = -ρ∇\cdot\pmb v \\ \frac{d\pmb r}{dt} = \pmb v \\ \end{cases} [/math]

получаем уравнение движения каждой частицы:

[math]\frac{d\pmb v_{a}}{dt} = - \sum_{b}m_{b}(\frac{P_{b}}{{ρ_{b}}^{2}} + \frac{P_{a}}{{ρ_{a}}^{2}} + П_{ab})∇_{a}W_{ab} \qquad \qquad \qquad \qquad \qquad a,b = 1...N,[/math]

где давление P в нашей задаче вычисляется по формуле:

[math]P = B (({\frac{ρ}{ρ_{0}}})^{\gamma} -1)[/math],

B - модуль упругости, [math]ρ_{0}[/math] - равновесная плотность.

А [math]П_{ab}[/math] - член, отвечающий за вязкость, который в нашем случае рассчитывается так:

[math]П_{ab} = -\frac{\alpha h c}{ρ_{a b}}(\frac{\pmb v_{ab} \cdot \pmb r_{ab}}{{r_{ab}}^{2}+0.01{h}^{2}}) [/math]

Поставленная задача

Рассмотреть двумерную задачу о разрушении плотины: как поведет себя столб жидкости при разрушении (исчезновении) одной из стенок, удерживающей этот столб.

Реализация

Текст программы на языке JavaScript:

Файл "4.js"

  1 // m: Вязкоупругий шар
  2 // Используется исходник Воробьёва Сергея от исходиника программы Цветкова Balls Версия 6.3 от 05.05.2014
  3 
  4 
  5 function MainBalls(canvas, slider_01, text_01, slider_02, text_02) {
  6 
  7     canvas.onselectstart = function () {return false;};     // запрет выделения canvas
  8 
  9     // Предварительные установки
 10 
 11     var context = canvas.getContext("2d");                  // на context происходит рисование
 12     canvas.oncontextmenu = function (e) {return false;};    // блокировка контекстного меню
 13 
 14     var Pi = 3.1415926;                 // число "пи"
 15 	var count = 0;						// счетчик в первом присваивании ro0
 16 	
 17 	var γ = 7;
 18 	var ro0 = 1;
 19 	var alpha = 1;
 20 	var speedOfLight = 300000000;
 21 	//var Bb = ro0 * speedOfLight* speedOfLight / 7;                          // модуль упругости по определению, но частицы разлетаются
 22 	var Bb = 10;											//был 1, потом 5
 23 	var Cab = 0.5;										//более красивое при ~2, но изначально было 5
 24 	
 25 	
 26     var m0 = 1;                         // масштаб массы
 27     var t0 = 1;                         // масштаб времени (период колебаний исходной системы)
 28     var a0 = 1;                         // масштаб расстояния (диаметр шара)
 29 
 30     var g0 = a0 / t0 / t0;              // масштаб ускорения (ускорение, при котором за t0 будет пройдено расстояние a0)
 31     var k0 = 2 * Pi / t0;               // масштаб частоты
 32     var C0 = m0 * k0 * k0;              // масштаб жесткости
 33     var B0 = 2 * m0 * k0;               // масштаб вязкости
 34 
 35     // *** Задание физических параметров ***
 36 
 37     var Ny = 50	;   //увеличить до 50-100, было 16, потом 22                     // число шаров, помещающихся по вертикали в окно (задает размер шара относительно размера окна)
 38     var m = 1* m0;                     // масса
 39     var CWall = 10 * C0;                // жесткость стен
 40     var CBall = 0.1 * CWall;            // жесткость между частицами
 41     var BVisc = 0.008 * B0;             // вязкость среды        было 0.008
 42     var BWall = 0.03 * B0;              // вязкость на стенках
 43     
 44     var r = 0.5 * a0;                   // радиус частицы в расчетных координатах
 45     var K = 0.7; //   ???                       // все силы, зависящие от радиуса, ограничиваются значением, реализующимся при r/a = K
 46     var a = 2 * r;                      // равновесное расстояние между частицами
 47     var aCut = 2.5 * r;                   // радиус обрезания
 48 	var alfa = 2;                     // коэффициент для хрупкого вз. Лен-Дж
 49 
 50     // *** Задание вычислительных параметров ***
 51 
 52     var fps = 20;                       // frames per second - число кадров в секунду (качечтво отображения)
 53     var spf = 100;                      // steps per frame   - число шагов интегрирования между кадрами (скорость расчета)
 54     var dt  = 0.04 * t0 / fps;          // шаг интегрирования (качество расчета)    //было 0.04*to/fps
 55 	var mg = 1 * m * g0;             // сила тяжести в начальный момент времени
 56 
 57     // Выполнение программы
 58 	var sqrt3 = Math.sqrt(3);
 59     var r2 = r * r;                     // ___в целях оптимизации___
 60     var a2 = a * a;                     // ___в целях оптимизации___
 61     var D = a2 * CBall / 72;            // энергия связи между частицами
 62     var LJCoeff = 12 * D / a2;          // коэффициент для расчета потенциала Л-Дж
 63     var bet = Math.pow(13 / 7, 1/6) * a;    // коэффициент для SLJ потенциала
 64     var bet2 = bet * bet;                     // ___в целях оптимизации___
 65     var SLJDenominator = 1 / (aCut * aCut - bet2);    // знаменатель для расчета SLJ потенциала
 66 	var sqrtkoef = Math.sqrt(alfa/(1+alfa));  //___в целях оптимизации___
 67 
 68     var Ka = K * a;                     // ___в целях оптимизации___
 69     var K2a2 = K * K * a2;              // ___в целях оптимизации___
 70 
 71     var dNd = null;                     // ссылка на захваченный курсором шар (drag & drop)
 72 	
 73 	var numberOfBorder = 0;
 74 
 75     this.setSlider_01 = function(c) {mg = c * m * g0;}; // функция для слайдера гравитации; 
 76 	this.setSlider_02 = function(c) {10;}; // функция для слайдера количества шаров по Х; 
 77 	
 78     // Настройка интерфейса
 79 
 80     slider_01.min = 0;               slider_01.max = 5;
 81     slider_01.step = 0.05;
 82     slider_01.value = mg / m / g0;          // начальное значение ползунка должно задаваться после min и max
 83     text_01.value = mg / m / g0;
 84 
 85 	
 86 //	slider_02.min = 1;
 87 //	slider_02.step = 1;
 88 //	slider_02.max = w-3;
 89 	
 90 
 91 	
 92     // Запуск новой системы
 93 
 94     // следующие переменные должны пересчитываться каждый раз, когда мы изменяем значение Ny
 95     var scale, w, h;
 96     var rScale13, rScaleShift;
 97     this.newSystem = function() {
 98         scale = canvas.height / Ny / a0;    // масштабный коэффициент для перехода от расчетных к экранным координатам
 99         w = canvas.width / scale;           // ширина окна в расчетных координатах
100         h = canvas.height / scale;          // высота окна в расчетных координатах
101 
102         rScale13 = r * scale * 1.3;         // ___в целях оптимизации___
103         rScaleShift = r * scale / 5;        // ___в целях оптимизации___
104 		
105 		this.setBorder();
106         this.setQuad(10,10);   //размером 10х10                 // сразу создаем конфигурацию
107 		this.setFence();
108     };
109 
110     // Работа с мышью
111 
112     var mx_, my_;                               // буфер позиции мыши (для расчета скорости при отпускании шара)
113 
114     canvas.onmousedown = function(e) {          // функция при нажатии клавиши мыши
115         var m = mouseCoords(e);                 // получаем расчетные координаты курсора мыши
116         // цикл в обратную сторону, чтобы захватывать шар, нарисованный "сверху"
117         // (т.к. цикл рисования идет в обычном порядке)
118         for (var i = balls.length - 1; i >= 0; i--) {
119             var b = balls[i];
120             var rx = b.x - m.x;
121             var ry = b.y - m.y;
122             var rLen2 = rx * rx + ry * ry;              // квадрат расстояния между курсором и центром шара
123             if (rLen2 <= r2) {                          // курсор нажал на шар
124                 if (e.which == 1) {                     // нажата левая клавиша мыши
125                     dNd = b;
126                     dNd.xPlus = dNd.x - m.x;            // сдвиг курсора относительно центра шара по x
127                     dNd.yPlus = dNd.y - m.y;            // сдвиг курсора относительно центра шара по y
128                     mx_ = m.x;    my_ = m.y;
129                     canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
130                 } else if (e.which == 3)                // нажата правая клавиша мыши
131                     balls.splice(i, 1);                 // удалить шар
132                 return;
133             }
134         }
135 
136         // если не вышли по return из цикла - нажатие было вне шара, добавляем новый
137         if (e.which == 1) {
138             dNd = addNewBall(m.x, m.y, true);   // добавляем шар и сразу захватываем его курсором
139             if (dNd == null) return;            // если шар не добавился (из за стен или других шаров) - возвращаемся
140             dNd.xPlus = 0;  dNd.yPlus = 0;      // держим шар по центру
141             mx_ = m.x;    my_ = m.y;
142             canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
143         }
144     };
145 
146     document.onmouseup = function(e) {          // функция при отпускании клавиши мыши
147         canvas.onmousemove = null;              // когда клавиша отпущена - функции перемещения нету
148         dNd = null;                             // когда клавиша отпущена - захваченного курсором шара нету
149     };
150 
151     function mouseMove(e) {                     // функция при перемещении мыши, работает только с зажатой ЛКМ
152         var m = mouseCoords(e);                 // получаем расчетные координаты курсора мыши
153         dNd.x = m.x + dNd.xPlus;
154         dNd.y = m.y + dNd.yPlus;
155         dNd.vx = 0.6 * (m.x - mx_) / dt / fps;   dNd.vy = 0.6 * (m.y - my_) / dt / fps;
156         mx_ = m.x;    my_ = m.y;
157     }
158 
159     function mouseCoords(e) {                   // функция возвращает расчетные координаты курсора мыши
160         var m = [];
161         var rect = canvas.getBoundingClientRect();
162         m.x = (e.clientX - rect.left) / scale;
163         m.y = (e.clientY - rect.top) / scale;
164         return m;
165     }
166 
167     // Работа с массивом
168 
169     var balls = [];                             // массив шаров
170     var addNewBall =  function(x, y, check) {
171         // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
172         if (check) {
173             if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
174             for (var i = 0; i < balls.length; i++) {
175                 var rx = balls[i].x - x;
176                 var ry = balls[i].y - y;
177                 var rLen2 = rx * rx + ry * ry;
178                 if (rLen2 < 4 * r2) return null;
179             }
180         }
181 
182         var b = [];
183 
184         b.x = x;                b.y = y;        // расчетные координаты шара
185         b.fx = 0;               b.fy = mg;      // сила, действующая на шар                  								    , вернул знак
186         b.vx = 0;               b.vy = 0;       // скорость 
187 		b.ro = 0;                         
188 		b.ro0 = ro0;
189 		b.ax = 0;
190 		b.ay = 0;
191 		b.m = m;
192 		
193 		b.latX = Math.floor(b.x / aCut);
194 		b.latY = Math.floor(b.y / aCut);
195 		
196 		balls[balls.length] = b;                // добавить элемент в конец массива
197         return b;
198     };
199   
200     this.setTriangularLattice = function(latX,latY) {            // задать на поле треугольную решетку (Серегин код, апгрейднутый мной)
201 //		countFence = 0;
202 		
203 	//        balls = [];
204        	for (var j = Math.floor(latY); j >= 0 ; j--)    
205           for (var i = 1; i < Math.floor(latX / r)-1  ; i++)
206               if ((i + j) % 2 == 0) addNewBall(4 * r + i*r, h - 4*r - r*sqrt3*j - r, false);
207 		  
208 	};
209     
210     this.setQuad = function(latX,latY) {               // квадратная конфигурация (мой код)  размера latX x latY
211 	    for (var j = Math.floor(latY); j > 0; j--)
212 		   for (var i = 0; i < Math.floor(latX); i++)
213 				addNewBall(r*2*(i+2) + r, h - 2*r*(j+1) -r ,false);
214 	
215     };
216 	
217 	this.setBorder = function(){                                                 //установка границ на пределе прорисовки
218 		balls = [];
219 		numberOfBorder = 0;
220 		
221 		this.deleteFence();
222 		for (var i = 0; i < 2; i++)
223 		{
224 			for (var j = 0; j < w ; j++){
225 				addNewBall(j + r, i*2*r + r,false);
226 				numberOfBorder++;			
227 			}
228 			for (var j = 0; j < w + 1; j++){
229 				addNewBall(j - r,h - 2 * r * i - r,false);
230 				numberOfBorder++;
231 			}
232 			for (var j = 2; j < h - 2; j++){
233 				addNewBall(2 * r * i + r, j + r, false);
234 				numberOfBorder++;
235 			}
236 			for (var j = 2; j < h - 2; j++){
237 				addNewBall(w - 2 * r * i - r / 2, j + r, false);
238 				numberOfBorder++;
239 			}
240 		}
241 	}	
242  	var countFence = 0;
243 	var countVarFence = true;
244 	this.setFence = function(){
245 		if (countVarFence){
246 			countFence = 0;
247 			for(var i = 0; i < 2; i++)
248 				for(var j = Math.round(3*h/4); j > 0 ; j--){
249 					addNewBall(26 * r + 2 * r * i , h - 3 * r - 2*r*j,false);
250 					countFence++;
251 			}
252 			countVarFence = false;
253 		}
254 
255 	}
256 	
257 	this.deleteFence = function(){
258 		b1 = balls;
259 		balls = [];
260 		for(var i = 0; i < b1.length - countFence; i++)
261 			balls[i] = b1[i];
262 		
263 		countFence = 0;
264 		countVarFence = true;
265 
266 	}
267 
268 
269   
270  // this.setEmpty = function() {balls = [];};                    //за ненадобностью, заменен на задание границ
271 
272     // Основной цикл программы
273 
274     function control() {
275         physics();
276         draw();
277     }
278 
279     // Расчетная часть программы
280 
281     function physics() {                            // то, что происходит каждый шаг времени
282         for (var s = 1; s <= spf; s++) {
283 
284 			for (var i = 0; i < balls.length; i++) {					//обнуление 
285                 balls[i].fy = 0;
286 				balls[i].fx = 0;
287 				balls[i].ro = 0;                         
288 				balls[i].ax = 0;
289 				balls[i].ay = 0;
290             }
291 
292 			impAverageX = 0;
293 			
294 
295 			for (var i = 0; i < balls.length; i++) {						//расчет плотностей
296 				var b1 = balls[i];
297                 for (var j = 0; j < balls.length; j++) {  
298 					var b2 = balls[j];	
299 					if((b2.latX == (b1.latX - 1)) | (b2.latX == (b1.latX)) | (b2.latX == (b1.latX + 1)) | (b2.latY == (b1.latY - 1)) | (b2.latY == (b1.latY)) | (b2.latY == (b1.latY + 1))){
300 	
301 						var rx = b1.x - b2.x;   
302 						var ry = b1.y - b2.y;         					// вектор смотрит на первый шар (на b1 из b2)
303 						var r2 = rx * rx + ry * ry;                         // квадрат расстояния между шарами
304 						var rLen = Math.sqrt(r2);                           //расстояние между шарами
305 						if (rLen <= aCut) {
306 		
307 							b1.ro += b2.m * (5 / (Pi * aCut * aCut)) * (1 + 3 * rLen / aCut) * (1 - rLen / aCut)* (1 - rLen / aCut)* (1 - rLen / aCut);  // ядро Люси
308 							//console.log(b1.ro);
309 						}
310 					}
311 				}
312 				balls[i] = b1;
313 
314 			}
315 			
316 /*			if (count != balls.length){														//изменение плотности при появлении нового шарика
317 
318 				var ro0 = balls[0].ro
319 				for(var i = 1; i < balls.length; i++){
320 					if (balls[i].ro < ro0){
321 						ro0 = balls[i].ro;
322 					}
323 							
324 				}
325 
326 				for (var i = 0; i < balls.length; i++){
327 					balls[i].ro0 = ro0;
328 				}
329 				count = balls.length;
330 
331 			}
332 */		
333 	
334 			for (var i = 0; i < balls.length; i++) {									    //расчет давлений для каждой частицы
335 				balls[i].p = Bb * (Math.pow(balls[i].ro/balls[i].ro0,γ) - 1); 
336 	//			balls[i].p = Bb*(balls[i].ro-ro0);                                      //с сайта кафедры
337  			}
338 
339 			
340 			for (var i = numberOfBorder; i < balls.length - countFence; i++) {									//расчет ускорений
341 				var b1 = balls[i];
342                 for (var j = 0; j < balls.length; j++) {                   
343                     var b2 = balls[j];	
344 					if((b2.latX == (b1.latX - 1)) | (b2.latX == (b1.latX)) | (b2.latX == (b1.latX + 1)) | (b2.latY == (b1.latY - 1)) | (b2.latY == (b1.latY)) | (b2.latY == (b1.latY + 1))){
345 
346 						var rx = b1.x - b2.x;   
347 						var ry = b1.y - b2.y;         					// вектор смотрит на первый шар (b1 из "нуля")
348 						var r2 = rx * rx + ry * ry;                         // квадрат расстояния между шарами
349 						var rLen = Math.sqrt(r2);                           //расстояние между шарами
350 						
351 						var vx = b1.vx - b2.vx;
352 						var vy = b1.vy - b2.vy;
353 
354 						if (rLen <= aCut) {
355 							var PabX = -alpha * aCut * Cab /((b1.ro + b2.ro)/2) * vx * rx/(rLen * rLen + 0.01 * aCut * aCut);
356 							var PabY = -alpha * aCut * Cab /((b1.ro + b2.ro)/2) * vy * ry/(rLen * rLen + 0.01 * aCut * aCut);
357 							
358 							b1.ax += b2.m * 60 / (Math.pow(aCut,6) * Pi) * (aCut - rLen) * (aCut - rLen) * (b2.p * rx / (b2.ro * b2.ro) + b1.p * rx / (b1.ro * b1.ro) + PabX * rx);
359 							b1.ay += b2.m * 60 / (Math.pow(aCut,6) * Pi) * (aCut - rLen) * (aCut - rLen) * (b2.p * ry / (b2.ro * b2.ro) + b1.p * ry / (b1.ro * b1.ro) + PabY * ry);
360 						}
361 					}
362 				}
363 			
364                 if (b1 == dNd) continue;  // если шар схвачен курсором - его вз. со стенами и перемещение не считаем
365 
366 //                if (b1.y + r > h) { b1.fy += -CWall * (b1.y + r - h) - BWall * b1.vy; }
367 //                if (b1.y - r < 0) { b1.fy += -CWall * (b1.y - r) - BWall * b1.vy;}
368 //                if (b1.x + r > w) { b1.fx += -CWall * (b1.x + r - w) - BWall * b1.vx; }
369 //                if (b1.x - r < 0) { b1.fx += -CWall * (b1.x - r) - BWall * b1.vx; }
370 				balls[i] = b1;
371 			}
372 		    for (var i = numberOfBorder; i < balls.length - countFence; i++){
373 				
374 				b1 = balls[i];
375 				b1.vx += (b1.fx + b1.ax) * dt;        
376 				b1.vy += (b1.fy + b1.ay + mg / m) * dt;
377 				b1.x += b1.vx * dt;             
378 				b1.y += b1.vy * dt;
379 				b1.latX = Math.floor(b1.x / aCut);
380 				b1.latY = Math.floor(b1.y / aCut);
381 				balls[i] = b1;
382 //				impAverageX += m * b1.vx;
383 				
384 			}
385 			//console.log(impAverageX/balls.length)
386 			//console.log(numberOfBorder);
387 				
388     }
389 }
390 
391     // Рисование
392   context.fillStyle = "#d3692e";
393     function draw() {
394         context.clearRect(0, 0, w * scale, h * scale);      // очистить экран
395         for (var i = 0; i < balls.length; i++){
396             var xS = balls[i].x * scale;           var yS = balls[i].y * scale;
397             context.beginPath();
398             context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
399             context.closePath();
400             context.fill();
401         }
402     }
403 
404     // Запуск системы
405     this.newSystem();
406     setInterval(control, 1000 / fps);
407     // след. функция обновляет информацию о количестве частиц на поле
408     setInterval(function(){document.getElementById('ballsNum').innerHTML = balls.length;}, 1000 / 20);
409 }

Возможности программы

В текущей версии программы можно рассматривать разрушение столба жидкости в изначально разных конфигурациях:

- квадратная решетка размером 10x38;

- треугольная решетка размером 11x38;

- пустое поле, в котором можно создать собственную конфигурацию (ЛКМ - создание частицы).

В любой момент можно поставить или убрать забор. Также есть возможность изменения гравитации от 0g до 5g.

Результаты

Рисунок 1. Начальная конфигурация
Рисунок 2. Успокоившийся столб жидкости
Рисунок 3. Отскок от границы
Рисунок 4. Осевшее на дно












































Ссылки по теме

  • J.J.Monaghan "Smoothed particle hydrodynamics", 1992;
  • W.G.Hoover, C.G.Hoover "SPAM-Based Recipes for Continuum Simulation".