Акустический журнал, 2020, T. 66, № 6, стр. 613-621

ВВЕДЕНИЕ

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

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

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

ПОСТАНОВКА ЗАДАЧИ И МОДОВОЕ ПРЕДСТАВЛЕНИЕ РЕШЕНИЯ

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

Рассмотрим трехмерное уравнение Гельмгольца

для акустического давления P(r, θ, z) в цилиндрических координатах $r,\theta ,z,$ где r – расстояние от оси симметрии волновода, θ – угловая (азимутальная) координата, а z – глубина (ось z направлена вниз).

Батиметрия в волноводе мелкого моря, показанном на рис. 1, описывается соотношением

где H₀ – глубина моря вдали от каньона (“невозмущенная” глубина), $\delta H$ – глубина каньона, V(r) – функция, описывающая форму поперечного сечения каньона (радиальный профиль возмущения глубины), такая что max(V(r)) = 1. Мы будем считать, что носителем функции V(r) является отрезок [r₁, r₂], т.е. что горизонтальной проекцией каньона является кольцо с границами r = r₁ и r = r₂.

На поверхности океана, как обычно, мы ставим условие Дирихле (условие мягкой границы)

На границе раздела вода–дно z = H(r, θ) поставим условия непрерывности акустического давления и нормальной компоненты колебательной скорости, которые имеют вид

где n – единичная нормаль к границе раздела. Отметим, что здесь и далее в работе используются обозначения: F_± ≡ F|_{z=H± 0} (т.е. нижние индексы + и – означают значения физических величин по разные стороны от границы раздела вода–дно).

Для того чтобы исключить появление непрерывного спектра, ограничим расчетную область снизу плоскостью z = H_max. На этой границе потребуем выполнения условия Дирихле (условия мягкой границы)

В целях обеспечения единственности решения краевой задачи для уравнения Гельмгольца (1) потребуем также выполнения условия излучения [1] при

Решение краевой задачи для уравнения Гельмгольца (1) с граничными условиями (3), (6) и условиями на границе раздела (4), (5) будем искать в виде модового разложения

где собственные функции ${{\phi }_{j}}\left( {z,r,\theta } \right)$ получаются из решения следующей акустической спектральной задачи: где скорость звука c = c(z), вообще говоря, зависит от глубины и, в частности, имеет скачок на интерфейсе (границе раздела вода–дно). Заметим, что в задаче Штурма–Лиувилля (8) условие на границе раздела содержит производную по z, в отличие от условия (5), в которое входит производная по нормали к интерфейсу [1].

После подстановки (7) в уравнение (1) получаем

В условиях данной задачи, очевидно, функции ${{\phi }_{m}}\left( {z,r} \right)$ зависят только от горизонтальной переменной r (как от параметра).

После применения оператора $\int {( \cdot )} \frac{1}{\rho }{{\phi }_{j}}(r,\theta ,z)dz$ (где ρ(z) – зависимость плотности от глубины) к последнему уравнению и упрощения полученных выражений с помощью условия нормировки собственных функций

воспользуемся адиабатическим приближением и отбросим члены, ответственные за взаимодействие мод, как это сделано в [1]. В результате мы получим отделенные друг от друга уравнения для модовых амплитуд ${{A}_{m}}(r,\theta ){\text{:}}$также известные как уравнения горизонтальной рефракции. С математической точки зрения они являются двумерными уравнениями Гельмгольца и, таким образом, обеспечивают редукцию размерности исходной математической модели. В дальнейшем будем рассматривать краевую задачу, получающуюся добавлением к уравнению (10) условия излучения при r = $\sqrt {{{x}^{2}} + {{y}^{2}}} \to \infty .$

СТРУКТУРА ПОЛЯ В ГОРИЗОНТАЛЬНОЙ ПЛОСКОСТИ: МЕТОД РАЗДЕЛЕНИЯ ПЕРЕМЕННЫХ

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

и, следуя общей схеме метода Фурье, будем искать его решение в виде

