Моделирование жидкости методом SPH
Материал из Department of Theoretical and Applied Mechanics
Версия от 11:25, 23 декабря 2016; Shvarevng (обсуждение | вклад)
Курсовые работы по МДС: 2016-2017 > Моделирование жидкости методом SPH
Курсовой проект по Механике дискретных сред
Исполнитель: Шварёв Николай
Группа: 09 (43604/1)
Семестр: осень 2016
Содержание
Описание метода
SPH (Smoothed Particle Hydrodynamics, Гидродинамика сглаженных частиц) - вычислительный лагранжев метод для симуляции жидкостей и газа, который используется в биологии, астрофизике, баллистике, вулканографии, океанографии.
Сутью метода SPH является разбиение жидкости на дискретные элементы-частицы. Любая физическая величина любой частицы может быть получена путём суммирования соответствующих величин всех частиц которые находятся в пределах двух сглаженных с помощью функции ядра длин.
Значение любой величины A на любом расстоянии r задается формулой:
где - масса частицы j, - значение величины A для частицы j, - плотность частицы j, W - функция ядра, h - радиус обрезания.
Таким образом, плотность любой частицы высчитывается по формуле:
В нашем случае роль функции ядра исполняет функция Люси:
где
Используя уравнение Эйлера, уравнения изменения плотности и положения
получаем уравнение движения каждой частицы:
где давление P в нашей задаче вычисляется по формуле:
,
B - модуль упругости, - равновесная плотность.
А - член, отвечающий за вязкость, который в нашем случае рассчитывается так:
Поставленная задача
Рассмотреть двумерную задачу о разрушении плотины: как поведет себя столб жидкости при разрушении (исчезновении) одной из стенок, удерживающей этот столб.
Реализация
Текст программы на языке 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.
Результаты
Ссылки по теме
- J.J.Monaghan "Smoothed particle hydrodynamics", 1992;
- W.G.Hoover, C.G.Hoover "SPAM-Based Recipes for Continuum Simulation".
