Математические методы моделирования аномалий ветровых волн на поверхности океана

Полный текст

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

К настоящему времени накоплен большой объем экспериментальных данных, посвященных наблюдению за морским волнением, в том числе за эффектами, связанными с воздействием на поверхность океана ветра, течений, а также наблюдению за эффектами взаимодействия поверхностных волн и течений. Эти экспериментальные данные получены средствами радиолокационного зондирования, исследованию собственного радиоизлучения взволнованной морской поверхности, лазерного зондирования, исследованию поверхности в видимом оптическом диапазоне [1$–$8]. Большое количество работ посвящено построению гидродинамических моделей для различных аномалий поверхности океана, в том числе с учетом неоднородной структуры в толще морской воды [1, 9$–$15].

Исследования поверхности океана средствами радиолокации. Воздействие течений приводят к ряду интересных и сложных явлений на поверхности моря. К одному из таких явлений относится возникновение зон деформации амплитуды поверхностного волнения за счет воздействия вышедших из глубины на поверхность воды течений.

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

Рис. 1. Радиолокационное изображение, полученное с помощью космического аппарата “Алмаз-1” в Северном море и батиметрическая схема участка съемки

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

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

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

Остановимся сначала на радиолокационных измерениях. Представим себе (разумеется, упрощенно) морскую поверхность в виде синусоидальной волны. Хорошо известно, что в случае падения плоской электромагнитной волны на такую поверхность при определенном соотношении между длиной падающей электромагнитной волны $\Lambda$ (или ее волновым числом ${\rm K}=2\pi /\Lambda$ ), углом падения $\alpha$ и периодом отражающей поверхности $\lambda$ (см. рис. 2) часть энергии падающей волны возвращается точно назад $–$ эффект резонансного рассеяния Брэгга$–$Вульфа. Это соотношение имеет вид

$\frac{\Lambda }{\lambda }=2\sin \alpha$ (1.1)

Рис. 2. Явление резонансного рассеивания на периодической поверхности

Оказывается (это проверено в эксперименте и с помощью численных расчетов), что сигнал, отраженный точно назад, примерно пропорционален амплитуде $h$ синусоиды, описывающей отражающую поверхность, если величина $h$ невелика (например, $h/\lambda <0.1$ ). Таким образом, по интенсивности обратного рассеяния радиоволн можно было бы судить об амплитуде поверхностных волн, длина которых связана с длиной электромагнитной волны соотношением (1.1).

Однако поверхность моря представляет собой, естественно, не одну синусоиду.

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

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

Данный результат является крайне важным для проблемы выявления аномалий поверхности, характеризующихся увеличением или, наоборот, уменьшением средней амплитуды в коротковолновой части спектра морского волнения. Во-первых, оказывается доступным простое средство определить среднюю амплитуду резонансных гармоник поверхности путем усреднения амплитуды пришедшего сигнала по достаточному количеству облученных поверхностей. Во-вторых, на основе результатов моделирования появляется возможность установить диапазон углов, для которых такая процедура корректна. Связь между усредненными амплитудами волн и отраженных сигналов была обнаружена ранее эмпирически. Считалось, что зависимость существует только лишь для скользящих углов падения. Численное моделирование показывает, что зависимость (не обязательно являющаяся прямой пропорциональностью) среднего отраженного сигнала от резонансной гармоники поверхности прослеживается вплоть до 30° (угол отсчитывается от вертикали).

На рис. 3 изображен радиолокационный кадр со стороной $L$, элемент разрешения локатора (а) и квадрат, по элементам которого проводится усреднение (l). Характерные масштабы: $L~10-20$ км, $a~2-4$ км, $\lambda ~200-400$ м.

Рис. 3. Обработка радиолокационного кадра с помощью процедуры усреднения.Сравнение различных масштабов задачи

Исследования поверхности океана средствами пассивной радиолокации. Важным и эффективным средством исследования аномалий поверхности моря является исследование собственного радиоизлучения поверхности. С помощью измерений собственного радиоизлучения можно получить интересные результаты о возмущениях (о которых было сказано выше при описании метода радиолокационного зондирования) поверхности океана на большой площади. Собственное радиоизлучение поверхности весьма чувствительно к небольшим флуктуациям поверхности. Эффект резонансной зависимости собственного радиоизлучения поверхности был впервые исследован в работе [4]. Было установлено, что на определенных углах наклона зондирования (около 30° с вертикалью) собственное инфракрасное радиоизлучение взволнованной поверхности максимально чувствительно к амплитуде миллиметровых и сантиметровых волн излучающей поверхности. Поэтому при исследовании поверхности средствами пассивной радиолокации крутые углы зондирования предпочтительнее пологих.

