О свойствах волн Кельвина на сетке Текст научной статьи по специальности «Математика»

Вычислительные технологии Том 13, № 1, 2008

О свойствах волн Кельвина на сетке*

С. В. Смирнов Институт автоматики и процессов управления ДВО РАН, Владивосток, Россия e-mail: smirnoff@iacp.dvo.ru

The effects of a space grid approximation on free Kelvin waves are investigated using the /-plane shallow-water equations with the Arakawa B-grids and the free-slip boundary condition on a wall.

Введение

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

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-01-96020) и Президиума РАН (программа № 14).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

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

1. Волны Кельвина

Запишем линейную систему уравнений мелкой воды [3] в декартовой системе координат (х, у, г) с направленной вверх осью г:

дьп* - ¡у* = -ддхг]*, (1)

дгь* + ¡п* = —дду п*, (2)

дгп* + Н (дхп* + ду у*) = 0. (3)

Здесь и далее £ — время; п* и у* — компоненты вектора скорости по направлениям х и у; д — ускорение свободного падения; п* — превышение уровня жидкости над его невозмущенным положением; д5Г — частная производная функции Г по переменной е. Пусть жидкость находится в бассейне с одной прямой вертикальной стенкой и глубиной Н > 0 при х < 0. Значения глубины Н и параметра Кориолиса ¡ полагаем постоянными. Пусть ¡ > 0. Запишем условие отсутствия потока через стенку:

п* |х=о = 0. (4)

Известно, что для системы (1)-(3) с краевым условием (4) решением типа волны, захваченной вращением у вертикальной стенки, является волна Кельвина [5]:

(у*,п*)=(УН О ПА ехр(¡Н + £ - гк*у), п* = 0, х < 0, (5)

где шк — безразмерная частота; к* — волновое число, к* = 2п/Л*; Л* — длина волны; Па — множитель. Амплитуда волны Кельвина экспоненциально убывает при удалении от берега. Это убывание характеризуется горизонтальным масштабом

Л* = (6)

который называется радиусом деформации Россби. Фазовые скорости Ск волн Кельвина совпадают с групповыми скоростями,

Ск = ¡Щт, шк = к^л/дН- (7)

Отметим, что в стационарном случае система (1)-(3) переходит в вырожденную систему уравнений геострофики:

-!'уд = -ддх%, ¡пд = -дду %, (8)

дхпд + ду Уд = 0; (9)

краевое условие (4) можно представить в следующем виде:

ду Пд 1х=о = 0- (10)

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

2. Разностные уравнения и краевые условия

Пусть система уравнений (1)-(3) решается разностным методом на сетке типа В (рис. 1, а), по классификации Аракавы [2]. Сеточные переменные размещаются на прямоугольных подсетках с координатами узлов

хт = Ах(т — М), уп = Ау п, Ах = г А, Ау = А, А > 0, где т и п — индексы узлов,

0 < г < 1.

(11) (12)

В узлах с целыми индексами расположены компоненты скорости и и V, в узлах с "полуцелыми" индексами — п. Будем использовать обозначение Fmn = ^(хт,уп). Запишем пространственную разностную аппроксимацию уравнений (1)-(3) [2]:

д4и — /V = — ,

д^ + /и = —д5уГ]х, дп + Н (^хйу + ¿уух) = 0, где разностные операторы определены следующим образом:

(х, у) = Ах-1 ^ (х + Ах/2, у) — F (х — Ах/2, у)), 5уF(х, у) = Ау-1 (F (х, у + Ау/2) — F (х, у — Ау/2)), ^(х, у) = 2 (F (х + Ах/2, у) + F (х — Ах/2, у)), —у 1

^(х, у) = 2 ^ (х, у + Ау/2) + F (х, у — Ау/2)).

(13)

(14)

(15)

(16)

(17)

(18)

(19)

На твердой границе находятся узлы п (рис. 1, б). В фиктивных узлах, расположенных на расстоянии Ах/2 от границы, положим

иМ,п = —иМ—1 ,п, ^М.п = ^М-1,1

(20)

■ V и,\ ■ V ,

■ и77 , ■

Ах

<-5—>

Рис. 1. Расположение узлов и границы: а — размещение узлов и, V, п на сетке В; б — граница в узлах п

