Моделирование поля скорости течения на свободной поверхности стратифицированной жидкости
- Авторы: Князьков Д.Ю.1
-
Учреждения:
- Институт проблем механики им. А.Ю. Ишлинского РАН
- Выпуск: Том 88, № 5 (2024)
- Страницы: 745-757
- Раздел: Статьи
- URL: https://journals.rcsi.science/0032-8235/article/view/280966
- DOI: https://doi.org/10.31857/S0032823524050074
- EDN: https://elibrary.ru/JPKACG
- ID: 280966
Цитировать
Полный текст
Аннотация
В работе моделируется поле скорости на свободной поверхности идеальной стратифицированной жидкости, порожденное вышедшими на поверхность внутренними гравитационными волнами. Написана компьютерная программа, позволяющая рассчитывать все компоненты поля скорости на поверхности. Показано, что результаты расчетов для вертикальной компоненты скорости согласуются с известными асимптотиками, полученными в приближении дальнего поля для случаев равномерного и прямолинейного движения точечного массового источника горизонтально (B. Voisin) или под фиксированным углом к горизонту (M.M. Scase и S.B. Dalziel) в равномерно стратифицированной жидкости.
Полный текст
Введение
Внутренние гравитационные волны существуют в толще океана из-за неравномерности плотности воды, вызванной неоднородностями солености и температуры, загрязнением. Они могут быть вызваны взаимодействием подводных течений с неровностями дна, вершинами подводных хребтов, подводными кабелями или трубопроводами, а также перемещением крупных рыб или косяков промысловой рыбы. Моделирование гравитационных волн важно при изучении океана, наблюдениях за загрязнениями, при прогнозировании погоды [1–3]. Гравитационные волны могут иметь амплитуду более ста метров, могут распространяться на километры от источника возмущений и могут быть обнаружены как непосредственно подводными датчиками [4, 5], которые измеряют колебания солености, так и с помощью радиометрии [1–3, 6–8]. В последнем случае регистрируется влияние гравитационных волн на ветровую рябь на свободной поверхности океана [9, 10].
Ранее решалась задача расчета поля скоростей в толще равномерно [9] или произвольно стратифицированной [11] жидкости, на поверхности идеальной однородной жидкости [9, 12]. Асимптотические и численные исследования распространения гравитационных волн в толще жидкости проводились для постоянной [13–15] или произвольной (по вертикали) стратификации [16]. В настоящей работе рассматривается задача моделирования поля скорости течений на свободной поверхности идеальной стратифицированной жидкости, сформировавшегося под действием внутренних гравитационных волн, распространяющихся в ее толще. Полученные поля могут быть далее использованы при расчете результата взаимодействия ветровой ряби с поверхностным течением [17], что важно для моделирования формы возмущения морской поверхности, при радиометрических исследованиях океана [1, 2, 10].
Постановка задачи
Рассмотрим задачу моделирования распространения внутренних гравитационных волн, порожденных движением точечного массового источника в стратифицированной идеальной жидкости (рис. 1). Изменения плотности жидкости по отношению к базовой стратификации малы, поэтому для записи уравнения движения можно воспользоваться приближением Буссинеска.
Рис. 1. Область Ω
Для нахождения вертикальной компоненты скорости жидкости vz в области Ω = [xmin, xmax] × [ymin, ymax] × [–H, H] ищется решение уравнения [18, 19]
(2.1)
с условиями свободной поверхности в плоскости z = H:
, (2.2)
и нулевыми условиями на оставшейся части границы области Ω:
(2.3)
где N (z) – частота плавучести. Стратификация может быть неравномерной, то есть частота плавучести N (z) может изменяться по глубине по произвольному закону, что соответствует естественной среде, где профиль стратификации может иметь сложную структуру и может меняться, например, в зависимости от сезона. В начальный момент времени
(2.4)
Функция f задает источник:
Заметим, что . В случае равномерного движения источника возмущения его положение определяется функцией r (t) = x0 + Ut, где U – скорость прямолинейного движения массового источника.
Можно выписать задачи для остальных компонент скорости [18], а именно, для компоненты vx:
(2.5)
с граничными условиями
(2.6)
(2.7)
начальным условием
, (2.8)
и для компоненты vy:
(2.9)
(2.10)
(2.11)
(2.12)
Для нахождения полей компонент скорости vx (x, t) и vy (x, t) во всей области Ω нужно решить задачу (2.5)–(2.8) и задачу (2.9)–(2.12), при этом в правых частях условий (2.6) и (2.10) используется функция vz, найденная ранее в результате решения задачи (2.1)–(2.4).
Моделирование поля скоростей на свободной поверхности
Пусть требуется найти горизонтальные компоненты поля скоростей жидкости vx, vy не во всей области Ω, а только на свободной поверхности z = H. Тогда сначала следует решить задачу (2.1)–(2.4), то есть найти вертикальную компоненту скорости vz (x, t), на всем временном отрезке [0, T], где T – полное время моделирования. Далее, найденное поле вертикальной компоненты скорости подставляется в правую часть соотношений (2.6) и (2.10) и для каждой точки свободной поверхности (xm, yn, H) решаются соответствующие задачи Коши для обыкновенных дифференциальных уравнений второго порядка (2.6) и (2.10), в результате чего определяются горизонтальные компоненты скорости vx (x, y, H, t) и vy (x, y, H, t), в точках свободной поверхности на всем интервале времени t ∈ [0, T].
3.1. Численная схема для нахождения вертикальной компоненты скорости жидкости vz в области Ω. Для решения краевой задачи (2.1)–(2.4) в области Ω используется неявная разностная схема на равномерной расчетной сетке, состоящей из Nx × Ny × Nz узлов. Обозначим
где m = 1, ..., Nx, n = 1, ..., Ny, k = 1, ..., Nz, s = 2, ..., T / q, h – шаг по пространству, q – шаг по времени. На каждом шаге по времени s решается система линейных алгебраических уравнений (СЛАУ)
(3.1)
для неизвестных , m = 1, ..., Nx, n = 1, ..., Ny, k = 1, ..., Nz. Для апроксимации вторых производных по времени и по пространственным переменным в (2.1) используются центральные разности:
Здесь Δxyz, Δxy – соответствующие разностные операторы. Таким образом, внутри области Ω, то есть для 1 < m < Nx, 1 < n < Ny, 1 < k < Nz имеем:
(3.2)
На свободной поверхности, то есть для k = Nz граничное условие (2.2) апроксимируется следующим образом:
, (3.3)
где обозначает апроксимацию центральной разностью второй производной по времени:
На остальной части границы области , то есть для m = 1, m = Nx, n = 1, n = Ny, k = 1, из условия (2.3) имеем:
(3.4)
Коэффициенты матрицы Q системы (3.1) определяются левыми частями соотношений (3.2)–(3.4), а вектор правой части rs – правыми частями этих соотношений.
3.2. Расчет горизонтальных компонент скорости vx и vy на свободной поверхности. Опишем подробнее процедуру нахождения горизонтальной x-компоненты скорости жидкости на свободной поверхности, vx (x, y, H, t), (x, y) ∈ [xmin, xmax] × [ymin, ymax], t ∈ [0, T]. Выпишем условие (2.6) в точке свободной поверхности (xm, yn), где x = xmin + h (m – 1), y = ymin + h (n – 1), m = 1, ..., Nx, n = 1, ..., Ny:
, (3.5)
где обозначено ζ (t) = (xm, yn, H, t). Зададим начальное условие
(3.6)
Функция vz правой части (3.5) найдена ранее в результате решения задачи (2.1)–(2.4). Для апроксимации второй производной по времени используется схема центральных разностей:
,
где ζs ≡ ζ (sq), ζ0 = ζ1 = 0. Решив задачу Коши (3.5), (3.6) для всех точек (xm, yn), m = 1, ..., Nx, n = 1, ..., Ny найдем поле vx на свободной поверхности. Аналогично (с использованием условия (2.9)) расчитывается поле горизонтальной компоненты скорости на свободной поверхности vy.
3.3. Компьютерная реализация. Описанный выше алгоритм был реализован в компьютерной программе на языке программирования С++. Разработанная программа позволяет моделировать распространение гравитационных волн от движущегося точечного источника возмущений, рассчитывать изменение во времени поля вертикальной компоненты скорости жидкости vz во всей области Ω, изменение во времени горизонтальных компонент скорости жидкости vx и vy на свободной поверхности жидкости. При решении краевой задачи (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),
,
где ψs + 1, j – j-е приближение к ψs + 1. Таким образом, на каждом шаге выполняется 5–10 итераций, в результате чего величина невязки e (s, j) падает на 1–2 порядка до значений около 10–2. Весь такой расчет (до выхода источника возмущений за пределы области Ω) занимает около суток. В программе реализовано периодическое сохранение текущих результатов, что позволяет продолжить счет в случае сбоя или когда требуемое время расчета превышает максимальное допустимое на суперкомпьютере время непрерывного счета. Таким образом находится vz во всей области Ω для t ∈ [0, T]. Далее решается набор задач (3.5), (3.6), и в результате получается искомое поле течений (vx, vy) на свободной поверхности z = H при t ∈ [0, T]. Для реализации элементарных операций с разрежеными матрицами и решения системы (3.1) используется пакет GNU Scientific Library (GSL) [21]. Исходный код программы доступен в интернете по адресу https://bitbucket.org/Jclash/inwaves.
Сравнение с асимптотиками
Для верификации правильности расчетов было проведено сравнение результатов расчета копмьютерной программы с известными ассимптотическими решениями. Пусть жидкость равномерно стратифицирована (N = const), а гравитационные волны возбуждаются точечным массовым источником, т.е. функция f правой части уравнений (2.1), (2.5), (2.9) задана как
Источник возмущения движется равномерного и прямолинено со скоростью V в неограниченной области, то есть его движение описывается функцией r (t) = x0 + Ut, где |U| = V. В этом случае известны ассимтотики дальнего поля для скорости жидкости v, когда источник движется горизонтально [13] или под фиксированным углом к горизонту [14].
4.1. Горизонтальное движение. В случае, когда источник движется горизонтально прямолинейно и с постоянной скоростью, скорость жидкости в дальнем поле r1 >> V / N для x > 0 может быть вычислена с использованием следующей формулы [13] в сферической системе координат (r, ξ, φ) с центром в точке текущего положения точечного источника r (t) (рис. 2):
(4.1)
Рис. 2. Схема расчета асимптотики поля скорости при горизонтальном движении массового источника из [13]
Здесь
,
x1 = x – r, r1 = (x1, y1, z1) – вектор из точки текущего положения источника O1 в точку M с координатами (x, y, z) в которой расчитывается поле скоростей, , ξ = arccosx1 / r1, , H (x) – функция Хевисайда,
Пример сравнения результата численного моделирования и асимптотики (4.1) показан на рис. 3.
Рис. 3. Моделирование горизонтального движения в идеальной равномерно стратифицированной жидкости, N = 0.8. Скорость движения V = 1 м/с, вертикальная составляющая скорости жидкости vz (x, 1, z) показана в вертикальном сечении y = 1. Время движения T = 125 с для (a). Показано численное решение (а) и аналитическая аппроксимация (б)
4.2. Движение под фиксированным углом к горизонту. Аналогичная асимптотика может быть выписана в случае, когда источник возмущения движется равномерно под фиксированным углом γ к горизонту, то есть U = (Vcosγ, 0, Vsinγ). Асимптотика дальнего поля для скорости в точке M в системе координат (r1, θ, φ) с центром в точке текущего положения источника O1 (рис. 4) имеет следущий вид [14]:
, (4.2)
где
,
,
, θ = arccos (x1 / r1), , α = 90° – γ, и тригонометрическии функции углов β, ω и η находятся из следующих соотношений:
(4.3)
Рис. 4. Схема расчета асимптотики поля скорости при движении массового источника под углом к горизонту из [14]
Поле скоростей (4.2) вычисляется как сумма волн, соответствующих действительным корням кубического уравнения (4.3), при которых выполнено двустороннее неравенство 0 < τs < T, где T – полное время движения,
Для вычисления асимптотик (4.1) и (4.2) была написана компьютерная программа на языке программирования Python. Было проведено сравнение между численными решениями задачи (2.1)–(2.4) в области Ω и приближениями (4.1) и (4.2). Пример такого сравнения показан на рис. 3 и рис. 5. Установлено, что результаты расчетов программы согласуются с этими асимптотиками на достаточном расстоянии от источника и линии движения, когда после начала движения прошло достаточно продолжительное время (то есть, когда картина возмущений, формирующихся вокруг движущегося источника, принимает установившийся характер).
Рис. 5. Моделирование распространения внутренних волн от массового источника, который движется равномерно и прямолинейно под фиксированным углом γ = 30° к горизонту, N = 0.8. Скорость движения V = 1 м/с, вертикальная составляющая скорости жидкости vz (x, 1, z) показана в вертикальном сечении y = 1. Время движения T = 139.3 с. Показаны численное решение (а) и аналитическая аппроксимация (б)
Пример расчета
Пусть область Ω = [–100, 100] × [–50, 50] × [–25, 25] с размерами 200.0 × 100.0 × 50.0 заполнена однородно стратифицированной (частота плавучести N = 0.8) идеальной жидкостью, глубина 2H = 50.0 м. Промоделируем распространение внутренних гравитационных волн в области Ω, в качестве источника возмущения возьмем горизонтально движущий точечный массовый источник. Такая вытянутая вдоль линии движения форма области выбрана для того, чтобы за время движения T картина течений на поверхности (движущаяся вслед за источником возмущения) успела установиться. Ширина области ymax – ymin выбрана достаточной для того, чтобы края области (на которых к тому же задано нефизическое условие vz = 0) не влияли на картину возмущения поверхности. Глубина области 2H выбрана достаточной, чтобы отраженное от дна (где задано условие непротекания vz = 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П. Результаты расчета вертикальной компоненты скорости vz (x, y, z, t) приведены на рис. 6. Для момента времени t = 175 с показаны (а) величина вертикальной компоненты vz (x, y, H) на свободной поверхности z = H и (б) величина вертикальной компоненты vz (x, 0, z) в вертикальном сечении y = 0, проходящем через линию движения источника возмущения.
Рис. 6. Величина вертикальной компоненты скорости vz на свободной поверхности (а) и в толще жидкости (б) в момент времени t = 175 с
Результаты решения задачи (2.1)–(2.4) загружались с вычислительного кластера на персональный компьютер, на котором уже решалась существенно менее требовательная к вычислительным ресурсам задача (3.5), (3.6). Таким образом определялись компоненты скорости vx и vy на свободной поверхности z = H. Эти компоненты показаны на рис. 7, а и б. Для исключения влияния переходных процессов на картину течения, в расчетах взаимодействия течения с ветровым волнением следует брать поля скоростей vx, vy на достаточном удалении от точки, находящейся над точкой, из которой начинается движения источника. Например, на рис. 8 показаны поля скоростей не на всей верхней границе области Ω, а только на учаcтке свободной поверхности [25, 125] × [–50, 50].
Рис. 7. Поле течений (vx, vy) на свободной поверхности 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
Россия, Москва
Список литературы
- Нестеров С.В., Шамаев А.С., Шамаев С.И. Алгоритмы, методы и средства компьютерной радиотомографии приповерхностного слоя Земли. М.: Научный мир, 1996. 296 с.
- Ulaby F.T., Long D.G. Microwave Radar and Radiometric Remote Sensing. Artech, 2015. 1116 p.
- Jackson C.R., da Silva J.C.B., Jeans G. et al. Nonlinear internal waves in synthetic aperture radar imagery // Oceanography. 2013. V. 26. № 2. P. 68–79.
- Baydulov V.G., Knyazkov D., Shamaev A.S. Motion of mass source in stratified fluid // J. Phys.: Conf. Ser. V. 2224. 2021 2nd Int. Symp. on Automation, Information and Computing (ISAIC 2021) December 03 – 06 2021 Online. P. 012038-1–8. 2022. https://doi.org/10.1088/1742-6596/2224/1/012038
- Байдулов В.Г. О решении обратной задачи движения источника в стратифицированной жидкости // Волны и вихри в сложных средах: 12-ая межд. конф. – школа молодых ученых; 01 – 03 декабря 2021 г. Сб. матер. школы. М.: ООО ИСПОпринт, 2021. С. 31–35.
- Ulaby F.T., Moore R.K., Fung A.K. Microwave Remote Sensing. Active and Passive. Massachusetts: Addison-Wesley Publishing Company, 1981. 456 p.
- Knyazkov D. Diffraction of Plane Wave at 3-dimensional Periodic Layer // AIP Conference Proceedings. 2018. V. 1978. P. 470075-1–4. https://doi.org/10.1063/1.5044145
- Булатов М.Г., Кравцов Ю.А., Лаврова О.Ю. и др. Физические механизмы формирования аэрокосмических радиолокационных изображений океана // УФН. 2003. Т. 173. № 1. С. 69–87.
- Князьков Д.Ю., Байдулов В.Г., Савин А.С., Шамаев А.С. Прямые и обратные задачи динамики поверхностного волнения, вызванного обтеканием подводного препятствия // ПММ. 2023. Т. 87. Вып. 3. С. 442–453. https://doi.org/10.31857/S0032823523030074
- Гавриков А.А., Князьков Д.Ю., Романова А.В. и др. Моделирование влияния волнения поверхности на спектр собственного излучения океана // Программные системы: теория и приложения. 2016. Т. 7. Вып. 2(29). С. 73–84.
- Knyazkov D., Shamaev A. Rectilinear motion of mass source in non-uniformly stratified fluid. AIP Conf. Proc. 2024. V. 3094(1). P. 500028-1–4. https://doi.org/10.1063/5.0210166
- Горелов А.М., Носов В.Н., Савин А.С., Савина Е.О. Метод расчета поверхностных возмущений над точечным источником и диполем // Изв. РАН. МЖГ. 2009. № 1. С. 203–207.
- Voisin B. Internal wave generation in uniformly stratified fluids. Part 2. Moving point sources // J. Fluid Mech. 1994. V. 261. P. 333–374.
- Scase M.M., Dalziel S.B. Internal wave fields and drag generated by a translating body in a stratified fluid // J. Fluid Mech. 2004. V. 498. P. 289–313.
- Матюшин П.В. Процесс формирования внутренних волн, инициированный началом движения тела в стратифицированной вязкой жидкости // Изв. РАН. МЖГ. 2019. № 3. С. 83–97.
- Bulatov V.V. Mathematical modeling of dynamics of internal gravity waves in the ocean with arbitrary distribution of buoyancy frequency // Fluid Dyn. 2023. V. 58 (Suppl 2). P. 274–285. https://doi.org/10.1134/S0015462823603169
- Зарубин Н.А., Шамаев А.С. Исследование взаимодействия поверхностных ветровых волн с течением // Морские интеллект. технол. 2023. Т. 3. № 4. С. 93–99.
- Bulatov V.V., Vladimirov Yu.V. Wave Dynamics of Stratified Mediums. M.: Наука, 2012. 584 p.
- Булатов В.В., Владимиров Ю.В. Дальние поля внутренних гравитационных волн от движущихся источников возмущений // Вестн. МГТУ им. Н.Э. Баумана. Сер. Естеств. науки. 2018. № 4. С. 73–89. https://doi.org/10.18698/1812-3368-2018-4-73-89
- Saad Y., Schultz M.H. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems // SIAM J. on Sci.&Statist. Comput. 1986. V. 7:3. P. 856–869.
- Galassi M., Davies J., Theiler J. et al. GNU Scientific Library Reference Manual (3rd Ed.). Network Theory Ltd, 2009. 592 p.
- Чашечкин Ю.Д., Гуменник Е.В., Сысоева Е.Я. Трансформация плотностного поля трехмерным телом, движущимся в непрерывно стратифицированной жидкости // ПМТФ. 1995. № 1. С. 20–32.
Дополнительные файлы
