Моделирование волновой динамики стратифицированных сред с учетом нелинейности, вязкости, вращения Земли и сжимаемости Текст научной статьи по специальности «Физика»

УДК 539.59:534.1

Моделирование волновой динамики стратифицированных сред с учетом нелинейности, вязкости, вращения земли и сжимаемости

© В.В. Булатов1, Ю.В. Владимиров1 1 ИПМех им. А.Ю. Ишлинского РАН, Москва 119526, Россия.

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

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

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

Постановка задачи. Система уравнений гидродинамики с учетом нелинейности, вязкости, сжимаемости среды и вращения Земли в адиабатическом приближении для уравнения состояния имеет вид [1, 3-5]

Ор ,ТТ

Р— + ^-=-р/и2 + У ш Ох

ёи2 Ор гтт

р~Г+~г = -Р/и1+у

ш оу

о

А 3и1 + — ёгуи

Ох

о

А 3и2 +—ё1уи

Оу

+ К

+ К

У'

ШЖ Ор

Р-+ + ЕР = у

Ш Ох

О

А3Ж +—&уи Ог

+ ;

(1)

1 ёр = ёр с2 Ш Ш

— + Жури = М, Ог

причем

Ш О тт О тт О тжг О

— =—+и1—+и2—+ж—.

Ш Ог Ох Оу Ох Здесь р, р — давление и плотность соответственно; / — параметр Кориолиса, / = 20 бшЗ; О — частота вращения Земли, 0 = = 7,27 • 10-5 с-1; 3 — географическая широта; у — коэффициент

вязкости, у = 10 2 см2/с; А3 = А +

О2

Л О2 О2 тт А = ——+ ——; Т — вектор

Ох2 ' Ох2 Оу2 скорости, и = (и1,и2, Ж); Гх, Гу, — плотности действующих на жидкость объемных сил; е — ускорение свободного падения; с — скорость звука в океане, с = 1,5 • 105 см/с2; М — плотность источников массы [6-8].

Далее будем рассматривать или безграничную по вертикали среду (ось г направлена вертикально вверх), или слой, ограниченный дном г = -Н и свободной поверхностью 2 = д(х, у, г). При наличии вязкости на дне необходимо задать условие прилипания: и1 = и2 = Ж = 0 при г = - Н (у, г), которое формировало бы соответствующий пограничный слой. Однако для внутренних гравитационных волн в океане характерны малые значения скоростей частиц (около 10 см/с и ниже) и соответственно малые значения градиентов

скорости [9-11]. Ввиду этого на дне (2 = -Н(х,у)) имеет место условие непротекания

W-U, —-U2 dH = 0. (2)

dx dy

Для горизонтального дна (H = const, z = -H) это условие упрощается:

W = 0. (3)

На свободной поверхности океана (z = с;(x, у, z)) задают два граничных условия: кинематическое и динамическое. Согласно кинематическому условию, у частиц жидкости на поверхности океана нормальная к поверхности компонента скорости должна совпадать со скоростью смещения поверхности:

W = д- + и, | + U2 d?£ = d- (Z = f). (4)

dt dx dy dt

Согласно динамическому условию, давление на поверхности океана должно совпадать с атмосферным, которое далее полагаем равным нулю:

p [x, y, x y, t)] = (5)

Систему уравнений (1) и граничные условия (2)-(5) линеаризуем относительно состояния покоя:

Z

и1 = и2 = W = 0; p = p°(z); p = p°( z) = -gjp°( Z )dz.

°

Для этого положим p = p° + p*, p = p° + p*; U, = Uf, U2 = U*, W = W* и запишем уравнения для p*, p*, U*, U*, W* (индекс «*» далее опускаем):

dU, dp Q F S p° — + = Qx + Fx = Sx; dt dx

dU2 dp Q F

p°^r =Qy + Fy = Sy;

dt dy dW dp

p° — + -Z- + gp = Qz + Fz=Sz; (6)

dt dz

dU, dU2 dW p° —L + —2 +-

dx dy dt

p+ w dP° = T dt dz

= R + M = Ф;

Здесь правые части уравнений Qx, Qy, Qz, Я, Т содержат члены,

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

Qx =

Ои1 ( + )

= -Р—— (Р + Р0)

Ог

тт Ои1 тт Ои1 тжгОи1

+ + Ж—т1 + /и2

Ох

Оу

Ог

+ у

Л ТТ ОБ

Ази1 + — Ох

Ои2 . -Р—--(Р + Р0)

Ог

и, и + и2

Ох

Qy =

Ои2 Ои2

Оу

-ж-

Ог

- ¡их

у

А тт ОБ

Ази2 + — Оу

ОЖ ,

-Р(Р + Р0) Ог

Qz =

„ ОЖ ОЖ ТТ,ОЖ

и1-+ и2-+ Ж-

Ох Оу Ог

+ у

А зЖ +

ОБ Ог

Я = --V

С 2

Ор тт Ор тт Ор „ , Ор

—+и1—+и2 —+Ж—

Ог Ох Оу Ог

-РБ;

Т = и ОР-и2 Р-Ж Р Л

Ох

Оу

Ог с2

Ор тт Ор тт Ор ТТЛ Ор

— + и1 — + и2 — + Ж—

