§ 13. Моделирование движения поверхности жидкости

Производная Лагранжа

Эта производная является частным случаем (с точки зрения физики) полной производной, а полная производная — это производная от функции по какому-то параметру (обычно, это время), которая учитывает зависимость ее аргументов от этого параметра.

Законы сохранения для жидкости

Общий закон сохранения для свойств жидкости (массы, импульса и энергии) можно сформулировать так: свойство в замкнутом объеме изменяется только за счет обмена с внешней средой (соседними объемами) и за счет источников/стоков внутри этого объема. В общем случае этот закон можно записать так.

Здесь $\phi$ — это свойство жидкости, $\vec\upsilon$ — скорость жидкости, $s$ — источники/стоки (внешние и внутренние силы). Используя теорему Остроградского—Гаусса, заменяем интеграл по площади на интеграл по объему и получаем

Теперь мы можем внести производную по времени под знак интеграла, предположив, что объем не изменяется по времени.

Поскольку это уравнение должно выполняться для любого замкнутого объема, то мы можем убрать знак интеграла. Тогда получаем

Наконец, раскрыв градиент, мы можем переписать уравнение через производную Лагранжа.

• Закон сохранения массы

  $\rho_t + \nabla\cdot\left(\rho\vec\upsilon\right) = 0, \qquad s \equiv 0, \qquad \phi \equiv m = \rho V.$

• Закон сохранения импульса

  $\frac{D}{Dt}\left(\rho\vec\upsilon\right) = \nabla\cdot\sigma + F_\text{внешние}, \qquad s \equiv \nabla\cdot\sigma + F_\text{внешние}, \qquad \phi \equiv m \vec\upsilon = \rho V \vec\upsilon.$

  Здесь $\sigma$ — тензор напряжений, который можно записать как $\sigma = -p I + \tau$. Здесь $p$ — давление, $I$ — единичная матрица, $\tau$ — модифицированный тензор напряжений.

• Закон сохранения энергии

  $\frac{\partial}{\partial t}\left(\rho\frac{1}{2}|\vec\upsilon|^2\right) + \nabla\cdot\left(\rho\frac{1}{2}|\vec\upsilon|^2\vec\upsilon\right) + s = 0, \qquad s \equiv \nabla\cdot\sigma + F_\text{внешние}, \qquad \phi \equiv m \frac{|\vec\upsilon|^2}{2} = \rho V \frac{|\vec\upsilon|^2}{2}.$

Эти три уравнения являются основой для решения основных задач гидродинамики. Конкретная задача получается путем добавления дополнительных уравнений, задающих начальные и граничные условия, а также состояние системы.

Волны на границе раздела двух жидкостей

Для описания таких волн используют кинематическое и динамическое граничные условия.

• Кинематическое граничное условие заключается в том, что частицы жидкости на границе двигаются с той же скоростью, что и сама граница: проекция скорости частицы жидкости на нормаль к поверхности равна проекции скорости поверхности на ту же нормаль.

  $\vec\upsilon_1\cdot\vec{n} = \vec\upsilon_2\cdot\vec{n}.$

• Динамическое граничное условие заключается в том, что давление на границе должно быть непрерывным: внутреннее давление одной жидкости равно внутреннему давлению другой.

  $p_1 = p_2.$

В этих условиях граница раздела двух жидкостей подвижная, это означает, что на каждом временном шаге решения находящиеся на границе узлы придется перемещать в соответствии со скоростями жидкостей. Если же граница проходит не по узлам, а сквозь ячейки, то придется изменять сетку, чтобы учесть граничные условия. Этот процесс ресурсоемкий для компьютера, поэтому для описания волн на границе двух жидкостей используют упрощенные модели.

Волны малых амплитуд

Высота таких волн должна быть много меньше их длины. Под это описание подходят океанские волны кроме волн-убийц и волн на мелководье. Поскольку такие волны не опрокидываются, то границу раздела двух жидкостей можно задать функцией, которая определяет координату $z$ поверхности:

Тогда граничные условия перепишутся следующим образом.

Итого, кинематическое граничное условие упростилось, а динамическое теперь можно использовать для подстановки атмосферного давления как константы в закон сохранения импульса.

Также в рамках теории волн малых амплитуд упрощаются и сами законы сохранения: движение предполагается безвихревым, а жидкость несжимаемой.

Используя равенство $\nabla\times(\nabla\phi) = 0$, представим вектор скорости как градиента некоторой скалярной функции — потенциала скорости $\phi$.

Теперь можно переписать законы сохранения и граничные условия, используя потенциал скорости.

• Закон сохранения массы

  $\nabla^2 \phi = 0.$

• Закон сохранения импульса (динамическое граничное условие)

  $\rho\frac{\partial}{\partial t}\nabla\phi = -\nabla p + F, \\ \rho\nabla\left(\frac{\partial}{\partial t}\phi + p\right) = F, \\ \rho\nabla\left(\frac{\partial}{\partial t}\phi + p + gz\right) = 0, \\ \frac{\partial}{\partial t}\phi + p + gz = 0.$

• Кинематическое граничное условие

  $\frac{D}{Dt} \zeta = \nabla\phi\cdot\vec{n}$

  или его линеаризованная форма

  $\frac{\partial\zeta}{\partial t} = \frac{\partial\phi}{\partial z}$

Линеаризованная система этих уравнений имеет аналитическое решение в виде ряда Фурье, слагаемые которого домножены на экспоненту.

Волны малых амплитуд — это наиболее простая и распространенная теория, которая используется для моделирования поведения волн на компьютере. Для этого достаточно сгенерировать случайные фазы для каждой суммируемой гармоники, а амплитуды вычислить из спектра. В итоге получится следующая формула для взволнованной поверхности.