Подставим это выражение в уравнение (11) и выполним разделение переменных, разделив полученное равенство на Q(r)ψ(θ) и умножив на r²:

Оставляя в (12) слагаемые, зависящие от r, в левой стороне уравнения, и перенося слагаемые, зависящие от θ, в правую сторону, получаем

Поскольку величина в левой части уравнения зависит только от r, а в правой – только от θ, обе они равны некоторой постоянной a² (постоянная разделения). Тогда получим следующее уравнение для азимутальной компоненты частного решения $\psi \left( \theta \right)$

Добавляя к данному уравнению граничные условия и условие периодичности при θ = ±π, получаем “азимутальную” задачу Штурма–Лиувилля (ШЛ), где a² есть спектральный параметр. Заметим, что на этом этапе обнаруживается принципиальное отличие описываемой здесь процедуры решения от метода построения поля шепчущей галереи в работе [7], где решалась задача ШЛ по радиальной координате.

Нетрудно заметить, что собственные значения задачи ШЛ для уравнения (14) с условиями периодичности суть $a_{m}^{2} = {{m}^{2}},$ m = 0, 1, 2,… Для каждого собственного значения m² существуют две линейно независимых собственных функции, а именно: $\psi _{m}^{0}(\theta )$ = $\frac{1}{{\sqrt \pi }}\cos (m\theta )$ и $\psi _{m}^{1}(\theta )$ = $\frac{1}{{\sqrt \pi }}\sin (m\theta )$ (значение верхнего индекса 0 относится к косинусам, а 1 – к синусам). Для удобства будем полагать, что θ_s = 0 (очевидно, это не ограничивает общность). В этом предположении решение A_m будет представлять собой суперпозицию $\left\{ {\psi _{m}^{0}} \right\},$ так как $\psi _{m}^{1}(0)$ = 0 при всех m.

Решение уравнения (10), таким образом, представимо в виде разложения

Подставляя суперпозицию гармоник (15) в уравнение (10), получаем

Пользуясь тем, что $\frac{{{{\partial }^{2}}{{\psi }_{m}}}}{{\partial {{\theta }^{2}}}}$ = –m²ψ_m, упростим последнее равенство, после чего применим к обеим его частям оператор $\int {( \cdot )} {{\psi }_{m}}(\theta )d\theta .$ Воспользовавшись ортогональностью системы собственных функций $\psi _{m}^{0}(\theta ),$ получим уравнения для радиальных компонент решения Q_j,m(r)

Вне отрезка [r₁, r₂] уравнение (16) имеет постоянные коэффициенты и фактически сводится к уравнению Бесселя. Из условия ограниченности решения при r = 0 мы делаем вывод, что при r ≤ r₁ оно должно с точностью до константы представлять собой функцию Бесселя

В то же время из условия излучения на бесконечности следует, что при r ≥ r₂ решение должно обращаться в функцию Ганкеля первого рода, т.е.

Теперь осталось лишь найти значения функции Q_j,m(r) на отрезке [r₁, r₂]. С этой целью составим нормальную систему, эквивалентную уравнению (16)

Решение этой системы может быть найдено методом двустороннего численного интегрирования. Сначала мы решаем ее на отрезке [r₁, r_s], начиная с точки r = r₁ и используя начальное условие

с тем чтобы найти величину решения и его производной Q^left = Q_jm(r_s), P^left = P_jm(r_s) слева от точки излучения r = r_s. После этого эта система решается на отрезке [r_s, r₂] с начальным условием при r = r₂. В результате получаем значения функции Q(r) и ее производной Q^right = Q_jm(r_s), P^right = = P_jm(r_s) справа от точки излучения r = r_s.

Наличие в правой части (16) дельта-функции (источника, сосредоточенного при r = r_s) эквивалентно тому, что решения системы (19) удовлетворяют следующим условиям

при r = r_s (непрерывность Q(r) и наличие у P(r) разрыва первого рода).