Радиометры космического базирования получают информацию об интенсивности радиоизлучения поверхности моря на площади в десятки тысяч квадратных километров. Описанные выше аномалии морского волнения могут быть обнаружены также во многих случаях и радиометрическими средствами. Математическое моделирование собственного радиоизлучения взволнованной поверхности осуществляется в работах [1, 6]. Интересный эксперимент по определению параметров излучающей поверхности в лабораторных условиях проделан под руководством Ю.Г. Трохимовского. Была установлена возможность определить характеристики формы волнистой поверхности воды, которая создана с помощью частично погруженной в воду капроновой сетки. Было проведено сравнение экспериментальных данных о собственном радиоизлучении волнистой поверхности воды и результатов математического моделирования данного излучения. Оно основано на использовании связи между мощностью собственного радиоизлучения в данном направлении и величиной поглощенной энергии когерентной плоской волны, падающей на поверхность в противоположном направлении [5].

Взаимодействие приповерхностных течений с поверхностными волнами. Вопросам количественного описания трансформаций ветровых поверхностных волн под воздействием возмущений, вышедших из глубины, посвящено также большое количество работ. Обратим здесь внимание на две группы исследований, связанных, во-первых, с выводом и анализом так называемого кинетического уравнения для плотности волнового действия [3] и численным моделированием на его основе и, во-вторых, посвященных выводу и анализу двумерного уравнения типа Шредингера и расчетам на его основе трансформаций поверхностных волн, удовлетворяющих так называемому “условию синхронизма” [9]. Приведем результаты эксперимента, выполненного на Черном море в 80-х гг. прошлого века, и обсудим его результаты в свете существующих теоретических моделей взаимодействия поверхностных волн с течениями. На рис. 4 представлена схема обтекания подводного препятствия, имеющего вид продолговатого хребта высотой около 10 м и длиной около 100 м, расположенного на глубине примерно 50 м.

Рис. 4. Схема обтекания подводного препятствия

На рис. 5 представлена картина трансформаций ветровой ряби волн с длиной около 4 см над препятствием, схема обтекания которого показана на рис. 4.

Рис. 5. Трансформация ветровой ряби над обтекаемыми неоднородностями дна

Эта картина получена с помощью радиолокатора, расположенного на вертолете, который совершил несколько пролетов над зоной трансформаций. Течение набегает на препятствие со скоростью около 3 м/с по направлению, указанному стрелкой. Трансформация ветровой ряби наблюдается внутри “восьмерки”. Для того, чтобы дать этим экспериментальным данным теоретическое объяснение, необходимо сначала сделать расчет поля скоростей течений, вышедших из глубины на поверхность, который может быть выполнен с помощью различных теоретических моделей.

Расчеты приводят к результатам, отраженным на рис. 6, из которого видно, что возмущение скорости локализовано строго над началом препятствия, а зона существенного возмущения скорости имеет небольшую протяженность $–$ примерно 10 м.

Рис. 6. Поле скоростей на поверхности воды непосредственно над движущимся телом

Вернемся к картине возмущения поверхностных волн, полученной с помощью радиолокатора вертолетного базирования, представленной на рис. 5. Кроме того, известна экспериментальная картина трансформации ветровых волн с групповыми скоростями $V_1=0.1$ м/с, $V_2=0.2$ м/с, $V_3=4.2$ м/с с длинами волн, соответственно ${\Lambda }_1=0.04$ м, ${\Lambda }_2=0.24$ м.

Зона трансформации ветровых волн в диапазонах ${\Lambda }_1$ и ${\Lambda }_2$ примерно соответствует “восьмерке” на рис. 5. Согласно теории резонансного рассеяния амплитуды первых двух кривых примерно пропорциональны амплитудам поверхностных волн длиной 2 и 12 см.

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

Кроме того, многочисленные экспериментальные данные говорят о том, что изменение спектра под действием течений имеют место и в длинноволновой части спектра, и для волн средней и малой длины. Ниже приводится модель, которая может претендовать на объяснение эффекта присутствия в спектре трансформаций ветрового волнения компонент различной длины, включая короткие волны (миллиметры$–$сантиметры).

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

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

2. Постановка задачи о деформации ветрового волнения под действием течений

Приведем в этом разделе для полноты изложения вывод основного уравнения, следуя [9]. Далее асимптотическими методами осуществляется вывод уравнения для приближенной модели трансформации спектра водной поверхности, ранее такой подход применялся в [17].

Будем рассматривать невязкую жидкость с безвихревым движением [16]. Введем потенциал скорости $\Phi$; из условия несжимаемости следует, что ${\nabla }^2\Phi =0$. Обозначим свободную поверхность $z=\zeta (x,y,t)$. Пусть ${\nabla }_{xy}=\left(\frac{\partial }{\partial x},\frac{\partial }{\partial y}\right)$ $–$ оператор горизонтального градиента. Кинематическое условие для рассматриваемого движения на свободной поверхности имеет вид:

