Журнал вычислительной математики и математической физики, 2021, T. 61, № 4, стр. 572-579

ВВЕДЕНИЕ

Важным механизмом возбуждения внутренних гравитационных волн (ВГВ) в природных (океан, атмосфера Земли) и искусственных стратифицированных средах является их генерация источниками возмущений различной физической природы: естественного (движущиеся возмущения атмосферного давления, обтекание неровностей рельефа океана, подветренные горы) и антропогенного (морские технологические конструкции, схлопывание области турбулентного перемешивания, подводные взрывы) характеров (см. [1]–[5]). В частности, касательное напряжение ветра, создаваемое движущимся ураганом, может формировать на океанической поверхности структуру в виде движущейся воронки с почти радиальной симметрией. Как показывают результаты мониторинга Мирового океана, движущиеся возмущения морской поверхности являются одним из основных природных механизмов генерации интенсивных ВГВ (см. [4], [5]). Аналитические результаты решений задач о генерации ВГВ представляются в самой общей интегральной форме, и в этом случае полученные интегральные представления требуют разработки численных и асимптотических методов их исследования. При математическом моделировании генерации ВГВ, возбуждаемых нелокальными источниками возмущений, наиболее распространенными являются два способа (см. [6], [7]). Первый заключается в численном решении системы уравнений гидродинамики, описывающей ВГВ, к недостаткам которого следует отнести ограниченность области пространства, в котором возможно численное решение задачи (см. [8]–[10]). При изучении дальнего распространения ВГВ прямые численные расчеты нецелесообразны, так как в дальней зоне волновые поля относительно малы по амплитуде и обычно их можно описать посредством линейных уравнений. Кроме того, эффекты вязкости, вращения среды и ее сжимаемости пренебрежимо малы и не сказываются на дальнем распространении ВГВ (см. [1]–[3], [6], [7]). Поэтому волновое поле в дальней зоне можно описать сравнительно простыми аналитическими формулами. Распространение диспергирующих ВГВ в стратифицированных средах и создаваемые ими волновые картины на больших расстояниях от источников возмущений (много больших его характерных размеров) практически не зависят от их формы и определяются только законом дисперсии и скоростью источника. Поэтому второй способ состоит в том, чтобы заменить функцию, описывающую форму нелокального источника, функцией, имеющей достаточно простое аналитическое представление (см. [11]–[14]). Целью настоящей работы является построение аналитических решений, описывающих поля ВГВ, возбуждаемых движущимся нелокальным модельным источником возмущений в слое стратифицированной среды конечной толщины.

1. ПОСТАНОВКА ЗАДАЧИ И ИНТЕГРАЛЬНЫЕ ФОРМЫ РЕШЕНИЙ

Рассматривается задача о полях ВГВ, возникающих при движении нелокального источника возмущений в слое невязкой стратифицированной среды толщины $H$. Источник движется на поверхности среды $z = 0$ с постоянной скоростью $V$ в горизонтальном направлении оси $x$, ось $z$ направлена вверх, форма нелокального источника описывается функцией $s(x,y)$. Рассматривается установившийся режим волновых колебаний. В движущейся системе координат и в приближении Буссинеска имеем следующее уравнение для малых возмущений вертикальной компоненты скорости $W(x,y,z)$ (см. [1], [6], [7]):

где ${{N}^{2}}(z)$ – квадрат частоты Брента–Вяйсяля (частоты плавучести), которая далее предполагается постоянной, N²(z) = N² = const, ${{\rho }_{0}}(z)$ – невозмущенная плотность среды по глубине, $g$ – ускорение свободного падения. Линеаризованное граничное условие на поверхности имеет вид (см. [2], [3], [6], [7])

На дне используется условие непротекания

В безразмерных координатах и переменных $x{\kern 1pt} * = \pi x{\text{/}}H$, $y{\kern 1pt} * = \pi y{\text{/}}H$, $z{\kern 1pt} * = \pi z{\text{/}}H$, $M = V{\text{/}}c$, $W{\kern 1pt} * = W{\text{/}}c$, $c = NH{\text{/}}\pi $ уравнение (1.1) и граничные условия (1.2), (1.3) перепишутся следующим образом (индекс “*” далее опускается):

