Моделирование поля скорости течения на свободной поверхности стратифицированной жидкости

Полный текст

1. Введение

Внутренние гравитационные волны существуют в толще океана из-за неравномерности плотности воды, вызванной неоднородностями солености и температуры, загрязнением. Они могут быть вызваны взаимодействием подводных течений с неровностями дна, вершинами подводных хребтов, подводными кабелями или трубопроводами, а также перемещением крупных рыб или косяков промысловой рыбы. Моделирование гравитационных волн важно при изучении океана, наблюдениях за загрязнениями, при прогнозировании погоды [1–3]. Гравитационные волны могут иметь амплитуду более ста метров, могут распространяться на километры от источника возмущений и могут быть обнаружены как непосредственно подводными датчиками [4, 5], которые измеряют колебания солености, так и с помощью радиометрии [1–3, 6–8]. В последнем случае регистрируется влияние гравитационных волн на ветровую рябь на свободной поверхности океана [9, 10].

Ранее решалась задача расчета поля скоростей в толще равномерно [9] или произвольно стратифицированной [11] жидкости, на поверхности идеальной однородной жидкости [9, 12]. Асимптотические и численные исследования распространения гравитационных волн в толще жидкости проводились для постоянной [13–15] или произвольной (по вертикали) стратификации [16]. В настоящей работе рассматривается задача моделирования поля скорости течений на свободной поверхности идеальной стратифицированной жидкости, сформировавшегося под действием внутренних гравитационных волн, распространяющихся в ее толще. Полученные поля могут быть далее использованы при расчете результата взаимодействия ветровой ряби с поверхностным течением [17], что важно для моделирования формы возмущения морской поверхности, при радиометрических исследованиях океана [1, 2, 10].

2. Постановка задачи

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

Рис. 1. Область Ω

Для нахождения вертикальной компоненты скорости жидкости v_z в области Ω = [x_min, x_max] × [y_min, y_max] × [–H, H] ищется решение уравнения [18, 19]

$\Delta \frac{{\partial }^2}{\partial t^2}{v}_z+N{\left(z\right)}^2{\Delta }_{xy}{v}_z=\frac{{\partial }^3}{\partial t^2\partial z}f(x-r(t\left)\right); x=(x, y, z)\in \Omega$ (2.1)

с условиями свободной поверхности в плоскости z = H:

$\frac{{\partial }^3}{\partial t^2\partial z}{v}_z=g{\Delta }_{xy}{v}_z; z=H$, (2.2)

и нулевыми условиями на оставшейся части границы области Ω:

$v_z=0; x\in \partial \Omega , z\ne H$ (2.3)

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

$v_z{|}_{t=0}=0$ (2.4)

Функция f задает источник:

$f(x, y, z)=\frac{B{A}^3}{\sqrt{{\pi }^3}}{e}^{-A^2(x^2+y^2+z^2)}$

Заметим, что ${\int }_{R^3}f\left(x\right)dx=B$. В случае равномерного движения источника возмущения его положение определяется функцией r (t) = x₀ + Ut, где U – скорость прямолинейного движения массового источника.

Можно выписать задачи для остальных компонент скорости [18], а именно, для компоненты v_x:

$\Delta \frac{{\partial }^2}{\partial t^2}{v}_x+N{\left(z\right)}^2{\Delta }_{xy}{v}_x=\left(\frac{\partial }{\partial t^2}+N{\left(z\right)}^2\right)\frac{{\partial }^2}{\partial t\partial x}f(x-x_0-Ut); x\in \Omega$ (2.5)

с граничными условиями

$\frac{{\partial }^2}{\partial t^2}{v}_x=-g\frac{\partial }{\partial x}{v}_z; z=H$ (2.6)

$v_x=0; x\in \partial \Omega , z\ne H$ (2.7)

начальным условием

$v_x{|}_{t=0}=0$, (2.8)

и для компоненты v_y:

$\Delta \frac{{\partial }^2}{\partial t^2}{v}_y+N{\left(z\right)}^2{\Delta }_{xy}{v}_y=\left(\frac{\partial }{\partial t^2}+N{\left(z\right)}^2\right)\frac{{\partial }^2}{\partial t\partial y}f(x-x_0-Ut); x\in \Omega$ (2.9)

$\frac{{\partial }^2}{\partial t^2}{v}_y=-g\frac{\partial }{\partial y}{v}_z; z=H$ (2.10)

$v_y=0; x\in \partial \Omega , z\ne H$ (2.11)

$v_y{|}_{t=0}=0$ (2.12)

Для нахождения полей компонент скорости v_x (x, t) и v_y (x, t) во всей области Ω нужно решить задачу (2.5)–(2.8) и задачу (2.9)–(2.12), при этом в правых частях условий (2.6) и (2.10) используется функция v_z, найденная ранее в результате решения задачи (2.1)–(2.4).

3. Моделирование поля скоростей на свободной поверхности

Пусть требуется найти горизонтальные компоненты поля скоростей жидкости v_x, v_y не во всей области Ω, а только на свободной поверхности z = H. Тогда сначала следует решить задачу (2.1)–(2.4), то есть найти вертикальную компоненту скорости v_z (x, t), на всем временном отрезке [0, T], где T – полное время моделирования. Далее, найденное поле вертикальной компоненты скорости подставляется в правую часть соотношений (2.6) и (2.10) и для каждой точки свободной поверхности (x_m, y_n, H) решаются соответствующие задачи Коши для обыкновенных дифференциальных уравнений второго порядка (2.6) и (2.10), в результате чего определяются горизонтальные компоненты скорости v_x (x, y, H, t) и v_y (x, y, H, t), в точках свободной поверхности на всем интервале времени t ∈ [0, T].

3.1. Численная схема для нахождения вертикальной компоненты скорости жидкости v_z в области Ω. Для решения краевой задачи (2.1)–(2.4) в области Ω используется неявная разностная схема на равномерной расчетной сетке, состоящей из N_x × N_y × N_z узлов. Обозначим

$$\begin{gathered} \psi \left(s\right)=v_z\left(qs\right) \\ {\psi }_{m,n,k}^s=v_z(x_{\min }+h(m-1), y_{\min }+h(n-1), -H+h(k-1), qs) \\ f(m,n,k,s)=f(x_{\min }+h(m-1), y_{\min }+h(n-1), -H+h(k-1), qs), \end{gathered}$$

где m = 1, ..., N_x, n = 1, ..., N_y, k = 1, ..., N_z, s = 2, ..., T / q, h – шаг по пространству, q – шаг по времени. На каждом шаге по времени s решается система линейных алгебраических уравнений (СЛАУ)

$Q{\psi }^{s+1}=r^s$ (3.1)

для неизвестных ${\psi }_{m,n,k}^{s+1}$, m = 1, ..., N_x, n = 1, ..., N_y, k = 1, ..., N_z. Для апроксимации вторых производных по времени и по пространственным переменным в (2.1) используются центральные разности:

$\frac{{\Delta }_{xyz}\psi (s+1)-2{\Delta }_{xyz}\psi \left(s\right)+{\Delta }_{xyz}\psi (s-1)}{q^2}+{\Delta }_{xy}\psi (m, n, k, s)=\frac{{\partial }^3}{\partial t^2\partial z}f(m, n, k, s)$

Здесь Δ_xyz, Δ_xy – соответствующие разностные операторы. Таким образом, внутри области Ω, то есть для 1 < m < N_x, 1 < n < N_y, 1 < k < N_z имеем:

$$\begin{gathered} \frac{1}{h^2{q}^2}({\psi }_{m-\mathrm{1,}n,k}^{s+1}+{\psi }_{m+\mathrm{1,}n,k}^{s+1}+{\psi }_{m,n-\mathrm{1,}k}^{s+1}+{\psi }_{m,n+\mathrm{1,}k}^{s+1}+{\psi }_{m,n,k-1}^{s+1}+{\psi }_{m,n,k+1}^{s+1}-6{\psi }_{m,n,k}^{s+1})= \\ =-\frac{-2{\Delta }_{xyz}\psi \left(s\right)+{\Delta }_{xyz}\psi (s-1)}{q^2}-{\Delta }_{xy}\psi \left(s\right)+\frac{{\partial }^3}{\partial t^2\partial z}f(m, n, k, s) \end{gathered}$$ (3.2)