$\frac{\partial \zeta }{\partial t}={\left.\frac{\partial \Phi }{\partial z}\right|}_{z=\zeta }-\left({\left.{\nabla }_{xy}\Phi \right|}_{z=\zeta },{\nabla }_{xy}\zeta \right)$

Динамическое условие на поверхности обусловлено требованием, чтобы перепад давления по обе стороны определялся только поверхностным натяжением. Учитывая известное уравнение Бернулли для движения, получим на поверхности следующее уравнение:

$\frac{p}{\rho }+g\zeta +\frac{\partial \Phi }{\partial t}+{\left.\frac{1}{2}{\left({\nabla }_{xy}\Phi \right)}^2\right|}_{z=\zeta }=0$,

где $p$ $–$ плотность воды, а $\rho$ $–$ давление. Пусть рассматривается модель жидкости бесконечной глубины, тогда, пренебрегая влиянием атмосферы, получим систему уравнений:

$\begin{array}{l}\Delta \Phi =Q\\ \frac{\partial \zeta }{\partial t}={\left.\frac{\partial \Phi }{\partial z}\right|}_{z=\zeta }-\left({\left.{\nabla }_{xy}\Phi \right|}_{z=\zeta },{\nabla }_{xy}\zeta \right)\\ g\zeta +\frac{\partial \Phi }{\partial t}+{\left.\frac{1}{2}{\left|{\nabla }_{xy}\Phi \right|}_{z=\zeta }\right|}^2+\frac{1}{2}{\left({\left.\frac{\partial \Phi }{\partial z}\right|}_{z=\zeta }\right)}^2=\mathrm{0,}\end{array}$ (2.1)

где $Q$ $–$ мощность источников течения. Представим потенциал $\Phi$ в виде суммы компоненты $F$, отвечающей за неоднородное течение, и компоненты $\Psi$, связанной с поверхностным течением. Смещение поверхности $\zeta$ также разделим на смещение $\eta$ для неоднородного и $\xi$ для течения на поверхности. Разложим компоненту $F$ в ряд Тейлора в точке $z=\eta$ :

 $F={\left.F\right|}_{z=\eta }+{\left.\frac{\partial F}{\partial z}\right|}_{z=\eta }\xi +\dots$

Учитывая, что для $F$ система (2.1) выполнена, получим, что первое уравнение системы переходит в $\Delta \Psi =0$.

Рассмотрим подробнее третье уравнение системы. Подставляя в него $\Phi =F+\Psi$, получим:

$\begin{array}{l}\frac{\partial \Psi }{\partial t}+\frac{\partial F}{\partial t}+{\left.\frac{{\partial }^{2}F}{\partial t\partial z}\right|}_{z=\eta }\xi +g\xi +g\eta +\frac{1}{2}{\left|{\nabla }_{xy}\left(\Psi +F+{\left.\frac{\partial F}{\partial z}\right|}_{z=\eta }\right)\right|}^2+\\ +\frac{1}{2}{\left(\frac{\partial \Psi }{\partial z}+\frac{\partial F}{\partial z}+{\left.\frac{{\partial }^{2}F}{\partial z^2}\right|}_{z=\eta }\xi \right)}^2=0\end{array}$

Здесь

$\begin{array}{l}\frac{1}{2}{\left|{\nabla }_{xy}\left(\Psi +F+{\left.\frac{\partial F}{\partial z}\right|}_{z=\eta }\right)\right|}^2=\frac{1}{2}\left({\left(\frac{\partial F}{\partial x}\right)}^2+{\left({\left.\frac{{\partial }^{2}F}{\partial x\partial z}\right|}_{z=\eta }\xi \right)}^2+{\left(\frac{\partial \Psi }{\partial x}\right)}^2+\right.\\ +2\frac{\partial F}{\partial x}\frac{\partial \Psi }{\partial x}+{\left.2\frac{\partial F}{\partial x}\frac{{\partial }^{2}F}{\partial x\partial z}\right|}_{z=\eta }\xi +{\left.2\frac{{\partial }^{2}F}{\partial x\partial z}\right|}_{z=\eta }\frac{\partial \Psi }{\partial x}\xi +{\left(\frac{\partial F}{\partial y}\right)}^2+{\left(\frac{\partial \Psi }{\partial y}\right)}^2+\\ +{\left({\left.\frac{{\partial }^{2}F}{\partial x\partial z}\right|}_{z=\eta }\xi \right)}^2+2\frac{\partial F}{\partial y}\frac{\partial \Psi }{\partial y}+\left.{\left.2\frac{\partial F}{\partial y}\frac{{\partial }^{2}F}{\partial y\partial z}\right|}_{z=\eta }\xi +{\left.2\frac{\partial \Psi }{\partial y}\frac{{\partial }^{2}F}{\partial y\partial z}\right|}_{z=\eta }\right),\end{array}$

будем считать амплитуду поверхностных волн малой, поэтому можно пренебречь квадратичными членами по $\Psi$ и $\xi$, тогда получим