Относительно функции $s(x,y)$ можно сделать следующие предположения: эта функция обладает цилиндрической симметрией $s(x,y) = s(\rho )$, $\rho = \sqrt {{{x}^{2}} + {{y}^{2}}} $ и $\int_{ - \infty }^\infty \rho s(\rho )d\rho < \infty $. Далее задача будет решаться методом Фурье. Для аналитического описания формы нелокального источника необходимо выбрать такую функцию $s(\rho )$ и двойное преобразование Фурье (с точностью до множителя $2\pi $, совпадающее с преобразованием Фурье–Бесселя), котороe имеет особые точки только в виде простых полюсов. В этом случае появляется возможность построить и исследовать аналитические представления решений для функции $W(x,y,z)$. Более сложная структура особых точек Фурье образа функции $s(\rho )$, например, наличие точек ветвления, не позволяет решить задачу аналитически. Рассмотрим семейство функций ${{g}_{m}}(r,a) = {{(r{\text{/}}2a)}^{m}}{{K}_{m}}(ar){\text{/}}\Gamma (m + 1)$, где $\mathop K\nolimits_m (\tau )$ – функция порядка МакДональда, $\Gamma (\tau )$ – гамма-функция Эйлера. Функции $\mathop g\nolimits_m (r,a)$ имеют непрерывные производные до порядка $2m - 1$ включительно. Преобразование Фурье–Бесселя для функции ${{g}_{m}}(r,a)$ имеет вид

где ${{J}_{0}}(kr)$ – функция Бесселя нулевого порядка (см. [15], [16]). Так как искомое решение $W(x,y,z)$ должно иметь непрерывные вторые производные, то далее в качестве функции, описывающей форму нелокального источника, будем использовать следующее выражение: $s(\rho ) = \mathop g\nolimits_2 (\rho ,a) = - {{s}_{0}}\mathop \rho \nolimits^2 \mathop a\nolimits^2 {{K}_{2}}(a\rho ){\text{/}}2$, где параметр ${{s}_{0}}$ – максимальная высота (вертикальный масштаб) нелокального источника ($s(0) = - {{s}_{0}}$), a – характерная ширина (горизонтальный масштаб) этого источника. С ростом значения параметра a характерная ширина источника уменьшается. Выбор этой модельной функции позволяет описать генерацию ВГВ движущейся областью возмущения атмосферного давления, которая, как правило, обладает достаточно выраженной радиальной симметрией (см. [5], [9], [11]).

Решение (1.4)–(1.6) будем искать в виде интеграла Фурье

Тогда для определения функции $\varphi (\mu ,\nu ,z)$ необходимо решить краевую задачу

Здесь $\mathop k\nolimits^2 = \mathop \mu \nolimits^2 + \mathop \nu \nolimits^2 $, $S(k)$ – преобразование Фурье–Бесселя от функции $s(\rho )$:

Решение (1.8)–(1.10) имеет вид

2. ПОСТРОЕНИЕ АНАЛИТИЧЕСКИХ РЕШЕНИЙ

Далее подставим (1.11) в (1.7) и сместим контур интегрирования на комплексной плоскости переменной $\mu $ вверх на величину $\varepsilon > 0$. Контур интегрирования необходимо сдвинуть для того, чтобы удовлетворить условию излучения, т.е. отсутствию волн уходящих вперед от источника возмущений. Разложение функции $\varphi (\mu ,\nu ,z)$, определяемой выражением (1.11), в ряд по собственным функциям краевой задачи (1.8) с нулевыми граничными условиями при $z = 0, - \pi $ имеет вид