Такие условия применяют для реализации "свободного скольжения" на стенке [4]. Всем узлам на границе припишем уравнение (15), которое с учетом (20) принимает вид

дг'Пы-i/2,n+i/2 + H (Syv - 2 _1>n+1/2 = 0. (21)

Таким образом, можно считать, что краевые условия описываются уравнением (21).

Решения системы (13)—(15) с краевым условием (21) будем искать в виде распространяющихся вдоль берега сеточных захваченных волн:

(u,v,n)m,ra = (uo, vo, no)Cm-M exp (¿/wí — ¿kn), m < M, (22)

ICI > 1, (23)

k = k*A, (24)

0 < k < п. (25)

Здесь k — безразмерное волновое число. Сеточную волну называют двухшаговой, когда Л* = 2A и, следовательно, k = п. Если решение вида (22) при неограниченном измельчении шага сетки стремится к решению (5), (7) исходной дифференциальной задачи, такое решение будем называть "сеточной" волной Кельвина. Подставив (22) в (13)—(15), получим:

k

¿uor^/wA — vorv7/A + no (C — 1) g cos - = 0, (26)

k

uov7/ A + zvoVC/wA — ¿no (1 + C) g sin - = 0, (27)

k k r-

uoH (C — 1) cos - — iVorH (1 + C) sin - + znorVC/wA = 0. (28)

Запишем три вспомогательных соотношения:

k k\ / k k\ (1 — C) cos - + rw (1 + C) sin - J uo + ¿ Í w (1 — C) cos - + r (1 + C) sin - J vo = 0, (29)

f k k\