Равенства (22) образуют систему линейных алгебраических уравнений, из которой легко определить коэффициенты α_jm и ${{\beta }_{{jm}}}$ и, таким образом, найти вид решения (16) на отрезке [r₁, r_s].

ПРИМЕР РАСЧЕТА ЗВУКОВОГО ПОЛЯ

Рассмотрим в качестве примера двухслойный волновод мелкого моря, глубина которого без учета каньона составляет H₀ = 90 м, скорости звука принимают значения c_w = 1500 м/с и c_b = 2000 м/с в водном слое и дне соответственно, а плотности в этих в этих средах равны ρ_w = 1 г/см³, ρ_b = 2 г/см³. Пусть точечный источник с частотой f = 50 Гц располагается на глубине z = z_s = 10 м и на удалении r_s = 5000 м от оси симметрии каньона. Глубину каньона в нашем примере мы примем равной δH = 10 м, а его радиальный профиль на отрезке [r₁, r₂] зададим формулой

где r₁ = 4800 м, r₂ = 5200 м. Задачу расчета волновых чисел и собственных функций в данном волноводе решаем методом конечных разностей [1] с шагом по глубине 0.25 м и нижней границей расчетной области $z = {{H}_{{{\text{max}}}}} = 500$ м (заметим, что значение этой величины достаточно велико для того, чтобы не сказываться на значениях волновых чисел мод дискретного спектра). Результаты расчета волновых чисел для различных значений глубины моря в нашем примере представлены в таблице 1.

На рис. 2 и 3 показаны графики уровней акустического поля при ${{N}_{\theta }} = 800$ и ${{N}_{\theta }} = 1200$ (в дБ отн. уровня в 1 м от источника). Видно, что первые 800 членов ряда (15) (рис. 2) описывают компоненты поля, для которых волновые векторы образуют достаточно большие углы с касательной к окружности вида r = const, в то время как для захваченных каньоном волн волновые числа не превышают некоторых критических значений (они будут вычислены в следующих разделах). Интерференционная картина поля, представленная на рис. 3, стабилизируется, начиная с значения N_θ = 900. Таким образом, фактически при расчете поля во всех случаях достаточно суммирования лишь конечного числа членов ряда (15), однако это число зависит от r_s. На рис. 3 хорошо видно, что часть акустической энергии источника фокусируется над каньоном, причем можно заметить, что поле в этой области имеет выраженную модовую структуру в горизонтальной плоскости.

Интeресно также отметить, что помимо фокусировки над каньоном наблюдается еще один эффект, связанный с горизонтальной рефракцией. Он заключается в формировании узкого пучка, распространяющегося в направлении, близком к касательному к окружности r = r_s (более точно, пучок направлен из точки x = 5 км, y = 0 в точку x = 4 км, y = 6 км). Уровень акустического поля в этом пучке заметно превышает уровни, формируемые источником в других направлениях на том же расстоянии от точки излучения. Забегая вперед, отметим, что такому пучку соответствуют горизонтальные лучи, которые выходят из источника под углами относительно касательной к r = r_s, близкими к нулю, но все же превышающими критический угол, при котором осуществляется их захват каньоном.

ГОРИЗОНТАЛЬНЫЕ ЛУЧИ И ОЦЕНКА ДОЛИ АКУСТИЧЕСКОЙ ЭНЕРГИИ, ФОКУСИРУЕМОЙ КАНЬОНОМ

Исследуем свойства решения краевой задачи для уравнения (10) с использованием теории горизонтальных лучей. Здесь и далее будем предполагать номер моды j фиксированным и опускать соответствующий нижний индекс во всех выражениях.

Модовые амплитуды A_m в лучевом приближении представляются в виде