Ряд (2.1) сходится к $\varphi (\mu ,\nu ,z)$ поточечно (по переменной $z$), но неравномерно, так как не выполняется граничное условие (1.9) при $z = 0$. По $m$ ряд (2.1) сходится условно и медленно, так как общий член ряда имеет порядок $1{\text{/}}m$. Далее для исследования функции $\phi (\mu ,\nu ,z)$, определяемой выражением (1.11), используем теорему Миттаг–Леффлера о разложении мероморфной (по переменной $kr$) функции для случая простых полюсов (см. [17]): Выражение (2.2), объединив в котором второе и третье слагаемые, можно представить в виде

Первое слагаемое в (2.3) – простейшая линейная функция, удовлетворяющая граничным условиям (1.9), (1.10), второе слагаемое – взято с обратным знаком разложение в ряд по тригонометрическим функциям, третье слагаемое – это ряд (2.1). Ряд (2.3) также может быть получен из ряда (2.2) с помощью метода Крылова улучшения сходимости тригонометрических рядов (см. [18]). Этот метод широко используется в задачах о генерации полей ВГВ источниками различной физической природы в стратифицированных средах (см. [6], [7]).

Для проведения интегрирования по переменной $\mu $ в (1.7) необходимо исследовать поведение полюсов функции $\varphi (\mu ,\nu ,z)$. Эта функция имеет следующие полюса по $\mu $. Две серии полюсов, $\mathop {\mu = \pm \mu }\nolimits_m (\nu )$, $\mathop {\lambda = \pm i\lambda }\nolimits_m (\nu )$, которые определяют соответствующие дисперсионные зависимости: $\mathop \mu \nolimits_m (\nu ) = {{\left( {0.5\left( {{{B}_{m}} - {{C}_{m}}} \right)} \right)}^{{1/2}}}$, $\mathop \lambda \nolimits_m (\nu ) = {{\left( {0.5\left( {{{B}_{m}} + {{C}_{m}}} \right)} \right)}^{{1/2}}}$, где $\mathop B\nolimits_m = {{({{({{m}^{2}} - {{M}^{{ - 2}}} + {{\nu }^{2}})}^{2}} + 4{{\nu }^{2}}{{M}^{{ - 2}}})}^{{1/2}}}$, ${{C}_{m}} = {{m}^{2}} - {{M}^{{ - 2}}} + {{\nu }^{2}}$. Также существуeт пара полюсов $\mu = \pm p(\nu )$, $p(\nu ) = {{({{a}^{2}} + {{\nu }^{2}})}^{{1/2}}}$, которые зависят только от фурье-образа функции $s(\rho )$, определяющей форму нелокального источника. Вычеты в точках $\mathop {\mu = \pm \mu }\nolimits_m (\nu )$, $\mathop {\lambda = \pm i\lambda }\nolimits_m (\nu )$ одинаковы для выражений (2.1), (2.3), вычеты в точках $ \pm ip(\nu )$ – различные для этих выражений.

Проведем почленное интегрирование по переменной $\mu $ в (1.7). Контур интегрирования необходимо замкнуть вниз при $x > 0$ и вверх при $x < 0$. На фиг. 1 штриховой линией изображен контур интегрирования и качественная картина расположения полюсов для первой волновой моды $n = 1$. Все дальнейшие вычисления будут приведены для следующих значений безразмерных параметров: $M = 3.05$, ${{S}_{0}} = 0.046$, $a = 2.44$. В размерных единицах типичная глубина океана равна $H \approx {{10}^{3}}$ м, тогда вертикальный масштаб нелокального источника будет составлять десятки метров, горизонтальный масштаб будет равен нескольким сотням метров и скорость движения источника – десятки сантиметров в секунду. Поэтому использованные для численных расчетов характерные значения параметров $M,\;{{S}_{0}},\;a$ полностью соответствуют наблюдаемым в Мировом океане пространственно-временным характеристикам областей возмущения морской поверхности от движущихся тайфунов, генерирующих интенсивные ВГВ (см. [4], [11]).

