Уравнения Буссинеска на вращающейся сфере Текст научной статьи по специальности «Математика»

Вычислительные технологии Том 16, № 6, 2011

Уравнения Буссинеска на вращающейся сфере*

З.И. Федотова, Г. С. Хакимзянов Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: zf@ict .nsc.ru, khak@ict .nsc.ru

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

Ключевые слова: поверхностные волны на воде, уравнения мелкой воды на сфере, уравнения Буссинеска, нелинейно-дисперсионные уравнения.

Введение

В работе [1] были представлены полные нелинейно-дисперсионные (НЛД) модели на поверхности вращающейся притягивающей сферы, предназначенные для математического моделирования продолжительного по времени распространения трансокеанических волн. Хотя ввиду уменьшения размерности задачи полная НЛД-модель по сравнению с трехмерными уравнениями Эйлера является более экономичной, для проведения практических расчетов вполне целесообразно ее дальнейшее упрощение. Это касается волновых процессов вдали от берега, где течение не всегда можно полагать полностью линейным, однако нелинейностью в дисперсионных членах можно пренебречь без потери точности моделирования. Такие слабо нелинейные дисперсионные модели на плоскости хорошо известны. Это модели Перегрина [2], Дорфмана-Яговдика [3], IIйогу [4] и др.

В настоящей работе из полной НЛД-модели выводится модель Буссинеска на сфере. Строго говоря, это слабо нелинейная слабо дисперсионная модель, описывающая распространение волн на поверхности тонкого слоя жидкости на вращающейся притягивающей сфере, полученная из полных уравнений НЛД-модели мелкой воды [1] в предположении, что выполнено соотношение а = 0(^2), где а и ^ — малые параметры, характеризующие меру нелинейности процесса и частотную дисперсию. Эти параметры появляются в полной НЛД-модели в процессе ее вывода при соответствующем масштабировании трехмерных уравнений Эйлера.

Путем перехода к линейным скоростям и линейным координатам установлена связь между полученной НЛД-моделью на сфере и НЛД-моделями на плоскости, частными случаями которых являются модели Перегрина и Дорфмана — Яговдика. Выполнена модификация полученной модели с целью построения на ее основе эффективных чис-

*Работа выполнена при финансовой поддержке РФФИ (гранты № 10-05-91052-НЦНИа и 09-05-00294а), Программы государственной поддержки научных школ РФ (НШ-6068.2010.9) и Проекта IV.31.2.1. программы фундаментальных исследований СО РАН.

ленных алгоритмов, опыт разработки которых был приобретен при решении аналогичных задач на плоскости [5, 6].

Настоящую статью можно рассматривать как вторую часть работы [1]. Для удобства сохранены основные обозначения.

1. Постановка задачи для уравнений Эйлера и их обезразмеривание

Введем сферическую систему координат OA0r с началом в центре Земли, вращающейся с постоянной угловой скоростью Q, вектор которой направлен к Северн ому полюсу; A — долгота, отсчитываемая к востоку от некоторого меридиана (0 < A < 2п), в = п/2 — — дополнение до широты ^ (—п/2 < ^ < п/2). На поверхности Земли, аппроксимируемой сферой радиуса R, рассмотрим слой жидкости, толщина которого по сравнению с радиусом R предполагается малой, поэтому величина g = |g|, где g — направленная к центру Земли сила ньютоновского притяжения, и плотность воды р принимаются постоянными по всему слою, р = 1,

Введенная сферическая система координат связана с неподвижной декартовой системой координат Ож^ж3 с осью Ож3, проходящей по оси вращения Земли, и координатной плоскостью Ож1ж2, совпадающей с экваториальной. Ориентация декартовой системы такова, что Северный полюс соответствует координате ж3 = R, Сферические и декартовы координаты связаны выражениями

ж1 = r cos(A + Ш) sin в, ж2 = r sin(A + Ш) sin в, ж3 = r cos в. (1)

Обозначим через J = —r2 sin в якобиан преобразования (1), а через g00, g«0, g«e (а, в = 1, 2, 3) — ковариантные компоненты метрического тензора этого преобразования:

g00 = 1 + ^V sin2 в, gio = Hr2 sin2 в, g20 = g30 = 0, (2)

gil = r2 sin2 в, g22 = r2, g33 = 1, gae = 0 при а = в. (3)

Тогда уравнения, описывающие течение идеальной несжимаемой жидкости, можно записать в следующем виде [1]:

(Jv1)A + (Jv2), + (Jw)r = 0, (4)

vt + (u ■ V)v + w vr + Vp = r, (5)

Wt + u ■ Vw + WWr + Pr = —g + r3, (6)