$\frac{1}{2}{\left|{\nabla }_{xy}\left(\Psi +F+{\left.\frac{\partial F}{\partial z}\right|}_{z=\eta }\right)\right|}^2={\left.\frac{1}{2}{\left({\nabla }_{xy}F\right)}^2+\left({\nabla }_{xy}F\frac{\partial {\nabla }_{xy}F}{\partial z}\right)\xi +\left({\nabla }_{xy}F{\nabla }_{xy}\Psi \right)\right|}_{z=\eta }$

Аналогично распишем для последнего слагаемого равенство:

$\begin{array}{l}\frac{1}{2}{\left({\left.\frac{\partial F}{\partial z}\right|}_{z=\eta }+{\left.\frac{{\partial }^{2}F}{\partial z^2}\right|}_{z=\eta }\xi +\frac{\partial \Psi }{\partial z}\right)}^2=\frac{1}{2}{\left({\left.\frac{\partial F}{\partial z}\right|}_{z=\eta }\right)}^2+\frac{1}{2}{\left({\left.\frac{{\partial }^{2}F}{\partial z^2}\right|}_{z=\eta }\xi \right)}^2+\\ +\frac{1}{2}{\left({\left.\frac{\partial \Psi }{\partial z}\right|}_{z=\eta }\right)}^2+{\left.\left(\frac{\partial F}{\partial z}\frac{\partial \Psi }{\partial z}\right)\right|}_{z=\eta }+{\left.\frac{\partial F}{\partial z}\frac{{\partial }^{2}F}{\partial z^2}\right|}_{z=\eta }\xi +\frac{{\partial }^{2}F}{\partial z^2}\frac{\partial \Psi }{\partial z}\xi \end{array}$

Учитывая, что $F$ $–$ стационарный поток, для которого система приобретает вид (2.1), получим следующее выражение для третьего уравнения в (2.1) относительно потенциала $\Psi$:

${\left.\frac{\partial \Psi }{\partial t}\right|}_{z=\eta }+g\xi +\left({\nabla }_{xy}F\frac{\partial {\nabla }_{xy}F}{\partial z}\right)\xi +{\left.\left({\nabla }_{xy}F{\nabla }_{xy}\Psi \right)\right|}_{z=\eta }+{\left.\left(\frac{\partial F}{\partial z}\frac{\partial \Psi }{\partial z}\right)\right|}_{z=\eta }+{\left.\frac{\partial F}{\partial z}\frac{{\partial }^{2}F}{\partial z^2}\right|}_{z=\eta }\xi =0$

Аналогично, расписав второе уравнение системы (2.1) через $F$, $\Psi$, $\xi$, $\eta$, получим уравнение:

$\frac{\partial \xi }{\partial t}+{\left.\frac{{\partial }^{2}F}{\partial z^2}\right|}_{z=h}\xi +{\left.\frac{\partial \Psi }{\partial z}\right|}_{z=\eta }-\left({\nabla }_{xy}F{\nabla }_{xy}\xi \right)-\left({\nabla }_{xy}\Psi {\nabla }_{xy}\eta \right)=0$

Тогда система (2.1) переходит в систему:

$\begin{array}{l}\Delta \Psi =0\\ {\left.\frac{\partial \Psi }{\partial t}\right|}_{z=\eta }+g\xi +\left({\nabla }_{xy}F\frac{\partial {\nabla }_{xy}F}{\partial z}\right)\xi +{\left.\left({\nabla }_{xy}F{\nabla }_{xy}\Psi \right)\right|}_{z=\eta }+{\left.\left(\frac{\partial F}{\partial z}\frac{\partial \Psi }{\partial z}\right)\right|}_{z=\eta }+{\left.\frac{\partial F}{\partial z}\frac{{\partial }^{2}F}{\partial z^2}\right|}_{z=\eta }\xi =0\\ \frac{\partial \xi }{\partial t}+{\left.\frac{{\partial }^{2}F}{\partial z^2}\right|}_{z=h}\xi +{\left.\frac{\partial \Psi }{\partial z}\right|}_{z=\eta }-\left({\nabla }_{xy}F{\nabla }_{xy}\xi \right)-\left({\nabla }_{xy}\Psi {\nabla }_{xy}\eta \right)=0\end{array}$ (2.2)

Теперь предположим, что наклон поверхности жидкости мал, т.е. $\left|{\nabla }_{xy}\eta \right|\ll 1$ и, обозначив ${\nabla }_{xy}F$ через $U$, запишем систему (2.2) в виде:

$\begin{array}{l}\Delta \Psi =0\\ \left(\frac{\partial }{\partial t}+\left(U,{\nabla }_{xy}\right)\right){\left.\Psi \right|}_{z=\eta }+g\xi =0\\ \left(\frac{\partial }{\partial t}+\left(U,{\nabla }_{xy}\right)\right)\xi -{\left.\frac{\partial \Psi }{\partial z}\right|}_{z=h}=0\end{array}$ (2.3)

