CC BY
22
Вычислительные технологии
Том 20, № 3, 2015
Численное моделирование распространения длинных поверхностных волн по вращающейся сфере в рамках полной нелинейно-дисперсионной модели
О. И. Гусев, Г. С. Хлкимзянов*
Институт вычислительных технологий СО РАН, Новосибирск, Россия *Контактный e-mail: khak@ict.nsc.ru
Для численного моделирования процесса распространения длинных поверхностных волн предложен алгоритм, основанный на расщеплении системы нелинейно-дисперсионных уравнений на вращающейся сфере на равномерно эллиптическое уравнение для дисперсионной составляющей давления и гиперболическую систему уравнений мелкой воды первого приближения с модифицированным ис-точниковым членом в правой части уравнения импульса. Алгоритм реализован в виде явной двухшаговой схемы предиктор-корректор, на каждом шаге которой поочередно решаются задачи, полученные в результате расщепления. На модельных задачах о распространении волн над ровным дном дана оценка важности учета эффектов, связанных со сферичностью Земли и ее вращением, зависимости дисперсионных эффектов от дальности распространения волн и размеров области начального возмущения свободной границы.
Ключевые слова: вращающаяся сфера, мелкая вода, длинные поверхностные волны, нелинейно-дисперсионные уравнения, численное моделирование, дисперсия, сила Кориолиса.
Введение
До последнего десятилетия моделирование процессов распространения длинных поверхностных волн в океане выполнялось в рамках бездисперсионной модели мелкой воды (NLSW-модели) с использованием различных программных систем. Так, в основе системы TUNAMI, часто применявшейся для моделирования волн цунами, лежит консервативная конечно-разностная схема leap-frog на разнесенной сетке, аппроксимирующая уравнения NLSW-модели с учетом реальной батиметрии [1, 2]. В системе MOST [3] используется метод расщепления по пространственным переменным, и этот программный комплекс также находит широкое применение [4, 5]. Вычислительный алгоритм программного комплекса MGC представляет собой модифицированную конечно-разностную схему Мак-Кормака, аппроксимирующую нелинейные уравнения мелкой воды в сферической системе координат [6, 7]. С помощью этого комплекса решались практические задачи распространения волн цунами, а его модификации использовались при решении задач наката волн на берег [8] и генерации волн подводными оползнями [9].
При наличии на поверхности воды коротких волн существенную роль начинает играть частотная дисперсия. Еще в 1996 г. в работе [10] было указано, что "... учет
© ИВТ СО РАН, 2015
нелинейности и дисперсии волн цунами является принципиальным и необходимо развитие соответствующих методов". Позднее на основе анализа результатов математического моделирования и натурных данных, в особенности данных о крупнейшем цунами 2004 г. в Индийском океане, сформировались более четкие представления о математических моделях и алгоритмах, способных описывать различные стадии катастрофических цунами с требуемой для практики точностью [11, 12]. Главный вывод заключался в том, что для приемлемого описания продолжительного явления и в больших по широтному и долготному направлениям акваториях требуются нелинейно-дисперсионные модели, учитывающие неровность дна и подвижность некоторых его фрагментов, криволи-нейность береговых линий, наличие островов и эффекты, связанные со сферичностью и вращением Земли.
К настоящему времени накоплен пока что небольшой опыт численного моделирования процессов распространения длинных поверхностных волн в океане с учетом дисперсии. В некоторых работах (например, в [13, 14]) сферичность Земли явно не учитывалась, а расчетной областью служила проекция океанической акватории на касательную плоскость. Фактически расчет выполнялся в плоской области с помощью плановой слабо нелинейной дисперсионной (ШКЬЮ) модели типа Буссинеска без учета силы Ко-риолиса. Показано, что учет дисперсии изменяет волновую картину, полученную с помощью бездисперсионной КЬБШ-модели.
В работе [15] изучалось трансокеаническое распространение волн цунами от огромного гипотетического оползня, который потенциально может произойти на острове Ла Пальма, входящем в состав Канарских островов. Также использовались ШКЬЮ-уравнения, но теперь уже записанные в сферической системе координат вращающейся Земли. Применяемая модель не позволяла учитывать подвижность дна, поэтому начальный этап генерации волн рассчитывался с помощью известной трехмерной модели гидродинамики и полученные поля служили в качестве начальных данных для ШКЬЮ-модели, по которой рассчитывалось длительное распространение волн. Слабо нелинейная модель не могла воспроизводить картину течения в прибрежной зоне из-за быстрого роста здесь нелинейных эффектов, поэтому расчеты выполнялись до выхода волн на прибрежный шельф. Показано, что при длительном распространении дисперсия оказывает заметное влияние на волновую картину. В отличие от КЬБШ-модели, которая дает выраженную одиночную волну, при использовании ШКЬЮ-модели возникает цуг волн, при этом амплитуда головной волны может стать меньше волн в дисперсионном хвосте [16]. Кроме того, на мелководье возможно появление ондулярных боров [17], которые не описываются КЬБШ-моделью [18].
Более глубокие исследования влияния дисперсии на процессы распространения волн цунами выполнены в недавней работе 2013 г. [19], в которой также используется ШКЬЮ-модель, но влияние дисперсии проанализировано сразу для нескольких реальных событий. Отмечается, что около побережья ШКЬЮ-модель дает неудовлетворительные результаты, поэтому для сквозного счета желательно использовать полностью нелинейную дисперсионную модель. В том же году опубликована работа [20], где выведены полные нелинейно-дисперсионные ^N1^) уравнения мелкой воды на вращающейся сфере, при этом в качестве скорости в ЕМЬЮ-модели взята, в отличие от [21], горизонтальная составляющая скорости трехмерного течения на некоторой поверхности, лежащей между дном и свободной границей. Предложенная полная модель имеет ряд существенных недостатков. В частности, для нее задача Коши может быть некорректной, поскольку она основывается на модели Nwogu [22], решение которой неустойчиво
по начальным данным при некоторых конфигурациях дна [23]. Отметим также, что авторы работы [20] пока что не представили результатов моделирования на основе полной модели, выполнив исследование влияния дисперсии и силы Кориолиса при распространении длинных волн над ровным дном лишь в рамках WNLD-модели, являющейся слабо нелинейным аналогом выведенной ими FNLD-модели. Для сферической геометрии конечно-разностная схема с TVD-ограничителями построена путем обобщения аналогичной схемы [24] для декартовых координат. Численный алгоритм для WNLD-модели реализован в программном комплексе FUNWAVE [25].
Полная FNLD-модель на сфере, в которой в качестве скорости берется осредненная по толщине слоя горизонтальная составляющая скорости трехмерного течения, была получена в 2010 г. [26]. В [27] показано, что FNLD-модель из [26] можно вывести без предположения о потенциальности исходного трехмерного течения. В отличие от [20], она записывается в виде законов баланса массы и импульса, а также имеет согласованное с трехмерной 3D-моделью Эйлера уравнение баланса полной энергии, что не только подтверждает физическую состоятельность модели, но и позволяет осуществлять дополнительный контроль при проведении расчетов. Кроме того, показано, что структура уравнений FNLD-модели [27] и такие положительные свойства, как консервативность и наличие уравнения баланса полной энергии, согласованного с 3D-моделью, сохраняются при поэтапном упрощении дисперсионной составляющей. Так, уравнения полученной слабо нелинейной дисперсионной модели [27] и классической (бездисперсионной) модели мелкой воды на сфере имеют тот же вид, что и уравнения FNLD-модели на сфере с тем лишь отличием, что меняется выражение для кинетической энергии и для WNLD-модели члены с давлением вычисляются с использованием более сложных формул, чем в NLSW-модели. Таким образом, обеспечивается единая форма записи уравнений всей иерархической цепочки моделей мелкой воды на сфере [28], что позволило сконструировать для них единый численный алгоритм.
В настоящей работе рассматривается численный алгоритм для моделирования процесса распространения длинных поверхностных волн в рамках полной FNLD-модели, учитывающей сферичность Земли и ее вращение. На основе анализа результатов решения модельных задач о распространении волн над ровным дном дана оценка важности учета эффектов, связанных со сферичностью Земли и ее вращением, зависимости дисперсионных эффектов от дальности распространения волн и размеров области начального возмущения свободной границы. Расчеты выполнены с помощью программы NLDSW_sphere [29], предназначенной для решения в общей постановке задач о распространении трансокеанических цунами в рамках FNLD-уравнений в сферических координатах с учетом подвижности дна, силы Кориолиса и центробежной силы, криволи-нейности береговых линий и наличия островов.
1. Постановка задачи
Уравнения полной FNLD-модели на вращающейся притягивающей сфере имеют следующий вид [26]:
(HR sin 6)t + (Ни) л + (Hv sin 9)e = 0, (1)
(HuR sin 9)t + ( Hu2 + g—j + (Huv sin d)e = gHh\ — Huv cos 9—
V 2 / л (2)
— fvHR sin 9 + — фh л,
Г ( Н 2\ 1 Н 2
(HvR sin 9)t + (Huv)x + [Hv2 + g—) sin 9 = gHhe sin 9 + g— cos 9+
L \ 2 J -10 2
(3)
+Hu2 cos 9 + /иЯЛ sin 9 + (рв — фhe) sin 9 + fl2HR2 sin2 9 cos 0,.
=0 при замене (16)
Здесь R — радиус сферы, вращающейся с постоянной угловой скоростью П вокруг оси Oz неподвижной декартовой системы координат Oxyz, координатная плоскость Оху которой совпадает с экваториальной плоскостью сферы. Для описания течения воды используется вращающаяся вместе со сферой система координат ОХвг, начало которой находится в центре сферы. Здесь Л — долгота, отсчитываемая в направлении вращения от некоторого меридиана (0 < Л < 27т), 9 = п/2 — р — дополнение до широты (—к/2 < р < тг/2); г — радиальная координата, отсчитываемая от центра сферы. Ньютоновская сила притяжения д, действующая на частицу жидкости единичной массы, считается направленной к центру Земли. Толщина слоя воды Н = r¡ + h > 0 предполагается малой по сравнению с радиусом сферы, поэтому величина д = |д| и плотность воды р принимаются постоянными во всем жидком слое, ограниченном снизу непроницаемым подвижным дном, а сверху — свободной поверхностью:
г = R — h(X,9,t), г = R + r](X,9,t). (4)
Через и и v обозначены физические компоненты вектора скорости (u = Re1 sin 9, v = Re2, с1 = X, с2 = в), f = 2П cos 9 — параметр Кориолиса, выраженный через дополнение к широте 9, при этом предполагается, что
9о < 9 < -к — 9о, (5)
где 90 = const > 0 — малый угол (полюсы исключаются из рассмотрения).
Величины р и ф, входящие в правые части уравнений движения (2), (3), представляют собой дисперсионные составляющие проинтегрированного по толщине слоя давления р в FNLD-модели и соответственно давления на дне р0:
Р = — V, ро = дН — ^. (6)
Дисперсионные добавки вычисляются по следующим формулам [30]:
Н 3 гг 2 тт 2
V = -3-Qi + -уQ2, Ф = ~2ГQi + HQ2, (7)
где Qi = D(V ■ c) — (V- c)2, Q2 = D2h, c =(c1,c2),
^ ^ ( d d \ ^ 1 d 2 d = at +c'V- V = {ах'вв)' c'V = Ж + CV
V- c = d + 1 (Jc2)e , J = —R2 sin 9. (8)
J
В более подробной записи выражения для Q1 и Q2 принимают следующий вид:
Q1 = (V ■ c)t + (и (V ■ c)A + v (V ■ c), sin 9) — (V ■ c)
2
= {DK)t + fíS^e ÍU (Dh)x + " {DK)e 81П 9) '
где
v ■ с = R—e (UA +(v Sin ), Dh = h + жЬ К +vhe Sin 0) •
FNLD-модель (1)-(3) называется полной [30] по той причине, что она выведена без предположения о малости амплитуды волн и в уравнениях сохранены все нелинейные члены, связанные с дисперсией. Ее можно использовать для расчета поверхностных волн, распространяющихся над неровным дном как по глубокой воде, так и в прибрежной зоне, при этом учет дисперсии волн может дать более точные результаты, чем бездисперсионная NLSW-модель. FNLD-модель позволяет моделировать также волны, генерируемые продолжительными по времени подвижками фрагментов дна, что расширяет круг задач, которые могут решаться в рамках известных WNLD-моделей мелкой воды на сфере [15, 16, 19]. Кроме того, эта модель обладает рядом дополнительных преимуществ, о которых говорилось во Введении.
Что касается воспроизведения дисперсии, то считается, что модели с усредненной скоростью, к которым принадлежит и рассматриваемая FNLD-модель (1)-(3), не обладают улучшенными дисперсионными свойствами, присущими линеаризованной модели из работы [22]. Однако надо заметить, что полученные на основе [22] нелинейные модели (например, [31]), в которых в качестве скорости берется горизонтальная составляющая вектора скорости трехмерного течения на некоторой вполне определенной поверхности, дают действительно лучшие дисперсионные свойства только в линейном одномерном случае. Для нелинейных уравнений это преимущество теряется, особенно при подвижном дне, что и подтверждается вычислительными экспериментами [32, 33], результаты которых показывают, что модель [31] не имеет преимуществ по дисперсионным свойствам перед нелинейно-дисперсионными моделями с осредненной скоростью.
Уравнения (1)-(3) имеют дивергентную форму записи своих левых частей. В рассматриваемом ниже численном алгоритме решения задач в рамках FNLD-модели будет использоваться также недивергентная форма этих уравнений
Ht + —i— [(Hu)x + (Hv sin в)в] = 0, (9)
R Sin и
1 , 1 1 1 Pa -^hAUv , ,
u + ñ^neUUA + Rvue + Kw= Шпё - rctg 0 - fv, (10)
vt + 1 uva + ^vЩ + 1gríe = 1 Pe + -ctg в + fu + Q2Ksinocosfl, (11)
R Sin и R R R H R v
=0 при замене (16)
в которой в явном виде выделены слагаемые правой части уравнений движения, связанные с силой Кориолиса и центробежной силой.
При пренебрежении дисперсионными членами, т. е. при p = ф = 0, из FNLD-модели получается бездисперсионная NLSW-модель мелкой воды на вращающейся притягивающей сфере [34], уравнения которой могут быть записаны в дивергентной форме вида (1)-(3) или недивергентной, аналогичной (9)—(11):
Ht + —Ц [(Hu)a + (Hv sin^] = 0, (12)
R Sin a
1 1 1 UV
ut + R^neUUA + RVUe + RSñe9VA = - R ctg ^ - fV, (13)
1 1 1 м2 vt + О • nuvа + ~БVVe + "т;9Ve = -77ctg в + fu + О2 К sin 9 cos 9j . (14)
К sin У К К К V
=0 при замене (16)
В начальный момент времени задавалось отклонение свободной границы от ее невозмущенного положения и значения компонент вектора скорости. Заметим, что для рассматриваемых в сферической геометрии моделей мелкой воды (9)—(11) и (12)-(14) невозмущенная свободная граница в состоянии покоя отличается от сферической [27] и описывается уравнением г = К + r¡00(9), где
щ0(9) = — О2К2 sin2 9 + const. (15)
29
В практических задачах высота волны и глубина отсчитываются от поверхности покоящейся воды. Поэтому будет более естественным задавать функции, описывающие дно и свободную границу, не в виде отклонений (4) от поверхности сферы радиуса К, а в виде отклонений h и f¡ от невозмущенной свободной границы, как это обычно и делается при решении уравнений мелкой воды в плоской геометрии. Для этого достаточно сделать следующие замены:
-h = 'По0 - h, r¡ = Щ0 + fj. (16)
Посмотрим, к каким изменениям в уравнениях приведут замены (16). Очевидно, что полная глубина не зависит от выбора поверхности отсчета: Н = r¡ + h = f¡ + h. Кроме того, r¡\ = f¡x и gr¡e = gщ+О2К2 sin 9 cos 9, поэтому, задавая глубину и свободную границу в виде отклонений от невозмущенной свободной границы, мы получаем ту же NLSW-систему (12)-(14), но только с отброшенным слагаемым О2К sin 9 cos 9 в уравнении (14) и с заменой h и на h и соответственно.
При подстановке выражений (16) в системы FNLD-уравнений (1)-(3) или (9)-(11) последние слагаемые в правых частях уравнений (3) и (11), связанные с центробежной силой, также сокращаются. Но в силу равенства D2h = D2h — D2r¡00 изменятся выражения для дисперсионных составляющих давления (7) в FNLD-модели, поскольку Q2 теперь не будет совпадать с D2h.
Кроме начальных условий уравнения дополняются также краевыми условиями, которые могут иметь различный вид на разных частях границы области решения. Криволинейная береговая линия аппроксимируется ломаной Г, звенья которой проходят по меридианам и параллелям. На линии Г ставится условие непротекания
u ■ n|r = 0, (17)
где n — вектор внешней нормали к соответствующему звену ломаной Г, u = (м, v)T. Таким образом, задача наката на береговой откос заменяется задачей наката на вертикальную стенку, установленную вдоль береговой линии на некоторой заданной глубине hw, причем дно акватории в окрестности береговой линии предварительно выравнивается так, что если реальная глубина меньше некоторого заданного значения hw, то она полагается равной hw.
2. Численный алгоритм
Система уравнений ЕЫЬВ-модели (1)-(3) не является системой типа Коши — Ковалевской в силу того, что уравнения движения (2), (3) содержат смешанные производные
третьего порядка по времени и пространственным переменным от компонент вектора скорости. Непосредственная аппроксимация этих производных приводит к сложной разностной задаче, не поддающейся исследованию. В плоском случае для аналогичной системы ЕЫЬВ-уравнений [35] плодотворным для конструирования численного алгоритма оказалось предварительное расщепление этой системы на скалярное уравнение эллиптического типа и систему уравнений гиперболического типа [36 - 39].
В настоящей работе такой же подход реализован для системы ЕЫЬВ-уравнений (1)-(3) на вращающейся сфере. В результате расщепления получается расширенная система уравнений, состоящая из эллиптического уравнения для дисперсионной составляющей р проинтегрированного по глубине давления р и гиперболической системы уравнений (1)-(3) с производными по времени первого порядка, которая отличается от классической модели мелкой воды только наличием в правой части уравнений движения дополнительных слагаемых, связанных с дисперсионными добавками р и ф. При таком расщеплении численный алгоритм целесообразно строить на поочередном решении на каждом шаге по времени эллиптического уравнения и гиперболической системы.
Вывод уравнения для р дан в Приложении, поэтому приведем здесь лишь окончательный его вид:
1 /Рл _ Ур ■ УН. \"| + Г/Фе. _
где к0 = коо + (Аю1)л + (^02)0,
Ур -УН.
Н
?) в
Не ) вт 9
- кор = Р,
:18)
6Н л
, 12(г - 3) „2 . л ,
коо =--К V, «01 = ц-2 • о,
Н 3г Н 2г втс/
к
о2
6 Не Н 2г
вт 9.
19)
Далее для краткости записи используются обозначения (16), поэтому в развернутой записи уравнения (18) надо всюду вместо Н подставлять выражение
Н = Н - 7]оо.
(20)
С учетом этого замечания формула для скалярного произведения У р ■ УН запишется
как
Ур^УН=К
рлНл ып29
+ ре (Не - Що,в
- П2 Ур ■ УН - ре— вт 9 есв 9.
(21)
Отличительной особенностью уравнения (18) является то, что ни левая, ни правая его части не содержат производных по времени от зависимых переменных Н, и и оно отличается от аналогичного уравнения в плоском случае [37] лишь наличием членов, обусловленных сферичностью и вращением Земли, в том числе и членов, входящих в слагаемые правой части Р = Р1 + Р2 + Р3:
Р1
1 ( + ЯН
9'1л +--Нл - а,1
-вШ 9
+
Я \
9Щ +--Не - а2 вш 9
(22)
Р2 = -6ЯК2 й1п 9, г = 4 + УН ■ УН,
Н
2
2 / \ 2
Рз = —^ [ил + (V вт в)е) - 2(илУе - ьлЩ) - 2(иь )л еtg 9 - (ь2 есв в)е вт их /
где
Я = (-дУП + а) ■УН +
1
К2 вт 9 \ вт 9
(
увт 9
Нлл + 2иу Нле + V Нее вт 9 + В
)
е
л
е
(U V \
RSñ> h* + r4,
U V
R sin и R
a = (ai,a2)T, a\ = — (2uv ctg0 + fvR)sinв, a2 = u2ctg0 + fuR. (23)
В слагаемом В сгруппированы члены, связанные с подвижностью дна. Для стационарного дна имеем В = 0.
Легко показать, что при Н > 0 и условии (5) уравнение (18) является равномерно эллиптическим. Для корректности задач, связанных с решением этого уравнения, важен знак коэффициента к0 при функции р. Так, например, одним из критериев того, что задача Дирихле для уравнения вида (18) с заданной правой частью F(х, t) не может иметь более одного решения, является неравенство к0 > 0 [40]. В случае ровного дна h = h0 = const это неравенство принимает следующий вид:
_ 12 R2 sin 9 ко = Н Зг
где
Q2 / 8Н \
г — 3--(н sin2 9 — т sin 2в +-cos 20)
2
Q4R2 . 2
г = 4 +--— sin2 2 в.
4д 2
> 0, (24)
Используя значения R = 6.38 • 106 м, П = 7.29 • 10-5 с-1, д = 9.81 м • с-2, а также предполагая, что градиент функции f¡ ограничен и выполнено условие (5), приходим к оценке к0 ^ 0, т. е. рассматриваемый критерий заведомо выполняется. В этом случае для нахождения численного решения уравнения (18) можно построить разностную схему с положительно определенным оператором. Для сильно неровного дна коэффициент к0 может стать отрицательным в некоторых подобластях, что обычно приводит к ухудшению свойств разностных схем и уменьшению скорости сходимости итерационных методов решения уравнения (18). Для одномерных задач анализ таких случаев выполнен в работе [41]. Как правило, избежать подобных нежелательных ситуаций удается путем сглаживания функции h, описывающей дно реальной акватории.
Разностное уравнение, аппроксимирующее (18) со вторым порядком, построено с помощью интегроинтерполяционного метода [42]. Для каждого внутреннего узла прямоугольной равномерной сетки с шагами h1 и h2 в направлении осей OX и Ов соответственно берется прямоугольник ABCD с вершинами в центрах соседних ячеек (рис. 1, а), и уравнение (18) интегрируется по этому прямоугольнику. После применения формулы Грина получается соотношение с контурными и двойными интегралами:
J Ф1 сШ - j Ф1 сШ + У Ф2dX - J Ф2(1Х - jj kapdXde = jj FdXdO, (25)
ВС AD DC AB ABCD ABCD
при этом двойной интеграл от функции (22) также заменяется контурным:
jI F1dXd9 = j fudd - У fudd + у Í12dX - У fndX, (26)
D B C A D D C A B
где
! 1 m v^-v^ \ 2 (ve v^-v^ \ .
Ф = ¡rn0 (Н — hx), ф = VH — 4sin(27)
fu = -n~fí (9V, + — fll) , Ы = (<ñ + 0he — a2) sin 0. sin
8
J 2-
D ^ 1 W\
N C -----1
E
S B 2
8 4
D с
1 A ! 0 \б .
J 2-
Jl
Jl
в 8 \4
D 1 W N 1
0 j
A 5 S B 2
Рис. 1. Контуры интегрирования и шаблоны разностных уравнений для р во внутреннем (а) в граничном (б) и угловом (в) узлах сетки
7
3
3
7
0
5
6
3
6
2
Для приближенного вычисления интегралов использовалась квадратурная формула трапеций с аппроксимацией функций и их первых производных со вторым порядком в центрах ячеек [36]. В результате для сеточной функции р получается 9-точечное неоднородное разностное уравнение в каждом внутреннем узле сетки.
Аналогичным образом могут быть получены разностные уравнения в граничных узлах, лежащих на непроницаемой стенке. Пусть, например, узел сетки Xj1,j2 = ( \j1, 0j2) принадлежит участку границы Г, проходящему по некоторой параллели (жирная линия на рис. 1, б), и в локальной окрестности этого узла область решения лежит выше рассматриваемой части границы. Если, например, узлы Xj1-1j2 и x.,1+1j-2 также лежат на Г, то контур интегрирования ABCD выбирается так, что две его вершины С и D совпадают с центрами соседних ячеек, имеющих общую вершину в точке x 1, 2, а две другие — A и B — с серединами сторон этих ячеек (рис. 1, б). В этом случае сторона A B контура интегрирования целиком лежит на Г, поэтому при вычислении интегралов по AB, входящих в интегральные соотношения (25), (26), необходимо использовать краевые условия для функции р.
Из условия (17), которое в рассматриваемом случае имеет вид v = 0, и на основе уравнения (11), которое предполагается выполненным вплоть до границы и в котором учтены замены (16), получаем следующее краевое условие:
R
О
H
1 (pe \ и
-9 Ve) + r ctg 0 + fu
0.
Дисперсионная составляющая давления на дне связана с искомыми величинами и р выражением (см. Приложение)
Ф = 1 (6Нг + НЯ + Ур ■ ун),
после подстановки которого в (29) краевое условие принимает вид
(29) H, u
(30)
(
pe v p vh
H
H
e
h e
6h e
R H2
p
1 f - Q,\ u2 n "
К ( 9Ve + ~he I - Rctg0 - fu
или (с учетом обозначений (19), (23), (27), (28))
(Ф2 - ko2'P - /12)
0.
(31)
АВ
Г
Г
Г
Краевое условие (31) используется при вычислении интегралов по стороне АВ. Аппроксимируя двойной интеграл из левой части (25) выражением
JJ ко(рс1\сШъ JJ коо^Х(Ю +( J когсШ^ когсШ + J ЫЛ ^ к02<1 Л)р(0) (32)
АВСИ АВСИ ВС АО ИС АВ
и полагая
У Ф2^Л « Ф2(0)НЬ У /12$Л « /12(0)^1, У ЫЛ « ко2(0),
А В А В А В
получаем, что в силу краевого условия (31) все разностные выражения, аппроксимирующие интегралы по АВ в соотношениях (25), (26), (32), взаимно сократятся. Таким образом, при использовании интегроинтерполяционного метода можно не аппроксимировать краевое условие вида (31), достаточно лишь учесть его при получении разностных уравнений в граничных узлах.
Для граничных узлов других типов используется такой же способ получения разностных уравнений с выбором в каждом отдельном случае соответствующего контура интегрирования. В качестве примера на рис. 1, в изображен контур интегрирования и 8-точечный шаблон для одного из восьми типов угловых узлов. Поскольку разностные уравнения выписываются для всех узлов расчетной сетки, число таких уравнений совпадает с общим количеством внутренних и граничных узлов.
Для решения системы разностных уравнений используется итерационный метод. В общем случае область решения имеет сложную форму границы и не будет выпуклой ни в одном из двух координатных направлений. Наличие островов, на границах которых ставятся условия непротекания вида (31), делает область многосвязной. Поэтому применение итерационных методов, основанных на прогонках, затруднительно. По этой причине в программном модуле КЬБ8Ш_8рЬеге [29] реализован метод последовательной верхней релаксации, более подходящий для многосвязных областей сложной конфигурации.
Для заданных функций р и ф система уравнений (1)-(3) при условии (5) и Н > 0 имеет строго гиперболический тип, поэтому для ее численного решения можно использовать весь арсенал методов, созданных для гиперболических систем. В одномерном случае для решения аналогичной системы использовалась схема предиктор-корректор [37, 41], надлежащий выбор схемного параметра в которой гарантирует выполнение ТУВ-свойства для линеаризованного аналога разностной схемы и сохранение монотонности численного решения скалярных уравнений [43]. Некоторые особенности применения этой схемы для двумерных задач отражены в [44]. В настоящей работе эта схема используется со схемным параметром, равным нулю, что повышает риск возникновения численных осцилляций на крутых фронтах волн при использовании КЬБШ-модели, но приводит к малой диссипации разностной схемы для ЕЫЬВ-модели. Кроме того, на шаге предиктор теперь вычисляются не потоки, а непосредственно величины Н, и, V. Связано это с тем, что именно эти величины, а не потоки требуются для вычисления правой части Р уравнения (18).
Кратко опишем алгоритм решения расширенной системы уравнений (см. также [36, 38]). Как уже говорилось, в начальный момент времени задаются возвышение г/ и вектор скорости и. Если дно предполагается подвижным, то при Ь = 0 должны быть
известными также величины к и кЭтих данных достаточно, чтобы определить в начальный момент времени дисперсионную составляющую проинтегрированного давления р путем численного решения уравнения (18) и дисперсионную составляющую давления на дне ф из конечно-разностного аналога формулы (30).
Пусть на п-м слое по времени величины Нп, ип, рп, фп известны. Далее используется метод предиктор-корректор, каждый из двух шагов которого состоит из двух этапов. На шаге предиктор вначале в центрах ячеек определяются величины Нп+1/2, ип+1/2 как решение явных разностных уравнений, аппроксимирующих недивергентную систему (9)—(11), с правыми частями, взятыми с п-го временного слоя. Затем решается разностное уравнение относительно рп+1/2, коэффициенты и правая часть которого вычисляются с использованием новых значений Нп+1/2, ип+1/2, и по формуле (30) определяется ффп+1/2. Найденные на шаге предиктор значения Нп+1/2, ип+1/2, рп+1/2 и фп+1/2 используются затем на шаге корректор для определения окончательных значений Нп+1 и ип+1 путем численного решения системы уравнений (1)-(3) с дивергентной формой левой части. В последнюю очередь вычисляются значения функций рп+1 и фп+1.
Для КЬБШ-уравнений этот алгоритм упрощается за счет того, что из него исключаются этапы вычисления дисперсионных добавок. В этом случае его можно назвать сферическим аналогом схемы Лакса — Вендроффа.
Важным свойством разработанного численного алгоритма является сохранение им состояния покоя в случае неподвижного (В = 0) и не очень искривленного ( к0 > 0) дна: если на п-м слое по времени жидкость покоилась (ип = 0), а свободная граница была невозмущенной (г/п = 0), то это состояние покоя сохранится и на ( п + 1)-м временном слое. Добиться выполнения этого свойства удалось за счет согласованной аппроксимации левых и правых частей уравнения для р и уравнений импульса.
3. Некоторые результаты численного моделирования
К настоящему времени сложился большой набор тестовых задач [45], принятых в качестве обязательных тестов для проверки достоверности результатов, получаемых с помощью численных алгоритмов расчета течений жидкости со свободной границей в одно-или двумерном приближении с использованием декартовых координат, а также для сравнения характеристик различных алгоритмов [46]. Для верификации методов расчета в рамках нелинейно-дисперсионных моделей на вращающейся сфере пока что не существует общепризнанных тестов, но работа по накоплению необходимого объема материалов с расчетными данными уже началась. К таким материалам можно отнести и результаты решения тестовой задачи, приведенной ниже.
Мы ограничимся здесь решением модельной задачи в простой области, простирающейся от 100 до 300° с запада на восток и от -60 до 65° с юга на север (рис. 2). Для наглядности здесь и далее вместо переменной используется географическая широта р = к/2 — 9. В рассматриваемой области содержится значительная часть Тихого океана, и она не включает полюса (см. условие (5)). Идеализация заключается в том, что глубина океана считается постоянной: к = 4 км. Кроме того, идеализирован источник волны цунами — начальное возмущение гауссовой формы:
Начальная амплитуда волны равнялась а0 = 5 м, центр возмущения располагался в точке ( Ао,ро) = (280, —40°), а параметр т выбирался равным т = 8 • 10-10, 8 • 10-11,
(33)
Рис. 2. Схема модельной акватории и изображения свободной поверхности (вид сверху), рассчитанные по РКЬБ-модели для источника Шз на момент времени £ = 6 ч (а) и 23 ч (б)
8 • 10-12 м-2, что соответствует эффективной протяженности начального возмущения, приблизительно равной ~ 107.3 км, ~ 339 км и ~ 1073 км соответственно.
В формуле (33)
р(Х, ф) = К • агеео^ еов(^) еов(^0) сов(А — Ло) + вт(^) вт(^0) ^ ,
т. е. р есть расстояние по дуге большого круга между точками сферы с координатами (Л, ф) и (А0, ^о), поэтому линии уровня поверхности (33) представляют собой концентрические окружности. Эффективная протяженность определяется как диаметр окружности, являющейся линией уровня а0/10. На наш взгляд, начальное возвышение (33) более подходит для обязательного теста, чем несимметричный источник из [20], поскольку при пренебрежении вращением Земли (П = 0) оно будет генерировать симметричные волны с линиями уровня в виде концентрических окружностей, а при учете вращения Земли использование (33) позволяет оценить при £ > 0 влияние силы Кориолиса, например, по отклонениям линий уровня свободной границы от окружностей.
На границах области решения ставились неотражающие краевые условия типа Зом-мерфельда, аналогичные применявшимся ранее в плоском случае [36]. Большинство расчетов выполнялось на двухминутной сетке с числом узлов 6001 х 3751. Вычисления проводились до момента £ = 30 ч физического времени распространения волн, при этом среднее время расчета на одном ядре процессора кластера НГУ составляло порядка 6 сут. Время расчета в рамках КЬБШ-модели было примерно в пять раз меньше. В некоторых задачах использовалась равномерная в плоскости ОХр сетка с шагом 40 угловых секунд.
3.1. Влияние центробежной силы
Для оценки влияния центробежной силы на процесс распространения волн проведены расчеты по ЕЫЬВ-модели и ее модифицированной версии ЕКЬВ0. Как уже говорилось, после замены (16) последние слагаемые в уравнениях (3), (11) исчезнут, однако выражения, связанные с центробежной силой и имеющие в качестве множителя величину П2, останутся в членах этих уравнений, содержащих производную . Кроме того, выражения с П2 присутствуют в левой части уравнения (18) ив слагаемых и Р2 правой его части. Модифицированная ГКЬВ0-модель получается занулением всех таких слагаемых. Таким образом, в ГКЬВ0-модели, в отличие от ЕМЬВ, влияние центробежной силы не учитывается.
Значения относительных разностей 5* в зависимости от эффективной протяженности источника, %
Мареограф wl w2
м1 —0.38 —0.19 —0.023
м2 0.27 0.036 0.006
Колебания свободной границы фиксировались двумя виртуальными мареографами, установленными в точках М1 = (200°, 0°) и М2 = (120°, 45°) (рис. 2). На мареограммах результаты расчетов по двум моделям практически неразличимы, поэтому приведем здесь лишь таблицу со значениями относительных разностей максимумов в мареограм-мах, рассчитанных по разным моделям:
. = шахм, (ЕЫЬБО) — шахм, (ЕЫЬБ) 2
* = шахм4 (ЕЫЬБО) , г= , '
Из таблицы видно, что с увеличением протяженности источника влияние центробежной силы ослабевает и значения относительной разности 8* во всех случаях не превышают половины процента. Из этого можно сделать вывод о слабом влиянии центробежной силы в рассматриваемом классе задач.
3.2. Влияние сферичности
Влияние сферичности оценивается здесь путем сравнения результатов расчетов на основе двух полных нелинейно-дисперсионных моделей: ЕЫЬЮ-модели на сфере и плановой ЕЫЬЮ-модели на плоскости [37, 38]. При моделировании поверхностных волн в рамках плановой модели начальные возмущения (33) имели такую же ширину Ш1, Ш2 и Ш3, как и в сферическом случае, а два виртуальных мареографа располагались на том же расстоянии от источника, что и рассмотренные выше М1 и М2. Для выделения чисто сферических эффектов при использовании ЕЫЬЮ-модели на сфере полагалось П = 0, т.е. эффекты, обусловленные вращением Земли, здесь не учитывались.
На рис. 3 показаны мареограммы, рассчитанные на основе указанных моделей. Видно, что влияние сферичности увеличивается с ростом размеров источника. Кроме того, во всех рассмотренных случаях амплитуды волн, зафиксированные мареографами, в расчетах на сфере больше, чем на плоскости. При этом в мареографе М1 различия еще не слишком большие, но с увеличением пройденного волной пути различие амплитуд волн возрастает (см. мареограммы М2) и может достигать сотен процентов. Причина столь сильного различия волновых картин кроется в том, что в плоской задаче при разбегании волн происходит монотонное уменьшение их амплитуд на лучах, исходящих из центра возмущения. В сферической же задаче криволинейные лучи, проходящие по дугам больших окружностей, вначале расходятся, но потом сходятся в точке, диаметрально противоположной центру возмущения ( А0,р0). Поэтому амплитуда волн падает, как и в плоском случае, на участках расходимости лучей, но затем возрастает при их сходимости. Несколько таких лучей, имеющих в плоскости координат О Ар криволинейную форму, изображено на рис. 2 в виде сплошных линий, исходящих из вершины ( А0, р0). Видно, что мареограф М2 расположен недалеко от точки сходимости лучей, поэтому в сферической задаче здесь получается волна существенно большей амплитуды, чем в плоской задаче.
Рис. 3. Мареограммы М1 (а- в) и М2 (г- е), полученные в расчетах на основе сферической РМЬБ-модели (!) и плановой модели (2) при эффективной протяженности начального возмущения (а, г), Ш2 (б, д) и Шз (в, е)
Следует заметить, что описанные выше особенности распространения волн по сферической поверхности должны проявляться по всей сфере, однако из-за ограниченности модельной области и использования неотражающих краевых условий образуются подобласти, в которые волны практически не доходят. Это происходит из-за того, что лучи, идущие к этим подобластям из центра начального возмущения, выходят за границы модельной области, при этом теряется информация о волнах, идущих вдоль этих лучей. Примером такой подобласти в рассмотренной задаче может быть окрестность угловой точки (100, —60°). Когда лучи вновь приходят в область решения, то неотражающие граничные условия не впускают волну, идущую вдоль этих лучей. Поэтому волновая картина становится несимметричной, как это видно из рис. 2, б, на котором изображен вид сверху на свободную поверхность с использованием шкалы оттенков серого цвета, отображающей величину отклонения поверхности воды от невозмущенного положения. Средняя часть волны имеет фронт, перпендикулярный лучам. Он хорошо воспроизводится в расчете, поскольку лучи, пришедшие в этот участок волны, не выходили за пределы области решения. На рис. 2, а изображена свободная поверхность в более ранний момент времени. Здесь практически все лучи пришли на фронт волны, не выходя из области решения, поэтому форма фронта слабо отличается от дуги окружности.
Отметим работу [47], в которой также указывается на важность учета сферичности. В этой работе в рамках WNLD-модели численно исследовалось цунами в Индийском океане в декабре 2004 г. Вычислительная область имела значительно меньшую протя-
женность (не более 40° в каждом из координатных направлений), чем в рассматриваемой нами модельной задаче, поэтому эффекты, связанные со сферичностью Земли, проявились не столь ярко: изменение максимальной амплитуды волн составило лишь 30 % в сравнении с плоским случаем.
3.3. Влияние силы Кориолиса
Рассмотрим теперь влияние вращения сферы на процесс распространения волн. Как уже было выяснено ранее, центробежная сила, обусловленная вращением Земли, оказывает ничтожно малое влияние на процесс распространения волн, поэтому определяющее значение здесь будет иметь сила Кориолиса, которая представлена слагаемыми в правых частях уравнений (10), (11) и (18). На рис. 4 показаны мареограммы М\ и М2, полученные в расчетах по полной FNLD-модели с учетом вращения Земли и при отсутствии вращения (О = 0) для источников разной эффективной протяженности. Видно, что в противовес сферичности учет вращения сферы понижает максимальные положительные амплитуды волн, зафиксированные выбранными мареографами, причем влияние вращения увеличивается с расстоянием и эффективной протяженностью начального возмущения. Таким образом, сила Кориолиса может внести значимый вклад в волновую картину только при распространении длинных волн на большие расстояния.
Расчеты с большими, чем у Земли, значениями угловой скорости О показали, что совместное действие центробежной силы и силы Кориолиса оказывает заметное влия-
г\, м 0.041
о
-0.04-0.08
1 -0.2
44000 48000 52000 с
л -1
-----2
1 /> 1 У
¥ V
-0.2-
44000 48000 52000 t■, с
[ ; 1 / Л 1 , \ —1
| !/ \ \ 1 '/ \ 1 I у \ 1 -----2
\и \ \ а \ 4 т \ V /1 \ * / 1 \ \ / 1 \ \ |_ 1 *
| Л 4 1 \ 4 1 V 1 1 \ * \ 4
44000 48000 52000 t■, с
Г), м 0.04 И—
о
-0.04-0.08-0.12-
88000
-1-1
96000 с
0.4-
-0.8
/\ —1
/ ¡\ -----2
7\\ V \ * 9 \ 1 | \ 1
1 \ 1 1 \ 1 1 \ \ / ' 1 1 / | \ 4 // 1 \ \ ь
1 \ \ / -1- -
-0.4-----------
88000
96000 с
Рис. 4. Мареограммы М\ (а- в) и М2 (г- е), полученные в расчетах на основе сферической РМЬБ-модели с учетом (!) и без учета (2) силы Кориолиса при эффективной протяженности начального возмущения (а, г), Ш2 (б, д) и (в, е)
а
в
г
е
ние на волновую картину. В частности, при больших значениях П появляются мощные остаточные вихревые поля в месте расположения начального возмущения (33) и замедляется скорость распространения волн [10].
Что касается более ранних исследований, которые проводились в рамках ШКЬЮ-моделей, то в [47] отмечается, что учет силы Кориолиса изменяет максимальную амплитуду волн вплоть до 15 %, в [15] — 1.5... 2.5 %, в [20] — до 5 %. В последней работе делается также вывод о том, что влияние силы Кориолиса увеличивается с ростом эффективной ширины источника, при этом сила Кориолиса пытается удержать часть возмущения свободной поверхности около его начального положения, т. е. способствует возникновению остаточного вихревого поля [48].
3.4. Дисперсионные эффекты
Для оценки вклада частотной дисперсии в процесс распространения волн будем сравнивать результаты расчетов по полной РКЬЮ-модели и бездисперсионной КЬБШ-модели с начальным возмущением (33). Как видно из рис. 5, при использовании РКЬЮ-модели за головной волной возникает цуг волн — дисперсионный хвост из коротких волн малой амплитуды, который не воспроизводится КЬБШ-моделью.
Расчеты на основе РКЬЮ-модели показывают, что дисперсионные эффекты сильнее проявляются для более компактных начальных возмущений. Чтобы в этом убедиться, достаточно сравнить поведение мареограмм, показанных сплошными линиями 1 на рис. 4. Видно, что возмущение Ш1 генерирует дисперсионный хвост, не наблюдаемый для более протяженных источников и
Еще одной особенностью волновых полей, полученных с учетом дисперсии, является то, что при длительном распространении волн максимум амплитуды может переместиться в дисперсионный хвост, т. е. амплитуда головной волны может стать меньше амплитуды следующих за ней волн из дисперсионного хвоста (сравни линии 1 на рис. 4, а и 4, г). Такой же эффект наблюдался при использовании ШКЬБ-моделей [15, 16, 19]. Более того, в дисперсионном хвосте номер волны с максимальной амплитудой может возрастать с расстоянием [10].
Рис. 5. Волновая картина, полученная в расчетах на основе моделей РКЬБ (а) и NLSW (б). Источник ^1, £ = 5 ч
Кроме того, влияние дисперсии увеличивается с ростом пройденного волной пути. Это также видно из рис. 4. Например, для возмущения Ш2 дисперсия на ближнем мареографе М1 еще не проявилась (рис. 4, б, линия 1), а на дальнем М2 стала чуть заметней (рис. 4, д, линия 1). Что касается источника то дисперсия генерируемых им волн начинает проявляться на расстояниях, меньших, чем удаление М1. Для того чтобы точнее определить расстояние, на котором должны проявиться дисперсионные эффекты, были использованы дополнительные виртуальные мареографы М^ (г = 3 ... 6) с координатами (А^,^) такими, что долгота всех мареографов совпадала с долготой центра начального возвышения А^ = Ао = 280°, а широта принимала разные значения: ^з = -35°, = -30°, = -20°, = 0. Таким образом, эти мареографы расположены значительно ближе к источнику начального возмущения, чем М1. Некоторые мареограммы представлены на рис. 6.
Из рис. 6, а - в видно, что для менее протяженного источника дисперсия проявилась в полной мере (первая волна из дисперсионного хвоста полностью отделилась от головной) на мареографе М5, установленном на расстоянии порядка 2200 км от центра начального возмущения. Видно также, что учет дисперсии ведет к понижению амплитуды головной волны по сравнению с МЬБШ-моделью и уменьшению скорости ее движения. На этот дисперсионный эффект указывалось в одной из первых работ, посвященных численному моделированию дисперсионных волн в плановой постановке [49].
Рис. 6. Колебания свободной границы, рассчитанные по моделям FNLD (1) и NLSW (2) для источников w1 (а- в) и (г- е) и зафиксированные мареографами м3 (а), м4 (б, г), м5 (в, д) и мб (е)
Что касается более протяженного источника W2, то дисперсия волн не проявляется даже на мареографе М6, удаленном на 4400 км (рис. 6, г - е). Таким образом, влияние учета дисперсии волн уменьшается с увеличением эффективной ширины источника (см. также [20]). Общим для FNLD- и NLSW-моделей является то, что на не очень удаленных мареографах период головной волны растет с увеличением расстояния, а ее амплитуда падает.
Для того чтобы приближенно определить расстояние Ldisp, на котором могут стать заметными дисперсионные эффекты для заданной волны длиной Л, распространяющейся над горизонтальным дном h0 = const, рассмотрим модельное RLW-уравнение (Regularised Long Wave) [18]
Vt + Co Г]х = v r]xxt,
которое можно получить, например, из линеаризованной системы FNLD-уравнений в одномерном приближении, приведенной в [37]. Здесь с0 = у/gh0, v = h§/6. Для этого уравнения дисперсионное соотношение ш(к) и фазовые скорости v и vu распространения заданной волны и волны, соответствующей волновому числу , определяются по формулам
(,) = сок = Co = Co
W(fcj=1 + z/к2 , V= 1 + »(2т:/Л)2, Щ = 1 + и(2т:/Лк)2 ,
где Лк = 2т:/к.
За время t заданная волна проходит путь Ldisp = vt, а отделившаяся от нее волна с меньшей длиной Лк < Л — меньший путь L& = Vkt = LdispVk/v < Ldisp. Несовпадение пройденных путей и есть проявление дисперсии. Далее предположим, что длина отделившейся волны связана с Л соотношением Л& = аЛ (а Е (0,1)) и что дисперсия проявляется в достаточной степени к тому моменту времени, когда короткая волна полностью отделяется от длинной, т. е. выполняется равенство Lд. = Ldisp — (1 + а)Л/2. Отсюда следует формула для определения Ldisp:
Ldisp(1 — v) = ЦтЛ
или
L
disp
( _*_\
1 + у(2т/Лк )2
1
o
1 + а —z—Л.
\ 1 + и(2п/Л)2 / С учетом равенства = аЛ получаем окончательный вид выражения для Ь^р
1 / 3Л3 \ 1 { Л3 \
^ = 2(1-0) V + Ц) м 2(1-0) \ + °'152а2Ц) ■ (34)
Графики зависимостей (34) представлены на рис. 7.
Отметим, что аналог формулы (34) был получен ранее [10] на основе анализа решений уравнения Кортевега — де Вриза. Величина названа здесь длиной дисперсии.
Для того чтобы применить формулу (34) на практике, необходимо знать значения а и Л. Для источника мареограф М5 зарегистрировал отделившуюся короткую волну с периодом, примерно в 3 раза меньшим, чем период головной волны (см. рис. 6, в). Поэтому с определенным приближением можно положить а = 0.33. Что касается длины головной волны, то она в подобных оценках [20], как правило, полагается равной
Рис. 7. Зависимость длины дисперсии от длины волны Л при ао = 0.33 и Ьо = 1 км (1); 2 км (2); 3 км (3); 4 км (4); 5 км (5); 6 км (б)
эффективной протяженности источника, хотя в общем случае отождествлять длину головной волны с размером очага неправомерно [10]. При сделанных предположениях получаем значение = 1058 км, т. е. формула (34) предсказывает, что на мареографе М5 дисперсионные эффекты обязательно проявятся. В этом смысле расчетные данные не противоречат приближенной аналитической оценке. Для источника Ш2 формула (34) дает значение Ь&вр = 30 139 км, если принять, что а = 0.33 и Л = 339 км. Таким образом, выделение короткой волны из головной не произойдет не только в месте установки мареографа М6, но и в точке установки самого удаленного мареографа М2.
Есть и другие критерии для определения важности учета дисперсионных эффектов. Один из наиболее употребимых — неравенство Ра < 4, при выполнении которого считается, что дисперсия волн станет заметной. Здесь Ра — параметр Калига, приведенный, например, в [20]:
*=(£Г ■ Е ■ (35)
где Ш — протяженность источника; Ь — расстояние, пройденное головной волной. Рисунку 6, в соответствуют значения Ш = 107.3 км и Ь = 2200 км, поэтому Ра ^ 6. Таким образом, согласно критерию, связанному с параметром (35), дисперсия в точке установки мареографа М5 еще не должна проявиться, что противоречит расчетным данным. Как указано в [19], критерий Ра < 4 слишком жесткий, поэтому предлагается использовать более адекватный критерий г > 0.1, где г — "нормализованное время дисперсии",
г = 6^0 (36)
£ — время распространения волн. Полагая £ = 12 000 с, получаем г ~ 0.18 > 0.1, т.е. для источника дисперсия на мареографе М5 должна проявиться, что согласуется с нашими расчетными данными (см. рис. 6, в).
На наш взгляд, использование длины дисперсии (34) для приближенной оценки важности учета дисперсии волн является более предпочтительным, чем парамет-
ров (35), (36), поскольку значения дают наглядное представление в количественном выражении о расстоянии, на котором должны сказаться дисперсионные эффекты.
Необходимо отметить более ранние работы [14, 47, 50], в которых дисперсия исследовалась на основе неполных ШКЬБ-моделей. В [14, 50] при исследовании распространения волн цунами в Индийском океане в 2004 г. показано, что в глубоководной западной части океана различие результатов, полученных в рамках ШКЬБ- и КЬБШ-моделей, достигает 20 %, в [47] — 60 %. Расчеты распространении цунами 11 марта 2011 г. подтвердили важность учета дисперсионных эффектов: различия максимальных амплитуд в Ш^Б- и КЬБШ-моделях достигали 60 % [20].
Заключение
В работе представлен численный алгоритм для моделирования процесса распространения длинных поверхностных волн в рамках полной нелинейно-дисперсионной модели, учитывающей сферичность Земли и ее вращение. Алгоритм основан на решении расширенной системы уравнений, состоящей из эллиптического уравнения для дисперсионной составляющей давления и гиперболической системы уравнений первого порядка, которые отличаются от уравнений классической модели мелкой воды лишь наличием дополнительных слагаемых, связанных с дисперсионными добавками. Алгоритм реализован в виде явной двухшаговой схемы предиктор-корректор, на каждом шаге которой поочередно решаются задачи, полученные в результате расщепления.
На модельной задаче о распространении волн над ровным дном показано, что в земных условиях центробежную силу можно не принимать во внимание, а учет сферичности поверхности, по которой распространяются волны, может иметь определяющее значение. Влияние силы Кориолиса увеличивается с ростом размеров начального возмущения и при распространении волн на большие расстояния, при этом учет вращения сферы понижает максимальные положительные амплитуды волн. Дисперсия также проявляется при длительном распространении волн и может приводить к существенному уменьшению максимальной амплитуды головной волны, которая распадается на цуг волн меньшей длины. Но, в отличие от силы Кориолиса, дисперсия сильнее проявляется при малых размерах начального возмущения, при этом учет дисперсии приводит к замедлению скорости распространения волн.
Отметим, что на частотную дисперсию может влиять не только эффективная протяженность источника, но и его форма. Поэтому в реальных задачах о цунами необходимо дополнительное исследование с использованием реальной формы очагов. В практических задачах на первый план могут выйти также особенности рельефа дна. Ответы на эти вопросы будут даны в следующей статье авторов. Кроме того, в настоящей работе мы совсем не коснулись вопроса, связанного с теоретической оценкой численной дисперсии используемой расчетной схемы. Этой теме также будет посвящена отдельная публикация. Расчеты в небольших областях показали, что для сетки с разрешением 40 угловых секунд численная дисперсия пренебрежимо мала по сравнению с физической. Для источников с размерами, меньшими рассмотренных, приходится использовать еще более подробные сетки, что приводит к значительному росту времени решения задач в рамках нелинейно-дисперсионной модели. Поэтому в дальнейшем планируется реализация параллельной версии алгоритма для многопроцессорных ЭВМ с соответствующей модификацией алгоритмических и программных модулей.
Приложение
Приведем вывод уравнения (18) для дисперсионной составляющей давления. Для этого будем использовать компактную форму записи уравнений FNLD-модели на сфере, приведенную в работе [27]:
Ht + V ■ (Не) = 0, vt + (е ■ V)v + V = г + HjVh. (37)
Здесь v = (vi, v2) — ковариантный вектор скорости,
V\ = (П + ci)R2 sin2 В, v2 = R2c2, (38)
r = (0, r2)т, r2 = (П + C1)2 R2 sin В cos B.
Перепишем уравнение движения (37) в виде, позволяющем выразить полную производную вектора скорости
Dv = -9V„+ + r, (39)
H
где
Dv = (Dvi,Dv2)т , Dvi = (vi)t + е ■ Vvi, Dv2 = (v2)t + е ■ V^
Далее, как и в плоском случае [37], выразим дисперсионную составляющую ф давления на дне через р — дисперсионную составляющую проинтегрированного давления. Из определений (7) получаем, что
ф=2H+hq>- (40)
где
Q2 = D2h = D(Dh) = (Dh)t + е ■ V(Dh) = (ht + е ■ Vh)t + е -V (ht + е ■ Vh) =
= htt + еt ■ Vh + 2е ■Vht + е -V (е ■ Vh) = В + ct ■ Vh + е -V (е ■ Vh) и через В обозначено слагаемое, наличие которого обусловлено подвижностью дна:
В = htt + 2е ■ Vht.
Преобразуем также последнее слагаемое в выражении для Q2:
е -V (е ■ Vh) = с1 (cihx + c2he)x + с2 (cihx + c2he)в =
= clc\ h\ + (ci)2hxx + éS^he + clc2 hxe + с2 cleh\ + Séhxe + c2 c^he + (c2)2hee = = (c14 + c2c¿) hx + (c^ + c2c¡) he + c1 (c1(hx)x + c2(hx)e) + c2 (c1(h)A + c2(h)e) =
= ((е ■ V)c) ■ Vh + е ■ ((е ■ V)Vh).
Следовательно,
Q2 = В + ct ■ Vh + ((е ■ V)c) ■ Vh + е ■ ((е ■ V)Vh),
или
Q2 = В + (Dc) -Vh + е ■ ((е ■ V)Vh), (41)
где
Dc = (Dc1,Dc2)t , Dc1 = c¡ + c1^ + c24, Dc2 = ct2 + c1c2 + C24 (42) (До) ■ Vh = (с! + c1 c\ + c2c¿) h\ + (c2 + c1c2 + c2c¡) h0, c ■ [(c ■ V)Vh] = (c1)2hAA + 2c1c2hxe + (c2)2h^. Используя выражение (41), перепишем формулу (40)
ф = + \ (Dc) ■ Vh + В + о ■ ((о ■ V)Vh)
(43)
и выразим Dc через Dv. Для этого введем ковариантные компоненты метрического тензора на сфере [51] д11 = R2 sin2 9, д22 = R2 и, согласно формулам (38), запишем ковариантный вектор скорости в виде v = QG + Qc, где
G = ( Y ) ■ » = ( Y ) .
Тогда
Dv = QDG + D (Qc) = c2ÜGe + QDc + cVG^ = QDc + c2 (Q + c1) Ge, (44)
при этом
„ ( 2R2 sin в cos в N G = ^ 0 ).
Следствием равенств (44) является выражение Dc = Q 1Dv — с2 (Q + с1) Q 1Go, которое с учетом формулы (39) принимает следующий вид:
Dc = Q ->( — gVr, + ) + C-ia, (45)
где
^ 0 \ „ = 1 = 1 ,2 = ± = ± 0 g22 ), 9 911 R2sin2e, 9 922 R2'
-1 =( 911 о \
V о g22 ) ,
^ = r — C2 (Q + c1) Ge = (Q + c1) R2 sin 0 cos ) .
Подставляя полученное выражение для Dc в формулу (43), запишем
, 3 р Н ÍVp ■Vh — ф |Vh|2 \ * = ¿ + 4 { * Н + Q)
где
V fi ■Vh = gnp\hx + g22pe he = ¿ ( —+ fie he
R2 \ sin2 в
|Vh|2 = Vh ■Vh = g^h2 + /2h2 = ^ (^ +
w л
R2\ sin2 в
Q = —gVri ■Vh + c ■ [(c ■ V) Vh] + В + a ■Vh, VvVh = g11 Vxhx + Лh, = -L f + ^h,)
л2 \ sin2 в )
а
УН = (П + с1) ^ в - 2с2Нл + (П + с1)кв яп2 в
Следовательно, для дисперсионной составляющей давления на дне справедливо следу■ ющее выражение:
1 /6 р
Ф = - + + Ур -Ун
где
г = 4+ |УН|2 ■
)
46)
Теперь можно приступить к выводу уравнения для функции р. Непосредственно из определений (7) следует, что
н 3 ц
Р = + 2Ф. (47)
Преобразуем в этом выражении член Q1. Из определения операторов полной производной В и дивергенции (8) получаем, что
В (У ■ с) = (У ■ с)4 + с1 (У ■ е)Л + с2 (У ■ о), =
У ■ с + с1
1 , ( Зс\
1Л +
З
+ с2
1, ( зг\
1Л +
З
У ■ о + (с1с1Л)Л - (сЛ)2 + (с2с1)Л - + с
[(З с2)л]
З
0 + с2
( З 2)
З
(З с?)
= ((4)л + (¿сЛ + с2 4)л +
[с 1( З с2)х]в 4 ( Зс2)
-(
З
А )2 _ ¿А
СЛ) Л о■
З
+
,2 (ЗС2), " З
2 (
З
Легко видеть (см. формулу (42)), что сумма выделенных фигурными скобками слагаемых равна (Вс1)х, поэтому
В (У ■ с) = (Вс1) л + +
( Зс2)в , ( Зс1сЛ)в 12 , ( Зс2с2), , (З|9 (с2)0,
З
- с^сЛ +
З
З
с2 ( ЗС2)в с2в ( ЗС2)в 12
Зв +--;--+ (Сх)
З2
З
СЛС6> ■
Слагаемые, выделенные в этом равенстве, в сумме дают величину
З (с2 + с1сЛ + с2с2)
( З В 2)
З З
а выражение в квадратных скобках преобразуется с учетом (8) к виду
(сЛ)2 +
(^ )2
поэтому В (У ■ с)
(Вс 1)х +
( ЗВ с 2)6 З
(У ■ с)2 - 2сЛ
- (У ■ с)2 + 2с Лс2 + 2 сЛ с2 Ц + 2с2 с2в Ц + (с2)2 ^ - 2с(
(З с2)6 З
З
2 2Зв с.
'2 \ 2
Зй,
З
З
. .1 г2
Л
1
е
Л
9
Наконец, учитывая, что Jee = — J, получаем формулу
„2\2
D (V ■ c) = V ■ (Dc) — (V ■ c)2 + 2 (с\с2 — c¿c2) + 2с2 (с\ + cg) ctg0 — (с2)
Следовательно, для Q 1 имеем выражение
Q = V ■ (Dc) — 2 (V ■ c)2 + 2 (c\c2e — el с2) + 2 с2 (c\ + c¡) ctg9 — (c2)2 .
(48)
(49)
Подставляя в (47) выражение (49) и используя формулы (45), (46) для До и ф, соответственно, получаем следующее уравнение относительно зависимой переменной р:
Н
3
V fi 6Vh
fi = i2{ Vi6-1(— gVr]+IT — ж fi — Vh7 —
Q (Vfi ■ Vh)Vh
Hr
+a
— 2 (V ■ c)2 + 2 (ele2, — c¿c2) + 2c2 (4 + c2) ctg9 — (c2)2 } +
+f (я + HQ + Vv VM.
(50)
Заметим, что умножение ковариантных векторов а, Vг], V р, ^к слева на матрицу О-1 переводит эти векторы в контравариантные, через которые и определяется дивергенция вектора. Например,
V ■ V'q
a
( Jg11 Ух) х + ( Jg22 щ )в j
( Jg11a !)a + ( Jg 22a2)
1
R2 sin0 \sin61
1
J R2 sin 9 V sin 9
Учитывая это замечание, а также равенство
6
Gn?+(г] в sin0)*)=Аг],
(^ + (a2 Ы).) .
6Vh V ■ ( — fi
)
\Н2rr ) Н2г перепишем уравнение (50) в компактной форме
Lfi = F,
Vfi■Vh + 6fiV
\H 2r)
(51)
(52)
где
Lfi = V
V fi (V fi ■ Vh) Vh
H
Hr
— 6 fi
2 r — 3 í Vh \ H3 r + ^ \H4-)
Р = д^ + V ■ - а) - ^ + 2 (V ■ о)2 - 2 (^2 - с1с2) - 2с2 (с* + <$) ctg0 + (с2)2
Используя формулу вида (51) для дивергенции ковариантного вектора, выпишем выражение для оператора Ь через частные производные:
Lfi
1
1 (fiл Vfi ■ Vh
R2 sin 9
sin 9 \ H
Hr
h;
+
{
, 2 r — 3
— 6^Н~з— +
1
R2 sin 9
1
h\
sin0 \H2r/\
H
+
Vfifi ■ Vh, Hr
he) sin9
A
После домножения обеих частей уравнения (52) на R2 sin б1 и перехода к физическим компонентам скорости u, v получаем следующее уравнение для р:
1 /px Vp •Vh
Vsin9\H
Hr
h
+
pe Vp •Vh.
H
H
?) si
he) sin 9
—6 p
LH3
2 (r — 3) R2 sin 9 / hx \ / he
H2r sin 9 ) x \H2r
e
F,
где F = Fi + F2 + F3,
Fi
1 ( ^—h
[9 Vx + — hx — ai
sin
+
sin^f g щ + —he — (¿2
(53)
F2 = — 6—R2 sin 9, 2 Hr '
F3
sin
ux + (v sin^),? — 2(uxve — vxue) — 2(uv)x ctg0 — (v2 cos 0)e,
Q = — gVr] • Vh +
1
R2 sin sin
u
hxx + 2uw hxe + v2hee sin 0 ) + В + a • Vh
in
(u V \
RÜTe hx +r4
ai = — 2uw cos 9 — fvR sin 9, a2 = u ctg 9 + fuR + П R sin 9 cos 9
ah
aihx R2 sin2 9
h
+ (u2ctg 9 + /uR) ^ + heП2 sin 9 cos 9,
f = 2П cos в — параметр Кориолиса и отдельно выделены слагаемые с множителем П2, наличие которых обусловлено действием центробежной силы.
Замечание. Положим П = 0 ив некоторой окрестности фиксированной точки ( А*, 9*) введем невырожденное преобразование координат
X = R (А — А*) sin9*, у = —R (9 — 9*), а также соответствующие компоненты вектора скорости
u = X = RA sin 9* = ua, a
sin
sin
v = y = — R9 = — v.
Предполагая, что рассматриваемая окрестность мала в направлении широты, т. е. мала величина 8 = 9 - (9*, переходя в уравнении (53) к новым переменным х, у, и, V и пренебрегая в преобразованном уравнении членами, имеющими порядок 0(8) или О (1/R), получаем известное уравнение на плоскости для дисперсионной составляющей р [37]. Это является косвенным подтверждением правильности вывода уравнения (53).
x
e
x
e
2
Благодарности. Исследование выполнено при поддержке Российского научного фонда (проект № 14-17-00219).
Список литературы / References
[1] Imamura, F. Simulation of wave-packet propagation along sloping beach by TUNAMI-code. Long-wave Runup Models / Eds H. Yeh, P. Liu, C. Synolakis. Singapore: World Scientific, 1996. P. 231-241.
[2] Зайцев А.И., Куркин А.А., Левин Б.В., Пелиновский Е.Н., Ялчинер А., Троицкая Ю.И., Ермаков С.А. Моделирование распространения катастрофического цунами (26 декабря 2004 г.) в Индийском океане // Доклады РАН. 2005. Т. 402, № 3. С. 388-392. Zaitsev, A.I., Kurkin, A.A., Levin, B.V., Pelinovsky, E.N., Yalciner, A., Troitskaya, Yu.I., Ermakov, S.A. Numerical simulation of catastrophic tsunami propagation in the Indian Ocean // Doklady Earth Sciences. 2005. Vol. 402, No. 4. P. 614-618.
[3] Titov, V.V., Synolakis, C.E. Numerical modeling of long wave runup using VTCS-3. Longwave Runup Models / Eds H. Yeh, P. Liu, C. Synolakis. Singapore: World Scientific, 1996. P. 242-248.
[4] Wei, Y., Bernard, E., Tang, L., Weiss, R., Titov, V., Moore, C., Spillane, M., Hopkins, M., Kanoglu, U. Real-time experimental forecast of the Peruvian tsunami of August 2007 for U.S. coastlines // Geophysical Research Letters. 2008. Vol. 35. L04609.
[5] Tang, L., Titov, V.V., Bernard, E., Wei, Y., Chamberlin, C., Newman, J.C., Mofjeld, H., Arcas, D., Eble, M., Moore, C., Uslu, B., Pells, C., Spillane, M.C., Wright, L.M., Gica, E. Direct energy estimation of the 2011 Japan tsunami using deep-ocean pressure measurements // Journal of Geophysical Research. 2012. Vol. 117. C08008.
[6] Shokin, Yu.I., Babailov, V.V., Beisel, S.A., Chubarov, L.B., Eletsky, S.V., Fedotova, Z.I., Gusyakov, V.K. Mathematical modeling in application to regional tsunami warning systems operations. Notes on Numerical Fluid Mechanics and Multidisciplinary Design. Vol. 101. Computational Science and High Performance Computing III. Berlin: Springer-Verlag, 2007. P. 52-69.
[7] Федотова З.И. О применении разностной схемы Мак-Кормака для задач длинноволновой гидродинамики // Вычисл. технологии. 2006. Т. 11. Специальный выпуск, посвященный 85-летию со дня рождения Н.Н. Яненко, ч. 2. С. 53-63.
Fedotova, Z.I. On application of the MacCormack difference scheme for problems of longwave hydrodynamics // Computational Technologies. 2006. Vol. 11. Special issue, devoted to N.N. Yanenko's 85-th anniversary, pt 2. P. 53-63. (in Russ.)
[8] Gusyakov, V.K., Fedotova, Z.I., Khakimzyanov, G.S., Chubarov, L.B., Shokin,
Yu.I. Some approaches to local modelling of tsunami wave runup on a coast // Russian Journal of Numerical Analysis and Mathematical Modelling. 2008. Vol. 23, No. 6. P. 551-565.
[9] Beisel, S.A., Chubarov, L.B., Dutykh, D., Khakimzyanov, G.S., Shokina, N.Yu.
Simulation of surface waves generated by an underwater landslide in a bounded reservoir // Russian Journal of Numerical Analysis and Mathematical Modelling. 2012. Vol. 27, No. 6. P. 539-558.
[10] Пелиновский Е.Н. Гидродинамика волн цунами. Н. Новгород: Ин-т прикл. физики РАН, 1996. 276 c.
Pelinovsky, E.N. Tsunami Wave Hydrodynamics. Nizhny Novgorod: Institute Applied Physics Press, 1996. 276 p. (in Russ.)
[11] Dalrymple, R.A., Grilli, S.T., Kirby, J.T., Watts, P. Tsunamis and challenges for accurate modeling // Oceanography. 2006. Vol. 19, No. 1. P. 142-151.
[12] Murty, T.S., Rao, A.D., Nirupama, N., Nistor, I. Numerical modelling concepts for tsunami warning systems // Current Science. 2006. Vol. 90, No. 8. P. 1073-1081.
[13] Glimsdal, S., Pedersen, G.K., Atakan, K., Harbitz, C.B., Langtangen, H.P., Lovholt, F. Propagation of the Dec. 26, 2004, Indian Ocean Tsunami: Effects of dispersion and source characteristics // International Journal of Fluid Mechanics Research. 2006. Vol. 33, No. 1. P. 15-43.
[14] Horrillo, J., Kowalik, Z., Shigihara, Y. Wave dispersion study in the Indian Ocean-tsunami of December 26, 2004 // Marine Geodesy. 2006. Vol. 29. P. 149-166.
[15] Lovholt, F., Pedersen, G., Gisler, G. Oceanic propagation of a potential tsunami from the La Palma Island // Journal of Geophysical Research. 2008. Vol. 113. C09026.
[16] Lovholt, F., Pedersen, G., Glimsdal, S. Coupling of dispersive tsunami propagation and shallow water coastal response // Open Oceanography Journal. 2010. Vol. 4. P. 71-82.
[17] Grue, J., Pelinovsky, E.N., Fructus, D., Talipova, T., Kharif, C. Formation of undular bores and solitary waves in the Strait of Malacca caused by the 26 December 2004 Indian Ocean tsunami // Journal of Geophysical Research. 2008. Vol. 113. C05008.
[18] Peregrine, D.H. Calculations of the development of an undular bore // Journal of Fluid Mechanics. 1966. Vol. 25, pt 2. P. 321-331.
[19] Glimsdal, S., Pedersen, G.K., Harbitz, C.B., Lovholt, F. Dispersion of tsunamis: does it really matter? // Natural Hazards and Earth System Science. 2013. Vol. 13. P. 1507-1526.
[20] Kirby, J.T., Shi, F., Tehranirad, B., Harris, J.C., Grilli, S.T. Dispersive tsunami waves in the ocean: Model equations and sensitivity to dispersion and Coriolis effects // Ocean Modelling. 2013. Vol. 62. P. 39-55.
[21] Grilli, S.T., Harris, J.C., Tajalli Bakhsh, T., Masterlark, T.L., Kyriakopoulos, C., Kirby, J.T., Shi, F. Numerical simulation of the 2011 Tohoku tsunami based on a new transient FEM co-seismic source: comparison to far- and near-field observations // Pure and Applied Geophysics. 2013. Vol. 170, No. 6-8. P. 1333-1359.
[22] Nwogu, O. Alternative form of Boussinesq equations for nearshore wave propagation // Journal of Waterway, Port, Coastal, and Ocean Engineering. 1993. Vol. 119, No. 6. P. 618-638.
[23] Lovholt, F., Pedersen, G. Instabilities of Boussinesq models in non-uniform depth // International Journal for Numerical Methods in Fluids. 2009. Vol. 61. P. 606-637.
[24] Shi, F., Kirby, J.T., Harris, J.C., Geiman, J.D., Grilli, S.T. A high-order adaptive time-stepping TVD solver for Boussinesq modelling of breaking waves and coastal inundation // Ocean Modelling. 2012. Vol. 43-44. P. 36-51.
[25] Shi, F., Kirby, J.T., Tehranirad, B. Tsunami benchmark results for spherical coordinate version of FUNWAVE-TVD (Version 2.0). Research Report No. CACR-12-02. Center for Applied Coastal Research, University of Delaware, 2012. 39 p.
[26] Fedotova, Z.I., Khakimzyanov, G.S. Nonlinear-dispersive shallow water equations on a rotating sphere // Russian Journal of Numerical Analysis and Mathematical Modelling. 2010. Vol. 25, No. 1. P. 15-26.
[27] Fedotova, Z.I., Khakimzyanov, G.S. Nonlinear-dispersive shallow water equations on a rotating sphere and conservation laws // Journal of Applied Mechanics and Technical Physics. 2014. Vol. 55, No. 3. P. 404-416.
[28] Шокин Ю.И., Федотова З.И., Хакимзянов Г.С. Иерархия моделей гидродинамики длинных поверхностных волн // Докл. РАН. 2015. Т. 462, № 2. С. 168-172.
Shokin, Yu.I., Fedotova, Z.I., Khakimzyanov, G.S. Hierarchy of nonlinear models of the hydrodynamics of long surface waves // Doklady Physics. 2015. Vol. 60, No. 5. P. 224-228.
[29] Гусев О.И. Модуль расчета поверхностных волн NLDSW_sphere. Роспатент. Свид. о гос. рег. прогр. для ЭВМ. № 2015616421 от 09.06.2015 г.
Gusev, O.I. Program code NLDSW_sphere for calculation of surface waves. Rospatent. Svidetel'stvo o gosudarstvennoy registratsii programmy dlya EVM. No. 2015616421 ot 09.06.2015. (in Russ.)
[30] Fedotova, Z., Khakimzyanov, G. Full nonlinear dispersion model of shallow water equations on a rotating sphere // Journal of Applied Mechanics and Technical Physics. 2011. Vol. 52, No. 6. P. 865-876.
[31] Lynett, P.J., Liu, P.L.-F. A numerical study of submarine-landslide-generated waves and run-up // Proceedings of the Royal Society of London. Series A. 2002. Vol. 458. P. 2885-2910.
[32] Chubarov, L.B., Eletskij, S.V., Fedotova, Z.I., Khakimzyanov, G.S. Simulation of surface waves generation by an underwater landslide // Russian Journal of Numerical Analysis and Mathematical Modelling. 2005. Vol. 20, No. 5. P. 425-437.
[33] Shokin, Yu.I., Fedotova, Z.I., Khakimzyanov, G.S., Chubarov, L.B., Beisel, S.A.
Modelling surfaces waves of generated by a moving landslide with allowance for vertical flow structure // Russian Journal of Numerical Analysis and Mathematical Modelling. 2007. Vol. 22, No. 1. P. 63-85.
[34] Черевко А.А., Чупахин А.П. Уравнения мелкой воды на вращающейся притягивающей сфере. 1. Вывод и общие свойства // ПМТФ. 2009. Т. 50, № 2. C. 24-36. Cherevko, A.A., Chupakhin, A.P. Equations of the shallow water model on a rotating attracting sphere. 1. Derivation and General Properties // Journal of Applied Mechanics and Technical Physics. 2009. Vol. 50, No. 2. P. 188-198.
[35] Fedotova, Z.I., Khakimzyanov, G.S. Shallow water equations on a movable bottom // Russian Journal of Numerical Analysis and Mathematical Modelling. 2009. Vol. 24, No. 1. P. 31-41.
[36] Гусев О.И. Алгоритм расчета поверхностных волн над подвижным дном в рамках плановой нелинейно-дисперсионной модели // Вычисл. технологии. 2014. Т. 19, № 6. С. 19-40. Gusev, O.I. Algorithm for surface waves calculation above a movable bottom within the frame of plane nonlinear dispersive model // Computational Technologies. 2014. Vol. 19, No. 6. P. 1940. (in Russ.)
[37] Гусев О.И., Шокина Н.Ю., Кутергин В.А., Хакимзянов Г.С. Моделирование поверхностных волн, генерируемых подводным оползнем в водохранилище // Вычисл. технологии. 2013. Т. 18, № 5. С. 74-90.
Gusev, O.I., Shokina, N.Yu., Kutergin, V.A., Khakimzyanov, G.S. Numerical modelling of surface waves generated by underwater landslide in a reservoir // Computational Technologies. 2013. Vol. 18, No. 5. P. 74-90. (in Russ.)
[38] Шокин Ю.И., Бейзель С.А., Гусев О.И., Хакимзянов Г.С., Чубаров Л.Б., Шокина Н.Ю. Численное исследование дисперсионных волн, возникающих при движении подводного оползня // Вестник ЮУрГУ. 2014. Т. 7, № 1. С. 121-133.
Shokin, Yu.I., Beisel, S.A., Gusev, O.I., Khakimzyanov, G.S., Chubarov, L.B., Shokina, N.Yu. Numerical modelling of dispersive waves generated by landslide motion // Bulletin of the South Ural State University. 2014. Vol. 7, No. 1. P. 121-133. (in Russ.)
[39] Khakimzyanov, G.S., Gusev, O.I., Beisel, S.A., Chubarov, L.B., Shokina, N.Yu.
Modelling of tsunami generated by submarine landslides in the Black Sea // Russian Journal of Numerical Analysis and Mathematical Modelling. 2015. Vol. 30, No. 4. (in print).
[40] Ладыженская О.А., Уральцева Н.Н. Линейные и квазилинейные уравнения эллиптического типа. М.: Наука, 1973. 576 с.
Ladyzhenskaya, O.A., Ural'tseva, N.N. Linear and Quasilinear Elliptic Equations. Moscow: Nauka, 1973. 576 p. (in Russ.)
[41] Гусев О.И. Об алгоритме расчета поверхностных волн в рамках нелинейно-дисперсионной модели на подвижном дне // Вычисл. технологии. 2012. Т. 17, № 5. С. 46-64. Gusev, O.I. On an algorithm for surface waves calculation within the framework of nonlinear dispersive model with a movable bottom // Computational Technologies. 2012. Vol. 17, No. 5. P. 46-64. (in Russ.)
[42] Хакимзянов Г.С., Шокин Ю.И., Барахнин В.Б., Шокина Н.Ю. Численное моделирование течений жидкости с поверхностными волнами. Новосибирск: Изд-во СО РАН, 2001. 394 с.
Khakimzyanov, G.S., Shokin, Yu.I., Barakhnin, V.B., Shokina, N.Yu. Numerical Simulation of Fluid Flows with Surface Waves. Novosibirsk: FUE "Publishing House SB RAS", 2001. 394 p. (in Russ.)
[43] Хакимзянов Г.С., Шокина Н.Ю. Метод адаптивных сеток для одномерных уравнений мелкой воды // Вычисл. технологии. 2013. Т. 18, № 3. С. 54-79.
Khakimzyanov, G.S., Shokina, N.Yu. Adaptive grid method for one-dimensional shallow water equations // Computational Technologies. 2013. Vol. 18, No. 3. P. 54-79. (in Russ.)
[44] Shokina, N.Yu. To the problem of construction of difference schemes on movable grids // Russian Journal of Numerical Analysis and Mathematical Modelling. 2012. Vol. 27, No. 6. P. 603-626.
[45] Synolakis, C.E., Bernard, E.N., Titov, V.V., Kanoglu, U., Gonzalez, F.I. Validation and verification of tsunami numerical models // Pure and Applied Geophysics. 2008. Vol. 165. P. 2197-2228.
[46] Horrillo, J., Grilli, S.T., Nicolsky, D., Roeber, V., Zhang, J. Performance benchmarking tsunami models for NTHMP's inundation mapping activities // Pure and Applied Geophysics. 2015. Vol. 170, No. 3-4. P. 1333-1359.
[47] Dao, M.H., Tkalich, P. Tsunami propagation modelling — a sensitivity study // Natural Hazards and Earth System Science. 2007. Vol. 7. P. 741-754.
[48] Носов М.А., Нурисламова Г.Н., Мошенцева А.В., Колесов С.В. Остаточные гидродинамические поля при генерации цунами землетрясением // Известия РАН. Физика атмосферы и океана. 2014. Т. 50, № 5. С. 591-603.
Nosov, M.A., Nurislamova, G.N., Moshenceva, A.V., Kolesov, S.V. Residual hydrodynamic fields after tsunami generation by an earthquake // Izvestiya RAN. Fizika Atmosfery i Okeana. 2014. Vol. 50, No. 5. P. 591-603. (in Russ.)
[49] Shokin, Yu.I., Chubarov, L.B. The numerical modelling of long wave propagation in the framework of non-linear dispersion models // Computers and Fluids. 1987. Vol. 15, No. 3. P. 229-249.
[50] Grilli, S.T., Ioualalen, M., Asavanant, J., Shi, F., Kirby, J.T., Watts, P. Source constraints and model simulation of the December 26, 2004, Indian Ocean tsunami // Journal of Waterway, Port, Coastal, and Ocean Engineering. 2007. Vol. 133. P. 414-428.
[51] Федотова З.И., Хакимзянов Г.С. Нелинейно-дисперсионные уравнения мелкой воды на вращающейся сфере // Вычисл. технологии. 2010. Т. 15, № 3. С. 135-145. Fedotova, Z.I., Khakimzyanov, G.S. Nonlinear dispersive equations of the shallow water on the rotating sphere // Computational Technologies. 2010. Vol. 15, No. 3. P. 135-145. (in Russ.)
Поступила в 'редакцию 27 апреля 2015 г.
32
Q. H. ryceB, r. C. XaKHM33HOB
Numerical simulation of long surface waves on a rotating sphere within the framework of the full nonlinear dispersive model
Gusev, Oleg I., Khakimzyanov, Gayaz S.*
Institute of Computational Technologies SB RAS, Novosibirsk, 630090, Russia * Corresponding author: Khakimzyanov, Gayaz S., e-mail: khak@ict.nsc.ru
The paper examines the sensitivity of long surface wave propagation to Earth sphericity, Coriolis force, centrifugal force and frequency dispersion.
For a numerical simulation, we propose the algorithm based on the partitioning of fully nonlinear dispersive equations on a rotating sphere based on a uniformly elliptic equation for the dispersion component of pressure and the hyperbolic system of shallow water equations for the momentum. The momentum equations yield the modified source term in the first approximation in the right hand side. These subproblems resulting from the partitioning are solved on each step of the explicit implemented two-step predictor-corrector scheme. The system of difference equations approximating the elliptic subproblem with the second order is constructed with use of integro-interpolation method and is solved by the SOR method.
A model domain includes the most part of the Pacific Ocean. The function, which defines distribution of depth from the still water surface was set to be constant. Gaussian perturbations of free surface with different effective width served as idealized sources of waves. Numerical results show that in terrestrial conditions centrifugal force can be neglected for all the considered sources, but other effects can change the wave pattern considerably. Thus, sphericity increases the maximum wave amplitude, but Coriolis force and dispersion decrease. Besides that, the influence of sphericity and Coriolis force increases with an effective source width, while the influence of dispersion decreases. Wave propagation distance enhances the influence of all these effects. The approximate formula for the quick assessment of the propagation distance sufficient for demonstration of the dispersion effects is derived and its good agreement with calculations is shown.
Keywords: rotating sphere, shallow water, long surface waves, nonlinear dispersive equations, numerical simulation, dispersion, Coriolis force.
Acknowledgements. The study was supported by the Russian Science Foundation (Project No. 14-17-00219).
Received 27 April 2015
© ICT SB RAS, 2015
