Simulation of the Flow Velocity Field on the Free Surface of a Stratified Fluid
- Авторлар: Knyazkov D.Y.1
-
Мекемелер:
- Ishlinsky Institute for Problems in Mechanics of the RAS
- Шығарылым: Том 88, № 5 (2024)
- Беттер: 745-757
- Бөлім: Articles
- URL: https://journals.rcsi.science/0032-8235/article/view/280966
- DOI: https://doi.org/10.31857/S0032823524050074
- EDN: https://elibrary.ru/JPKACG
- ID: 280966
Дәйексөз келтіру
Толық мәтін
Аннотация
The paper considers the problem of simulation of the velocity field on the free surface of an ideal stratified fluid generated by internal gravitational waves that reached the surface. The buoyancy frequency may vary with depth. The computer program has been written that allows calculating all components of the velocity field on the surface. It is shown that the calculation results for the vertical velocity component are consistent with the known asymptotics obtained in the far-field approximation for the cases of uniform and rectilinear motion of a point mass source horizontally (by B. Voisin) or at a fixed angle to the horizon (by M.M. Scase and S.B. Dalziel) in a uniformly stratified fluid.
Негізгі сөздер
Толық мәтін
Введение
Внутренние гравитационные волны существуют в толще океана из-за неравномерности плотности воды, вызванной неоднородностями солености и температуры, загрязнением. Они могут быть вызваны взаимодействием подводных течений с неровностями дна, вершинами подводных хребтов, подводными кабелями или трубопроводами, а также перемещением крупных рыб или косяков промысловой рыбы. Моделирование гравитационных волн важно при изучении океана, наблюдениях за загрязнениями, при прогнозировании погоды [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/.
Авторлар туралы
D. Knyazkov
Ishlinsky Institute for Problems in Mechanics of the RAS
Хат алмасуға жауапты Автор.
Email: knyaz@ipmnet.ru
Ресей, Moscow
Әдебиет тізімі
- Nesterov S.V., Shamaev A.S., Shamaev S.I. Methods, Algorithms and Tools of Aerospace Computer Tomography of Near-Surface Layer of the Earth. Moscow: Nauch. Mir, 1996. 272 p. (in Russian)
- 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 // Oceanogr., 2013, vol. 26, no. 2, pp. 68–79.
- Baydulov V.G., Knyazkov D., Shamaev A.S. Motion of mass source in stratified fluid // J. Phys.: Conf. Ser., 2021, vol. 2224, 2nd Int. Symp. on Automation, Information and Computing (ISAIC 2021), December 03–06 2021 Online. pp. 012038-1–8. 2022. https://doi.org/10.1088/1742-6596/2224/1/012038
- Baydulov V.G. On the solution of the inverse problem of the motion of a source in a stratified fluid // in: Proc. 12th Int. Conf. – School of Young Sci. Waves and Vortices in Complex Media. Moscow, Dec. 01–03, 2021. Moscow: ISPOPrint, 2021. pp. 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 Conf. Proc., 2018, vol. 1978, 470075-1–4. https://doi.org/10.1063/1.5044145
- Bulatov M.G., Kravtsov Yu.A. Lavrova O.Yu. et al. Physical mechanisms of aerospace radar imaging of the ocean // Phys. Usp., 2003, vol. 46, no. 1, pp. 63–79.
- Knyazkov D.Y., Baydulov V.G., Savin A.S., Shamaev A.S. Direct and inverse problems of the dynamics of surface waves caused by the flow around an underwater obstacle // Fluid Dyn., 2023, vol. 58, pp. 1725–1733. https://doi.org/10.1134/S0015462823603030
- Gavrikov A., Knyazkov D., Romanova A., Chernik V., Shamaev A. Simulation of influence of the surface disturbance on the ocean self-radiation spectrum // Progr. Syst.: Theory&Appl., 2016, vol. 7, iss. 2, pp. 73–84. https://doi.org/10.25209/2079-3316-2016-7-2-73-84
- Knyazkov D., Shamaev A. Rectilinear motion of mass source in non-uniformly stratified fluid. AIP Conf. Proc., 2024, vol. 3094(1), pp. 500028-1–4. https://doi.org/10.1063/5.0210166
- Gorelov, A.M., Nosov, V.N., Savin, A.S. et al. Method of calculating surface disturbances over a point source and a dipole // Fluid Dyn., 2009, vol. 44, no. 1, pp. 170–174. https://doi.org/10.1134/S0015462809010177
- Voisin B. Internal wave generation in uniformly stratified fluids. Part 2. Moving point sources // J. Fluid Mech., 1994, vol. 261, pp. 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, vol. 498, pp. 289–313.
- Matyushin P.V. Process of the formation of internal waves initiated by the start of motion of a body in a stratified viscous fluid // Fluid Dyn., 2019, vol. 54, no. 3, pp. 374–388. https://doi.org/10.1134/S0015462819020095
- Bulatov V.V. Mathematical modeling of dynamics of internal gravity waves in the ocean with arbitrary distribution of buoyancy frequency // Fluid Dyn., 2023, vol. 58 (Suppl 2), pp. 274–285. https://doi.org/https://doi.org/10.1134/S0015462823603169
- Zarubin N.A., Shamaev A.S. Investigation of the model of interaction of wind waves with the sea current // Marine Intellect. Technol., 2023, vol. 62, pp. 93–98. (in Russian) https://doi.org/10.37220/MIТ. 2023.62.4.070
- Bulatov V.V., Vladimirov Yu.V. Wave Dynamics of Stratified Mediums. Moscow: Nauka, 2012. 584 p.
- Bulatov V.V., Vladimirov Yu.V. Far fields of internal gravitational waves from moving perturbance sources // Herald of the Bauman Moscow State Tech. Univ., Nat. Sci., 2018, no. 4, pp. 73–89. (in Russian) 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, vol. 7:3, pp. 856–869.
- Galassi M., Davies J., Theiler J. et al. GNU Scientific Library Reference Manual (3rd Ed.). Network Theory Ltd, 2009. 592 p.
- Chashechkin Y., Gumennik E., Sysoeva E. Transformation of a density field by a three-dimensional body moving in a continuously stratified fluid // J. Appl. Mech.&Tech. Phys., 1995, vol. 36, pp. 19–29.
Қосымша файлдар