С помощью преобразования Фурье по $x$ и $y$ решим первое уравнение системы (2.3) с учетом краевого условия, в качестве которого выступает третье уравнение. Избавляясь от $\Psi$, подставим решение во второе уравнение и получим уравнение для поверхностного смещения $\xi (x,y,t)$:

$\left(\frac{\partial }{\partial t}+\left(U,{\nabla }_{xy}\right)\right)A\left(\frac{\partial }{\partial t}+\left(U,{\nabla }_{xy}\right)\right)\xi +g\xi =0$, (2.4)

где $A={\left(-{\Delta }_{xy}\right)}^{-1/2}$, $t>0$, $x\in \left(-\infty ,+\infty \right)$, $y\in \left(-\infty ,+\infty \right)$. Степень оператора Лапласа определяется по формуле:

${\left(-{\Delta }_{xy}\right)}^{-1/2}f(x,y)=\frac{1}{4{\pi }^2}{{\displaystyle \int }}^{\text{}}\exp \left(i{k}_{1}x\right)\exp \left(i{k}_{2}y\right)\frac{\widehat{f}\left(k_1,k_2\right)}{k}d{k}_{1}d{k}_2$,

где

$\widehat{f}\left(k_1,k_2\right)={{\displaystyle \int }}^{\text{}}\exp \left(-i{k}_{1}x\right)\exp \left(-i{k}_{2}y\right)f\left(x,y\right)dxdy$, $k=\sqrt{k_1^2+k_2^2}$

3. Асимптотический анализ решений уравнения динамики морской поверхности

Пусть вышедшее на поверхность течение имеет вид

$\overrightarrow{U}=\epsilon \overrightarrow{u}\left(x,y,t\right)+U_0{\overrightarrow{e}}_1,$

где $U_0{\overrightarrow{e}}_1$ $–$ вектор скорости основного потока, $\overrightarrow{u}\left(x,y,t\right)=\left(u_1\left(x,t\right),u_2\left(y,t\right)\right)$ $–$ вектор скорости возмущения основного потока, соразмерный по модулю со скоростью основного потока, $\epsilon$ $–$ малый параметр. Будем считать, что при $t<0$ возмущение $u$ основного потока отсутствует, далее при $t>0$ оно возникает, развиваясь от нуля до некоторого стационарного состояния $u^{*}\left(x,y\right)=\left(u_1\left(x\right),u_2\left(y\right)\right)$.

Будем искать решение задачи (2.4) в виде

$\xi ={\xi }_0+\epsilon {\xi }_1+\mathrm{...}$

Подставим $\xi$ в исходное уравнение (2.4) и приравняем коэффициенты при одинаковых степенях $\epsilon$, получим уравнения для ${\xi }_0$ и ${\xi }_1$, нулевого и первого приближений соответственно.

Нулевое приближение

Имеем

$\frac{{\partial }^2{\xi }_0}{\partial t^2}+2U_0\frac{{\partial }^2{\xi }_0}{\partial t\partial x}+U_0^2\frac{{\partial }^2{\xi }_0}{\partial x^2}=-g{A}^{-1}{\xi }_0$,     (3.1)

где $A^{-1}\xi (x,y,t)=\sqrt{-\Delta }\xi (x,y,t)=\frac{1}{4{\pi }^2}\exp {{\displaystyle \int }}^{\text{}}\left(i{k}_{1}x\right)\exp \left(i{k}_{2}y\right)\widehat{\xi }\left(k_1,k_2,t\right)kd{k}_{1}d{k}_2$.

Пусть решение уравнения (3.1) имеет вид

${\xi }_0=e^{i\Omega x}{e}^{i{\gamma }_{1}x}{e}^{i{\gamma }_{2}y}$ (3.2)

Подставим (3.2) в (3.1) и получим, что решение в виде (3.2) существует при $\Omega =-U_0{\gamma }_1\pm \sqrt[4]{g^2\left({\gamma }_1^2+{\gamma }_2^2\right)}$.

Первое приближение

Для первого приближения получаем уравнение

$\frac{{\partial }^2{\xi }_1}{\partial t^2}+2U_0\frac{{\partial }^2{\xi }_1}{\partial t\partial x}+U_0^2\frac{{\partial }^2{\xi }_1}{\partial x^2}=-g{A}^{-1}{\xi }_1+F\left(x,y,t\right)$, (3.3)

где

