Моделирование жидкости методом SPH — различия между версиями

Материал из Department of Theoretical and Applied Mechanics
Перейти к: навигация, поиск
(→‎Реализация)
 
(не показано 16 промежуточных версий 2 участников)
Строка 1: Строка 1:
−
[[Курсовые работы по МДС: 2016-2017]] > '''Моделирование жидкости методом SPH''' <HR>
+
[[Курсовые работы по ВМДС: 2016-2017]] > '''Моделирование жидкости методом SPH''' <HR>
  
 
'''''Курсовой проект по [[Механика дискретных сред|Механике дискретных сред]]'''''
 
'''''Курсовой проект по [[Механика дискретных сред|Механике дискретных сред]]'''''
Строка 61: Строка 61:
 
== Реализация ==
 
== Реализация ==
  
−
{{#widget:Iframe |url=http://tm.spbstu.ru/htmlets/Shvarev/SPH/SPH.html |width=800 |height=1200 |border=0 }}
+
{{#widget:Iframe |url=http://tm.spbstu.ru/htmlets/Shvarev/SPH/SPH.html |width=800 |height=800 |border=0 }}
 +
 
 +
<div class="mw-collapsible mw-collapsed">
 +
'''Текст программы на языке JavaScript:''' <div class="mw-collapsible-content">
 +
Файл '''"4.js"'''
 +
<syntaxhighlight lang="javascript" line start="1" enclose="div">
 +
// m: Вязкоупругий шар
 +
// Используется исходник Воробьёва Сергея от исходиника программы Цветкова Balls Версия 6.3 от 05.05.2014
 +
 
 +
 
 +
function MainBalls(canvas, slider_01, text_01, slider_02, text_02) {
 +
 
 +
    canvas.onselectstart = function () {return false;};    // запрет выделения canvas
 +
 
 +
    // Предварительные установки
 +
 
 +
    var context = canvas.getContext("2d");                  // на context происходит рисование
 +
    canvas.oncontextmenu = function (e) {return false;};    // блокировка контекстного меню
 +
 
 +
    var Pi = 3.1415926;                // число "пи"
 +
var count = 0; // счетчик в первом присваивании ro0
 +
 +
var γ = 7;
 +
var ro0 = 1;
 +
var alpha = 1;
 +
var speedOfLight = 300000000;
 +
//var Bb = ro0 * speedOfLight* speedOfLight / 7;                          // модуль упругости по определению, но частицы разлетаются
 +
var Bb = 10; //был 1, потом 5
 +
var Cab = 0.5; //более красивое при ~2, но изначально было 5
 +
 +
 +
    var m0 = 1;                        // масштаб массы
 +
    var t0 = 1;                        // масштаб времени (период колебаний исходной системы)
 +
    var a0 = 1;                        // масштаб расстояния (диаметр шара)
 +
 
 +
    var g0 = a0 / t0 / t0;              // масштаб ускорения (ускорение, при котором за t0 будет пройдено расстояние a0)
 +
    var k0 = 2 * Pi / t0;              // масштаб частоты
 +
    var C0 = m0 * k0 * k0;              // масштаб жесткости
 +
    var B0 = 2 * m0 * k0;              // масштаб вязкости
 +
 
 +
    // *** Задание физических параметров ***
 +
 
 +
    var Ny = 50 ;  //увеличить до 50-100, было 16, потом 22                    // число шаров, помещающихся по вертикали в окно (задает размер шара относительно размера окна)
 +
    var m = 1* m0;                    // масса
 +
    var CWall = 10 * C0;                // жесткость стен
 +
    var CBall = 0.1 * CWall;            // жесткость между частицами
 +
    var BVisc = 0.008 * B0;            // вязкость среды        было 0.008
 +
    var BWall = 0.03 * B0;              // вязкость на стенках
 +
   
 +
    var r = 0.5 * a0;                  // радиус частицы в расчетных координатах
 +
    var K = 0.7; //  ???                      // все силы, зависящие от радиуса, ограничиваются значением, реализующимся при r/a = K
 +
    var a = 2 * r;                      // равновесное расстояние между частицами
 +
    var aCut = 2.5 * r;                  // радиус обрезания
 +
var alfa = 2;                    // коэффициент для хрупкого вз. Лен-Дж
 +
 
 +
    // *** Задание вычислительных параметров ***
 +
 
 +
    var fps = 20;                      // frames per second - число кадров в секунду (качечтво отображения)
 +
    var spf = 100;                      // steps per frame  - число шагов интегрирования между кадрами (скорость расчета)
 +
    var dt  = 0.04 * t0 / fps;          // шаг интегрирования (качество расчета)    //было 0.04*to/fps
 +
var mg = 1 * m * g0;            // сила тяжести в начальный момент времени
 +
 
 +
    // Выполнение программы
 +
var sqrt3 = Math.sqrt(3);
 +
    var r2 = r * r;                    // ___в целях оптимизации___
 +
    var a2 = a * a;                    // ___в целях оптимизации___
 +
    var D = a2 * CBall / 72;            // энергия связи между частицами
 +
    var LJCoeff = 12 * D / a2;          // коэффициент для расчета потенциала Л-Дж
 +
    var bet = Math.pow(13 / 7, 1/6) * a;    // коэффициент для SLJ потенциала
 +
    var bet2 = bet * bet;                    // ___в целях оптимизации___
 +
    var SLJDenominator = 1 / (aCut * aCut - bet2);    // знаменатель для расчета SLJ потенциала
 +
var sqrtkoef = Math.sqrt(alfa/(1+alfa));  //___в целях оптимизации___
 +
 
 +
    var Ka = K * a;                    // ___в целях оптимизации___
 +
    var K2a2 = K * K * a2;              // ___в целях оптимизации___
 +
 
 +
    var dNd = null;                    // ссылка на захваченный курсором шар (drag & drop)
 +
 +
var numberOfBorder = 0;
 +
 
 +
    this.setSlider_01 = function(c) {mg = c * m * g0;}; // функция для слайдера гравитации;
 +
this.setSlider_02 = function(c) {10;}; // функция для слайдера количества шаров по Х;
 +
 +
    // Настройка интерфейса
 +
 
 +
    slider_01.min = 0;              slider_01.max = 5;
 +
    slider_01.step = 0.05;
 +
    slider_01.value = mg / m / g0;          // начальное значение ползунка должно задаваться после min и max
 +
    text_01.value = mg / m / g0;
 +
 
 +
 +
// slider_02.min = 1;
 +
// slider_02.step = 1;
 +
// slider_02.max = w-3;
 +
 +
 
 +
 +
    // Запуск новой системы
 +
 
 +
    // следующие переменные должны пересчитываться каждый раз, когда мы изменяем значение Ny
 +
    var scale, w, h;
 +
    var rScale13, rScaleShift;
 +
    this.newSystem = function() {
 +
        scale = canvas.height / Ny / a0;    // масштабный коэффициент для перехода от расчетных к экранным координатам
 +
        w = canvas.width / scale;          // ширина окна в расчетных координатах
 +
        h = canvas.height / scale;          // высота окна в расчетных координатах
 +
 
 +
        rScale13 = r * scale * 1.3;        // ___в целях оптимизации___
 +
        rScaleShift = r * scale / 5;        // ___в целях оптимизации___
 +
 +
this.setBorder();
 +
        this.setQuad(10,10);  //размером 10х10                // сразу создаем конфигурацию
 +
this.setFence();
 +
    };
 +
 
 +
    // Работа с мышью
 +
 
 +
    var mx_, my_;                              // буфер позиции мыши (для расчета скорости при отпускании шара)
 +
 
 +
    canvas.onmousedown = function(e) {          // функция при нажатии клавиши мыши
 +
        var m = mouseCoords(e);                // получаем расчетные координаты курсора мыши
 +
        // цикл в обратную сторону, чтобы захватывать шар, нарисованный "сверху"
 +
        // (т.к. цикл рисования идет в обычном порядке)
 +
        for (var i = balls.length - 1; i >= 0; i--) {
 +
            var b = balls[i];
 +
            var rx = b.x - m.x;
 +
            var ry = b.y - m.y;
 +
            var rLen2 = rx * rx + ry * ry;              // квадрат расстояния между курсором и центром шара
 +
            if (rLen2 <= r2) {                          // курсор нажал на шар
 +
                if (e.which == 1) {                    // нажата левая клавиша мыши
 +
                    dNd = b;
 +
                    dNd.xPlus = dNd.x - m.x;            // сдвиг курсора относительно центра шара по x
 +
                    dNd.yPlus = dNd.y - m.y;            // сдвиг курсора относительно центра шара по y
 +
                    mx_ = m.x;    my_ = m.y;
 +
                    canvas.onmousemove = mouseMove;    // пока клавиша нажата - работает функция перемещения
 +
                } else if (e.which == 3)                // нажата правая клавиша мыши
 +
                    balls.splice(i, 1);                // удалить шар
 +
                return;
 +
            }
 +
        }
 +
 
 +
        // если не вышли по return из цикла - нажатие было вне шара, добавляем новый
 +
        if (e.which == 1) {
 +
            dNd = addNewBall(m.x, m.y, true);  // добавляем шар и сразу захватываем его курсором
 +
            if (dNd == null) return;            // если шар не добавился (из за стен или других шаров) - возвращаемся
 +
            dNd.xPlus = 0;  dNd.yPlus = 0;      // держим шар по центру
 +
            mx_ = m.x;    my_ = m.y;
 +
            canvas.onmousemove = mouseMove;    // пока клавиша нажата - работает функция перемещения
 +
        }
 +
    };
 +
 
 +
    document.onmouseup = function(e) {          // функция при отпускании клавиши мыши
 +
        canvas.onmousemove = null;              // когда клавиша отпущена - функции перемещения нету
 +
        dNd = null;                            // когда клавиша отпущена - захваченного курсором шара нету
 +
    };
 +
 
 +
    function mouseMove(e) {                    // функция при перемещении мыши, работает только с зажатой ЛКМ
 +
        var m = mouseCoords(e);                // получаем расчетные координаты курсора мыши
 +
        dNd.x = m.x + dNd.xPlus;
 +
        dNd.y = m.y + dNd.yPlus;
 +
        dNd.vx = 0.6 * (m.x - mx_) / dt / fps;  dNd.vy = 0.6 * (m.y - my_) / dt / fps;
 +
        mx_ = m.x;    my_ = m.y;
 +
    }
 +
 
 +
    function mouseCoords(e) {                  // функция возвращает расчетные координаты курсора мыши
 +
        var m = [];
 +
        var rect = canvas.getBoundingClientRect();
 +
        m.x = (e.clientX - rect.left) / scale;
 +
        m.y = (e.clientY - rect.top) / scale;
 +
        return m;
 +
    }
 +
 
 +
    // Работа с массивом
 +
 
 +
    var balls = [];                            // массив шаров
 +
    var addNewBall =  function(x, y, check) {
 +
        // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
 +
        if (check) {
 +
            if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
 +
            for (var i = 0; i < balls.length; i++) {
 +
                var rx = balls[i].x - x;
 +
                var ry = balls[i].y - y;
 +
                var rLen2 = rx * rx + ry * ry;
 +
                if (rLen2 < 4 * r2) return null;
 +
            }
 +
        }
 +
 
 +
        var b = [];
 +
 
 +
        b.x = x;                b.y = y;        // расчетные координаты шара
 +
        b.fx = 0;              b.fy = mg;      // сила, действующая на шар                      , вернул знак
 +
        b.vx = 0;              b.vy = 0;      // скорость
 +
b.ro = 0;                       
 +
b.ro0 = ro0;
 +
b.ax = 0;
 +
b.ay = 0;
 +
b.m = m;
 +
 +
b.latX = Math.floor(b.x / aCut);
 +
b.latY = Math.floor(b.y / aCut);
 +
 +
balls[balls.length] = b;                // добавить элемент в конец массива
 +
        return b;
 +
    };
 +
 
 +
    this.setTriangularLattice = function(latX,latY) {            // задать на поле треугольную решетку (Серегин код, апгрейднутый мной)
 +
// countFence = 0;
 +
 +
//        balls = [];
 +
      for (var j = Math.floor(latY); j >= 0 ; j--)   
 +
          for (var i = 1; i < Math.floor(latX / r)-1  ; i++)
 +
              if ((i + j) % 2 == 0) addNewBall(4 * r + i*r, h - 4*r - r*sqrt3*j - r, false);
 +
 
 +
};
 +
   
 +
    this.setQuad = function(latX,latY) {              // квадратная конфигурация (мой код)  размера latX x latY
 +
    for (var j = Math.floor(latY); j > 0; j--)
 +
  for (var i = 0; i < Math.floor(latX); i++)
 +
addNewBall(r*2*(i+2) + r, h - 2*r*(j+1) -r ,false);
 +
 +
    };
 +
 +
this.setBorder = function(){                                                //установка границ на пределе прорисовки
 +
balls = [];
 +
numberOfBorder = 0;
 +
 +
this.deleteFence();
 +
for (var i = 0; i < 2; i++)
 +
{
 +
for (var j = 0; j < w ; j++){
 +
addNewBall(j + r, i*2*r + r,false);
 +
numberOfBorder++;
 +
}
 +
for (var j = 0; j < w + 1; j++){
 +
addNewBall(j - r,h - 2 * r * i - r,false);
 +
numberOfBorder++;
 +
}
 +
for (var j = 2; j < h - 2; j++){
 +
addNewBall(2 * r * i + r, j + r, false);
 +
numberOfBorder++;
 +
}
 +
for (var j = 2; j < h - 2; j++){
 +
addNewBall(w - 2 * r * i - r / 2, j + r, false);
 +
numberOfBorder++;
 +
}
 +
}
 +
}
 +
var countFence = 0;
 +
var countVarFence = true;
 +
this.setFence = function(){
 +
if (countVarFence){
 +
countFence = 0;
 +
for(var i = 0; i < 2; i++)
 +
for(var j = Math.round(3*h/4); j > 0 ; j--){
 +
addNewBall(26 * r + 2 * r * i , h - 3 * r - 2*r*j,false);
 +
countFence++;
 +
}
 +
countVarFence = false;
 +
}
 +
 
 +
}
 +
 +
this.deleteFence = function(){
 +
b1 = balls;
 +
balls = [];
 +
for(var i = 0; i < b1.length - countFence; i++)
 +
balls[i] = b1[i];
 +
 +
countFence = 0;
 +
countVarFence = true;
 +
 
 +
}
 +
 
 +
 
 +
 
 +
// this.setEmpty = function() {balls = [];};                    //за ненадобностью, заменен на задание границ
 +
 
 +
    // Основной цикл программы
 +
 
 +
    function control() {
 +
        physics();
 +
        draw();
 +
    }
 +
 
 +
    // Расчетная часть программы
 +
 
 +
    function physics() {                            // то, что происходит каждый шаг времени
 +
        for (var s = 1; s <= spf; s++) {
 +
 
 +
for (var i = 0; i < balls.length; i++) { //обнуление
 +
                balls[i].fy = 0;
 +
balls[i].fx = 0;
 +
balls[i].ro = 0;                       
 +
balls[i].ax = 0;
 +
balls[i].ay = 0;
 +
            }
 +
 
 +
impAverageX = 0;
 +
 +
 
 +
for (var i = 0; i < balls.length; i++) { //расчет плотностей
 +
var b1 = balls[i];
 +
                for (var j = 0; j < balls.length; j++) { 
 +
var b2 = balls[j];
 +
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))){
 +
 +
var rx = b1.x - b2.x; 
 +
var ry = b1.y - b2.y;        // вектор смотрит на первый шар (на b1 из b2)
 +
var r2 = rx * rx + ry * ry;                        // квадрат расстояния между шарами
 +
var rLen = Math.sqrt(r2);                          //расстояние между шарами
 +
if (rLen <= aCut) {
 +
 +
b1.ro += b2.m * (5 / (Pi * aCut * aCut)) * (1 + 3 * rLen / aCut) * (1 - rLen / aCut)* (1 - rLen / aCut)* (1 - rLen / aCut);  // ядро Люси
 +
//console.log(b1.ro);
 +
}
 +
}
 +
}
 +
balls[i] = b1;
 +
 
 +
}
 +
 +
/* if (count != balls.length){ //изменение плотности при появлении нового шарика
 +
 
 +
var ro0 = balls[0].ro
 +
for(var i = 1; i < balls.length; i++){
 +
if (balls[i].ro < ro0){
 +
ro0 = balls[i].ro;
 +
}
 +
 +
}
 +
 
 +
for (var i = 0; i < balls.length; i++){
 +
balls[i].ro0 = ro0;
 +
}
 +
count = balls.length;
 +
 
 +
}
 +
*/
 +
 +
for (var i = 0; i < balls.length; i++) {     //расчет давлений для каждой частицы
 +
balls[i].p = Bb * (Math.pow(balls[i].ro/balls[i].ro0,γ) - 1);
 +
// balls[i].p = Bb*(balls[i].ro-ro0);                                      //с сайта кафедры
 +
}
 +
 
 +
 +
for (var i = numberOfBorder; i < balls.length - countFence; i++) { //расчет ускорений
 +
var b1 = balls[i];
 +
                for (var j = 0; j < balls.length; j++) {                 
 +
                    var b2 = balls[j];
 +
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))){
 +
 
 +
var rx = b1.x - b2.x; 
 +
var ry = b1.y - b2.y;        // вектор смотрит на первый шар (b1 из "нуля")
 +
var r2 = rx * rx + ry * ry;                        // квадрат расстояния между шарами
 +
var rLen = Math.sqrt(r2);                          //расстояние между шарами
 +
 +
var vx = b1.vx - b2.vx;
 +
var vy = b1.vy - b2.vy;
 +
 
 +
if (rLen <= aCut) {
 +
var PabX = -alpha * aCut * Cab /((b1.ro + b2.ro)/2) * vx * rx/(rLen * rLen + 0.01 * aCut * aCut);
 +
var PabY = -alpha * aCut * Cab /((b1.ro + b2.ro)/2) * vy * ry/(rLen * rLen + 0.01 * aCut * aCut);
 +
 +
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);
 +
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);
 +
}
 +
}
 +
}
 +
 +
                if (b1 == dNd) continue;  // если шар схвачен курсором - его вз. со стенами и перемещение не считаем
 +
 
 +