Здесь $u_n,v_n$ — волновые числа, $\omega_n$ — угловая частота, $\epsilon_n$ — случайная фаза, $S$ — частотный спектр морских волн. Спектр можно получить с помощью аналитических аппроксимаций или в дискретном виде из метеонаблюдений. Модель получила название Лонге—Хиггинса.

Трохоидальные волны Герстнера

Эти волны являются более точным описанием взволнованной поверхности. Модель основана на описании движения отдельных частиц. Раньше мы рассматривали конечные объемы, через границы которых входили и выходили частицы. Такая система координат называется эйлеровой. Этой системе противопоставляются лагранжевы координаты, которые задают положение отдельных частиц, а не конечного объема. Волны Герстнера являются решением системы уравнений из предыдущего раздела, записанного в такой системе координат.

Здесь $a,b,c$ — координаты неподвижных узлов сетки, вокруг которых вращаются частицы, $x,y,z$ — координаты частиц.

Волны Герстнера отличаются наличием вихрей, а значит, нельзя использовать потенциал скорости для вычисления силы давления под взволнованной поверхностью. Этим ограничением часто пренебрегают, поскольку альтернативного аналитического метода просто не существует. Также эти волны используют для спецэффектов в компьютерных играх и фильмах, поскольку волны выглядят реалистично, а формулы легко переписать с помощью преобразований Фурье, которые эффективно вычисляются на процессоре и видеокарте.

Гидродинамика сглаженных частиц

Это один из так называемых бессеточных методов, то есть методов, которые используют частицы для моделирования жидкости, а не конечные объемы. Взаимодействие частиц происходит в заданном радиусе от конкретно взятой частицы. Параметры частицы (плотность, объем, давление, скорость и т.д) определяются как взвешенная сумма параметров окружающих частиц, а в качестве весовой функции используюется гауссиан или аналогичная сглаживающая функция. Раньше мы аппроксимировали производные линейной комбинацией значений функции в окружающих узлах или ячейках сетки, а теперь мы аппроксимируем их взвешенной суммой значений функции в окружающих частицах.

Это основное уравнение для гидродинамики сглаженных частиц, которое определяет некоторое свойства в виде интеграла со сглаживающим ядром (взвешенной суммы). Здесь $Q$ — физическая величина (плотность, скорость, давление и т.п.), $\vec{x}$ — координата частицы, $h$ — радиус взаимодействия, $W$ — сглаживающее ядро (весовая функция).

Вычисление производных так же делается по этим формулам.

Для этого нам надо знать только производную весовой функции.

Для плотности уравнение принимает следующий дискретный вид.

Здесь индексы $i,j$ обозначают конкретные частицы.

Теперь перепишем известные нам уравнения с помощью этой формулы.

• Закон сохранения массы примет вид

  $\frac{d}{dt}\rho_i = \sum\limits_j \rho_j V_j (\vec\upsilon_i - \vec\upsilon_j) \nabla W_{ij}.$

• Закон сохранения импульса примет вид

  $\rho \frac{d}{dt}\vec\upsilon = -\nabla p + F,\\ \frac{d}{dt}(\rho_i\vec\upsilon_i) = -\sum\limits_j m_i m_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) \nabla W_{ij} + f_i.$

В качестве сглаживающего ядра выбирают следующие функции.

• Гауссиан

  $W(\vec{x},h) = \left(\frac{1}{h\sqrt{\pi}}\right)^d \exp\left(-\frac{|\vec{x}|^2}{h^2}\right).$

  Здесь $d$ — количество измерений.

• Полиномиальную аппроксимацию

  $W(\vec{x},h) = C \begin{cases} (2 - q)^3 - 4(1 - q)^3, & 0 \leq q \lt 1, \\ (2 - q)^3, & 1 \leq q \lt 2, \\ 0 & q \geq 2. \\ \end{cases}, \qquad q = \frac{|\vec{x}|}{h}, \qquad C = \begin{cases} \frac{1}{6h}, & d = 1, \\ \frac{5}{15\pi h^2}, & d = 2, \\ \frac{1}{4\pi h^3}, & d = 3. \\ \end{cases}$

Для вычисления давления используют аппроксимацию уравнения состояния идеального газа.

Задания

Волны Герстнера1 балл

Сгенерируйте трехмерную взволнованную морскую поверхность, сформированную тремя волнами с произвольно выбранными параметрами (амплитудой, направлением, частотой, длиной). Для этого создайте равномерную сетку на поверхности, и в каждом узле сетки поставьте начальную координату частицы. Затем считайте новые координаты частиц по формулам.

Волны малых амплитуд2 балла

Используйте модель Лонге—Хиггинса для моделирования трехмерной взволнованной морской поверхности. Для этого аппроксимируйте спектр с помощью модели Пирсона—Московица.

Здесь $S$ — частотный спектр морского волнения, $Q$ — угловое распределение энергии, $\omega=2\pi f$ — угловая частота, $f$ — частота волн, $\theta$ — направление волн, $g$ — ускорение свободного падения, $\alpha=8.1\times 10^{-3}, \beta=0.74$ — коэффициенты, $\omega_0=g/U$ — характерная угловая частота, $U$ — средняя скорость ветра.

Перейти между частотой и направлением и волновыми числами можно, используя дисперсионное соотношение для морских волн.

Случайные фазы распределены равномерно на промежутке $[0,2\pi]$.

Сглаженные частицы3 балла

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

git clone https://courses.igankevich.com/numerical-methods/sph.git
cd sph
make
cd build
./sph

Видео

2022

Доска

2022