a Ü 3/1 2\T

где w — радиальная, или вертикальная , составляющая скорости v , u = (v1,v2)T — вектор "горизонтальной" составляющей скорости, V = (d/dA, д/дв), r = (r1,r2)T,

(П + t;1)2 ддп (П + у1)2 ддп , (v2f дд22

Г1 = 0' = Гз = 2 дг~ ~2 дт7' (7)

v = (v1, v2)T, При этом ковариантные компоненты скорости vY (7 = 1, 2, 3) выражаются через контравариантные по формулам

VY = gY0 + gie^,

которые е учетом равенств (2), (3) приводят к выражениям

VI = 910 + 9иУ1, У2 = 922У2, Уз = V3 = т. (8)

Система уравнений (4)-(6) дополняется начальными и краевыми условиями. Не рассматривая постановку этих условий [1], отметим лишь, что поскольку было сделано предположение о шарообразной форме Земли, то естественным следствием этого является условие сферичности поверхности воды в состоянии покоя. Динамическое условие р = 0 на такой поверхности означает, что вектор Ур должен иметь только радиальную составляющую. Отсюда следует, что центробежной силой можно пренебречь, что и было сделано в [1]. При таком допущении первая формула в (2) примет вид д00 = 1, Для вывода уравнений мелкой воды на сфере введем характерные масштабы и в уравнениях (4)-(6) перейдем к безразмерным величинам. Пусть Ь и к0 — характерные масштабы в горизонтальном и вертикальном направлениях соответственно, а0 — характерная амплитуда волны, а = а0/к0, С величиной Ь свяжем характерный горизонтальный масштаб А0, измеренный в радианах:

А -Ь-£ К ц

где е = к0/К ^ 1, ^ = к0/Ь. Масштаб времени ¿0 определим как

Ь

¿0

^/дКо'

Тогда характерный масштаб угловой скорости распространения волны будет определяться формулой

А0 е л/дК~о л/Фо

Шо = — =--— = —Б—>

¿0 ^ Ь К

а за масштаб "горизонтальной" угловой скорости частиц жидкости в волне естественно принять величину аш0.

В соответствии с введенными масштабами определим безразмерные переменные:

-\-e~h_s— Н_г1-г7Ч п

А = —, 0 = —, п=—, 8 = —, Н = —, г] =—, ¿= —, \1 —

Ао Ао ко ко ко ао ¿о а^о

я V13 и> V

* = -' п ™= -=,

а^о аКл/дк0 а^удк0 дко

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

2. Полная НЛД-модель на вращающейся сфере

Для вывода уравнений мелкой воды определим искомые величины: Н = п + к — полная глубина слоя жидкости и с = с(А,9,1) — вектор скорости в приближенной модели. Здесь, как и в работе [1], с = и, где компоненты вектора и согласно условиям вывода модели не зависят от радиальной координаты.

Основное предположение, используемое в [1] при выводе НЛД-уравнений мелкой воды на сфере, состоит в том, что в тонком по сравнению с радиусом Земли слое воды изменением компонент метрического тензора по глубине можно пренебречь и в формулах (5), (6) использовать приближенные величины

#10 = ^Д2 sin2 0, gil = R2 sin2 0, g22

R2.

Тогда и для якобиана можно брать его приближенное значение

J = - R2 sin 0.

(9)

(10)

Дополнительное исследование, проведенное авторами, показывает, что такое упрощение приводит к исчезновению из уравнений Эйлера слагаемых порядка е. Вывод НЛД-уравнений осуществляется по схеме, примененной в [1], с той лишь разницей, что уравнения Эйлера использованы в обезразмеренном виде и, следовательно, при некоторых слагаемых появляются коэффициенты, зависящие от параметров а и р, С учетом этой поправки полученная НЛД-модель на вращающейся сфере примет следующий вид (для упрощения черточки над безразмерными переменными опущены):

vt + a(c ■ V)v + Vn — r = p2

Ht + aV ■ (Hc) = 0,

(и) (12)

Здесь

H = an + h, c =(c1,c2)T

v

(v1,v2)T, v1 =(fi + c1) sin2 0, v2 = c2

£

(?"i) ?"2)T) = r2 = a- (m + с1) c1 sin 6 eos 0,

v-c = c1;