На свободной поверхности, то есть для k = N_z граничное условие (2.2) апроксимируется следующим образом:

$\frac{1}{2h}\left(3\frac{{\partial }^2}{\partial t^2}{\psi }_{m,n,k}-4\frac{{\partial }^2}{\partial t^2}{\psi }_{m,n,k-1}+1\frac{{\partial }^2}{\partial t^2}{\psi }_{m,n,k-2}\right)=g{\Delta }_{xy}\psi \left(s\right)$, (3.3)

где $\frac{{\partial }^2}{\partial t^2}{\psi }_{m,n,k}$ обозначает апроксимацию центральной разностью второй производной по времени:

$\frac{{\partial }^2}{\partial t^2}{\psi }_{m,n,k}\equiv ({\psi }_{m,n,k}^{s+1}-2{\psi }_{m,n,k}^s+{\psi }_{m,n,k}^{s-1})/q^2$

На остальной части границы области $\partial \Omega$, то есть для m = 1, m = N_x, n = 1, n = N_y, k = 1, из условия (2.3) имеем:

${\psi }_{m,n,k}^{s+1}=0$ (3.4)

Коэффициенты матрицы Q системы (3.1) определяются левыми частями соотношений (3.2)–(3.4), а вектор правой части r^s – правыми частями этих соотношений.

3.2. Расчет горизонтальных компонент скорости v_x  и v_y на свободной поверхности. Опишем подробнее процедуру нахождения горизонтальной x-компоненты скорости жидкости на свободной поверхности, v_x (x, y, H, t), (x, y) ∈ [x_min, x_max] × [y_min, y_max], t ∈ [0, T]. Выпишем условие (2.6) в точке свободной поверхности (x_m, y_n), где x = x_min + h (m – 1), y = y_min + h (n – 1), m = 1, ..., N_x, n = 1, ..., N_y:

$\frac{{\partial }^2}{\partial t^2}\zeta \left(t\right)=-g\frac{\partial }{\partial x}{v}_z(x_m, y_n, H, t)$, (3.5)

где обозначено ζ (t) = (x_m, y_n, H, t). Зададим начальное условие

$\zeta \left(0\right)=\mathrm{0,} \frac{\partial }{\partial t}\zeta \left(0\right)=0$ (3.6)

Функция v_z правой части (3.5) найдена ранее в результате решения задачи (2.1)–(2.4). Для апроксимации второй производной по времени используется схема центральных разностей:

${\zeta }^{s+1}=2{\zeta }^s-{\zeta }^{s-1}-q^2{v}_z(x_m, y_n, H, (s+1\left)q\right); s=\mathrm{2,} \dots , T/q$,

где ζ^s ≡ ζ (sq), ζ⁰ = ζ¹ = 0. Решив задачу Коши (3.5), (3.6) для всех точек (x_m, y_n), m = 1, ..., N_x, n = 1, ..., N_y найдем поле v_x на свободной поверхности. Аналогично (с использованием условия (2.9)) расчитывается поле горизонтальной компоненты скорости на свободной поверхности v_y.