VC/Ar (1 — w2) vo + ( (1 — C ) cos - + rw (1 + C ) sin - J gno = 0, (30)

¿ ((1 — C2) R2 sin k + - rC w) rvo + (R2 (1 — C)2 (1 + cos k) + - C w2r2) uo = 0. (31)

Уравнения (29) и (30) получены путем исключения переменных no и uo из (26) и (27), уравнение (31) — исключением no из (26) и (28).

Систему уравнений (26)-(28) представим в матричном виде. Здесь и далее столбец переменных — (uo, vo, no)T. Приравняв нулю детерминант матрицы коэффициентов, получим

w ((w2 — 1) r2 — 101 + 1 (c + 1) ft) =0, (32)

где R — безразмерный радиус деформации Россби

R* VgH , Л

R =A= /A ' (33)

ft = - R2 (1 + r2 + cos k — r2 cos k) , (34)

ft = - R2 (1 — r2 + cos k + r2 cos k) . (35)

3. Сеточные захваченные волны

В уравнение (21) подставим выражения для функций и п из (22):

k k г 2u0H cos 2 + 2iv0Hr sin ^ — ¿v£/wr Дп0 = 0. (36)

Систему уравнений (26), (28) и (36) представим в матричном виде и приравняем нулю детерминант матрицы коэффициентов. Получим

(rw sin k — 2 R2 sin2 k + (1 — cos k) r2w2) £ — (1 — cos k) r2w2 + 2 R2 sin2 k + rw sin k = 0.

(37)

Систему уравнений (27), (28) и (36) представим в матричном виде и приравняем нулю детерминант матрицы коэффициентов. Получим

((rw — 2 R2 sin k) (1 — cos k) + w2 sin k) £ — (rw + 2 R2 sin k) (1 — cos k)+w2 sin k = 0. (38)

Уравнение (37) умножим на 1 — cos k и вычтем из него уравнение (38), умноженное на sin k. Из полученного уравнения выразим £:

4 R2rw sin k + 8 R4 sin2 k — w2^i £ =-^-• (39)

w2e2

Исключив из (37) и (38) переменную £, получим

sin k (r2w4 — (А + r2) w2 + 4 R4 sin2 k) = 0. (40)

Пусть выполняется уравнение

r2w4 — (в1 + r2) w2 + 4 R4 sin2 k = 0, (41)

решениями которого являются

wi = Vei + r2 — вз, (42)

r2

w2 =--^ л/в1 + r2 + вз, (43)

r

л/2

шз,4 = л/в1 + r2 Т вз, (44)

rv 2

вз = \/(в1 + r2)2 - 16 r2R4 sin2 k. (45)

В (44) и далее верхний знак в выражении в правой части соответствует первому значению индекса в левой части, нижний знак — второму значению индекса. Подставив выражения для ш из (42)-(44) в (39), получим выражения для £, которые можно привести к следующему виду:

£i = вт (r2 + вз + rWA + r2 + вз) , (46)

где

6 = (r2 - вз - ГV2V01 + r2 - £з) , (47)

Сз,4 = -1 (г2 ± вз Т г + г2 ± вз) . (48)

Запишем некоторые вспомогательные соотношения, которые следуют из определений (34), (35) и (45) при условиях (12) и (25):

в > 0, в12 = в22 + 16 R4r2 sin2 k, вз = r4 + 2 г2в1 + в2- (49)

С применением (49) получаем следующие результаты:

6£з = 1, && = 1, Ы > 1, (50)

в3 - г2 V Г4 + 2 Г2 |&| + в22 - Г2

"Р2Т-

|6|>^" >"-ñT¡- = 1, ^2 |k=n =1- (51)

Из (50) и (51) следует, что для решений 3 и 4 условие (23) не выполняется и требуется рассмотреть только решение 1, определяемое соотношениями (42) и (46), и решение 2, определяемое (43) и (47), при 0 < к < п. Отметим, что уравнение (32) выполняется при подстановке решения 1 или решения 2. Можно показать, что функция ^(к) достигает максимума в единственной точке к = ко:

dk

2 R2 (1 - г2) , n

= 0, cos ko =----, -. (52)

k=ko 2 R2r2 + r2 + 2 R2 + Ы 4 R2r2 + r2 + 4 R2

Графики функций w(k), w(k)/wK(k) и £(k) для решений 1 и 2 при R =1, r = 1/3 представлены на рис. 2, а, б и в соответственно. Графики функций uo/(ivo), Re (no/vo) x yjg/H и Im (no/vo) \/g/H представлены на рис. 3, а, б и в соответственно. Цифра 1 или 2 у кривой соответствует номеру решения. Значения функций, представленных на рис. 3, вычислены с применением (29) и (30).

Л

г>

Рис. 2. Графики w(k), ш(к)/шК(k) и {(k) при R = 1, r = 1/3

Рис. 3. Графики функций u0/(iv0), Re (n0/v0) \/g/H и Im (n0/v0) д/g/H при R = 1, r = 1/3

У решений £1 и £2 есть особенность при в2 = 0. На рис. 2, в этому случаю соответствует вертикальная штриховая линия. Обратимся к исходным уравнениям и предположим, что при в2 = 0 решение существует в следующем виде:

(им-i,n,VM-i,n) = (ui,vi)exp (i/wt — ikn),

Пм-i/2,n+i/2 = ni exp (i/wt — ik(n + 1/2)) , (53)

Um,n = vm,ra = nm+i/2,n+i/2 = 0, m < M — 1.

Подставив выражения из (53) в (13) и (14) при m = M — 1, в (15) — при m = M — 3/2

и в (21) соответственно, получим систему линейных уравнений для ui,vi и пь

k

iu/wrA — v/rA + nig cos 2 = 0, (54)

k

uiA / + iviA/w — inig sin 2 = 0, (55)

kk

ui cos 2 — ivir sin 2 = 0, (56)

kk

— 2ui H cos 2 — 2vi iHr sin 2 + ini/w r A = 0. (57)

Запишем условие линейной независимости уравнений (54)—(56):

1 — r2 + (1 + r2) cosk = 0 (в2 = 0). (58)

Отметим, что, подставив в (35) выражение для cos k из (58), получим в2 = 0. Условие линейной независимости уравнений (54), (56) и (57) можно представить в виде квадратного уравнения для w:

w2 (1 + r2) + w (1 + r2) — 4 R2 = 0, в2 = 0. (59)

С применением (58) преобразуем уравнение (41) к следующему виду:

(w2 (1 + r2) + w (1 + r2) — 4 R2)(w2 (1 + r2) — w (1 + r2) — 4 R2) = 0, в2 = 0. (60)

Очевидно, что решения уравнения (59) совпадают с двумя решениями (60). Этими двумя решениями являются Ш11в2=0 и 1в2=о:

. , \/г2 + 8 Я2 + 1 т VI + Г2 + 16 Д2^Г+г2 , VI + Г2 + 16 Я2 1

Ш12|д п = ±-. =- = ±-, ---. (61)

1,21 в2=0 V^VГTr2 ^VГTr2 2 ^ ;

Покажем, что решение 1 при неограниченном измельчении шага сетки стремится к решению (5), (7). В (29) и (30) подставим выражения (42) и (46) для ш и В полученные уравнения подставим выражения (24) и (33) для к и Я. Полагаем, что к* не зависит от А. Находим пределы

Т Uo

lim —

Д^о v0

Далее находим пределы

= 0, lim — t t Д^о v0

(62)

11ш Ш1 = ^к*, 11ш ег-м = 11ш еХ™/(гА) = ехр Г, (63)

А-0 / ' А-0^1 А—0

где ш1 и определены в (42) и (46) и предполагается, что хт = х и к* не зависят от А. Из (62) и (63) следует, что при А ^ 0 решение 1 переходит к пределу, совпадающему

с (5), (7).

При неограниченном измельчении шага сетки в направлении по нормали к границе решение ш1 стремится к шс — известному решению для чисто гравитационной волны на одномерной сетке с "разнесенными" узлами:

11шШ1 = шс, 11шег-м = 11ш£Гт/(гД) = ехр ( л/2/х- ) , (64)

г-0 ^(Л1 г-^1 у^Я (1 + 0С8 к))

2 VgЯ . к*А (65)

шс = 2^/ 81П~У • (65)

При исследовании зависимости решения 1 от сеточных параметров получены неравенства

дШ1 < 0, (66) 0<к<0

5 А

дг

< 0. (67)

0<fc<0

Кратко изложим доказательство неравенства (66). В (42) последовательно подставим выражения для в3 из (45), для в1 из (34), для k и R из (24) и (33), продифференцируем по А и представим в следующем виде:

дА- = Ai (8R4r2 sink (kcos k - 2 sink) - (в1 + R2ksink (1 - г2)) (вз - в1 - г2)) , (68)

где

A1 =-^ . = > 0, 0 < k < п. (69)

гА вз V^v/e1 + г2 - вз

Производную при к = п, А = Д1 определим как предел следующего отношения: дш1

д А

= lim ^lU=Al'fc*=n/AAl ^1 A=A2,fc*=n/Ai = _2-R n < 0, (70)

k=n A2 >Ai 0 А1 _ А2 rAV4 Rl2 + 1

где R1 = \JgH/(f Ai). Частная производная в (70) вычисляется при фиксированной длине волны Л* = 2n/k* = 2A1 > 2A2.

Получим верхнюю оценку для правой части (68). Отметим, что в3 входит в правую часть уравнения (68) с отрицательным сомножителем. Поэтому при оценке правой части (68) в3 можно заменить нижней оценкой. Запишем вспомогательные соотношения:

вз > Г2 + |в2| , о < k < п; (71)

k

1 - cos k - ^ sin k > 0, 0 < k < п. (72)

Рассмотрим два случая.

1. Пусть в2 > 0. При выводе неравенства заменим в3 на г2+в2 + C, где C — некоторая положительная константа, C > 0. С учетом (34), (35) и (72) получим

д^1

< _41 _ cosk _ 2 sink) _ AiC (01 + R2ksink (1 _ r2)) < 0. (73)

в2>0 V 2 /

дА

в2>0 4 2. Пусть в2 < 0. Теперь заменим вз на r2 _ в2 и с учетом (34), (35) и (25) получим

< 2 в2А1 R2 (2 + 2 cos k + k sin k) < 0. (74)

dw1

дА в2<0

Из (70), (73) и (74) следует (66) и можно сделать вывод, что решение монотонно стремится к Шк с уменьшением А.

Теперь докажем неравенство (67). В (42) последовательно подставим выражения для вз из (45) и для в1 из (34). Результат продифференцируем по г и представим в следующем виде:

дШ1 = А2 (г2 + в2 - вз) , (75)

где

= У-2«2 (1+со.к) > о. (76)

■г2вз\/в1 + г2 - вз

Воспользовавшись вспомогательным неравенством (71), получим (67). Из (67) следует, что ш1 монотонно стремится к Шс при уменьшении параметра г. Отметим, что ш1 = 0 при к = п, поэтому

дШ1 = 0. (77)

к=п

dr

Рассмотрим случай, когда к = п. Представим в матричном виде систему уравнений (27), (28) и (36) при к = п. Приравняв нулю детерминант матрицы коэффициентов, получим

ш/2Д3 г2Я (£ - 1) = 0. (78)

Уравнение (78) выполняется только при ш = 0. Подставим ш = 0 и к = п в уравнения (26)-(28) и (36). Из четырех полученных уравнений три уравнения совпадают. Результат можно записать в следующем виде:

(Л/ё/и0А - (£ + 1) П0) =0, ^4=^ = 0, ш|к=п = 0. (79)

V /

Из (79) следует, что решения существуют при любом £, удовлетворяющем условию (25). Например, можно положить £ = £1|к=п < 1. По-видимому, такие стационарные решения следует трактовать как сеточные двухшаговые захваченные волны, поскольку на границе не выполняется разностный аналог условия (10).

Заключение

В работе построены решения типа захваченных волн для системы дифференциально-разностных уравнений, представляющих собой пространственную разностную аппроксимацию линейных уравнений мелкой воды, и исследованы зависимости полученных решений от параметров сетки. Показано, что решение 1, определяемое соотношениями (42) и (46), при неограниченном дроблении шага сетки монотонно стремится к решению типа волны Кельвина исходной дифференциальной задачи. При неограниченном измельчении шага сетки только в направлении по нормали к границе, частота решения 1 монотонно стремится к частоте чисто гравитационной волны на одномерной сетке с "разнесенными" узлами. Решение 2, определяемое (43) и (47), характеризуется отрицательной фазовой скоростью и, следовательно, является чисто вычислительным решением. По-видимому, его необходимо отфильтровывать или подавлять. Важная особенность решения 2 — высокая частота во всем диапазоне длин волн. Решение такого вида отсутствует в работе [4], где предполагается геострофическое равновесие для продольной компоненты скорости.

У решения 1 поперечная компонента не равна нулю (см. рис. 3, а), отсутствует единый для всех длин волн масштаб экспоненциального убывания амплитуды при удалении от стенки, а при £1 < —1 от узла к узлу изменяется даже знак амплитуды (рис. 3, в). Поскольку решение 1 не обладает отличительными свойствами волны Кельвина, для него целесообразно применять термин "сеточная" волна Кельвина. Предположим теперь, что при х = х1 некоторая совокупность "сеточных" волн Кельвина образует волновой пакет. При х = х2 = х1 структура волнового пакета нарушается, поскольку амплитуда каждой волны в этой совокупности изменилась индивидуально — как функция длины волны. По-видимому, рассматривать распространение волновых пакетов можно только для совокупностей "сеточных" волн Кельвина с относительно близкими значениями £.

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

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

[1] Марчук Г.И., Дымников В.П., Залесный В.Б. Математические модели в геофизической гидродинамике и численные методы их реализации. Л.: Гидрометеоиздат, 1987. 296 с.

[2] Мезингер Ф., Аракава А. Численные методы, используемые в атмосферных моделях: пер. с англ. Л.: Гидрометеоиздат, 1979. 136 с.

[3] ВольцингЕр Н.Е., Пясковский Р.В. Теория мелкой воды. Океанологические задачи и численные методы. Л.: Гидрометеоиздат, 1977. 208 с.

[4] Hsieh W.W., Dayey M.K., WAjsowics R.S. The free Kelvin wave in finite-difference numerical models // J. of Phys. Oceanogr. 1983. Vol. 13, N 8. P. 1383-1397.

[5] Гилл А. Динамика атмосферы и океана: пер. с англ. В 2 т. М.: Мир, 1986. 815 с.

[6] Bryan K. A numerical method for the study of the circulation of the world ocean //J. Comput. Phys. 1969. Vol. 4. P. 347-376.

[7] ДиАнский Н.А., Багно А.В., Залесный В.Б. Сигма-модель глобальной циркуляции океана и ее чувствительность к вариациям напряжения трения ветра // Изв. РАН. Физика атмосферы и океана. 2002. Т. 38, № 4. С. 537-556.

Поступила в редакцию 1 июня 2004 г., в переработанном виде — 22 августа 2007 г.