//                if (b1.y + r > h) { b1.fy += -CWall * (b1.y + r - h) - BWall * b1.vy; }
 +
//                if (b1.y - r < 0) { b1.fy += -CWall * (b1.y - r) - BWall * b1.vy;}
 +
//                if (b1.x + r > w) { b1.fx += -CWall * (b1.x + r - w) - BWall * b1.vx; }
 +
//                if (b1.x - r < 0) { b1.fx += -CWall * (b1.x - r) - BWall * b1.vx; }
 +
balls[i] = b1;
 +
}
 +
    for (var i = numberOfBorder; i < balls.length - countFence; i++){
 +
 +
b1 = balls[i];
 +
b1.vx += (b1.fx + b1.ax) * dt;       
 +
b1.vy += (b1.fy + b1.ay + mg / m) * dt;
 +
b1.x += b1.vx * dt;           
 +
b1.y += b1.vy * dt;
 +
b1.latX = Math.floor(b1.x / aCut);
 +
b1.latY = Math.floor(b1.y / aCut);
 +
balls[i] = b1;
 +
// impAverageX += m * b1.vx;
 +
 +
}
 +
//console.log(impAverageX/balls.length)
 +
//console.log(numberOfBorder);
 +
 +
    }
 +
}
 +
 
 +
    // Рисование
 +
  context.fillStyle = "#d3692e";
 +
    function draw() {
 +
        context.clearRect(0, 0, w * scale, h * scale);      // очистить экран
 +
        for (var i = 0; i < balls.length; i++){
 +
            var xS = balls[i].x * scale;          var yS = balls[i].y * scale;
 +
            context.beginPath();
 +
            context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
 +
            context.closePath();
 +
            context.fill();
 +
        }
 +
    }
 +
 
 +
    // Запуск системы
 +
    this.newSystem();
 +
    setInterval(control, 1000 / fps);
 +
    // след. функция обновляет информацию о количестве частиц на поле
 +
    setInterval(function(){document.getElementById('ballsNum').innerHTML = balls.length;}, 1000 / 20);
 +
}
 +
 
 +