где $\bar {k}$ – отсчетное значение волнового числа (удобно принять в качестве $\bar {k}$ волновое число моды в точке излучения), A⁽⁰⁾(r, θ) – амплитуда горизонтального луча, а S(r, θ) – эйконал, удовлетворяющий уравнению Гамильтона–Якоби где ${{\nabla }_{ \bot }} = \left( {\frac{\partial }{{\partial r}},\frac{1}{r}\frac{\partial }{{\partial \theta }}} \right)$ – горизонтальная компонента градиента, n(r) – индекс рефракции, соответствующий выбранной вертикальной моде, который определяется формулой n(r) = k(r)/$\bar {k}.$

Для нахождения траектории лучей введем в рассмотрение касательные к ним векторы e = (e_r, e_θ) = = ∇_⊥S, где e_r и ${{e}_{{\theta }}}$ суть радиальная и тангенциальная компоненты данного вектора.

Система уравнений Гамильтона для лучевых траекторий имеет вид [8]

где τ – длина дуги вдоль луча. Из данного уравнения следует, что d(re_θ)/dτ = 0, поэтому re_θ = const = = σ. Из последнего соотношения легко получить закон Снеллиуса для радиально-слоистой среды [8, 7 ]: где χ – это угол между вектором e и его радиальной компонентой. Из (24) и (25) следует, что часть горизонтального луча между двумя точками заворота определяется формулой

В первую очередь нас будут интересовать лучи, захваченные волноводом в горизонтальной плоскости, который формируется в связи с наличием подводного каньона (см. рис. 4). Для захваченных лучей характерно наличие точек заворота r = $r_{1}^{t}$ и r = $r_{2}^{t},$ лежащих внутри кольца r₁ ≤ r ≤ r₂. Так как в точках заворота sin(χ) = 1, то для луча, выходящего из источника под углом χ = χ_s, при r = r_s из закона Снеллиуса получается следующее равенство