$\begin{array}{l}F=\Omega {\gamma }_{1}A\left(u_1{\xi }_0\right)+\Omega {\gamma }_{2}A\left(u_2{\xi }_0\right)+\Omega {\gamma }_1{u}_{1}A{\xi }_0+\\ +\Omega {\gamma }_2{u}_{2}A{\xi }_0+U_0{\gamma }_1^2{u}_{1}A{\xi }_0+U_0{\gamma }_1{\gamma }_2{u}_{2}A{\xi }_0-\\ -i{U}_0{\gamma }_{1}A\left(\frac{\partial u_1}{\partial x}{\xi }_0+u_1\frac{\partial {\xi }_0}{\partial x}\right)-i{U}_0{\gamma }_{2}A\left(\frac{\partial u_2}{\partial x}{\xi }_0+u_2\frac{\partial {\xi }_0}{\partial x}\right)\end{array}$

Пусть нам известно Фурье преобразование компонент вектора скорости, ${\widehat{u}}_1=v_1$ и ${\widehat{u}}_2=v_2$. Сделаем Фурье-преобразование от правой и левой частей уравнения (3.3) и будем искать его решение в виде ${\widehat{\xi }}_1=\widehat{\eta }(k_1,k_2)e^{i\Omega t}$.

Итак,

$\begin{array}{l}\widehat{F}=\frac{\Omega {\gamma }_1{\left(u_1{\xi }_0\right)}^{\wedge }+\Omega {\gamma }_2{\left(u_2{\xi }_0\right)}^{\wedge }}{k}+\frac{1}{\gamma }\left(\Omega {\gamma }_1{v}_1*\delta \left(k_1-{\gamma }_1,k_2-{\gamma }_2\right)+\right.\\ \left.+\Omega {\gamma }_2{v}_2*\delta \left(k_1-{\gamma }_1,k_2-{\gamma }_2\right)\right)e^{i\Omega t}+\frac{1}{\gamma }\left(U_0{\gamma }_1^2{v}_1*\delta \left(k_1-{\gamma }_1,k_2-{\gamma }_2\right)+\right.\\ \left.+U_0{\gamma }_1{\gamma }_2{v}_2*\delta \left(k_1-{\gamma }_1,k_2-{\gamma }_2\right)\right)e^{i\Omega t}-i{U}_0{\gamma }_1\frac{i{k}_1{v}_1*{\widehat{\xi }}_0+i{\gamma }_1{v}_1*{\widehat{\xi }}_0}{4{\pi }^{2}k}-\\ -U_0{\gamma }_2\frac{i{k}_1{v}_2*{\widehat{\xi }}_0+i{\gamma }_1{v}_2*{\widehat{\xi }}_0}{4{\pi }^{2}k},\end{array}$

где символом $*$ $–$ обозначен оператор свертки двух функций.

Обозначим ${\tilde{v}}_{\mathrm{1,2}}=v_{\mathrm{1,2}}(k_1-{\gamma }_1,k_2-{\gamma }_2)$. Следовательно,

$\begin{array}{l}{e}^{-i\Omega t}{\left(A^{-1}F\right)}^{\wedge }=e^{-i\widehat{\Omega }t}\cdot \widehat{F}\cdot \left|\overrightarrow{k}\right|=\Omega {\gamma }_1{\tilde{v}}_1+\Omega {\gamma }_2{\tilde{v}}_2+\\ +\frac{k}{\gamma }(\Omega {\gamma }_1{\tilde{v}}_1+\Omega {\gamma }_2{\tilde{v}}_2+U_0{\gamma }_1^2{\tilde{v}}_1+U_0{\gamma }_1{\gamma }_2{\tilde{v}}_2)+U_0{\gamma }_1\left(k_1-{\gamma }_1\right){\tilde{v}}_1+U_0{\gamma }_1^2{\tilde{v}}_1+\\ +U_0{\gamma }_2\left(k_1-{\gamma }_1\right){\tilde{v}}_2+U_0{\gamma }_1{\gamma }_2{\tilde{v}}_2=\\ =\left(\Omega {\gamma }_1\left(1+\frac{k}{\gamma }\right)+U_0{\gamma }_1\left(k_1+{\gamma }_1\right)\right){\tilde{v}}_1+\left(\Omega {\gamma }_2\left(1+\frac{k}{\gamma }\right)+U_0{\gamma }_2\left(k_1+{\gamma }_1\right)\right){\tilde{v}}_2\end{array}$

Таким образом, получаем Фурье-преобразование уравнения (3.3):

$\frac{\widehat{{\partial }^2{\xi }_1}}{\partial t^2}+2U_0+U_0^2\frac{\widehat{{\partial }^2{\xi }_1}}{\partial x^2}=-gk{\widehat{\xi }}_1+{\left(A^{-1}F\right)}^{\wedge }$ (3.4)

Далее, подставляем ${\widehat{\xi }}_1=\widehat{\eta }(k_1,k_2)e^{i\Omega t}$ в уравнение (3.4) и находим

$\widehat{\eta }=\frac{1}{k-{\left(\Omega +k_1{U}_0\right)}^2}\left.\left(\Omega {\gamma }_1\left(1+\frac{k}{\gamma }\right)+U_0{\gamma }_1\left(k_1+{\gamma }_1\right)\right){\tilde{v}}_1\right)+\left.\left(\Omega {\gamma }_2\left(1+\frac{k}{\gamma }\right)+U_0{\gamma }_2\left(k_1+{\gamma }_1\right)\right){\tilde{v}}_2\right)$ (3.5)

