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

[[Курсовые работы по ВМДС: 2016-2017]] > '''Моделирование жидкости методом SPH''' <HR>

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

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

'''Группа:''' [[Группа 09|09]] (43604/1)

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

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

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

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

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

<big><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></big>

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

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

<big><math> W(r < h) = \frac{5}{π{b}^{2}}(1+3\frac{r}{h}){(1-\frac{r}{h})}^{3},                </math></big>

где <math>r = |\pmb r - \pmb r_{j}|</math>

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

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

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

<big><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></big>

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

<big><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></big>

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

<big><math>P = B (({\frac{ρ}{ρ_{0}}})^{\gamma} -1)</math></big>,

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

А <math>П_{ab}</math> - член, отвечающий за вязкость

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

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

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

== Реализация ==

{{#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: Механика дискретных сред]]