3.3. Компьютерная реализация. Описанный выше алгоритм был реализован в компьютерной программе на языке программирования С++. Разработанная программа позволяет моделировать распространение гравитационных волн от движущегося точечного источника возмущений, рассчитывать изменение во времени поля вертикальной компоненты скорости жидкости v_z во всей области Ω, изменение во времени горизонтальных компонент скорости жидкости v_x и v_y на свободной поверхности жидкости. При решении краевой задачи (2.1)–(2.4) на каждом шаге s по времени решается СЛАУ (3.1) с матрицей Q определенной соотношениями (3.2)–(3.4), задающими численную схему. Матрица Q разреженная: каждая ее строка содержит всего 7 ненулевых элементов (по количеству используемых узлов сетки в шаблоне численного метода), поэтому для ее хранения используется специальный сжатый формат. Это позволяет минимизировать требования к объему оперативной памяти и вести расчет на сетках до 220 × 220 × 220 на персональном компьютере или до 400 × 400 × 400 на одном вычислительном узле суперкомпьютера. Для решения системы (3.1) используется обобщенный метод минимальной невязки (Generalized minimal residual method, GMRES). Этот итерационый проекционный метод решения СЛАУ был предложен в [20]. Расчет одного шага по времени в типовом расчете на сетке 250 × 250 × 250 занимает примерно 30 минут на одном вычислительном узле суперкомпьютера МСЦ РАН (Intel(R) Xeon(R) CPU E5530@2.40GHz, 8Gb ОЗУ). На s-м шаге по времени выполняется несколько итераций решения системы (3.1), при этом контролируется величина невязки e (s, j),

$e(s, j)=\frac{\left|\right|Q{\psi }^{s+\mathrm{1,}j}-r^s\left|\right|}{\left|\right|r^s\left|\right|}$,

где ψ^{s + 1, j} – j-е приближение к ψ^{s + 1}. Таким образом, на каждом шаге выполняется 5–10 итераций, в результате чего величина невязки e (s, j) падает на 1–2 порядка до значений около 10^–2. Весь такой расчет (до выхода источника возмущений за пределы области Ω) занимает около суток. В программе реализовано периодическое сохранение текущих результатов, что позволяет продолжить счет в случае сбоя или когда требуемое время расчета превышает максимальное допустимое на суперкомпьютере время непрерывного счета. Таким образом находится v_z во всей области Ω для t ∈ [0, T]. Далее решается набор задач (3.5), (3.6), и в результате получается искомое поле течений (v_x, v_y) на свободной поверхности z = H при t ∈ [0, T]. Для реализации элементарных операций с разрежеными матрицами и решения системы (3.1) используется пакет GNU Scientific Library (GSL) [21]. Исходный код программы доступен в интернете по адресу https://bitbucket.org/Jclash/inwaves.

4. Сравнение с асимптотиками

Для верификации правильности расчетов было проведено сравнение результатов расчета копмьютерной программы с известными ассимптотическими решениями. Пусть жидкость равномерно стратифицирована (N = const), а гравитационные волны возбуждаются точечным массовым источником, т.е. функция f правой части уравнений (2.1), (2.5), (2.9) задана как

$f(x, y, z)=Q\delta \left(x\right)\delta \left(y\right)\delta \left(z\right)$

Источник возмущения движется равномерного и прямолинено со скоростью V в неограниченной области, то есть его движение описывается функцией r (t) = x₀ + Ut, где |U| = V. В этом случае известны ассимтотики дальнего поля для скорости жидкости v, когда источник движется горизонтально [13] или под фиксированным углом к горизонту [14].

4.1. Горизонтальное движение. В случае, когда источник движется горизонтально прямолинейно и с постоянной скоростью, скорость жидкости в дальнем поле r₁ >> V / N для x > 0 может быть вычислена с использованием следующей формулы [13] в сферической системе координат (r, ξ, φ) с центром в точке текущего положения точечного источника r (t) (рис. 2):

$v(x, y, z)=H\left(x_1\right)\frac{NQ}{2\pi V{r}_1}\frac{\sqrt{{\sin }^2\xi +{\cos }^2\xi {\cos }^2\phi }}{\sin \xi }\frac{W}{\left|W\right|}\cos \left(\frac{N}{V}{r}_1\left|\sin \phi \right|\right)$ (4.1)

Рис. 2. Схема расчета асимптотики поля скорости при горизонтальном движении массового источника из [13]

Здесь

$W=(-r_{1}tg\xi \sin \xi , r_1\sin \xi \cos \phi , r_1\sin \xi \sin \phi )$,