4. Явное выражение решения задачи о трансформации спектра взволнованной поверхности под действием течений

Общее решение задачи о трансформации спектра принимает вид:

$\widehat{\xi }={\widehat{\xi }}_0+\epsilon {\widehat{\xi }}_1=4{\pi }^2\delta \left(k_1-{\gamma }_1,k_2-{\gamma }_2\right)e^{i\Omega t}+\epsilon \hat{\eta }{e}^{i\Omega t}$ (4.1)

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

При этом знаменатель в формуле (3.5) обращается в ноль на некоторой линии в плоскости $k_1,k_2$. Для этих значений формула (3.5) теряет смысл. Связано это с тем, что в принятой линейной модели данные значения $k_1$ и $k_2$ являются резонансными значениями, при которых амплитуда соответствующей компоненты спектра возмущенной волны неограниченно возрастает при стремлении времени к бесконечности. Но при больших амплитудах не будет выполнено предположение о малости градиентов поверхности, которое лежит в основе вывода уравнений рассматриваемой модели. Поэтому деформации спектров будут достоверными только в некоторой окрестности кривой, на которой знаменатель в формуле (3.5) не обращается в нуль.

Для исследования возможности наблюдения деформации поверхности моря с помощью неконтактных измерительных приборов удобно именно представление деформированной поверхности в спектральной форме, так как для моделирования флуктуаций сигналов, полученных от радиометров и радиолокаторов удобно именно такое представление. Формулы для расчетов возмущений сигналов от радиометров и радиолокаторов приведены, например, в [2]. Имеется, как уже было отмечено выше, большой экспериментальный материал, подтверждающий связь возмущений океанских течений за счет неоднородности дна с радиолокационными и радиометрическими изображениями поверхности. Очевидно, между двумя изображениями на рис. 1 существует связь, но ведь электромагнитные волны в соленую воду не проникают! Это явление давно вызывало большой интерес и, несомненно, требовало теоретического объяснения. В данной работе рассматривается связь вышедших на поверхность течений и деформаций поверхностного волнового поля в спектральных и координатных представлениях. При этом локализованное в некоторой области поверхности течение порождает также локализованную в некоторой другой области флуктуацию поверхностного волнения, что вытекает из свойств преобразования Фурье и наблюдается также в экспериментах, о которых упоминалось выше.

Кроме того, спектр возмущения поверхности может содержать более широкий набор частот, чем спектр вышедшего на поверхность поля скоростей течения, поскольку в формуле (3.5) присутствует сдвиг спектра поля скорости возмущающего течения. Подобный эффект также наблюдается в эксперименте, когда крупномасштабное течение, которое в спектре не содержит коротких по длине волны компонент, вызывает изменение по всему спектру деформаций поверхности. Так, конвекционные течения, возникающие в океане в результате перегрева верхних слоев океана, влияют на коротковолновую (миллиметры- сантиметры) составляющую спектра поверхности океана. Изменения в этой части спектра можно зафиксировать в инфракрасном диапазоне собственного радиоизлучения поверхности океана даже с помощью радиометров космического базирования. Такого рода измерения используются, в частности, для предсказания возникновения тайфунов в зоне перегрева океанской поверхности.

5. Переход от спектрального представления возмущенного профиля поверхности моря к координатному

Естественный вопрос $–$ как получить форму возмущенной поверхности в исходных координатах $\left(x,t\right)$. Спектральное представление содержит в знаменателе множитель, который обращается в ноль на целой линии, поэтому непосредственное применение формулы обращения двумерного преобразования Фурье затруднено. Ниже предлагается способ решения уравнения вида (3.3) в явном виде с помощью формулы, эффективно реализуемой при численных расчетах.

Рассмотрим уравнение

$\ddot{\xi }+2{{\dot{\xi }}^{\prime }}_{\widehat{x}}+{{\xi }^{″}}_{\widehat{x}\widehat{x}}+\sqrt{-{\Delta }_{\widehat{x}\widehat{y}}}\xi =f\left(\widehat{x},\widehat{y},\widehat{t}\right)$

Будем считать, что при $t<t_0<0$ $f\equiv 0$ и $\xi \equiv 0$. Используя замену переменных $\widehat{x}=x+t$, $\widehat{y}=y$, $\widehat{t}=t$, получаем уравнение в переменных $\left(x,y,t\right)$

$\ddot{\xi }+\sqrt{-{\Delta }_{xy}}\xi =f\left(x,y,t\right)$

Положив $f=\delta \left(x\right)\delta \left(y\right)$ и используя преобразование Фурье по переменным $\left(x,y\right)$, найдем фундаментальное решение этого уравнения, через решение обыкновенного дифференциального уравнения