Ог Ох Оу Ог

_ Ои1 Ои2 ОЖ где Б = —1 + —2 +-.

Ох Оу Ог

Внутренние гравитационные волны имеют на поверхности (г = 0) весьма малые возмущения д, поэтому, линеаризуя на ней граничные условия, имеем

Ж = Од/Ог; р(х, у, 0, г) = д(х, у, г)ЕР0(0).

Отсюда можно получить

Ор

"О7

- Жр(0) = 0.

(7)

Исключая из системы уравнений (6) переменные и1, и2, Р и р, находим для вертикальной компоненты Ж скорости обычное уравнение внутренних гравитационных волн с некоторой ненулевой правой частью [1, 6-8]:

d 2

-Чг A3W + N 2( z )AW = Z, dt

(8)

где N 2( z) — частота Брента—Вяйсяля, а

Z

N 2( z) д 3W 1 d2

g dzdt2 р0 dzdt

dSx dSy

dx dy

дФ

dt

AT.

Po dt ро

Граничное условие на дне (z = -H) сохраняет вид (2) или (3). Применяя к выражению (7) оператор A и выражая Apt через W, получаем условие при z = 0:

d3W - A W = J_ dzdt2 р0

d 2Ф d2 Sx d2Sy

dt2 dtdx dtdy

(9)

В уравнении (8) и граничном условии (9) правые части — это сумма слагаемых, зависящих от сторонних источников (массовых сил Е и плотности источников массы M), и малых поправок, учитывающих вязкость, сжимаемость и вращение среды, а также поправок, обусловленных использованием приближения Буссинеска и имеющих порядок N2/g « 1.

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

d2

-Чг A3W + N 2( z ) AW = sP(W), dt2

(10)

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

Исследуем вопрос о том, как влияет поправка sP на распространение одной отдельной волновой моды внутренних гравитационных волн. С целью упрощения выкладок рассмотрим плоский случай, т. е. U2 = 0. Обозначим U1 = U. Тогда задачу сформируем следующим образом. Пусть поправка sP «включается» при t = 0:

dt2

2 i я2

d2W | d2W dz2 dx2

Л

d2W

+ N 2( z ) d?r = sX(t) P(W ),

(11)

причем x(t) = 0 при t < 0 и x(t) = 1 при t > 0, а P(W) определяется выражением

P =

д2

dxdt

W

(д 2U д 2U Л _(д 2W d2W Л

dz2 дх2

- U

дz2 дх2

+Р

( д2и д 2W Л дtдz дtдx

др ди-дрдЖ_ дz дt дх дt

+g 4 г и р w р.

дх2 \ дх дz J

Здесь р — возмущение плотности, нормированное на некоторое типичное значение невозмущенной плотности р0.

Для дальнейшего анализа при t < 0 можно использовать решение уравнения (11) c нулевой правой частью в виде собственной волновой моды:

W0 = A<n (z, к) cos \an (к)t - кх],

где А — константа; к — спектральный параметр; <pn (z, к), a (к), — нормированные собственные функции и дисперсионные кривые основной вертикальной спектральной задачи внутренних гравитационных волн соответственно [1, 6-8]:

д2<п (z, к) + к 2 г N2(z) -1

д 2

Л

а2(к )

<п (z, к) = 0.

(12)

Отметим, что при г = 0 и г = -Н функции срп (г,к) = 0.

Тогда соответствующие нулевому приближению Жо горизонтальная скорость ио и возмущение плотности ро имеют вид

U0 = A <(z к) sin\an(к)t - кх];

Р0 =

к дz AN 2( z)

(13)

an(к)g

<рп (z, к) sin\an (к)t - кх].

При I > 0 решение (11) будем искать в виде ряда по степеням малого параметра а:

Ж = Жо + +е2Ж2 +... Тогда для определения функции Ж1 получаем уравнение

д2 (д2щ д2щ Л лг2/ чд2щ

1 _|___ _1_ А/2' 1

дt2

дz2 дх2

+ N 2( z = P(Wo).

(14)

Решение уравнения (14), совпадающее с Ж0 при t < 0, непрерывно вместе со своей производной по t. Поэтому при t = 0 функция Щ и ее производная по t обращаются в нуль: Ж1 =дШ1/дt = 0. Поскольку в реальных океанических условиях возбуждаются, как правило, только первые моды, в дальнейшем для определенности будем рассматривать первую волновую моду. Тогда правая часть уравнения (14) примет вид

P(Wo) =

= A22k бш [ (к )г - 2кх]

N 2( 2 )^(к)

gк

т z, к)

д 1д>1( 2, к) ^ дт( 2, к)

д22

д2

2

2 N 2( 2 )кюх(к) 2, к )]2

g

а>1(к) ( дт(2, к) д2т(2, к)

д2

д2 2

z, к )

д 3т( 2, к)

д2 3

(15)

®1(к)т(2, к) дт(2, к) дN 2(2)

2, к)Л2 дN2(2) I _

gk д2 д2 а>1 (к) ^ д2

_ Ф(2, к) в1п(2®1 (к)t - 2кх).

д2

Решение уравнения (14) будем искать в виде ряда по собственным функциям задачи (12):