x₁ = x – r, r₁ = (x₁, y₁, z₁) – вектор из точки текущего положения источника O₁ в точку M с координатами (x, y, z) в которой расчитывается поле скоростей, $r_1=\sqrt{x_1^2+y_1^2+z_1^2}$, ξ = arccosx₁ / r₁, $\phi =\arccos \left(y_1/\sqrt{y_1^2+z_1^2}\right)$, H (x) – функция Хевисайда,

$H\left(x\right)=\left\{\begin{array}{ll}\mathrm{1,}& x<0\\ 1/2,& x=0\\ \mathrm{0,}& x>0\end{array}\right.$

Пример сравнения результата численного моделирования и асимптотики (4.1) показан на рис. 3.

Рис. 3. Моделирование горизонтального движения в идеальной равномерно стратифицированной жидкости, N = 0.8. Скорость движения V = 1 м/с, вертикальная составляющая скорости жидкости v_z (x, 1, z) показана в вертикальном сечении y = 1. Время движения T = 125 с для (a). Показано численное решение (а) и аналитическая аппроксимация (б)

4.2. Движение под фиксированным углом к горизонту. Аналогичная асимптотика может быть выписана в случае, когда источник возмущения движется равномерно под фиксированным углом γ к горизонту, то есть U = (Vcosγ, 0, Vsinγ). Асимптотика дальнего поля для скорости в точке M в системе координат (r₁, θ, φ) с центром в точке текущего положения источника O₁ (рис. 4) имеет следущий вид [14]:

$v=\frac{NQ}{2\pi c_g}\frac{{(R_x^2+R_y^2)}^{1/2}}{R^2}\frac{R}{R}\frac{\cos \left(N\right|R_z|/c_g+\pi H\left(A\right)/2)}{|A{|}^{1/2}}$, (4.2)

где

$c_g=V\cdot \left(\frac{R}{R}-\frac{R}{R_z}{e}_z\right)\frac{R}{R}, A=-{\left(\frac{V}{c_g}\times \frac{R}{R}\right)}^2+2\frac{V}{c_g}\cdot \frac{R}{R_z}{e}_z$,

$R=r_1\left(\begin{array}{l}\sin \theta \cos \varphi {\cos }^2\alpha +\sin \alpha \sin \beta tg\eta -\sin \alpha \cos \alpha \cos \theta \\ \sin \theta \sin \varphi \\ \cos \theta {\sin }^2\alpha +\cos \alpha \sin \beta tg\eta -\sin \alpha \cos \alpha \sin \theta \cos \varphi \end{array}\right)$,

$\overrightarrow{OM}=(x_1, y_1, z_1)$, θ = arccos (x₁ / r₁), $\varphi =\arccos \left(x_1/\sqrt{y_1^2+x_1^2}\right)$, α = 90° – γ, и тригонометрическии функции углов β, ω и η находятся из следующих соотношений:

$$\begin{gathered} \cot \beta =\frac{\sin \theta \cos \varphi \sin \alpha +\cos \theta \cos \alpha }{\sqrt{{\sin }^2\theta {\sin }^2\varphi +{(\cos \theta \sin \alpha -\sin \theta \cos \varphi \cos \alpha )}^2}} \\ \cos \omega =\frac{(\cos \theta \sin \alpha -\sin \theta \cos \varphi \cos \alpha )\sin \alpha }{\sqrt{{\sin }^2\theta {\sin }^2\varphi +{(\cos \theta \sin \alpha -\sin \theta \cos \varphi \cos \alpha )}^2}} \\ t{g}^3\eta +\left(\frac{\cos \omega }{\cos \alpha }ctg\beta +2\right)tg\eta +\frac{\cos \omega }{\cos \alpha }-ctg\beta =0 \end{gathered}$$ (4.3)

Рис. 4. Схема расчета асимптотики поля скорости при движении массового источника под углом к горизонту из [14]

Поле скоростей (4.2) вычисляется как сумма волн, соответствующих действительным корням кубического уравнения (4.3), при которых выполнено двустороннее неравенство 0 < τ_s < T, где T – полное время движения,

${\tau }_s=T-\frac{r\sin \beta }{V{\cos }^2\eta }\frac{\cos \alpha tg\eta +\cos \omega }{\cos \omega tg\eta -\cos \alpha }$

Для вычисления асимптотик (4.1) и (4.2) была написана компьютерная программа на языке программирования Python. Было проведено сравнение между численными решениями задачи (2.1)–(2.4) в области Ω и приближениями (4.1) и (4.2). Пример такого сравнения показан на рис. 3 и рис. 5. Установлено, что результаты расчетов программы согласуются с этими асимптотиками на достаточном расстоянии от источника и линии движения, когда после начала движения прошло достаточно продолжительное время (то есть, когда картина возмущений, формирующихся вокруг движущегося источника, принимает установившийся характер).

Рис. 5. Моделирование распространения внутренних волн от массового источника, который движется равномерно и прямолинейно под фиксированным углом γ = 30° к горизонту, N = 0.8. Скорость движения V = 1 м/с, вертикальная составляющая скорости жидкости v_z (x, 1, z) показана в вертикальном сечении y = 1. Время движения T = 139.3 с. Показаны численное решение (а) и аналитическая аппроксимация (б)

5. Пример расчета

Пусть область Ω = [–100, 100] × [–50, 50] × [–25, 25] с размерами 200.0 × 100.0 × 50.0 заполнена однородно стратифицированной (частота плавучести N = 0.8) идеальной жидкостью, глубина 2H = 50.0 м. Промоделируем распространение внутренних гравитационных волн в области Ω, в качестве источника возмущения возьмем горизонтально движущий точечный массовый источник. Такая вытянутая вдоль линии движения форма области выбрана для того, чтобы за время движения T картина течений на поверхности (движущаяся вслед за источником возмущения) успела установиться. Ширина области y_max – y_min выбрана достаточной для того, чтобы края области (на которых к тому же задано нефизическое условие v_z = 0) не влияли на картину возмущения поверхности. Глубина области 2H выбрана достаточной, чтобы отраженное от дна (где задано условие непротекания v_z = 0, см. (2.3)) возмущение достаточно ослабло и не влияло на картину возмущений вблизи поверхности. Источник движется из точки с координатами (–75, 0, 20) со скоростью U = (1, 0, 0), то есть, движение происходит вдоль оси Ox на глубине 5 метров.

При решении задачи (2.1)–(2.4) разностным методом (3.1)–(3.4) использовалась расчетная сетка размером 378 × 186 × 63. Шаг по пространству брался равным h = 0.79, шаг по времени q = 0.079. Расчет проводился до момента T = 225 с, когда источник возмущения достигал границы области Ω. Один расчет занимал около двух суток на одном вычислительном узле суперомпьютера МВС-10П. Результаты расчета вертикальной компоненты скорости v_z (x, y, z, t) приведены на рис. 6. Для момента времени t = 175 с показаны (а) величина вертикальной компоненты v_z (x, y, H) на свободной поверхности z = H и (б) величина вертикальной компоненты v_z (x, 0, z) в вертикальном сечении y = 0, проходящем через линию движения источника возмущения.

Рис. 6. Величина вертикальной компоненты скорости vz на свободной поверхности (а) и в толще жидкости (б) в момент времени t = 175 с

Результаты решения задачи (2.1)–(2.4) загружались с вычислительного кластера на персональный компьютер, на котором уже решалась существенно менее требовательная к вычислительным ресурсам задача (3.5), (3.6). Таким образом определялись компоненты скорости v_x и v_y на свободной поверхности z = H. Эти компоненты показаны на рис. 7, а и б. Для исключения влияния переходных процессов на картину течения, в расчетах взаимодействия течения с ветровым волнением следует брать поля скоростей v_x, v_y на достаточном удалении от точки, находящейся над точкой, из которой начинается движения источника. Например, на рис. 8 показаны поля скоростей не на всей верхней границе области Ω, а только на учаcтке свободной поверхности [25, 125] × [–50, 50].

Рис. 7. Поле течений (v_x, v_y) на свободной поверхности z = H в момент времени t = 175 с: показаны x-компонента (а) и y-компонента (б)

Заключение

В настоящей работе решается задача моделирования распространения гравитационных волн от движущегося в идеальной стратифицированной жидкости массового источника. Частота плавучести может изменяться с высотой. Была разработана компьютерная программа, позволяющая рассчитывать изменение во времени поля скорости на свободной поверхности, вертикальной компоненты скорости в толще жидкости. Достоверность расчетов верифицирована сравнением с асимптотическими приближениями дальнего поля для случаев, когда жидкость равномерно стратифицирована, а источник движется равномерно и прямолинейно горизонтально [13] или под углом к горизонту [14]. Результаты расчетов согласуются с экспериментом по обтеканию сферы потоком стратифицированной жидкости [22].

Для отладки численных методов и для сравнения с известными асимптотиками считалось, что гравитационные волны возбуждаются точечным массовым источником. Однако в разработанном методе моделирования поверхностных течений возможно задать более сложные возмущающие воздействия: с помощью осциллирующего источника, нескольких движущихся источников. Например, представляет интерес задача моделирования возмущения от движущейся стаи промысловых рыб. С помощью разработанных компьютерных программ можно моделировать возмущение от такого сложного объекта: в толще жидкости расчет вести по асимптотическим формулам (4.1) или (4.2), а вблизи свободной поверхности использовать для расчета алгоритм, предложенный в разделе 3.

Результаты расчетов, приведенные выше, были получены как на персональных компьютерах, так и с использованием вычислительных кластеров Межведомственного суперкомпьютерного центра РАН (МСЦ РАН), г. Москва. Автор выражает глубокую признательность руководству и сотрудникам МСЦ РАН, предоставившим возможность и техническую поддержку этих расчетов.

Исследование выполнено за счет гранта Российского научного фонда № 24-61-00025, https://rscf.ru/project/24-61-00025/.

Об авторах

Д. Ю. Князьков

Автор, ответственный за переписку.
Email: knyaz@ipmnet.ru
Россия, Москва

Список литературы

Дополнительные файлы

Доп. файлы

Действие

1. JATS XML

Скачать

2. Рис. 1. Область Ω

Скачать (18KB)

Метаданные

3. Рис. 2. Схема расчета асимптотики поля скорости при горизонтальном движении массового источника из [13]

Скачать (5KB)

Метаданные

4. Рис. 3. Моделирование горизонтального движения в идеальной равномерно стратифицированной жидкости, N = 0.8. Скорость движения V = 1 м/с, вертикальная составляющая скорости жидкости vz (x, 1, z) показана в вертикальном сечении y = 1. Время движения T = 125 с для (a). Показано численное решение (а) и аналитическая аппроксимация (б)

Скачать (23KB)

Метаданные

5. Рис. 4. Схема расчета асимптотики поля скорости при движении массового источника под углом к горизонту из [14]

Скачать (6KB)

Метаданные

6. Рис. 5. Моделирование распространения внутренних волн от массового источника, который движется равномерно и прямолинейно под фиксированным углом γ = 30° к горизонту, N = 0.8. Скорость движения V = 1 м/с, вертикальная составляющая скорости жидкости vz (x, 1, z) показана в вертикальном сечении y = 1. Время движения T = 139.3 с. Показаны численное решение (а) и аналитическая аппроксимация (б)

Скачать (20KB)

Метаданные

7. Рис. 6. Величина вертикальной компоненты скорости vz на свободной поверхности (а) и в толще жидкости (б) в момент времени t = 175 c

Скачать (21KB)

Метаданные

8. Рис. 7. Поле течений (vx, vy) на свободной поверхности z = H в момент времени t = 175 c: показаны x-компонента (а) и y-компонента (б)

Скачать (20KB)

Метаданные