(Jofh Jn

J0 = — sin I

,2 „ ^ 2г ~ d

а от

аргумент 0 в тригонометрических функциях не масштабирован. Система уравнений

(11), (12) служит для определения полной глубины Н и контрвариантных компонент св (в = 1, 2) вектора скор ости с.

Модель гидродинамики поверхностных волн, описываемая системой уравнений (11),

(12), принадлежит к классу так называемых полных НЛД-моделей второго гидродинамического приближения. Из нее путем тех или иных упрощений могут быть получены многочисленные варианты промежуточных НЛД-моделей, "расположенных" между полученной полной НЛД-моделью и классической моделью мелкой воды первого приближения.

r

3. Формальный вывод уравнений Буссинеска

В классической теории Буссинеска предполагается а = 0(р2), Это означает, что оставаясь в пределах точности НЛД-модели, членами порядка ар2 можно пренебречь. Такая процедура существенно упрощает уравнения НЛД-модели,

Обозначим через Q выражение в квадратных скобках правой части уравнения (12) и выделим в Q члены, не зависящие от а, обозначив остальные как 0(a), Нетрудно видеть, что выражение для R1 можно записать в виде

Ri = (V- c)t + 0(a). (13)

Обратив внимание на то, что из уравнения неразрывности (11) вытекает соотношение

ht = 0(a), (14)

получим следующее выражение для R2:

R2 = -htt + ct-Vh + 0(a). (15)

a

Применяя эти подстановки, дифференцирование и приведение подобных, получим

Q = ^V (V • hct) - у V (V • c)t + AVhtt + o(a). (16)

Таким образом, отбрасывая члены порядка 0(а^2) и возвращаясь к размерным переменным, уравнения НЛД-модели в приближении Буссипеска, вытекающие из (11), (12), (16), можно записать в следующем виде:

Ht + V ■ (Hc) = 0, (17)

h h2 h vt + (с • V)v + gVr] - г = -V (V • hct) - —V (V • c)t + -Vhtt, (18)

где