Ж1 = Бт [(к > - 2кх ] йт (2,2к).

(16)

г =1

Функцию Ф( 2, к) в правой части уравнения (15) можно представить в виде

Ф ( 2, к ) = N 2( 2 )£ от ( 2,2к ),

г=1

(17)

где Ог =| Ф(2, кт (2,2к)й2.

- я

Подставляя выражения (16), (17) в уравнение (14), получаем

.2/

d = . Cp2(lk )

4k 2(4<(k) - of (2k)) Учитывая начальные условия, имеем

Wi =

=Z df( z,2k )

г=1

sin[k)t-kx]- o(2k)sin2(2kx) + 2o1(k)cos2(2kx)

О (2k )

x sin [2o1 (k)t - 2kx] - [o (k) - О (2k)] cos [2o1 (k)t - 2kx]

(18)

Из выражения (18) следует, что наибольший вклад в W1 дает слагаемое с множителем d1. Сравним поправку W1 с невозмущенным решением W0. Для этого заменим первое слагаемое в уравнении (18) резонансным членом (считая o1(2k) = 2o1(k)), чем только увеличим значение W1 :

W; = a°1(22k ) f( z,2k )t cos [[(2k )t -2kx ], (19)

8k

и оценим время, за которое W* станет сравнимо с W0. Для расчета коэффициента a1 можно принять, что ((z,k) «((z,2k), что справедливо для достаточно малых значений k. В результате получаем

„1 = 2A2 J N2(z)«(z,^{(^dz -

-J N2(z)(fz,k))3dz + Щ- J N2(z)(z,k))22)dz-

g -H o1( k ) -H Ôz

-3o(k )k 2 J (( z, k ) )2 dz J.

(20)

Чтобы воспользоваться формулами (19), (20), рассмотрим слой стратифицированной среды с глубиной H = 100 м при N ( z ) = const =

= 0,01 с-1. Тогда

V2sin(^z /H )

" = n4H

При к = 0,02 м-1 длина волны X ~ 300 м и A = 10-4м3/2с-2, что соответствует амплитуде скорости невозмущенной волны «14 см/с. Для случая N (z) = const последние два слагаемых в выражении (20) равны нулю, а первые два интеграла можно вычислить аналитически, и тогда a1 =-1,4 • 10-13. Время, за которое величина W* составит около 5 % от значения W0, будет t = 4к2A • 0,05/[с1®1(2к)] « «1,5 • 107 с. Полученная оценка показывает, что интервал времени, в течение которого нелинейные слагаемые становятся сравнимыми по порядку величин с линейными, существенно превышает типичные периоды колебаний внутренних гравитационных волн в океане, составляющие десятки минут [1, 10, 11].

Оценка влияния вращения Земли. Рассмотрим теперь моделирование дальних полей внутренних гравитационных волн, возбужденных в стратифицированной вращающейся среде, например, расположенным в точке x = y = 0, z = z0 импульсным ^-образным источником [2, 3]. С учетом приближений Буссинеска и «твердой крышки» вертикальная компонента скорости W удовлетворяет уравнению

d2 „ d2W

A3W + f2 ^ + N 2( z) AW = E(x, y, z, z0, t)

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

W = 0 (z = 0, - H); W = 0 при t < 0.

Вид правой части E (x, y, z, z0, t) зависит от характера источника. Например, в случае источника массы

E = 5{ x)5{ y) ,

dz dt2

а в случае источника силы, направленной по вертикали, „ dS(t) .(d2S(x) . . d2S(y)

E = Ч^5(z - z0) —¿Г-5(y) + 5(x)^rrL

dt у dx dy j

Рассмотрим теперь функцию Грина, т. е. решение задачи с правой частью E = 5(x)S(y)S(z - z0)S(t). Для этого воспользуемся методом разделения переменных (методом Фурье). Положим при Im а> > 0

О, и,У, г, го) =

1 да да да

=-з | ехр(/'®^)& | | ехр [¡(¿их + у у)] ](¿, х, у, г,г0 )ёхёу.

Функция <р(а, и,У, г, г0) удовлетворяет уравнению

(«2 _ /2)К + к2[N2(2)_»2

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

<р = 0 (г = 0, - Н),

где к2 = и2 +у2.

Если приближение «твердой крышки» при г = 0 заменить приближением свободной поверхности, то граничное условие < = 0 примет вид

д<(0) = gk2

& ®2 _ /2

Решение будем искать в виде линейной комбинации собственных функций <р„ (г, к) граничной задачи

<) + к 2 <( г, к) = 0,

дг2 к) _ /2 ^ ' ' '

или при < = <п, <п (0) = 0

д< (0) = gk2

<(0).

<п (0).

дг ®1(к) _/•

Эта задача имеет на отрезке [_Н, 0] полную систему собственных функций срп (г, к), причем для этих функций выполняется условие ортогональности

0

g<n (0) < (0) +| (N2(г) _ /2) <п (г) <т (= Зтп,

_ Н

где 5тп — символ Кронекера, и в приближении «твердой крышки» внеинтегральный член исчезает. Чтобы выразить решение через собственные функции <п (г, к), заметим, что в силу полноты и ортонор-мированности этих функций

S(z - zo) = [[2(z) - f2 ] ¿ <(z) <(zo).

n

от

Положив < =^ pn<n ( z, k )<n ( z0, k ), получаем

n

—2(k) - f2

pn ='

(2л)3 k2 [—(k) - f2 ] и в результате имеем

— (k ) - f2 )< ( z , k )<n ( zo, k )

(2 л)3k2 (a2 -—(k))

Проведем теперь обратное преобразование Фурье. Примем

от+i 0

g(/,v,t,z,zo) = j exp(-—)<(u,v—z,zo)d—

-OT+Ю

Подставляя вместо < полученный ряд и проводя почленное интегрирование, находим

i от B

g(/ v t, z, zo) = —ц Z k— "(k) (—«(k) - f 2 ) Vn (z, k) < (zo, k),

где Bn = exp [— (k )t ] - exp [-— (k )t ].

В результате представление для функции W принимает вид

от

W = £Wn ( X, y, z, zo, t),

n

где

OT+i o OT+Ю

Wn (X, y, z, zo, t) = —2 í í exp [i(/x + vy)]x

-от+i o -от+i o

Bn [(k) - f2 ] ( J ) ( J )d d

x —k—k— < ( z, k ) < ( zo, k )dvd /и.

k —n (k)

После перехода к полярным координатам / = k cos у, v = k sin у, x = r cosy, y = r sin у и интегрирования по у получаем, что Wn бу-

дут функциями r = ijx2 + y2 :

n

г г B

Wn(r,z,zo,t) = — F Jo(kr) -—(a2 - f 2)yn(z,к)yn{zo,k)dk.

&.1T ■> fr ml frl

4ж

к an (к )

Поскольку рассматривается асимптотика Wn при больших r (дальнее поле), функцию Бесселя J0(kr) можно заменить на ее

асимптотику при кг » 1, т. е. J0(kr) «^2/(nkr) • cos(kr-ж/4). В результате имеем

Wn =

exp(ix / 4)

4я32у[2г

i

exp [i(kr - ant)] g(k2 )dk exp [i(kr + at)] g(k2 )dk

Vk

i

4k

(21)

Здесь g(k2) = an (k)—f_ ^ (z, k(z0, k) — четная функция по k, а k an (k )

функция an (k) продолжена в область k < 0 нечетным образом, точка ветвления k-1/2 обходится в верхней полуплоскости.

Дисперсионная кривая функции an (k) есть монотонно возрастающая функция по k. Поэтому фазовая функция второго интеграла в выражении (21) не имеет стационарных точек и этот интеграл при t, r ^<х>, t / r = const экспоненциально мал.

Рассмотрим первый интеграл в формуле (21). Положим /,ч dan

cn (k) =- — групповая скорость n-й волновой моды, cn =

dk

= max cn (k). Ясно, что при r > cnt первый интеграл также экспонен-

k

циально мал. При r < cnt можно применить обычный метод стационарной фазы. Обозначим через Xj = Xj (r/t) стационарные точки первого интеграла, т. е. решения уравнения r/1 = cn(k). Вычисляя при t, r ^-<x>, r/1 = const асимптотику первого интеграла методом стационарной фазы, получаем

Wn ( z, z0, r, t )i

exp(i^ / 4) 442n

(

exp г

, ч ж . dc„ Xjr-an (Xj)t --sign ^^

k=Xj

dc„ (k)

дk

k=Xj

где

Т ®1(Х])" I ( ) ( ) Т =——7 Р" (^Р" ()-

X] Юп {%] )

При г « спг стационарные точки первого интеграла сливаются и асимптотика интеграла оказывается зависящей от того, где достигается максимум групповой скорости сп (к). Если I Ф 0, то для типичных профилей N2(г) максимум будет в единственной точке х* > 0. В окрестности этой точки фазовая функция допускает разложение

кг-тп(к)г = Хг -ю*г + (к-х*)(г -с*г)--(к- X)3 д2с"(х ) г +...

Ох

2

где ю" = ю„ (х*), с" = Сп (X ).

Аппроксимируя фазовую функцию кубическим полиномом и заменяя внеэкспоненциальный множитель под знаком интеграла его выражением при к = х*, находим

Ж

ехр \\ / 4 + ¡(х*г -ф)] юШ2 -12

4^3/2>/2г

х*3/Ю

Р 0?, х* )Рп( го , х ) X ,

где

X =| ехр <I

(к-х)(г - сЩ) --(к-х)3 Г

*)3(д2сп(х*) 6 ~ " ) дх2

>dk.

Данный интеграл выражается через функцию Эйри Лг(^), что дает следующую асимптотику волнового поля в окрестности волнового фронта:

Жп (г, г, го, г)■■

ехр(я / 4)

2^2

яг

1

—г 2

д 2 с п (х)

дх2

^ Ю2 - / 2

*3/2 *

х 3/2ю

Р ( г, х )р ( го, х У,

где

У = ехр\(х*г -фг )]]

(г - с^)/

-х д2сп(х*) 2 дх2

1/3"

Проведенные по этим формулам численные расчеты показывают, что для типичных длин внутренних гравитационных волн, характерных для реального океана вдали от полярных зон, учет влияния вращения Земли не вносит существенную поправку в дальние волновые поля и его необходимо делать только для исследования волновых возмущений в высоких широтах Мирового океана [10,11].

Влияние вязкости. Оценим теперь влияние эффектов вязкости на волновую динамику стратифицированных сред. При исследовании внутренних гравитационных волн, например монохроматических с частотой а, с учетом вязкости v необходимо решать дисперсионные соотношения вида [1, 2, 9]

k4 + k П (2к2 + i (N2 sin2 Р - а2 )) v) + 2ik( k„N2 sin p cos p/av + +ik( (N2 cos2 p-a2 )/av + k( = 0,

где kn, kT — нормальные и касательные составляющие волнового вектора соответственно, Р — угол наклона жесткой границы стратифицированного слоя.

Корни дисперсионного уравнения, описывающего пограничные течения, определяют толщину d пограничного слоя:

-,1/2

d « 2v sin ©/[N (sin2 p - sin2 ©)]

где sin2 0 = а /N, зависящую от масштаба Стокса / а и геометрии задачи. Известно, в частности, что при Р = 0 корни дисперсионного соотношения не имеют особенностей. В лабораторных условиях при наклоне дна Р = 15°, N « 1,0 с"1 и а/ N « 0,5 толщина пограничного слоя составляет порядка 2 мм при глубине бассейна 40 см. В условиях реального океана пограничные слои могут иметь толщину до нескольких метров [4, 10, 11]. И в том и в другом случае толщина пограничных слоев составляет несколько процентов от длины полуволны внутренних гравитационных волн для мод с малыми номерами, а поскольку в океане реально возбуждается несколько низших мод, вязкость не оказывает существенного влияния на волновую динамику природных стратифицированных сред [1, 3, 5].

Эффекты сжимаемости. При изучении внутренних гравитационных волн с учетом сжимаемости стратифицированной среды необходимо рассматривать частоту Брента—Вяйсяля N2(z) = - g

в

Ро dz

модифицированном виде:

n2( z) = -gdp - 4 = N2 - 4.

Ро dz c2 c2

где с — скорость звука.

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

исходной системе уравнений гидродинамики (1) уравнение несжима-

1 dp = dp с2 dt dt

шению порядка соответствующего дисперсионного уравнения

емости div U = 0 заменяется на —— -—¡-, что приводит к повы-

а4 - а2( N2 + k 2c2) + к 2c2 N2 sin2 © = 0,

где sin © = а / N.

Однако для длин волн, удовлетворяющих условию N2 А/g << 1 (автоматически это означает / < c2¡g = 200 км), дисперсионное уравнение распадается на два соответствующих звуковым и внутренним волнам уравнения

Nck2

а = kc; а =■

к

где к2 — горизонтальное волновое число. Следовательно, уравнение внутренних волн имеет тот же вид (см. выражение (12)), только вместо N(г) при слабой стратификации необходимо использовать Nc (г) [2, 9].

Рассмотрим основные особенности волновой динамики импульсов акустико-гравитационных волн в стратифицированных средах с учетом их сжимаемости. Акустико-гравитационные волны генерируются при различных атмосферных процессах и играют важную роль в динамике атмосферы Земли и других планет. Мощными импульсными источниками этих волн являются промышленные и ядерные взрывы, извержения вулканов и т. д. Учет сжимаемости важен и при изучении внутренних волн в океане, так как позволяет снять некоторые парадоксы установления поля внутренних гравитационных волн. Действительно, в сжимаемой среде такие волны в определенной степени «похожи» на обычные акустические, правда, в сильно анизотропной среде [4, 9, 10]. Будем исходить из линеаризованной системы уравнений гидродинамики стратифицированной сжимаемой среды

^ ди 1 дР = г .

А) = /у. (22)

дг ду

дЖ др А0 — + — + ЕЧ = /г; дг дг

дg дЖ -Г- + —Р0 +Ро дt дг

(ди 1 + дЦ2 + дЖ ] д$о

дх ду дг

дл

-,-2

д^р -роgW \ + ро

ди 1 | ди2 | дЖ

дх ду дг

Л

= Р(М *,

где g — возмущение плотности; /х, /у, / — компоненты внешней

объемной силы; 5о — источник возмущений плотности; М* — источник объемной скорости.

Характер источника определяется постановкой задачи. Например, взрыв является источником объемной скорости и источником возмущений плотности вследствие нагрева среды. Характерными источниками объемной скорости в атмосфере являются так называемые термики — перегретые области, поднимающиеся вверх и приводящие слои атмосферы в движение. Источники внешней силы возникают при нелинейных взаимодействиях волн различной природы, а также при обтекании тел (области срыва потока при обтекании являются источниками внешних объемных сил) [3, 9].

Систему (22) можно свести к одному уравнению, например, относительно вертикальной компоненты Ж скорости:

ЦЖ = П.

(23)

Здесь

Ц =■

2

д12

дt2

М

дг дг2

2

д2 д:

2

кдх2 ду2 у

(

дt2

+ М - ^С 2

- ,-.-1

П = р

дt2

2

д2 д:

2

дх2 ду2

дt

- Ш-РМ.

(24)

+ ^

-2

Р

д 2М д (д/у д/х

дг2 дt

X

ду дх

М — величина, характеризующая скорость изменения невозмущен-

др

ной плотности с высотой, м = -р —.

дг

Величина N = (мg - g2/с2 ^ определяет частоту Брента—Вяй-

сяля в сжимаемой стратифицированной среде [6-9]. При изучении внутренних волн традиционным является приближение Буссинеска, т. е. в системе (22) пренебрегают градиентом плотности в тех членах,

которые не умножаются на g [1, 6-8]. В этом случае оператор Ь\ сводится к оператору

L2 - „ 2 dt2

д2 (а2 а2 Л (д2 д2Л

dt2 dz2

дх2 ду2

а2 2 22 +Mg - gc ,

Если формально устремить скорость звука с к бесконечности, то получим оператор внутренних волн в несжимаемой среде:

Ьз = -

д2 ( д д2 Л ( д2 д2 Л

дt2

M—+

дz дz2

дх2 ду2

(д2 5t2 + Mg

Л

При использовании приближения Буссинеска для несжимаемой среды оператор Ь3 сводится к оператору

l4 = -

д2 ( д2

дt2

+ + д

дх2 дУ2 дz2

Mg

( д2 + д

J

2

дх ду

Если среда однородная, но сжимаемая, то в однородном поле силы тяжести она обязательно имеет экспоненциальную стратификацию ло = g|c2 [2, 9]. При 0 < л < ло стратифицированная среда оказывается неустойчивой несмотря на убывание плотности с высотой. Действительно, тогда N2 = g (¡-¡о )< 0. Отсюда следует, что в решениях появятся экспоненциально растущие со временем слагаемые. Волновой оператор для естественно стратифицированной среды имеет вид

L5 -

д2_ (_д 2 д12

д2 д2

д2

е2дг дх2 ду2 дz2

gc

-2

д

2

дz 2

Наконец, в пренебрежении гравитационным полем (g = 0) оператор Ь1 превращается в дважды продифференцированный волновой оператор

L6=дё

д2 ( д2

д2 д2

д

2

v c 2дt2 дх 2 ду 2 дz 2 j

Рассмотрим далее гармонический источник возмущений частоты а в среде с экспоненциальной стратификацией р — р0 exp (-mz) (M - const) и постоянной скоростью звука. В этом случае оператор L\ имеет постоянные коэффициенты и его фундаментальное решение G удовлетворяет уравнению

Цв = 5( х)5( у)5( г) 5^).

Тогда решение задачи с конкретным источником (/х, /у, /, 5о, М )

выражается в виде свертки в с функцией П. В частности, в случае точечного мгновенного источника объемной скорости имеем

П = 5'(г )5( х )5( у )5'( г),

т. е. вертикальная компонента Ж скорости выражается через функцию в следующим образом:

д3

Ж =в дЛ2 дг

С помощью подходящей замены переменных удобно избавиться от нечетной производной д/дг в уравнении (23). Действительно, имеем

Ц = Ц7ехр (-/г/2).

Здесь

¿7 =

д

2 (

дt2

--2

дt2

■ + М2/4 -

д

2 Л ( Л2

дг 2

д2 д:

|2 \ ( д 2

дх2 ду 2у

дt2

■ + N2

Л

а ехр(-/г /2) — оператор умножения на функцию.

Проведем одностороннее преобразование Фурье. Спектр фундаментального решения Г запишем в виде

Г(ю, х, у, г ) = ехр (/иг/2) Го(ю, х, у, г ),

где Г0(ю, х, у, г) — фундаментальное решение оператора Гельмголь-ца в анизотропной среде [2, 3, 9]. Тогда

¿8Го = -5( х )5( у )5( г ), (25)

где

д2 д

2

д 2

Ц = (ю2 - N2 + + ю2^-г + ю4с-2-02 М2/4.

1 71 дх2 ду2 дг 2 ^ '

дг 2

-1/2

В координатах х1 = х, у1 = у, г1 = ю (ю2 - N2 ) 2 уравнение (25) принимает вид

(Д1 + ю2 (2 - N2 )-1 ((с-2 - и2/4)) Го =

-1 (ю2 - N2)-1/2 5(хО 5(у1) 5(п), (26)

= -ю

где

а2 а2 а2

Ai =—- +—- + -

ах2 ¿yf arf

2

Используя известное фундаментальное решение обычного оператора Гельмгольца, из уравнения (26) получаем

Го = ~(4жсо)~1 (о2 - N2 )"1/2 R-1 exp (ВД).

Здесь

Ri = R (i - N2 о"2 sin2 0)1/2; R = ( X2 + y2 + z2 )1/2; sin0 = z/R; k2 = о2 (о2 - N2 )-1 (о2с-2 - ¡л2/4) .

Таким образом,

Г(о, х, y, z) = A0 (о) R-1 exp [ik0 (о) R + ¡ z/ 2], (27)

где

Aq(o) =-^

у1(о2 - N2 ) (о2 - Nqq )

ko (о) = с

- ^,-1

(о2 - N:2 )(о2 - N0);

а2 - ы2 N1 = /и2е2 ¡4; Ы2в = N2sin2 в.

Формула (27) описывает поле гармонического источника в сжимаемой стратифицированной среде. Обсудим физический смысл входящих в нее величин. Функция кв(а) — волновое число плоской акустико-гравитационной волны, распространяющейся под углом в к горизонту. Действительно, поле плоской волны н'о = ехр х х \_-iat + ¡кв(а)(хcosв + 2sinв)] удовлетворяет уравнению Ь^о = 0. Зависимость волнового числа кв(а) от угла в обусловлена анизотропией среды. Акустико-гравитационные волны обладают диспер-

а а

сией, а их фазовая скорость Cf =■

ko(o) (о2 - N2)(о2 - NQ)

зависит от частоты (рис. 1).

о

- N 2

Дисперсионные кривые, показанные на рис. 2, имеют две ветви: акустическую (при ш > Ыл) и гравитационную (при 0 < ш < И). Фазовые и групповые скорости (акустическая ветвь), стремятся к скорости звука с при ш ^го. С ростом со волны, описываемые акустической ветвью, переходят в обычные звуковые. Нижней критической частотой этих волн является Ил — наименьшая частота распространения акустических волн. При частотах от N до Ыл (верхний диапазон запирания) волновое число кв(ш) чисто мнимое. В этом диапазоне нет распространяющихся плоских волн, а только экспоненциально затухающие при удалении от источника возмущений.

Рис. 1. Дисперсионные кривые Рис. 2. Фазовые (сплошные ли-о(к) акустико-гравитационных нии) и групповые (пунктирные волн линии) скорости акустико-грави-

тационных волн при в = 0 (/ ) и в > 0 (2)

Кривые, изображенные на рис. 2, при 0 < ш < N описывают гравитационную ветвь акустико-гравитационных волн. Верхняя критическая частота этой ветви равна N, нижняя — зависит от угла в. Если волна распространяется строго горизонтально (в = 0), то нижняя критическая частота отсутствует. Фазовая и групповая скорости горизонтально распространяющейся волны при ш ^ 0 стремятся к значению С0 = Ыс/Ыл < с. Если в Ф 0, то групповая скорость имеет максимум на некоторой частоте ш*(в), отличной от нуля.

Картины волновых фронтов поля, возбуждаемого гармоническим источником, сильно различаются в случаях ш > Ыл и ш < N (при N < ш < NA поле синфазно и экспоненциально убывает при удалении от источника). При ш > N л фронты имеют вид эллипсоидов вращения (рис. 3, а), причем с ростом частоты их эксцентриситет убывает, а в пределе при ш ^го они превращаются в сферы, как для обычных

звуковых волн. При о ^ Na эллипсы становятся все более сплюснутыми в вертикальном направлении. При о < N фронты представляют собой семейство однополостных гиперболоидов вращения с общей асимптотой под углом во = arcsin(o /N), зависящим от частоты о (рис. 3, б). При \в\ > во поле синфазно и экспоненциально убывает при удалении от источника. Внутри конуса (при \\ < во) поле осциллирующее, гребни волн «удаляются» от конуса \\ = во. На конусе поле гармонического точечного источника имеет корневую особен-

(4л)"1

ность Ав(о) =-, v ' . Если размеры ис-

NR^(о2 - N2 ) (sin2 во - sin2 в)

точника конечные, то особенность не возникает, но поле в основном сосредоточено в узкой области вблизи этого конуса, имеющей поперечные размеры, близкие к размерам источника.

Рис. 3. Волновые фронты для акустической (а) и гравитационной (б) ветвей акустико-гравитационных волн

Обсудим физический смысл величин Ne, N, Na. Частота Na определяется произведением малого параметра / на большой параметр с. Оценки показывают, что как для океана, так и для атмосферы частота Na близка к частоте N Брента—Вяйсяля: для атмосферы N « 1,9 • 10-2c-1, Na « 2,1 • 10-2 c-1, для океана N ~10-2c-1, NA ~ ~ 1,5 • 10"2 c-1 [2, 5, 9]. Из этих оценок следует, что при анализе дисперсии ке(ю) на частотах, близких к N, необходимо учитывать сжимаемость среды. Отметим, что всегда NA > N, NA - N2 = (/icj2 - g/с) =

= с2 (/- 2/0 )2/4. Таким образом, только при / = 2/ отсутствует

верхний диапазон запирания и получается единая дисперсионная кривая (рис. 4, а). При / = // частота Брента—Вяйсяля обращается в нуль и пропадает нижняя ветвь дисперсионной кривой (рис. 4, б), в отличие от предыдущего случая эта кривая не зависит от в.

а б

Рис. 4. Дисперсионные кривые О (к) для р = 2f0 (а) и для р = р0(б)

При рассмотрении низкочастотных внутренних гравитационных волн в атмосфере обычно учитывают вращение Земли [3, 5, 9]. Однако это приводит к изменению лишь параметра N в в уравнении (27):

Ng = N2 sin2 в + JZ cos2 в,

где Jz — z-компонента вектора угловой скорости вращения Земли.

Следовательно, вращение Земли заметно сказывается на поле внутренних волн лишь при очень малых в (в — Jz/N < 10 3 рад). Влияние вращения Земли приводит к тому, что зависимость поля акустико-гравитационных волн от R и t на горизонтальной плоскости и вблизи нее такая же, как для поля на конусе с углом в « Jz/N, рассчитанном без учета вращения. Поэтому вращение Земли можно учесть, например, изменив угловую зависимость Ng.

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

, , ч и \а>2 — Ng

кв(о) = 2 J —2-2.

2 V N2 — о

Очевидно, что при таком законе дисперсии пропадает акустическая ветвь дисперсионной кривой, а гравитационная ветвь с критиче-

скими частотами Ыд и N качественно (при а «N количественно) не меняется. Физически абсурдным также является приближение Буссинеска в сжимаемой среде [2, 9]. В этом случае Ыа = 0 и, следовательно, Ыа < N, что невозможно в устойчивых средах.

Дисперсионная кривая имеет две ветви (рис. 5) и качественно совпадает с дисперсионной кривой, изображенной на рис. 1, однако верхняя и нижняя границы диапазона запирания искажены. Поверхности равных фаз при а < N принимают форму двухполостных гиперболоидов вращения (рис. 6). Отметим, что в сжимаемой среде ( = ¿и0) волновое число кд не зависит от угла в (изотропный случай):

кв(а) = с "У а2 - N1 ,

где NA = (2с), и поэтому анизотропия среды проявляется лишь в появлении множителя ехр [gz|(2с2 )] в уравнении (27).

Рис. 5. Гравитационная ветвь дис- Рис. 6. Гравитационная ветвь вол-

персионных кривых 0)(к) аку- новых фронтов акустико-грави-

стико-гравитационных волн в при- тационных волн в приближении

ближении Буссинеска в сжимаемой Буссинеска в сжимаемой страти-

стратифицированной среде фицированной среде

В заключение отметим, что при изучении внутренних гравитационных волн в океане обычно используют модель несжимаемой среды и приближение Буссинеска. В этом приближении уравнение (26) сводится не к уравнению Гельмгольца, а к уравнению Лапласа, так что волнового распространения нет и среда движется синфазно, как при колебании тел в несжимаемой однородной среде [2]. Однако при а < N между обла-

стями вт и |в > вт наблюдается сдвиг фаз на ж/2. Еще одно отличие

от случая нестратифицированной жидкости — корневая особенность амплитуды на границе между областями (в = вт), поэтому для перехода к несжимаемой среде в приближении Буссинеска достаточно в уравнении (27) принять kв(ю) = 0.

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

Работа выполнена при поддержке РФФИ (проекты № 11-01-00335а, № 13-05-00151а).

ЛИТЕРАТУРА

1. Miropol'skiy Yu.Z., Shishkina O.V. Dynamics of internal gravity waves in the ocean. Boston, Kluwer Academic Publishers, 2001, 421 p.

2. Габов С.А. Новые задачи математической теории волн. Москва, Наука, 1998, 448 с.

3. Pedlosky J. Waves in the Ocean and Atmosphere: Introduction to Wave Dynamics. Berlin-Heidelberg, Springer, 2010, 260 p.

4. Коняев К.В., Сабинин К.Д. Волны внутри океана. Санкт-Петербург, Гидро-метеоиздат, 1992, 272 с.

5. Sutherland B.R. Internal gravity waves. Cambridge, Cambridge University Press, 2010, 394 p.

6. Булатов В.В., Владимиров Ю.В. Динамика негармонических волновых пакетов в стратифицированных средах. Москва, Наука, 2010, 470 с.

7. Bulatov V.V., Vladimirov Yu.V. Wave Dynamics of Stratified Mediums. Москва, Наука, 2012, 584 с.

8. Булатов В.В., Владимиров Ю.В. Волновые движения стратифицированных сред: теория и приложения. Saarbrucken, Palmarium Academic Publishing, 2012, 577 с.

9. Бреховских Л.М., Годин О.А. Акустика неоднородных сред. Т. 1: Основы теории отражения и распространения звука. Москва, Наука, 2007, 443 с.; Т. 2: Звуковые поля в слоистых и трехмерно-неоднородных средах. Москва, Наука, 2009, 426 с.

10. Поверхностные и внутренние волны в арктических морях. Санкт-Петербург, Гидрометеоиздат, 2002, 364 с.

11. Morozov E.G. Internal Tides. Global Field of Internal Tides and Mixing Caused by Internal Tides. Waves in geophysical fluids. Ed. by J. Grue and K. Trulsen. Wein, New York, Springer, 2006, pp. 271-332.

Статья поступила в редакцию 21.02.2013

Ссылку на эту статью просим оформлять следующим образом: В.В. Булатов, Ю.В. Владимиров. Моделирование волновой динамики стратифицированных сред с учетом нелинейности, вязкости, вращения земли и сжимаемости. Инженерный журнал: наука и инновации, 2013, вып. 2. URL: http://engjournal.ru/catalog/appmath/hidden/606.html

Булатов Виталий Васильевич — старший научный сотрудник ИПМех им. А.Ю. Ишлинского РАН; д-р физ.-мат. наук, д-р экон. наук, профессор, автор более 250 научных работ, в том числе 12 монографий; сфера научных интересов: математическое моделирование, гидрофизика, волновая динамика, экономическая теория. e-mail: bulatov@index-xx.ru

Владимиров Юрий Владимирович — старший научный сотрудник ИПМех РАН им. А.Ю. Ишлинского; канд. физ.-мат. наук; автор более 100 научных работ, в том числе 5 монографий; сфера научных интересов: математическое моделирование, гидрофизика, волновая динамика, прикладная математика. e-mail: vladimyura@yandex.ru