где l ∈ {1, 2. Для получения оценки величины сектора лучей, захватываемых каньоном, удобно ввести в рассмотрение угол α = π/2 – χ между горизонтальным лучом и касательной к окружности r = const в данной точке.

Рассмотрев график функции rk_j(r), представленный (для j = 4) на рис. 5в, можно получить из последнего соотношения следующее неравенство для захваченных каньоном лучей

Каньоном захватываются только лучи с углами выхода из источника α_s, для которых одновременно выполняются условия α_s < α_{s, 1} и α_s < α_{s, 2}, где по определению

Отметим, что из рис. 5 видно, что первое неравенство является заведомо менее ограничительным, чем второе. Таким образом, долю D акустической энергии источника, переносимой данной модой и захваченной каньоном, можно оценить по формуле

Действительно, D совпадает с отношением углового сектора, заполняемого захваченными лучами (очевидно, что он равен 4α_{s, 2}), к общей угловой мере 2π семейства лучей, выходящих из источника. Оценка этой величины в зависимости от радиальной координаты r_s точки излучения показана на рис. 5а (все значения параметров – такие же, как в предыдущем разделе). Видно, что в этом случае в толще воды непосредственно над каньоном фокусируется около 5% энергии, переносимой четвертой водной модой.

СВЯЗЬ ГОРИЗОНТАЛЬНЫХ ЛУЧЕЙ С ЧЛЕНАМИ РЯДА В РАЗЛОЖЕНИИ МОДОВЫХ АМПЛИТУД

Пусть r_max – точка, в которой функция rk(r) имеет максимум (см. рис. 5). Тогда из закона Снеллиуса (25) следует, что самый большой угол α между некоторым захваченным лучом и касательной к окружности в данной точке будет именно при r = r_max. Для этого угла α(r_m) будет выполняться пара неравенств

Далее, волны вида cos(nθ)Q_n(r) по азимутальной переменной являются стоячими, причем

где l = rθ – длина дуги окружности. Поэтому можно считать, что горизонтальное волновое число k(r) складывается из нормальной k_r и тангенциальной k_𝓵 компонент, причем k_𝓵 = $\frac{n}{r}.$ Очевидно, что угол α между горизонтальным волновым вектором и касательной к окружности и тангенциальная компонента волнового числа связаны соотношением

Таким образом, n = rk(r)cos α(r) = r_maxk(r_max)  × × cos α(r_max). Иными словами, лучевой инвариант Снеллиуса σ в волновой интерпретации и есть число n, определяющее зависимость компоненты решения от угловой переменной.

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

В рассматриваемом нами примере максимальное значение функции rk₄(r) на интервале [r₁, r₂] достигается в точке r_max = 5082 м и составляет r_maxk₄(r_max) = 887.9959 (см. рис. 5в). Нижняя граница номеров n, соответствующих захваченным лучам, есть k_minr₂ = 874.6199. Таким образом, на основании этих данных неравенство (31) для рассматриваемой задачи выглядит следующим образом

На рис. 6 показаны графики функций Q_4,m(r). Видно, что энергия соответствующих компонент решения Q_4,m(r)ψ_m(θ) локализована в горизонтальной плоскости в относительно узком кольце, являющемся проекцией “внешней” половины каньона. Контурный график акустического поля, формируемого этими компонентами, отдельно показан на рис. 7.

ЗАКЛЮЧЕНИЕ И ВЫВОДЫ

В настоящей работе изложена методика расчета звуковых полей в волноводах мелкого моря с симметрией относительно вращений в горизонтальной плоскости в случае, когда источник не находится в центре симметрии. Расчеты проведены в адиабатическом приближении, т.е. пренебрегается взаимодействием мод. В дальнейшем изложенная здесь методика будет расширена на случаи, когда взаимодействием мод пренебречь нельзя, в частности, на волноводы с такими вариациями глубины, которые приводят к отсечке одной или нескольких водных мод. С теоретической точки зрения такое обобщение почти очевидно, и вместо нормальной системы (19) просто получим связанную систему из 2N связанных обыкновенных дифференциальных уравнений, а условие согласования (22) сведется к решению системы линейных уравнений с 2N неизвестными. Такая процедура построения Q_m(r), однако, приводит к некоторым вычислительным сложностям (в частности, матрица в аналоге системы (22) оказывается плохо обусловленной) и нуждается в регуляризации, которая выходит за рамки настоящей работы.

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

Заметим также, что изложенный здесь подход в некоторых аспектах схож с методами расчета звукового поля, рассеянного подводной горой, предложенными в работах [9, 10]. Однако, в отличие от их авторов, мы уделяем особое внимание установлению связи между отдельными элементами разложения модовых амплитуд по гармоникам вида Q(r)ψ(θ) с физическими свойствами описываемых ими волн. Мы также устанавливаем связь этих компонент звукового поля некоторой моды с соответствующими ей горизонтальными лучами, что, в частности, позволяет нам получить простую лучевую оценку доли акустической энергии, сфокусированной в водном слое непосредственно над подводным каньоном. Несмотря на очевидную идеализацию в рассмотренной нами постановке задачи, выражающуюся в замкнутости каньона, лучевая оценка (29) остается справедливой и в случае каньона с открытыми концами (например, если его проекция на плоскость (x, y) представляет собой половину кольца).

Подчеркнем еще раз связь нашей работы с недавней статьей [7]. Глядя на рис. 5в, внимательный читатель легко заметит, что внутренняя граница каньона r = r₁ фактически не принимает участия в фокусировке звукового поля, и его локализация происходит в кольце между r = r₂ и второй точкой пересечения штрихпунктирной горизонтальной линии с графиком функции rk(r). По-видимому, в большинстве случаев, которые можно представить, эта точка будет находиться правее границы r = r₁, и, таким образом, захваченные каньоном лучи в таких случаях будут фактически лучами шепчущей галереи [7].

Работа выполнена в рамках темы госзадания ТОИ ДВО РАН “Математическое моделирование и анализ динамических процессов в океане” (№ 0271-2019-0001) и поддержана грантами РФФИ 18-05-00057_а и 18-35-20081_мол_а_вед. Авторы также признательны д-ру А.В. Цветковой (ИПМех РАН) и рецензентам работы за ряд предложенных поправок, которые позволили существенно ее улучшить.