c =(c1,c2)T, v = (Vi,V2)T, vi =(n + c^ R2 sin2 0, V2 = R2c2, (19)

r =(ri,r2)T, ri = 0, г2 = (2П + c1) c1R2 sin0 cos 0. (20)

Далее будет установлена связь с известными моделями на плоскости.

4. Уравнения Буссинеска для линейной скорости

При численном решении практических задач о распространении волн на поверхности сферы удобнее использовать линейные компоненты и и V вектор а с:

u = Rc1 sin 0, v = Rc2.

(21)

Из формул (17), (18) следует, что для зависимых переменных Н, и, V уравнения НЛД-модели будут выглядеть как

Ht +

1

R sin 0

(Hu) л + (Hv sin 0)g

0,

(22)

uux vv± I 9Vx R sind R R sind

ctg 9 ~R~

(2QR sin 0 + u) v+

h

2R sin 0

((V • Ma " \ (V • c)tA + htt^j

(23)

+

vt

uvx We gr]e rtglíonTj л i \ i R^Te + -R + -R = —{2QR sm 0 +

h_

2R

^(V • hct)e - ^ (V • c)te + fc^)

(24)

где

V • hct

1

({hut)\ + (/щ sin 6% ), V • с = д ^ (u\ + sin 6% ).

Rsin 9V 4 ' Rsin 9

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

Ht

1

(Hu)\ + (Hv sin 6)e

Ut

R sin 6

uu\ vue дщ uv ctg 6

vt

R sin 6 R R sin 6 R

uv\ vve grqe u2 ctg 6 Rsm9 + ~R+ ~

0,

- fv, fu,

R R

где f = 2П cos 9.

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

Vt +

1

R sin 6

ut

(hu)x + (hv sin 6)e

, gvx

= -ht,

vt

R sin 6

gm R

= -fv, fu.

5. Модель Буссинеска для линейных координат

Чтобы установить связь между полученными уравнениями Буссинеска на вращающейся сфере и известными слабо нелинейными дисперсионными моделями на плоскости, перейдем к координатам на поверхности Земли, Для этого в некоторой окрестности точки (Ао,0о) рассмотрим преобразование координат:

x = R sin 60 (А - А0), y = -R (6 - 60).

(25)

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

u = X = Rc1 sin 60, v = y = -Rc2

(26)

Из формул (25), (26) следует, что

d

u c2 = --, — = Rsm90 — , — = -R—, R sin 60' R' dA dX d6 dy}

Vc = V-U-

v ctg 6 R '

__ y du dv T T . . V-U = — + —, U = (u,v). dx dy

i

c

Дополнительно потребуются следующие выражения:

, , ( n vtx ctg в (hvt)x ctgd

(V • c)tx = (V • U)te---—, (V • hct)x = (V • hUt)x -

R

(V ■ c)ty = (V ■ U)ty -

R

Vty Ctg в

Vt

(V ■ hct)y = (V ■ hUt)y -

ty R R2 sin2 в' (hvt)y ctg в hvt

y t y R R2 sin2 в

Переходя в системе уравнений Буссинеска (17), (18) к производным по x,y и скорости с компонентами (26) и учитывая приведенные выше выражения для дисперсионных членов, получим следующую систему уравнений:

Hv ctg в

Ht + (Hu)* + (Hv)y

Du + a2 gnx = afv +

R

2uv ctg в Л

+a2

h

' , /г (fe^ctgg fe^ctgg' (V • hVt)x ~ ~ (V • U)te + /te----+-——

(27)

(28)

Dv

gny

1 1 u2ctg в --fu--ö-^--\

h 2

h

(V • hVt)y - - (V • U)ty + htty

a a*

hvty ctg 9

Ж

где

a

sin в0

D

d_

dt

d

u-

R

(hvt)y ctg в

R dy'

2hvt

3R2 sin2 в

(29)

1+O(í). Поскольку раеематри-

sin в dt dx

Предполагая, что величина ö = в — в0 мала, получим a ваемые НЛД-уравнения предназначены для моделирования волн цунами, то намеренно исключим из рассмотрения приполярные области Земли, что гарантирует ограничен-

ctg в u, v, H

производных, В силу этих условий, пренебрегая в уравнениях (27)-(29) членами, имеющими порядок O(ö) или O (1/R), получим плановые НЛД-уравнения на плоскости. Пренебрегая силой Кориолиса, запишем их в виде

Ut + (U -V) U + gVn

Ht + V ■ HU = 0,

h h2 h -V (V • hVt)-yV(V-U)t + -Vhtt.

(30)

(31)

Данные уравнения легко выводятся из НЛД-модели на плоскости [7], если к полученной в этой работе полной плановой модели применить процедуру перехода к слабо нелинейной модели типа модели Буссинеска по аналогии с настоящей работой. Отметим, что в случае неподвижного дна уравнения (30), (31) переходят в известную модель Перегрина [2]:

Н + V ■ Ни = 0,

Ut + (U -V) U + gVn

hd_

2 dt

V

(V . Ш) - ^ (V • U)

Кроме того, если предположить, что форма донной поверхности описывается формулой Л,(х,у,£) = Лю(х,у) + аЛ,(х,у,£), то из уравнений (30), (31) следует одна из моделей, описанных в работе [3].

f

2

6. Модифицированные уравнения Буссинеска

Предлагаемый численный алгоритм основан на модифицированной системе уравнений Буссинеска, которая легко расщепляется на две системы, одна из которых содержит уравнения первого порядка, другая — уравнения второго порядка эллиптического типа. Модификации подвергается уравнение движения (18), Учитывая равенства а = 0( и (13)—(15), получим, что в уравнении (12)

0 = £ (V • с*) + ^У (V • с*) + ^У (с, • Ш) + £ ^ + 0(/12). (32)

2 3 2 2 а

Преобразуем первые три слагаемых в правой части (32), Из уравнений движения (12) следуют равенства

^М + Па = 0(^2), ь2,г + Пв = 0(^2). (33)

Из них после перекрестного дифференцирования и вычитания получаются соотношение

Щ,ьв - Щ,ЬА = 0(^2) (34)

и равенства

4 = = --Й + 0(^2), с2 = г*>4 = + О(^). (35)

"о "о

Следовательно, для первого члена правой части (32) имеем выражение

^ (V • с4) = (с1х + = V2?? + (V), (36)

где через V2 обозначен оператор Лапласа

= ^ + (37)

"о "0

Третье слагаемое в правой части (32) с учетом равенств (35) примет вид

^У (с, • УЬ) = + ^ + ¿V). (38)

2 ^ у 2 V

Рассмотрим, наконец, второе слагаемое в правой части равенства (32):

1x2 V ГУ • с ) - 1x2 V /V + 1М2Л - ]г\(+ С учетом выражений (34) и (35) получим следующие цепочки равенств:

, , У1,4АА . ("о^2,4А)в УМАА . ("оУМв )в , гл! 2\ ^-,2 , 2\

• с0а = +-т—- = —¡2~ +-т—- + ) = V г;м + ),

"о "о "о "о

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

V (V ■ е4) = (V2(v - уо))4 - С + 0(А (39)

где

е = «1.с2)т о = о, <, = » (:$*)+„„ (.у »о=(Пя0пав). («)

С учетом равенств (14), (36), (38), (39) перепишем выражение (32):

где

Ч = - 2 V«. V2, - ¥< - V + ЧЛ) + - ^ + 0(„2).

Поэтому уравнение движения (12) переходит в следующее модифицированное уравнение модели Буееинеека:

V* + а(с • + У/7 - г = /х2 - + + 0(/х4). (41)

Отбрасывая в нем члены порядка 0(р4) и переходя к размерным величинам, получим модифицированную модель Буееинеека, у которой уравнение неразрывности по-прежнему имеет вид (17), а уравнения движения записываются в виде

V* + (с • + дУг] - г = - + Ч, (42)

где

q = (Qi,Q2)T = -g'-VhV2r] - g'-v ( ''л7Л2„ + - a'-^C + тУК,

2 \i?2sin26 R2 )~9 3C+2 T Л „ „ 1 , Л . х ( HR2 sin2 0

С = «1,<2) , Ci = 0, (2 = - р2 . 2д (2?7лл ctge + щ), Vo-. п

R2 sin2 0 0

V2

2 _ Ухх (щ sin 0), V ^ i?2sin20+ R2 sin в • l4dj

Остальные величины в (42) определяются согласно формулам (19), (20),

Таким образом, уравнения (17), (42) представляют новый вариант уравнений Буееинеека на сфере. Отметим, что преобразования выполнены на основе следствий из НЛД-уравнений без привлечения каких бы то ни было дополнительных ограничений на характер течения,

uv

получим следующий вид уравнений:

(HR sin 0)t + (Hu^ + (Hv sin 0)e = 0, (44)

d (h,2 \

R sin 6 щ + uu\ + vu$ sin 9 + gr¡\ = — (2ГШ sin в + u) v eos в + — ( —V2 (и i? sin 0) J + qi,

d / h2 R sin 0 \

i? sin 6 vt+uv\+vvo sin 6+grie sin в = (2ПЛ sin в + u) u eos в + — ( -V2v ) +q2 siné1.

di \ 3 J

Для построения численного алгоритма удобно рассматривать уравнения движения в форме

Решение системы (44)-(47) можно получить с помощью двухслойной разностной схемы, состоящей из двух шагов: сначала на каждом временном слое решаются уравнения первого порядка (44)-(46) и находятся величины H, U, V, а затем для определения компонент скорости u, v используются эллиптические уравнения (47) с вычисленными на первом шаге правыми частями.

Заключение

В настоящей работе получены уравнения Буссинеска на вращающейся притягивающей сфере и выполнена их модификация с ориентацией на создание вычислительных алгоритмов, которые можно было бы применить для моделирования распространения волн цунами. Преобразование модели было выполнено благодаря соотношениям (14), (34), (35), Это позволило получить при расщеплении уравнений движения оператор Лапласа, выполняющий роль регуляризатора численного решения [5, 6], Отметим, что как при выводе модели, так и при ее модификации не вводилось предположение о потенциальности течения.

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

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

[1] Федотова З.И., Хакимзянов Г.С. Нелинейно-дисперсионные уравнения мелкой воды на вращающейся сфере // Вычисл. технологии. 2010. Т. 15, № 3. С. 135-145.

[2] Peregrine D.H. Long waves on a beach // J. Fluid. Mech. 1967. Vol. 27, pt. 4. P. 815-827.

[3] Дорфман A.A., Яговдик Г.И. Уравнения приближенной нелинейно-дисперсионной теории длинных гравитационных волн, возбуждаемых перемещениями дна и распространяющихся в бассейне переменной глубины // Численные методы механики сплошной среды. 1977. Т. 8, № 1. С. 36-48.

[4] Nwogu О. Alternative form of Boussinesq equations for near shore wave propagation //J. Waterway, Port, Coastal and Ocean Engng. 1993. Vol. 119. P. 618-638.

[51 Fedotova Z.I., Pashkova V.Yu. Methods of construction and the analysis of difference schemes for nonlinear dispersive models of wave hydrodynamics // Russ. J. Numer. Anal. Math. Modelling. 1997. Vol. 12, No. 2. P. 127-149.

Ut = -uu\ — vue sin 9 — дщ — (2QR sin 9 + u) v cos 9 + q\, V = —uv\ — vve sin 9 — д'Цв sin 9 + (2QR sin 9 + u) u cos 9 + q2 sin 9,

(45)

(46)

(47)

[61 Chubarov L.B., Fedotova Z.I., Shkuropatskii D.A. Investigation of computational models of long surface waves in the problem of interaction of a solitary wave with a conic island // Ibid. 1998. Vol. 13, No. 4. P. 289-306.

[71 Федотова З.И., Хакимзянов Г.С. Нелинейно-дисперсионные уравнения мелкой воды на нестационарном дне // Вычисл. технологии. 2008. Т. 13, № 4. С. 114-126.

Поступила в редакцию 17 марта 2011 г.