</syntaxhighlight>
 +
</div>
  
 
== Возможности программы ==
 
== Возможности программы ==
  
−
== Анализ ==
+
В текущей версии программы можно рассматривать разрушение столба жидкости в изначально разных конфигурациях:
 +
 
 +
- квадратная решетка размером 10x38;
 +
 
 +
- треугольная решетка размером 11x38;
 +
 
 +
- пустое поле, в котором можно создать собственную конфигурацию (ЛКМ - создание частицы).
 +
 
 +
В любой момент можно поставить или убрать забор. Также есть возможность изменения гравитации от 0g до 5g.
 +
 
 +
== Результаты ==
 +
 
 +
[[Файл:SPH_Meow_1.png|thumb|Рисунок 1. Начальная конфигурация|слева|450px]]
 +
[[Файл:SPH_Meow_2.png|thumb|Рисунок 2. Успокоившийся столб жидкости|справа|450px]]
 +
 
 +
[[Файл:SPH_Meow_3.png|thumb|Рисунок 3. Отскок от границы|слева|450px]]
 +
[[Файл:SPH_Meow_4.png|thumb|Рисунок 4. Осевшее на дно|справа|450px]]
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
 +
 
  
 
== Ссылки по теме ==
 
== Ссылки по теме ==
 +
 +
* J.J.Monaghan "Smoothed particle hydrodynamics", 1992;
 +
 +
* W.G.Hoover, C.G.Hoover "SPAM-Based Recipes for Continuum Simulation".
  
 
[[Category: Механика дискретных сред]]
 
[[Category: Механика дискретных сред]]

Текущая версия на 15:18, 10 января 2017

Курсовые работы по ВМДС: 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".
Источник — «http://tm.spbstu.ru/?title=Моделирование_жидкости_методом_SPH&oldid=45740»