Действительные полюса ${{\mu }_{m}}(\nu )$ растут с увеличением номера моды $m$ и ограничены: ${{\mu }_{1}}(\nu ) < {{\mu }_{2}}(\nu ) < ... < N{\text{/}}V$. Мнимые полюса ${{\lambda }_{m}}(\nu )$ неограниченно растут с увеличением номера моды $m$ и ${{\lambda }_{m}}(\nu ) \approx m$ при больших значениях $m$. Для данных параметров задачи значение $p(\nu )$ находится в интервале ${{\lambda }_{2}}(\nu ) < p(\nu ) < {{\lambda }_{3}}(\nu )$. При параметрах, определяющих иные вертикальные и горизонтальные масштабы нелокального источника возмущений, значения $p(\nu )$ находятся в других интервалах, и вклад от этих полюсов может заметно влиять на структуру полного волнового поля особенно в ближней зоне.

Обозначим далее

При $x > 0$ имеем При $x < 0$ имеем

Сумма вычетов в действительных полюсах $\mu = \pm {{\mu }_{m}}(\nu )$ равна

Вычеты в точках $\lambda = \pm i{{\lambda }_{m}}(\nu )$ равны

Аналитические выражения для вычетов в точках $ \pm ip(\nu )$ в силу их громоздкости не приводятся. На фиг. 2 представлены результаты расчетов главного значения интеграла (2.4) при $m = 1$, $\nu = 0.9$, $z = - \pi {\text{/}}2$ (сплошная линия), вклада действительных вычетов (штриховая линия) и вклада мнимых вычетов (штрихпунктирная линия). Из приведенных результатов следует, что при больших значениях $x$ основной вклад вносят действительные вычеты, при малых значениях $x$ – мнимые вычеты. Тогда решение $W(x,y,z)$ задачи (1.1)–(1.3) можно представить в виде суммы волновых мод:

При больших значениях $x,y$ (вдали от нелокального источника возмущений) интегралы (2.5) можно вычислить методом стационарной фазы (см. [6], [7], [14]). Решение для вертикальной компоненты скорости ВГВ можно также представить в другом виде. Возьмем обратное преобразование Фурье от (2.3). Тогда, исходя из аналитического выражения для функции $s(\rho )$ и учитывая соотношение

имеем Отсюда можно получить где слагаемые ${{w}_{m}}(x,y,z)$ вычисляются аналогично ${{W}_{m}}(x,y,z)$, за исключением вычетов в точках $ \pm ip(\nu )$. На фиг. 3, 4 изображены графики (полигоны) зависимости полного решения $W(x,y,z)$ от значений $n$. Параметр $n$ определяет общее число суммируемых волновых мод в выражениях (2.5), (2.6). Численные расчеты проводились при значениях $x = 1$, $y = 0.5$, $z = - 0.55$ для двух способов представления решения, определяемых соответственно выражениями (2.5) (фиг. 3) и (2.6) (фиг. 4). Из представленных результатов следует, в частности, что точность порядка 1% достигается при $n = 4$ для второго метода (фиг. 4), в то время как эта же точность для первого метода достигается при $n = 40$ (фиг. 3). На фиг. 5 приведены зависимости $W(x,y,z)$ от вертикальной координаты $z$ при $x = 1$, $y = 0.5$: сплошная линия – расчеты по формулам (2.5) при $n = 23$, пунктир – расчеты по формулам (2.6) при $n = 4$. На фиг. 6 изображены зависимости $W(x,y,z)$ от горизонтальной координаты $x$, рассчитанные по формулам (2.6) при $z = - 0.6$, $y = 2.2$, $n = 4$. Как показывают численные расчеты, при больших значениях $x,y$ основной вклад в полное поле $W(x,y,z)$ вносят слагаемые, отвечающие действительным полюсам, вклад от мнимых полюсов, описывающих форму нелокального источника, экспоненциально мал. Вблизи нелокального источника возмущений основной вклад в полное поле $W(x,y,z)$ вносят слагаемые, отвечающие мнимым полюсам, и в ближней зоне представление решения в виде (2.6) обеспечивает более быструю сходимость ряда к искомому решению.

ЗАКЛЮЧЕНИЕ

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