$\ddot{\widehat{\xi }}+\sqrt{k^2+l^2}\widehat{\xi }=0$

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

$\widehat{\xi }\left(0\right)=0$, $\dot{\widehat{\xi }}\left(0\right)=1$

Получаем $\widehat{\xi }\left(t,\rho \right)=\theta \left(t\right)\sin \left(\sqrt{\rho }t\right)/\sqrt{\rho }$, где $\rho =\sqrt{k^2+l^2}$, а $\theta \left(t\right)$ $–$ функция Хевисайда.

Обращая образы Фурье, получаем

$\xi \left(x,y,t\right)=\frac{1}{4{\pi }^2}{\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle \underset{{ℝ}^2}{\int }{\displaystyle \underset{{ℝ}^2}{\int }\frac{e^{ik\alpha +il\beta }}{\sqrt{\rho }}\sin \left(\sqrt{\rho }t\right)}dkdlf\left(x-\alpha ,y-\beta ,t-\tau \right)d\alpha d\beta d\tau }}$

Последний интеграл сходится, если функция $f$ достаточно быстро убывает на бесконечности.

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

Далее приводится пример расчета поля скоростей, вызванного движущимся источником [14] в координатном и спектральном представлении (см. рис. 7, 8). С помощью него можно осуществить моделирование возмущения поверхности также в спектральном и координатном представлении. Для этого нужно воспользоваться формулой (4.1). Именно, нужно получить картину поверхностного поля скоростей в спектральном представлении, сдвинуть эту картину на постоянный волновой вектор возмущаемой ветровой волны и, далее, умножить на передаточную функцию (которая, заметим при этом, имеет особенности на целой линии). Затем осуществляется переход к исходным координатам $\left(x,t\right)$, описанным выше методом.

Рис. 7. Карты компонент $v_x$, $v_y$ поля скорости течений на поверхности жидкости в координатном представлении

Рис. 8. Вещественные составляющие Фурье-образов компонент скорости течения на поверхности жидкости:

компонента $\mathrm{Re}{\widehat{v}}_x$ а: в области [$–$2,2]×[$–$2,2] и б: $–$ в области [$–$0.25,0.25]×[$–$0.25,0.25] вблизи начала координат (изометрическая проекция);

компонента $\mathrm{Re}{\widehat{v}}_y$ в: в области [$–$2,2]×[$–$2,2] и г: $–$ в области [$–$0.25,0.25]×[$–$0.25,0.25] вблизи начала координат (изометрическая проекция)

Поле скоростей на поверхности, изображенное на рис. 7, рассчитано с применением сеточного метода для уравнения Буссинеска с краевым условием свободной поверхности на верхней границе жидкости. Установлено хорошее соответствие между исходными полями скоростей и результатами применения прямого и обратного преобразования Фурье к ним (в спектральной области вблизи начала координат, см. рис. 8), что свидетельствует о достаточной точности численных методов при использовании преобразования Фурье.

Важный вопрос, который можно исследовать в рамках описанной модели $–$ влияние неоднородной стратификации жидкости на картину возмущений поля ветрового волнения.

Изложенные в настоящей работе методы позволяют осуществить «сквозное» моделирование:

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

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

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

Об авторах

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

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

А. С. Шамаев

Email: sham@rambler.ru
Россия, Москва

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

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

Доп. файлы

Действие

1. JATS XML

Скачать

2. Рис. 1. Радиолокационное изображение, полученное с помощью космического аппарата “Алмаз-1” в Северном море и батиметрическая схема участка съемки

Скачать (432KB)

Метаданные

3. Рис. 2. Явление резонансного рассеивания на периодической поверхности

Скачать (81KB)

Метаданные

4. Рис. 3. Обработка радиолокационного кадра с помощью процедуры усреднения.Сравнение различных масштабов задачи

Скачать (125KB)

Метаданные

5. Рис. 4. Схема обтекания подводного препятствия

Скачать (56KB)

Метаданные

6. Рис. 5. Трансформация ветровой ряби над обтекаемыми неоднородностями дна

Скачать (95KB)

Метаданные

7. Рис. 6. Поле скоростей на поверхности воды непосредственно над движущимся телом

Скачать (44KB)

Метаданные

8. Рис. 7. Карты компонент vx , vy поля скорости течений на поверхности жидкости в координатном представлении

Скачать (120KB)

Метаданные

9. Рис. 8. Вещественные составляющие Фурье-образов компонент скорости течения на поверхности жидкости: компонента  а: в области [2,2]×[2,2] и б:  в области [0.25,0.25]×[0.25,0.25] вблизи начала координат (изометрическая проекция); компонента  в: в области [2,2]×[2,2] и г:  в области [0.25,0.25]×[0.25,0.25] вблизи начала координат (изометрическая проекция)

Скачать (284KB)

Метаданные