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

Курсовые работы по МДС: 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]