CC BY
24
Вычислительные технологии
Том 19, № 6, 2014
Алгоритм расчёта поверхностных волн над подвижным дном в рамках
__u u __u >к
плановой нелинеино-дисперсионнои модели
О. И. Гусев
Разработан численный алгоритм, основанный на разбиении на эллиптическую и гиперболическую части плановой нелинейно-дисперсионной модели на подвижном дне. Полученные решения сравнивались с расчётами по модели мелкой воды и экспериментальными данными на задачах о генерации волн оползнем и их распространении.
Ключевые слова: подводный оползень, поверхностные волны, цунами, уравнения мелкой воды, нелинейно-дисперсионные уравнения, численное моделирование, конечно-разностная схема.
Введение
Оползни могут представлять опасность для прибрежных объектов, располагающихся как вблизи океана (см., например, [1]), так и на берегах крупных внутренних водоёмов [2,3]. В последнее время проводится много работ по изучению процесса образования волн оползнем средствами математического моделирования (например, [4]), выполнены также экспериментальные исследования [5-7], которые могут служить для верификации моделей.
Для конкретного водоёма параметры оползня можно оценить не точно, а только в некотором диапазоне. Это порождает большое количество расчётов с перебором всех возможных вариантов, поэтому трёхмерные гидродинамические модели [8,9], требующие значительных вычислительных ресурсов, вряд ли подойдут для этих целей. Тем не менее расчёты по ним могут служить неким эталоном, так как они хорошо согласуются с экспериментальными данными [5,6]. Среди приближённых моделей теории мелкой воды можно выделить полные нелинейно-дисперсионные НЛД-модели [10-14], полученные без предположения о малости амплитуды волн, НЛД-модели типа модели Буссинеска [15,16] и бездисперсионные нелинейные НЛ-модели первого приближения. В работе [17] было показано, что полные НЛД-модели дают результаты более близкие к таковым трёхмерной модели потенциальных течений, чем модели типа модели Буссинеска и НЛ-модели.
В НЛД-уравнениях содержатся производные высокого порядка от искомых функций, прямая аппроксимация которых приводит к сложной разностной задаче, трудно поддающейся исследованию. Поэтому исходные системы уравнений разбиваются на части, решаемые поочерёдно. Для некоторых моделей типа модели Буссинеска в одномерном случае такой подход описан, например, в работе [18]. Статья [19] посвящена выводу
1 Институт вычислительных технологий СО РАН, Новосибирск, Россия
Контактный e-mail: [email protected]
Исследование выполнено при поддержке Российского научного фонда (проект № 14-17-00219).
полной НЛД-модели Вея на неподвижном дне и описанию разностной схемы для неё типа предиктор-корректор четвёртого порядка аппроксимации Адамса — Башфорта — Моултона. Похожий алгоритм для модели Лью — Лайнетта применялся в [14] на подвижном дне. Существенным недостатком таких алгоритмов является невозможность применения большого опыта, накопленного исследователями при численном решении гиперболических систем уравнений. В [9] описано разбиение полных НЛД-уравнений Железняка — Пелиновского на стационарном дне, в результате которого получаются части эллиптического и гиперболического типа, что позволило построить для этих частей хорошо изученные численные алгоритмы. Для подвижного дна подобное разбиение было реализовано только в одномерном случае [20-22].
В настоящей работе представлен численный алгоритм для решения полных плановых НЛД-уравнений Железняка — Пелиновского на подвижном дне, основанный на разбиении системы на гиперболическую и эллиптическую части. Гиперболическая задача в нашем случае отличается от таковой классической модели мелкой воды только правой частью, поэтому для неё используется зарекомендовавшая себя схема типа предиктор-корректор [23]. Подробно описан алгоритм для эллиптической части, основанный на интегро-интерполяционном методе. Предложены граничные условия свободного прохода волн типа волн Зоммерфельда для расширенной системы. Проведены сравнения расчётов по НЛД- и НЛ-моделям на задаче о распаде начального возмущения над ровным дном, а также с экспериментальными данными [7] в задаче о сходе двумерного твёрдого оползня. На основе сравнений обсуждаются важность учёта дисперсии в рассмотренных задачах и границы применимости моделей. Для модельной акватории в виде части реки исследуется влияние учёта пространственной неоднородности оползня. В задаче о взаимодействии волны с островом проводится сравнения расчётов по НЛД-модели с экспериментальными данными [24], показана возможность аппроксимации границ острова ломаной кривой для качественного совпадения волновых картин.
1. Постановка задачи
Пусть слой несжимаемой жидкости ограничен снизу подвижным дном г = -к(х,у,1), а сверху — свободной границей г = п(х,у, £), где £ — время, х, у, г — координаты точки в декартовой системе координат Охуг, ось Ог которой направлена вертикально, а координатная плоскость Оху совпадает с невозмущённой свободной поверхностью. В моделях мелкой воды искомыми величинами являются полная глубина слоя жидкости Н = п + Ь > 0 и вектор и = (и, у), связанный некоторым образом с вектором скорости трёхмерного течения. Если в качестве и использовать осреднённую по глубине горизонтальную составляющую, то НЛД-уравнения мелкой воды второго приближения принимают следующий вид [25]:
где р означает давление в НЛД-модели, проинтегрированное по глубине слоя жидкости, р0 — давление на дне:
Н + V ■ (Ни) = 0,
(1)
. . Ур Ро уу 7
(2)
д — ускорение свободного падения, р, ф — дисперсионные составляющие р и р0:
Н3 И2 Н2
У = -3- Я! + - Я2, ф = — + -я 2,
д
Я1 = и) - (V- и)2, Я2 = Б2й, Б = + и ■ V.
Форма записи уравнения движения (2) в терминах скорость — давление позволяет заметить аналогию с неоднородным уравнением движения для течения идеального газа. В газовой же динамике (а также в гидродинамике вязкой несжимаемой жидкости) для численного решения часто используется подход с выделением задачи определения давления в качестве одной из "простых задач". В НЛД-модели гидростатическая составляющая проинтегрированного давления вычисляется через полную глубину слоя жидкости, поэтому "простую задачу" сформулируем только для дисперсионной составляющей р, уравнение для которой имеет вид [21]
V•( V(Из ^-.(Ий) = * (3)
где
л „ ( „ RVh\ 6R ,2
F = V ■ i gVn + ^^ ) - Hr + 2(V ' u)2 - 2(uxVy - UyVx),
\2
r = 4 + (Vh)2, R = -gVn ■ Vh + u ■ ((u ■ V)Vh) + B, B = htt + 2(u ■ Vht).
Уравнение (3) является равномерно эллиптическим, не содержит производных по времени от компонент вектора скорости и полной глубины, и если коэффициент при p положителен (например, при h = h0 = const), то для нахождения численного решения можно построить разностные схемы с положительно определёнными операторами. Для сильно неровного дна коэффициент при p в некоторых точках может стать отрицательным. Анализ таких случаев выполнен в работе [20].
Введение новой зависимой переменной p позволяет записать уравнение движения (2) НЛД-модели в виде
/ ^ ^ Vp - фVh
ut + (u ■ V)u + gVn = У H -, (4)
при этом дисперсионная составляющая давления на дне связана с функциями p, H и u выражением
ф = + HR + Vp -Vh^. (5)
Таким образом, возникает вторая "простая задача", которая состоит в решении системы уравнений гиперболического типа (2), (4), аналогичной системе бездисперсионных уравнений мелкой воды и отличающейся от последней лишь правой частью. В настоящей работе для её численного решения используется явная схема предиктор-корректор [23], хорошо зарекомендовавшая себя при исследовании волновых процессов в рамках бездисперсионной модели мелкой воды. При построении численного алгоритма также будет использоваться запись уравнения движения (4) в дивергентном виде
H2
(Hu)t + V ■ (Huu) + gV — = gHVh + Vp - фVh. (6)
Для контроля вычислений можно применять закон изменения полной энергии Е в НЛД-модели [25]
(НЕ\ + V ■ \пИ(Е + Р)] = -рокг, (7)
где
Н
^ и ■ и И2 Н_ (Бк)2 д(Н - 2к)
Е = + и)2 + — (V ■ и)Бк + +
2 о 2 2 2
В случае стационарного дна (кг = 0) уравнение изменения энергии (7) принимает дивергентный вид и выражает собой закон сохранения полной энергии.
Представленные уравнения дополняются краевыми условиями, которые могут иметь различный вид на разных частях границы Г области решения П. В том случае, когда поверхностные волны рассматриваются в бассейне с вертикальными непроницаемыми стенками, на границе Г ставится условие непротекания
и ■ п|г = 0, (9)
где п — вектор внешней нормали к Г. Тогда краевое условие для функции р на прямолинейном участке границы Г примет вид
др дк HVк ^р + 6р тт( дп Ядк\
и ( д^т + -тг . (10)
дп дп Иг \ дп г дп ^
Оно получается после умножения уравнения движения (4) скалярно на п и учёта граничного условия (9) и выражения (5) для ф.
2. Алгоритм численного решения
В настоящем разделе приводятся подробное описание конечно-разностной схемы, предназначенной для решения задачи (3), (10) для р, и некоторые замечания об алгоритме решения всей задачи.
Декомпозиция исходной системы НЛД-уравнений на две "простые задачи" позволяет предложить следующий алгоритм расчёта на текущем слое по времени. На шаге предиктор сначала определяется функция рп из конечно-разностного уравнения, аппроксимирующего (3) с коэффициентами, вычисленными по известным значениям Нп и ип с п-го слоя по времени. Затем решается гиперболическая система уравнений (1), (4) для нахождения величин И * и и*, при этом в правой части уравнения движения используются известные значения рп, фп, Ип и ип. Далее вновь численно решается уравнение (3), в коэффициентах которого берутся величины, вычисленные на предикторе. Найденные на этом этапе значения р* и ф* используются затем на шаге корректор для определения окончательных значений Нп+1 и ип+1 путём численного решения системы уравнений гиперболического типа с дивергентной формой левой части (1), (6).
2.1. Аппроксимация уравнения для дисперсионной составляющей давления
В практических задачах, связанных с изучением волн цунами, область течения ограничена береговой линией весьма сложной формы. Использование регулярных криволинейных сеток, адаптирующихся к геометрии такой "изрезанной" границы, связано с большими алгоритмическими трудностями, в частности, с необходимостью применения многоблочных криволинейных сеток и их адаптацией к каждой конкретной задаче.
Более универсальными, хотя возможно и менее точными, могут оказаться алгоритмы расчёта на прямоугольных равномерных сетках.
Предположим, что область решения П содержится в некотором прямоугольнике [0,Ь1 ] х [0,Ь2], покрытом прямоугольной равномерной сеткой с узлами х (] = (]1 ,]2), 31 = 0,..., N, ]2 = 0,..., N) и шагами к1 = Ь1 /И1, к2 = Ь2/И2 в направлениях осей Ох и О у соответственно. Далее предполагается, что все криволинейные участки границы Г аппроксимированы ломаными, звенья которых параллельны осям координат и проходят по координатным линиям сетки. Для модифицированной границы и области, ограниченной этой модифицированной границей, будем использовать прежние обозначения Г и П соответственно. Расчётная сетка П^ = П^ У Г^ состоит из внутренних х Е П^ С П и граничных х Е Г^ С Г узлов. Множество узлов сетки разбивается в соответствии с их типом на непересекающиеся классы. Внутренним узлам приписан тип 0, узлы х Е Г^ могут иметь тип с 1-го по 12-й в зависимости от того, как два звена ломаной, представляющей границу Г, стыкуются в данном конкретном узле. Например, если оба звена параллельны оси Оу и внешняя нормаль к ним направлена против оси Ох, то узлы имеют тип 1, если по направлению этой оси — тип 3. Для угловых узлов, в которых стыкуются два звена, одно из которых параллельно оси Ох, а другое — оси Оу, и при этом звенья образуют между собой прямой внутренний (по отношению к П) угол, типы узлов могут быть равными 5, 6, 7, 8, а для внешних прямых углов — соответственно 9,10,11,12 (рис. 1, б).
Для сокращения изложения перепишем уравнение (3) в следующем виде:
дФ1 дФ2 7
"Б- + —--ко р = Е,
дх
ду
где введены обозначения
1 др др Ф1 = ¿11^ + к12^~,
дх ду
2 др др Ф = ¿12^" + ¿22^-,
дх ду
(11)
(12)
а
Рис. 1. Область течения ^ с границей Г, состоящей из ломаных, звенья которых параллельны осям координат (а), и типы узлов расчётной сетки (б)
"11
4 + h2y ,
12
ко = 6( ^ + ) +12koo, koi
dx
dy
hx hy
hx
H 2r
k22 = ko2
4 + hX
Hr
hy
H 2r
oo
r — 3 H3r '
;i3) '14)
fi = gnx +--hx,
r
^ dfi д ¡2 6R F = f + f — 6R + 2L,
от dy Hr R
¡2 = gny + ~hy' L = (V • u)2 — (uxVy — UyVx)
'15)
Разностную схему для задачи (3), (10) построим с помощью интегро-интерполяци-онного метода [9]. Пусть Xj Е П^ — внутренний узел расчётной сетки. Тогда он является общей вершиной четырёх соседних ячеек. Обозначим через ABCD контур прямоугольника с вершинами в центрах этих ячеек (рис. 2, а). Интегрируя уравнение (11) по выбранному контуру и применяя формулу Грина, получаем
У ФМу — J ФМу + J Ф20,х — J Ф2dx = JJ k0ipdxdy + JJ Fdxdy. (16)
BC AD DC AB ABCD ABCD
Для приближённого вычисления интегралов по сторонам контура будем использовать квадратурную формулу трапеций. Для примера выпишем результат для первого интеграла из (16):
ФЫу ^ —-Н2
2
BC
= Т (k»(B)|(B) + ki2(B)|(B) + kii(C)|(C) + ki2(C)|(C)
у I kii(B)
V3 + V6 — Vo — V2 2hi
+ ki2(B)
V3 + Vo — V6 — V2
i 1 ín\^3 + V7 — Vo — V4 , (n,
+kii(C)-—--+ ki2(C)
2h1
2h2
V4 + V7 — Vo — V3 2h2
'17)
8
Í2—
D 1 W
N C -----T
E
S B 2
J 2
8 4
D TC
i л! 1 0 Б .
J 2—»
Jl
D 1 W N 1
0 j
A 5 S B 2
f
а
в
4
8
4
7
3
3
7
ü
5
6
6
3
Рис. 2. Контур интегрирования и шаблон разностного уравнения для р во внутреннем узле (а) и в граничном узле типа 2 (б) и 11 (в)
Здесь использована стандартная локальная нумерация узлов (см. рис. 2, а), соседних с х. Предполагается, что коэффициенты кц, к12, к22, коо, ко1, ко2, а также выражения Я и г вычисляются в центрах ячеек. При этом в центрах ячеек первые производные по х и у, встречающиеся в указанных коэффициентах и в правой части уравнения (11), аппроксимируются по формулам, аналогичным разностным производным от функции р в (17), например,
Ои.^* и3 + и7 — ио — Щ Щ + и7 — щ — и3
аи(С) м--■ оу,(С) м-2-2-• (18)
Теперь выпишем подобные аппроксимации для двойных интегралов в (11), где также применяется формула трапеций:
Jу кор(х(у = JJ ЛБОВ ЛБОВ
. дк01 дк02 . . 10, 6 ( + + 12коо
ох оу
р(х(у
6ро \ J ко1(1у — J ко1(1у + J ко2(1х — J ко2(1х I + 12ро JJ kооdxdy ; ЗО ЛВ ВО ЛБ / ЛБОВ
„ Зко\(Б) + ког(С)и МА) + ко1р) , ЫС) + ,
6Ро ( -2-—2--2-—2 + -2--1 —
ЫА) + ко2(В)+ Зро коо(А) + коо(В) + коо(С) + коо(^)
-1 -2, (19)
Я(х(у = [ /1(у — [ /1(у + [ /2(х — [ /2(х — 6 [( ~я(х(у+
Ь(х(у « ^ (/1(В) + /1(С)) — ^ (/1(А) + /1р)) +
ЛБОВ
+ у (/2Р) + /2(С)) — -1 (/2(А) + /2(В)) —
—3№ (НГг (А) + НГг (В) + НГ (С) + НГ (°)) +2 Ц ь(Х(У. (20)
ЛБОВ
Последний интеграл вычисляется по-разному в зависимости от того, перед каким этапом решения гиперболической части — предиктором или корректором — определяется дисперсионная составляющая давления. Перед предиктором для вычисления разностных производных от компонент скорости, заданных в целых узлах, используются формулы вида (18), с помощью которых выражение для Ь определяется в центрах ячеек сетки и полагается
цЬ(х(у«^ (Ь(А>+ь(в>+ь(с )+шу (21)
ЛБОВ
При повторном вычислении функции р, т. е. перед выполнением корректора, используются компоненты скорости и полная глубина, определённые в центрах ячеек сетки
Коэффициенты разностного уравнения (23)
Тип узла а1 а2 аз а4 а5 ав а7 а8
0 в Л + в В - вл - вБ в б + во -во - во 7л 12б 1Ь 12п
1 0 -вБ в Б + во -во 0 12б "(о 0
2 во 0 во -во - во 0 0 11
3 в Л + во - вл 0 -во 7л 0 0 То
4 вл - вл - вБ вв 0 ^л ^Б 0 0
5 0 0 во -во 0 0 0
6 во 0 0 -во 0 0 0 11
7 вл -вл 0 0 7л 0 0 0
8 0 -вв вв 0 0 ^Б 0 0
9 во -вБ в Б + во -во - во 0 ^о 11
10 в л + во -вл во -во - во 7л 0 Тс То
11 вл + во -вл - вБ вв -во ^л ^Б 0 То
12 вл -вл - вБ в Б + во -во ^л 12Б то 0
после исполнения шага предиктор для гиперболической части расщеплённой системы. В этом случае применяется следующая аппроксимация:
// ЬЫу и (и(Б)+ "(С) - и(А) - "(Д) + V(D> + "(С> - - ^" 2
ЛБОВ
,'и(Б) + и(С) - и(А) - и(Б) ю(Б) + у(С) - и(А) - и(Б)
—Д1Д2
2/ 2^2 и(Р) + и(С) - и(А) - и(Б) у(Б)+ у(С) - у(А) - у(Р) N (22)
2К2 2/ )' ()
Таким образом, в каждом внутреннем узле х, Е П^ для сеточной функции р получается 9-точечное разностное уравнение вида
= Р, (23)
' г=о 3
с правой частью, определяемой аппроксимацией (20), при этом
8
ао = - ^ аг,
г=1
а для вычисления коэффициентов аг (г = 1,. ■■ , 8) используются формулы из первой строки таблицы, записанные с использованием обозначений
а к2 1 1 к2 1 , , Ь 2 /2 1 . Ы 1
в = 4/21к11 - 4/2к22, 7 = 421к11 + 4/2к22 + 2к12, 7 = 4/21к11 + 4/2к22 - 2
2.2. Граничные условия на непроницаемых границах
Покажем особенности вычисления коэффициентов и правой части уравнения (23) для граничных узлов, при этом будем предполагать, что дно вблизи суши является неподвижным. Для примера выберем узел х, Е Г^ типа 2, т. е. рассмотрим случай, когда два
звена ломаной, представляющей границу Г, стыкуются в рассматриваемом узле, имеющем локальный номер 0, и оба лежат на одной и той же горизонтальной координатной линии сетки с номером j2, причём область течения располагается выше этого участка границы (см. рис. 2, б). Выберем контур интегрирования — прямоугольник ABCD — так, чтобы две его вершины A и B находились в серединах смежных звеньев границы, а вершины C и D — в центрах приграничных ячеек с общей вершиной Xj. Интегрируя уравнение (11) по данному контуру, получим то же соотношение (16), однако аппрокси-мационные формулы для некоторых интегралов будут отличаться от тех, что использовались при вычислении интегралов для внутренних узлов. Так, для интегралов из левой части (16) применим следующую аппроксимацию:
/ ФЧу - / ФЫу + / ф2<Ъ - / ф2Л &
BC AD DC AB
2ф'(С) - hf+ -2
уФ'(С) - f <í'(D) + f (ф2(С) + <I>2(D)) - <иФ2(0), (24)
где величины Ф1, Ф2 в центрах ячеек вычисляются аналогично (17). Что касается величины Ф2(0), которую надо определить в граничном узле х, то, как будет показано ниже, никаких аппроксимаций для неё строить не нужно, поскольку она в силу заданного краевого условия (10) взаимно сокращается с другими членами уравнения.
Приведём теперь аналоги формул (19), (20) в том же граничном узле х:
JJ kopdxdy & 6po I J koidy - J koidy + J ko2dx - J ko2dx | +
ABCD \BC AD DC AB
±10 /7, л л a U ính kof(C) + kof(D)
+12po / koodxdy & 6po ( koi(C)y - koi(D)— +-^-hl-
ABCD
-ko2(0)hi) + 12pkoo(C} + koo(D) hhf = 3vo([koi(C) - koi(D)]hf+
+ [kof(C) + kof(D)]hi + [koo(C) + koo(D)]hihf) - 6yokof(0)hi, (25)
I[ Fdxdy - /1(С)^ - МП)^ + /2(С} + /2(Д) /1 - /1/2(0)-
ЛБОВ
- 3№ (ТТг (С) + !ТГ Р)) + 2 Ц Ldxdy. (26)
ЛБОВ
Краевое условие (10) для узлов типа 2 принимает следующий вид:
, + hypy 6 \ / R .
Py- ---+ —И = H [9Vy + ,
поэтому с учётом обозначений (12)—(15) его можно переписать как
Ф2 = 6ko2P + ¡2. (27)
Подставив формулы (24)-(26) в уравнение (16) и приняв во внимание краевое условие (27), убеждаемся, что подчеркнутые слагаемые в формулах (24)-(26) взаимно сокращаются. Таким образом, при использовании интегро-интерполяционного метода не требуется аппроксимировать краевое условие (27) конечными разностями, поскольку оно учитывается в граничных узлах в точной форме при вычислении интегралов по сторонам контура интегрирования, лежащим на границе расчётной области.
Применение интегро-интерполяционного метода для граничных узлов типа 2 приводит к разностному уравнению на 6-точечном шаблоне, изображённом на рис. 2, б. Для рассматриваемого граничного узла шаблон не содержит узлов, имеющих номера 5, 2 и 6 в шаблоне внутреннего узла (см. рис. 2, а). Тем не менее формально разностное уравнение для граничного узла типа 2 можно записать в виде 9-точечного уравнения (23), положив коэффициенты а5, а2, а6 равными нулю. Формулы для вычисления коэффициентов уравнения (23) в граничных узлах типа 2 приведены в третьей строке вышеприведённой таблицы.
Аналогично получаются разностные уравнения в граничных узлах других типов, при этом должны выбираться контуры интегрирования, соответствующие типам узлов. Пример контура интегрирования для узлов типа 11 показан штриховой линией на рис. 2, в. Для этих узлов разностное уравнение является 8-точечным и не содержит в своем шаблоне узла с номером 7 из 9-точечного шаблона внутреннего узла. В таблицу сведены формулы для вычисления коэффициентов уравнения (23) для узлов сетки всех возможных типов.
Отметим, что система сеточных уравнений вида (23) содержит уравнения для всех узлов х Е П^. Матрица этой системы является симметричной. Выбирая аналогично [26] гильбертово пространство сеточных функций, можно показать, что разностный оператор, соответствующий системе уравнений (23), является самосопряжённым, причём и в случае неровного подвижного дна. В [26] аналогичный результат получен только для ровного горизонтального дна. Кроме того, в случае прямоугольной области П можно показать, что оператор разностной задачи для дисперсионной составляющей давления является положительно определённым, если всюду в области выполняется неравенство ко > 0. Аналогичный результат для ровного горизонтального дна был установлен ранее в работе [26].
Полученная система уравнений (23) решается итерационным методом. В общем случае область П имеет сложную форму границы (см. рис. 2, б). Кроме того, она может быть многосвязной (см. ниже рис. 9, а). Поэтому для итерационного решения системы уравнений (23) был использован метод последовательной верхней релаксации, легко реализуемый для областей сложной конфигурации.
2.3. Неотражающие граничные условия
Рассмотрим постановку неотражающих условий на примере узла х на прямолинейной левой границе (^ = (0,]2), ]2 = 1, • • • ,Ы2 — 1). В расчётах по классической модели мелкой воды зачастую используются условия типа условий Зоммерфельда, которые в данном случае можно записать как
П — ^/-Пх = 0, щ — л/-их = 0, ух = 0.
При этом считается, что в окрестности границы глубина постоянна. Можно показать, что при явной аппроксимации указанных условий с первым порядком:
Нп+1(0-2) - Нп(0-2) - /щ-уНп(1-2) - Нп(0-2) = 0, т ' Н1 '
ип+1(0,-2) - иП(0,-2) - /Щ-)иП(1-2) - иП(0,-2) = 0 т ' Н1 '
^п+1(0,-2) = ^+1(1-2), (29)
где т — временной шаг, абсолютная амплитуда отражённых волн в расчётах по модели мелкой воды на рассмотренных задачах не превышает 0.1 % от амплитуды набегающих.
Уравнения вида (28) можно использовать и для гиперболической части НЛД-модели. Дополнительно необходима постановка краевого условия для эллиптического уравнения (3). Одним из возможных вариантов такого условия для левой границы является аналогичное (28) уравнение
щ - = 0. (30)
Аппроксимировать его будем также при помощи явной схемы первого порядка. На (п + 1)-м временном шаге перед предиктором, когда с предыдущего шага известны функции Нп,ип и р*, значение рп на левой границе находится из уравнения
¥>"(0,-2) - Р*(0-2) - УЩ-) Р*(1-2) - Р*(0-2) = 0. Перед корректором известны Н*,и* и рп, поэтому р* вычисляется из соотношения
¥*(0-2) - ¥п(0,-2) - ущ-)¥"(1,-2) - ¥"(0,-2) = 0.
т/2 Ь,1
При такой реализации краевых условий в рассмотренных расчётах по НЛД-модели абсолютная амплитуда отражённых волн составляет менее 0.1% от амплитуды набегающих, что в большом классе решаемых задач незначительно влияет на волновую картину.
3. Тестовые задачи
В настоящем разделе описаны четыре модельные задачи, на которых исследуется влияние учёта дисперсии волн и пространственной неоднородности оползня, а также соответствие полученных результатов некоторым экспериментальным данным.
3.1. Распространение волн над ровным дном
Первая модельная задача — распространение волн над ровным дном, вызванных начальным возмущением свободной поверхности вида
П(х, у, 0) = аое—[(х-х°)2+(у-уо)2], (31)
где а0 — высота возмущения, х0, у0 — координаты центра возмущения, т — параметр, влияющий на эффективную протяжённость. В качестве акватории выбран квадратный
бассейн со сторонами Ь1 = Ь2 = 10 км и глубиной = 100 м. На всех границах ставились неотражающие краевые условия. Начальное возмущение располагалось по центру бассейна, х0 = у0 = 5 км, и имело высоту а0 = 10 м. Параметр и принимал значения и = 2.5 • 10-6, 10-5 и 5 • 10-5 м-2, что соответствует эффективной протяжённости, приблизительно равной 3000, 1500 и 750 м (рис. 3, а).
На рис. 3, б представлен профиль свободной поверхности, вычисленной по НЛД-модели при и = 10-5 м-2, в момент времени £ = 140 с. Заметим, что картина волн симметрична относительно центра области (более точно это было проверено при сравнении соответствующих мареограмм), следовательно, построенный численный алгоритм сохранил изначальную симметрию задачи. Также видно, что за основной волной образуется цуг второстепенных волн меньшей амплитуды. Такая картина характерна только для моделей, учитывающих дисперсию волн, и не раз наблюдалась в экспериментах. При этом, как отмечалось, например, в работе [27], влияние дисперсии зависит от эффективной протяжённости начального возмущения. Для проверки этого свойства на рис. 4 приведены мареограммы в точке (0 км, 5 км), полученные в расчётах по НЛД-
а б
ю
Г], М
/1
//:' \ 2 .....3
/ : \
' 1 • \
/ /; \
У ' 1
2000
4000
6000
Я, м
8000
Рис. 3. а — сечение у = уо начальной свободной поверхности при т = 2.5
10-6м-2
-5 м-2
(1);
10-5 м-2 (2); 5 • 10 в расчёте по НЛД-модели при т = 10
(3); б — профиль свободной поверхности в момент времени £ = 140 с
-1П-5 м-2
300 с
300 и С
Рис. 4. Мареограммы в точке (0 км, 5 км), полученные в расчётах по НЛД- и НЛ-моделям для
разной эффективной ширины начального возмущения: т = 2.5 • 10
-6 м-2
(а), т = 10
-5 м-2
т = 5 • 10 5 м 2 (в). Сплошные линии — результаты расчётов по НЛД-модели, пунктирные по модели мелкой воды
(б),
а
и НЛ-моделям для т = 2.5 • 10-6, 10-5 и 5 • 10-5 м-2. Как и предполагалось, при очень протяжённых начальных возмущениях дисперсия волн не оказывает значительного влияния за рассмотренное время распространения, но при более коротких волнах для адекватного воспроизведения картины течения учёт нелинейно-дисперсионных свойств может быть решающим.
3.2. Движение твёрдого оползня по прямолинейному склону
Проведём сравнение полученных результатов с экспериментальными данными из работы [7]. Серия экспериментов была выполнена в бассейне, наполненном до глубины к0 = 1.5 м, длиной Ь1 = 30 ми шириной Ь2 = 3.7 м. В принятых в настоящей статье обозначениях на левой стенке при х = 0 м невозмущённая свободная поверхность сопрягалась с плоским склоном, образующим с ней угол 15°. В качестве модели оползня использовалось твёрдое тело из алюминия гладкой формы массой 16 кг, толщиной Т = 0.082 м, длиной Ь = 0.395 м и шириной т = 0.68 м. При таких параметрах объём оползня V = 6.57• 10-3 м3, а его плотность р = 2.45• 103 кг/м3. В начальный момент времени оползень находился на заглублении d — расстояние от невозмущённой свободной поверхности до вершины профиля оползня, которое определяло абсциссу его начального положения х°. Ордината же всегда полагалась равной у° = 1.85 м. В экспериментах менялся параметр d, здесь будет рассмотрено только d = 0.061 и 0.189 м — два предельных случая, для которых приведены графики мареограмм. Для таких заглублений абсцисса начального положения составляла соответственно х0 ~ 0.522 и ^ 1 м. Четыре мареографа располагались следующим образом: первый находился над вершиной оползня (х0, у0), остальные устанавливались одинаково для всех экспериментов и имели координаты (1.469 м, 2.2 м), (1.929 м, 1.85 м) и (1.929 м, 2.385 м).
Для описания движения оползня по криволинейному склону выделим в функции к(х,у,1) подвижную часть:
где кы — функция стационарного основания дна, по которому движется оползень с профилем hsi(x,y,t). Оползень ассоциируется с материальной точкой — центром его носителя (xc, yc), скользящей по дну согласно закону несвободного движения по криволинейной поверхности. На оползень действовали силы тяжести, плавучести, трения о дно и сопротивления воды. Подробный вывод и описание используемого закона можно найти в [28], здесь лишь отметим, что определяющими параметрами в нём, помимо размеров оползня, являются отношение плотности оползня к плотности воды y, коэффициенты гидродинамического сопротивления Cd, присоединённой массы воды Cw и трения о дно Cf = tan в*, где в* — угол трения. Для численного моделирования использовалась следующая форма оползня:
h(x,y,t) = -hbt(x,y) - hsi(x, y,t),
(32)
иначе,
где толщина Т была той же, что в эксперименте, а длины Ьх и Ьу вдоль осей Ох и Оу выбирались так, чтобы сохранялись объём оползня и соотношение Ьх/Ьу = Ь/ш, Ьх = 0.431 м, Ьу = 0.743 м. В [7] приводятся следующие вычисленные значения параметров закона движения оползня по плоскому склону:
При этом трением оползня о дно склона пренебрегалось. Для визуального сходства графиков движения оползня вдоль склона в настоящей работе вводилось "искусственное" трение: 9* = 1.7° при d = 0.061 ми 9* = 1.2° при d = 0.189 м.
На рис. 5 представлены записи двух мареографов для d = 0.061 м, полученные в эксперименте и расчётах по НЛД- и НЛ-моделям. На первой мареограмме (рис. 5, а) видно, что обе модели качественно хорошо согласуются с экспериментом, но НЛД-модель значительно лучше описывает амплитуду основной впадины. Из записей второго мареографа (рис. 5, б) следует, что НЛ-модель удовлетворительно описывает только первое возвышение и впадину, идущую за ним, а далее сильно упрощает картину течения. Напротив, НЛД-модель неплохо согласуется с экспериментом и после первой впадины, особенно при описании волны самой большой амплитуды, следующей сразу за ней. Однако стоит отметить некоторое преувеличение амплитуды второстепенных волн и их длины с течением времени. Последнее можно объяснить тем, что оползень в это время уже двигался по более глубокой части склона, где модели мелкой воды могут терять адекватность.
Для проверки данного предположения рассмотрим записи третьего мареографа при d = 0.061 м, располагающегося в более глубокой части склона (рис. 6, а). Видно, что несмотря на качественную схожесть волновых картин амплитуды волн в расчёте по НЛД-модели значительно отличаются от амплитуд в эксперименте. В случае d = 0.189 м, когда оползень стартует из более глубокого места, даже согласно записям первого мареографа (рис. 6, б), располагающегося над начальным положением, отличия расчётов НЛД-модели от эксперимента достаточно велики. Было проверено, что для остальных мареографов наблюдалась такая же картина. Из этих сравнений можно сделать
Cw = 0.607, Cd = 0.473 при d = 0.61, Cw = 0.582, Cd = 0.353 при d = 0.189.
'w
a
б
7], CM
-2
0
-2
0
0.4
0.8
1.2
0
2
3 t, с
Рис. 5. Записи первого (а) и второго (б) мареографов при d = 0.061 м: 1 — расчёт по НЛД-модели, 2 — эксперимент, 3 — расчёт по НЛ-модели
а
б
7], СМ
Г}, СМ
0.4
0
О
-2 -
1 2
-0.4
-4 Н—"—"—"—|—г
О 1
2
3 Ас 4
-0.8
О
0,4
0,8
1,2 1,6
Рис. 6. Записи первого мареографа при d = 0.061 м (а) и третьего мареографа при d = 0.189 м (б): 1 — расчёт по НЛД-модели; 2 — эксперимент
вывод, что в случае, когда оползень стартовал и двигался по верхней части склона, НЛД-модель хорошо описывала рассмотренный эксперимент, в отличие от НЛ-модели воспроизводя не только впадину над оползнем, но и цуг волн после неё. Однако при движении по глубокой части НЛД-модель может значительно завышать амплитуду. В одномерных задачах такой эффект был описан, например, в работе [14], где использовалась НЛД-модель Лью — Лайнетта, и в статьях [20-22], где была применена настоящая модель. При этом отмечается, что определяющую роль играет не заглубление оползня, а его отношение к собственной длине.
3.3. Сход подводного оползня в модельном участке реки
Сход подводных оползней может быть опасен не только в океанах, но и в крупных внутренних водоёмах [3]. В настоящем разделе рассмотрим задачу о сходе подводного оползня в модельном участке реки с прямолинейными границами. Форма основания дна, по которому двигался оползень, была однородна в направлении О у и имела вид дуги параболы:
где Н^^ = 10 м — заглубление на левой и правой стенках, £ = 500 м — точка вершины параболы, глубина в которой Н^ = 100 м. Таким образом, длина области в направлении Ох получилась = 1000 м, такое же значение было выбрано для направления Оу: Ь2 = 1000 м. Слева и справа модельный водоём ограничен вертикальными стенками, а на нижней и верхней границах ставились условия свободного прохода волн.
Форма оползня Н^ задавалась по формуле (33), в которой принимались следующие значения параметров: (х0, у0) = (109.1 м, 500 м), Ьх = 200 м, Ьу = 400 или 800 м. Абсцисса начального положения х0 была выбрана так, что начальное заглубление центра носителя оползня ¿0 = 45 м. В предыдущих исследованиях на одномерных задачах [21] показано, что при таком соотношении заглубления к длине оползня Ьх представленная НЛД-модель даёт очень близкие к модели потенциальных течений результаты, что может указывать на её адекватность в данном классе задач.
(34)
В законе движения оползня выбирались значения параметров, принятые в [21], за основу:
7 = 2, С = 1, Ст = 1, в* = 5°.
Их влияние на образование волн при исследовании в рамках НЛ-модели можно найти в работе [29]. Из-за ограниченности объёма настоящей статьи эти параметры будут фиксированы.
При сравнении с расчётами по одномерной модели покажем влияние учёта неоднородности формы оползня в двух направлениях. Одномерный расчёт в этом случае соответствует бесконечно широкому оползню. Рассмотрим мареограммы на стенках слева (0 м, 500 м) и справа (1000 м, 500 м) от центра оползня (рис. 7). Отметим, что законы движения оползня для представленных расчётов совпадают.
Сравнивая двумерные расчёты при разной ширине оползня, можно заметить, что с увеличением ширины в два раза амплитуда большинства волн также возрастает приблизительно в два раза. Этот эффект можно было ожидать в связи с увеличением в два раза объёма и, следовательно, массы оползня. Исключением стали только первые волны на левой стенке, которые по амплитуде среди всех трёх расчётов различаются не столь сильно. Причиной этому является начальное расположение оползня, близкое к левой стенке, так как волны накатились на неё за короткое время и не успели значительно распространиться в направлении Оу. С учётом последнего эффект "двумерности" в данном случае был меньшим, чем при описании остальных волн, где можно заметить большое различие в амплитудах. При этом качественно волновые картины достаточно похожи.
Для объектов, расположенных на левом и правом берегах, важной может быть информация о распределении максимальных за всё время расчёта заплесков. В одномерном расчёте, очевидно, максимальный заплеск будет одинаковым вдоль всей стенки, но в двумерных случаях картина оказывается нетривиальной. На рис. 8 представлены максимальные заплески вдоль левой и правой стенок за время £ = 200 с при ширине
Рис. 7. Рассчитанные мареограммы на стенках слева (а) и справа (б) от центра оползня: 1 — by = 400 м, 2 — by = 800 м, 3 — одномерный расчёт
0.5
^шах'М
0.4
0.3
0.2
1 -| 0.8 -
0.6-
0.4
250
500
750 v 1000 Л, м
250
500
750 v лд 1000 х, м
Рис. 8. Максимальная свободная поверхность за всё время расчёта при by = 400 м (а) и by = 800 м (б): 1 — на левой, 2 — на правой стенке
оползня by = 400 и 800 м. Из графиков видно, что максимальные заплески на правой, дальней от начального положения оползня, стенке значительно больше, чем на левой. При этом различаются и формы этих заплесков: справа максимум находится по центру стенки, напротив центра оползня, и при удалении от центра заплеск почти монотонно убывает, слева же максимумы расположены симметрично на небольшом отдалении относительно центра, при этом имеются ещё два симметричных локальных максимумов, близких по амплитуде к первым.
3.4. Моделирование взаимодействия уединённой волны с островом цилиндрической формы
Приведём результаты численного моделирования процесса взаимодействия уединённой волны с островом цилиндрической формы. В начальный момент времени волна располагалась так, что её гребень был параллелен оси ординат (рис. 9), т. е. все точки гребня имели одинаковую абсциссу х0. Начальное возвышение и скорость течения при £ = 0 с задавались по формулам
П0) = ao ch
V 2hoUo
(x - Xo)
u(x,y, 0) = Uo:
П(x,y, 0)
v(x, y, 0) = 0,
(35)
(36)
Н0 + п(х, у, 0)'
где и0 = у/д(Н0 + а0), Н0 — толщина слоя воды над горизонтальным дном, а0 — амплитуда начальной волны. На рис. 9, а штриховой линией показано положение х = х0 гребня уединённой волны в начальный момент времени.
Расчёты проводились для бассейна длиной ^ = 35 м, шириной = 30 ми глубиной Н0 = 32 см. Остров радиуса Яы = 2.32 м располагался в точке х* = 22.96 м, у* = 13.8 м в центре усечённого прямого кругового конуса с параметрами Я0 = 3.6 ми Я1 = 2.576 м — радиусы основания и верхней части соответственно, х = аге1ап(1/4) — угол между образующей конуса и плоскостью я = — Н0, Ны = 6.4 см — заглубление верхней части цилиндра (рис. 9, б). Дно в этом случае описывалось уравнением
—Ны при Яы < р*(х,у) < Я1,
я = —Н(х,у )=^ —Н0 + 1ап х(Я0 — р*) при Я1 < р* (х,у) < Я0, (37)
-Н0 при р*(х,у) > Я0,
а
Рис. 9. Схема области решения (а) с расположенным в ней препятствием (б) в виде кругового цилиндра, установленного на усечённый конус; в — аппроксимация границы сечения цилиндрического препятствия с помощью ломаной
I
41 Г|, см
2 -
О -2
28
с 40
32
36
Г, с
40
II
4 п Л > см
2-
28
32
г
32
*,с 40
36
Г, с 40
6
3
о
-3
6
3
о
-3
- Л, см 6 3 Л, см А
ЧгТу. 0 -3: ' \ • • • «
..... V .....
28
Л, см
28
32
32
36
40
к --ГГГТ^^ у
•ч»—Г» • 4
• » » •
36
£ с 40
28
Л, см
28
32
г
32
36
40
36
г, с
40
Рис. 10. Мареограммы, измеренные волномерами Вб (а), Вд (б), (в), В22 (г) в эксперименте [24] (маркеры) и полученные в расчётах по НЛД-модели (линии) для амплитуды набегающей волны ад = 1.44 см (I) и ад = 2.9 см (II)
а
в
а
в
Рис. 11. Распределение максимальных высот заплесков на круговой конус в эксперименте [24] (маркеры) и на круговой цилиндр в расчётах по НЛД-модели (линии) для амплитуды набегающей волны ад = 1.44 см (а), ад = 2.9 см (б)
где р*(х,у) = д/(х — х*)2 + (у — у*)2. Граница цилиндра аппроксимировалась ломаной, как показано на рис. 9, в. В начальный момент времени гребень волны располагался при х0 = 10 м, а в качестве её амплитуды выбирались два значения: а0 = 1.44 см и а0 = 2.9 см. Такие параметры задачи были выбраны для сравнения с экспериментальными данными, полученными в [24]. Следует отметить, что в эксперименте весь остров представлялся усечённым конусом, возвышающимся над водой.
В расчётах и экспериментах картина взаимодействия волны с островом фиксировалась 27-ю волномерами. По аналогии с [30,31] в настоящей работе рассматривались четыре из них (см. рис. 9, б) с координатами В6 = (19.36 м, 13.8 м), В9 = (20.36 м, 13.8 м), В16 = (12.96 м, 11.22 м) и В22 = (15.56 м, 13.8 м). Сравнение расчётов по НЛД-модели с экспериментальными данными в указанных волномерах для амплитуд а0 = 1.44 см иа0 = 2.9 см представлено на рис. 10, I, II соответственно. При сравнении представленных данных видно их неплохое согласование, но вместе с тем в эксперименте волномеры записывали более сложную картину при откате волны. Последнее связано с тем, что в отличие от экспериментов в расчётах накат осуществлялся на вертикальную преграду. При этом наибольшие отличия фиксировались волномером В22, который находился около дальней от начального положения волны точки острова.
Картину максимальных заплесков на остров можно наблюдать из рис. 11. Здесь угол в по оси абсцисс отсчитывался от точки острова с наименьшей ординатой (см. рис. 9, б). Результаты сравнения заплесков показали качественное совпадение представленных данных, при этом и в экспериментах, и в расчётах для разных начальных амплитуд волн различаются положения максимального заплеска на остров: при а0 = 1.44 см он находится около в = 270°, т.е. в ближней к начальному положению волны точке, а при а0 = 2.9 см — около в = 90°. В некоторых проведённых ранее исследованиях в рамках модели мелкой воды [31,32] такого эффекта не наблюдалось.
Заключение
В работе представлен численный алгоритм для решения полных плановых НЛД-урав-нений Железняка — Пелиновского на подвижном дне, основанный на расщеплении исходной системы на гиперболическую и эллиптическую части. Подробно описан алгоритм для эллиптической части, построенный при помощи интегро-интерполяционно-го метода. Для гиперболической части использовалась схема типа предиктор-корректор [23], хорошо зарекомендовавшая себя при решении классических уравнений мелкой воды.
На задаче о распаде начального возмущения с различной эффективной протяжённостью над ровным дном сравнивались расчёты по НЛД- и НЛ-моделям. Показано, что важность учёта дисперсии значительно зависит от эффективной протяжённости начального возмущения и может иметь решающую роль для адекватного описания волновой картины. При сравнении расчётов по указанным моделям с экспериментальными данными [7] о сходе подводного оползня продемонстрировано, что НЛД-модель точнее описывает амплитуды и качественную картину волн, полученные в эксперименте. Однако даже эта модель теряет адекватность, когда оползень движется по глубокой части исследуемого водоёма. Проведены расчёты в модельном участке реки с параболической формой дна, ограниченном около берегов вертикальными стенками. При сравнении ма-реограмм для берегов, полученных по одномерной и двумерной моделям при разной ширине оползня, сделаны выводы, что амплитуды большей части образованных волн возрастают пропорционально ширине оползня, а одномерная модель эти амплитуды сильно завышает. Качественно же картины волн во всех рассмотренных случаях получились похожими. Представлены графики распределения максимальных за всё время расчёта заплесков на берега. Показано, что на противоположном от начального положения оползня берегу заплески были значительно больше, при этом на разных берегах различаются расположения глобальных максимумов. В задаче о взаимодействии волны с коническим островом [24] границы острова аппроксимировались ломаной кривой, а верхняя часть конуса заменялась затопленным цилиндром. При таких условиях наблюдалось хорошее качественное соответствие полученных волновых картин с экспериментом, в том числе — совпадение расположения максимальных заплесков на остров. Можно ожидать, что применение алгоритмов наката волн с подвижной линией уреза улучшит амплитудное соответствие, особенно в волнах понижения при откате.
Список литературы
[1] Tappin D.R., Watts P., Grilli S.T. The Papua New Guinea tsunami of 17 July 1998: Anatomy of a catastrophic event // Natural Hazards Earth Syst. Sci. 2008. Vol. 8. P. 243-266.
[2] Ward S.N., Day S. The 1963 landslide and flood at vaiont reservoir Italy. A tsunami ball simulation // Ital. J. Geosci. 2011. Vol. 130, No. 1. P. 16-26.
[3] Диденкулова И.И., Пелиновский Е.Н. Цунамиподобные явления в российских внутренних водоёмах // Фундамент. и прикл. гидрофизика. 2009. № 3(5). С. 52-96.
[4] Lynett P.J., Liu P.L.-F. A numerical study of the run-up generated by three-dimensional landslides // J. Geophys. Res. 2005. Vol. 110. C03006, doi:10.1029/2004JC002443.
[5] Елецкий С.В., Майоров Ю.Б., Максимов В.В. и др. Моделирование генерации поверхностных волн перемещением фрагмента дна по береговому склону // Совместный выпуск журн. "Вычисл. технологии" и "Вестник КазНУ им. аль-Фараби". 2004. Т. 9, ч. II. С. 194-206.
[6] Grilli S.T., Watts P. Tsunami generation by submarine mass failure. I: Modeling, experimental validation, and sensitivity analyses // J. of Waterway Port Coastal and Ocean Eng. 2005. Vol. 131, No. 6. P. 283-297.
[7] Enet F., Grilli S.T. Experimental study of tsunami generation by three-dimensional rigid underwater landslides // Ibid. 2007. Vol. 133, No. 6. P. 442-454.
[8] Watts P., Grilli S.T. Modeling of waves generated by a moving submerged body. Applications to underwater landslides // Eng. Anal. with Boundary Elements. 1999. Vol. 23. P. 645-656.
[9] Хлкимзянов Г.С., Шокин Ю.И., Барахнин В.Б., Шокинл Н.Ю. Численное моделирование течений жидкости с поверхностными волнами. Новосибирск: Изд-во СО РАН, 2001. 394 с.
[10] Green A.E., Naghdi P.M. A derivation of equations for wave propagation in water of variable depth //J. Fluid Mech. 1976. Vol. 78, pt 2. P. 237-246.
[11] Железняк М.И., ПЕЛМНОвский Е.Н. Физико-математические модели наката цунами на берег // Накат цунами на берег. Горький: ИПФ АН СССР, 1985. С. 8-33.
[12] Базденков С.В., Морозов Н.Н., ПогуцЕ О.П. Дисперсионные эффекты в двумерной гидродинамике // Докл. АН СССР. 1987. Т. 293, № 4. С. 818-822.
[13] Алешков Ю.З. Течения и волны в океане. СПб.: Изд-во С.-Петербургского ун-та, 1996. 226 с.
[14] Lynett P.J., Liu P.L.-F. A numerical study of submarine-landslide-generated waves and run-up // Proc. Royal Soc. of London. A. 2002. Vol. 458. P. 2885-2910.
[15] Peregrine D.H. Long waves on a beach // J. Fluid Mech. 1967. Vol. 27, pt 4. P. 815-827.
[16] Дорфман А.А., Яговдик Г.И. Уравнения приближённой нелинейно-дисперсионной теории длинных гравитационных волн, возбуждаемых перемещениями дна и распространяющихся в бассейне переменной глубины // Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ, ИТПМ. 1977. Т. 8, № 1. С. 36-48.
[17] Шокин Ю.И., Чубаров Л.Б. О подходах к численному моделированию оползневого механизма генерации волн цунами // Вычисл. технологии. 2006. Т. 11. Спец. выпуск. Ч. 2. С. 100-111.
[18] Fedotova Z.I., Pashkova V.Yu. Methods of construction and the analysis of difference schemes for nonlinear dispersive models of wave hydrodynamics // Rus. J. Numer. Anal. Math. Modelling. 1997. Vol. 12, No. 2. P. 127-149.
[19] Wei G., Kirby J.T., Grilli S.T., Subramanya R. A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves //J. Fluid Mech. 1995. Vol. 294. P. 71-92.
[20] Гусев О.И. Об алгоритме расчёта поверхностных волн в рамках нелинейно-дисперсионной модели на подвижном дне // Вычисл. технологии. 2012. Т. 17, № 5. С. 46-64.
[21] Гусев О.И., ШокинА Н.Ю., Кутергин В.А., ХАкимзянов Г.С. Моделирование поверхностных волн, генерируемых подводным оползнем в водохранилище // Там же. 2013. Т. 18, № 5. С. 74-90.
[22] Шокин Ю.И., Бейзель С.А., Гусев О.И. и др. Численное исследование дисперсионных волн, возникающих при движении подводного оползня // Вестник ЮУрГУ. 2014. Т. 7, № 1. С. 121-133.
[23] Шокин Ю.И., ХАкимзянов Г.С. Схема предиктор-корректор, сохраняющая гидравлический скачок // Вычисл. технологии. 2006. Т. 11. Спец. выпуск. Ч. 2. С. 92-99.
[24] Briggs M.J., Synolakis C.E., Harkins G.S., Green D.R. Laboratory experiments of tsunami runup on circular island // Pure and Appl. Geoph. 1995. Vol. 144, No. 3/4. P. 569-593.
[25] Федотова З.И., ХАкимзянов Г.С. Анализ условий вывода НЛД-уравнений // Вычисл. технологии. 2012. Т. 17, № 5. С. 94-108.
[26] Barakhnin V.B., Khakimzyanov G.S. On the algorithm for one nonlinear dispersive shallow-water model // Rus. J. Numer. Anal. Math. Modelling. 1997. Vol. 12, No. 4. P. 293-317.
[27] Kirby J.T., Shi F., Tehranirad B. et al. Dispersive tsunami waves in the ocean: Model equations and sensitivity to dispersion and Coriolis effects // Ocean Modelling. 2013. Vol. 62. P. 39-55.
[28] Бейзель С.А., ХАкимзянов Г.С. Численное моделирование поверхностных волн, возникающих при движении подводного оползня по неровному дну // Вычисл. технологии. 2010. Т. 15, № 1. С. 105-119.
[29] Beisel S.A., Chubarov L.B., Khakimzyanov G.S. Simulation of surface waves generated by an underwater landslide moving over an uneven slope // Rus. J. of Numer. Anal. Math. Modelling. 2011. Vol. 26, No. 1. P. 17-38.
[30] Chubarov L.B., Fedotova Z.I., Shkuropatskii D.A. Investigation of computational models of long surface waves in the problem of interaction of a solitary wave with a conic island // Ibid. 1998. Vol. 13, No. 4. P. 289-306.
[31] Chubarov L.B., Fedotova Z.I. An Effective high accuracy method for tsunami runup numerical modeling // Submarine Landslides and Tsunamis: Proc. of the NATO Advanced Research Workshop on Underwater Ground Failures on Tsunami Generation. Modeling, Risk and Mitigation. Istanbul, Turkey. Dordrecht: Kluwer, 2003. P. 203-216.
[32] Liu P.L.-F., Cho Y.-S., Briggs M.J. et al. Runup of solitary waves on a circular island // J. Fluid Mech. 1995. Vol. 302. P. 259-285.
Поступила в редакцию 29 мая 2014 г.
Algorithm for surface waves calculation above a movable bottom within the frame of plane nonlinear dispersive model
Gusev Oleg I.1
The influence of the effects due to frequency dispersion on tsunami wave patterns is investigated. The numerical algorithm developed is based on the partitioning of the fully nonlinear dispersive shallow water equations with a movable bottom into elliptic and hyperbolic subproblems, which are solved alternately at each time step. At that, equations of the hyperbolic subproblem differ from the classic shallow water system in the right side only. Therefore, well-examined finite-difference schemas are implemented for both subproblems. In this paper, we describe in detail the integro-interpolation method for the elliptic subproblem. The hyperbolic subproblem is solved by the explicit schema of predictor-corrector type. For open bounds of the computational region Zommerfeld type conditions are proposed for both subproblems.
Obtained numerical solutions are compared with computations based on shallow water model and experimental data. Regarding the problem of disintegration of initial disturbance above the flat bottom, it is checked out that as the source width is decreased so dispersion is significantly increased. During the simulations of the landslide-generated tsunamis, we detect the loss of adequacy of the model when the ratio of the landslide length to its depth is small. In other cases the model reproduce wave pattern better than the classic
1 Institute of Computational Technologies SB RAS, Novosibirsk, 630090, Russia
Corresponding author: Gusev Oleg I., e-mail: [email protected]
Aknowlegements: Research was supported by RFBR (Russian Foundation for Basic Research) grant No. 14-17-00219.
shallow water model. In addition, significant dependence of the surface wave amplitudes versus the landslide width is demonstrated. Good agreement of numerical solutions with experimental data obtained in the problem of interaction of solitary wave with the conical island shows the applicability of the model in complex multi linked domains. Keywords: underwater landslide, surface waves, tsunami, shallow water equations, nonlinear dispersive equations, numerical simulation, finite-difference schema.
Received 29 May 2014
