Об алгоритме расчёта поверхностных волн в рамках нелинейно-дисперсионной модели на подвижном дне Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 17, № 5, 2012

Об алгоритме расчёта поверхностных волн

в рамках нелинейно-дисперсионной модели

>к

на подвижном дне*

О. И. Гусев

Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: gusev_oleg_igor@mail.ru

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

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

Введение

В последнее время широкое распространение получила гипотеза о том, что аномально большие для сейсмического происхождения волны цунами образованы с участием схода оползней [1, 2]. Численному моделированию данного явления посвящено много работ, при этом применяются различные модели для описания движения оползней и процесса образования и дальнейшего распространения поверхностных волн [3, 4]. Имеются также попытки экспериментального изучения поверхностных волн, образованных движением твёрдой модели оползня по плоскому подводному склону [5 - 7].

В существующих оползневых моделях можно выделить несколько основных подходов: оползень представляется как твёрдое тело, скользящее по плоскому откосу [5-8], или совокупность таких тел [9], как жидкость, отличающаяся от воды большими плотностью и вязкостью [10, 11], как некая упругопластическая среда. Для детального моделирования необходимо учитывать неровности дна: исследование некоторых случаев в [12] показало сильную зависимость образованных волновых режимов от геометрии подводного склона. В работе [13] предложен закон движения оползня по неровному склону с учётом сил тяжести, плавучести, трения и сопротивления воды. Этот закон движения применяется и в настоящем исследовании.

Моделирование поверхностных волн в рамках полных моделей гидродинамики требует большого времени расчёта, поэтому широкое применение получили приближённые модели. Классическая нелинейная модель мелкой воды (НЛ-модель) не воспроизводит дисперсию и не отражает неоднородность процесса в вертикальном направлении, поэтому для детального моделирования желательно использовать нелинейно-

* Работа выполнена при финансовой поддержке РФФИ (гранты № 10-05-91052-НЦНИ и 12-01-00721).

дисперсионные (НЛД)-модели. В работах [5, 14] применялись модели типа модели Бус-синеска [15], при выводе которых предполагалась малость амплитуды волн. При использовании НЛД-уравнений, свободных от этого ограничения (назовем их полными НЛД-моделями), точность воспроизведения волнового поля повышается, что показано в работе [16] путём сравнения с расчётами по полной гидродинамической модели, которая с высокой точностью воспроизводит как амплитудный, так и частотный характер волнового режима. Наиболее известными полными НЛД-моделями являются уравнения Грина — Нагди [17], модели Железняка — Пелиновского [18, 19], Алешкова [20] и Лью — Лайнетта [21].

В работе [19] получена модель Железняка — Пелиновского в предположении, что дно является неподвижным. НЛД-уравнения отличаются от классических уравнений мелкой воды наличием производных высокого порядка от искомых функций, что значительно усложняет численную реализацию модели. Один из способов решения таких систем — расщепление на эллиптическую и гиперболическую части, которые решаются на каждом шаге по времени. В [22] представлен алгоритм, основанный на одном из возможных способов расщепления системы из [19].

Обобщение модели Железняка — Пелиновского на случай подвижного дна проведено в [23]. Однако численной реализации этой модели до сих пор не было. Настоящая работа посвящена восполнению этого пробела для случая одной пространственной переменной. Проведено расщепление системы НЛД-уравнений на эллиптическую и гиперболическую части аналогично [22], описан конечно-разностный алгоритм решения этих частей. На ряде модельных задач численное решение сравнивалось с аналитическими и экспериментальными данными [24, 25], результатами, полученными по НЛ-модели и полной гидродинамической модели. Исходя из результатов сравнений установлена важность учёта дисперсионных свойств в исследуемом классе задач, а также оцениваются параметры оползня, при которых разработанная НЛД-модель даёт удовлетворительное соответствие с полной гидродинамической моделью.

1. Расщепление нелинейно-дисперсионных уравнений на две части

Настоящий раздел посвящён алгоритму расщепления нелинейно-дисперсионных уравнений Железняка — Пелиновского на две части и выводу некоторых замечаний о корректности поставленной задачи.

1.1. Исходные нелинейно-дисперсионные уравнения

Пусть декартова система координат Oxz расположена таким образом, что ось Ox совпадает с невозмущённой свободной поверхностью воды, а ось Oz направлена вертикально вверх. Свободная поверхность задаётся уравнением z = n(x,t), а подвижное дно уравнением z = -h(x,t), где t — время. Жидкость считается несжимаемой.

Подробный вывод нелинейно-дисперсионных уравнений Железняка — Пелиновского в случае подвижного дна приведён в [23], здесь лишь отметим, что при выводе существенным является предположение о потенциальности трёхмерного течения.

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

Н + (Ни)х = 0,

1 (Н3 Н2 \ / Н

Щ + + Пх = Н ^ "У + ^ - "х "2+ Д2

д д

где = их + иихх — иХ, Я2 = Б2", Б = ——Ь и—, и — осреднённая по глубине

д£ дх

горизонтальная компонента скорости, Н = п + " — полная глубина. Обезразмеривание проводилось по формулам

_ х - Г^ т " _ п -¡т Н _ и х = ^ 1 = tVlr, " = 7Г' п = 7-, Н = 7-, и =

V "о "о "о "о

где д — ускорение свободного падения, "0 — характерная глубина. Для упрощения записи в системе (1) и далее черта над безразмерными величинами опущена.

Отметим, что первое уравнение системы (1), уравнение неразрывности, имеет такой же вид, как и в классической нелинейной теории мелкой воды, а второе — уравнение движения, отличается правой частью.

Описанная система уравнений дополняется начальными и краевыми условиями. Будем рассматривать слой жидкости, ограниченный слева и справа вертикальными непроницаемыми стенками, расположенными в точках х = 0 и х = Ь. Краевые условия в этих точках принимают вид

и(0^) = 0, и(_М) = 0, t > 0. (2)

Предполагается, что в начальный момент времени заданы скорость жидкости и форма ее свободной границы:

и(х, 0) = и0(х), п(х, 0) = п0(х), х € [0,Ь]. (3)

1.2. Расщепление на эллиптическую и гиперболическую части

Обозначим ускорение частиц жидкости через d = Би = щ + иих. Тогда Д1 и Я2 можно переписать в виде

Д1 = dx — 2иХ, Я2 = + и2"хх + В, (4)

где В = Ни + 2и"х4 — слагаемое, отвечающее за подвижность дна. Введём новые переменные

Н3 Н2 Н2

У = -3- Я1 + ^ ^2, Ф = — Д1 + НЯ2. (5)

р 12^ — 6фН р —+ 4фН Л1 =-—-, ^2 =-—-, (6)

а уравнение движения в новых переменных примет вид

d = —н--пх. (7)

В результате подстановки выражений (6) и (7) в (4) получим систему уравнений для ^ и Ф

Ух — Ф"х \ р. 2 _ 12^ — 6фН

Тогда

Н пх ) 2их Н3 '

/ ж

^Х — ф"х ^ 2> „ 6У 4ф

-Н--пП "ж + Щ "хх + в = — — +

Из второго уравнения получаем

ф = ^ + пи + 6Р,

г Иг

где г = 4 + к2х, П

/Щ> /^хх I Пх /

. Подставляя это выражение в первое уравнение

системы, получаем одно уравнение относительно одной неизвестной функции р

(к\Рх)х - кор = Р,

в котором

к\ = тт", к

иг

Р9к 01 , 107

--+12коо,

дх

ко1 =

Кх

и2г

коо =

г_з 6

р = (КхП+Пх)х - ип+2'

Уравнение (9) — обыкновенное дифференциальное уравнение второго порядка. Дополним уравнение (9) краевыми условиями. Они получаются при подстановке граничного условия (2) и выражения (8) для ф в уравнение (7) и после элементарных преобразований принимают следующий вид:

4

иг

Рх

6Кх

иг

р = НХК + Пх |х=0; х=Ь ■

:ю)

Решив уравнение (9), по формуле (8) определяем функцию ф и далее решаем систему гиперболического типа

и + (И')х = о,

щ + иих + Пх которая в дивергентной записи имеет вид

Рх - фКх

и

д и ^

дЬ дх

:12)

где и — вектор решения, ^ — вектор потоков

{(и) =

и

и ии

ии И и2 + И2/2

О

о

Рх - фКх + ИКх

1.3. Некоторые замечания о корректности задачи (9), (10)

Предполагая, что полная глубина И всюду положительна, а дно описывается функцией с ограниченной производной |Кх| < С, получаем, что коэффициент к\ в уравнении (9) удовлетворяет оценке

1 4 1

0 <яТГС £ ^ £ И

Выполнение аналога данного условия в двумерном случае дает гарантию того, что уравнение для р является равномерно эллиптическим. В связи с этим будем называть ОДУ (9) "эллиптическим". Рассмотрим теперь коэффициент к0 в уравнении (9):

ко = 6

Кх \ 2 г - 3 Жг)^ + И3 г

6

И 3г2

КххИ(4 - Н1) + 2(4 + /х)(1 - КхПх)

:13)

Для корректности рассматриваемой краевой задачи для ^ достаточно, чтобы k0 был положительным [26]. Проверим выполнение этого условия на некоторых простых примерах.

Если дно горизонтальное (h = h0), то

3

ko = H >

Для плоского откоса (h = ax + ß) имеем

12

k0 = H3(4 + a2)(1 - аПх). Поэтому k0 > 0 при выполнении условия

anx < 1, (14)

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

Для корректности постановки граничных условий потребуем, чтобы коэффициент при ^ имел вполне определённый знак [26, 27], а именно, чтобы он был неотрицательным в граничной точке x = 0 и неположительным при x = L. Эти условия эквивалентны неравенствам

hx(°,t) > °, hx(L,t) < 0, (15)

означающим, что с приближением к границам глубина бассейна уменьшается.

2. Конечно-разностный алгоритм

В настоящем разделе подробно описана конечно-разностная схема, предназначенная для решения задачи для а также даны некоторые замечания об алгоритме решения всей задачи.

2.1. Разностные уравнения эллиптической части во внутренних узлах

Пусть расчётная область П = [0, Ь] покрыта равномерной сеткой П с узлами х = 0,... , N) и шагом к = Ь/Ж. Разностную схему для задачи (9), (10) построим интегро-интерполяционным методом.

Для получения разностных уравнений во внутренних узлах х (^ = 1,... , N — 1) рассмотрим элементарную ячейку сетки ^ = [ж^-1/2, с центром х, где =

X ± к/2, и проинтегрируем уравнение (9) по этой ячейке:

х .7 + 1/2 х.+1/2 х .7 + 1/2

У ¿ж — ^У (ко1)х ^¿ж — 12 У коо^ж =

Х 7 —1/2 х 7 —1/2 х 7 —1/2

xj + 1/2 xj + 1/2 xj + 1/2

R

„ „ „ H

xj-1/2 xj-1/2 xj-1/2

2 uXdx + (hxR + nx)x dx — 6 / — dx. (16)

В зависимости от выбора аппроксимационных формул для входящих в равенство (16) производных и интегралов будут получаться те или иные разностные уравнения. Здесь будут использованы следующие аппроксимации:

х3 + 1/2

(к1Рх)х -X ~ (к1 Рх)^-+1/2 - (к1Рх)^-1/2 ,

х3-1/2

х 3 + 1/2

(ко1)х Р-х ~ (ко1,^+1/2 - ко1,^-1/^ Ру

1/2

Г^'+1/2

кооР(х ~ К

Г^'+1/2

к00,у+1 /2 + к00,у-1/2

Ру, 2 и2х(х ~ К

и2 + и2 их,У+1/2 + их,У-1/2

х 3 1/2

х 3 1/2

х3 + 1/2

(КхП + Пх)х —х ~ (КхП + Пх)у+1/2 - (КхП + Пх)у-1/2,

х3-1/2

х3 +1/2

6 / И-х ~ зл^+Иг^) •

3 И \ИУ±1/2 И 3-1/2]

где

х3-1/2

_ Ру±1 - Рз Рх,3±1/2 = К , к1,3±1/2

к

К

01,3+1/2 7772

х,3±1/2

и

х,3±1/2

(И 2г)у±1/2

из±1 - из

коо

■оо,з±1/2 = 7773

(Иг)3±1/2 ! г3±1/2 - 3

(И3г)у±1/2 ;

к

г 3+1/2 = 4 + Кх,3±1/2,

К(х3±1) - К(х3)

х,3±1/2

К

К

п = Пз±1 - П п

Цх,з±1/2 = -Т-, п3±1/2

К

К

К3±2 - К3±1 - + К3-1

, = . 2 К2 :

;,3±1/2 1 К±1 - 2К(ху±1/2) + Кз

'и/ Кхх I Пх К^х

3+1/2

К2

В3±1/2 = (Кй + 2иКхг)з±1/2 , Кít,j±1/2

^ = 1,...,^ - 2, ^ = 0, ^ = N - 1,

Кй,3±1 + Кй,3

К

КП+1 - 2КП + КГ1

3 ' 3

П > 1,

хг,3±1/2

Ъп±1 — Ъп

х,3±1/2 х,3±1/2

Т

т — шаг по времени. Поскольку величины п,И,и всегда определены в целочисленных узлах х3, то применяются осреднения

П3+1/2

П±1 + П3 2 :

И3

3+1/2

И3±1 + И3

и3+1/2

и3±1 + и3 2

4

г

2

2

2

При использовании выписанных аппроксимаций интегральное соотношение (16) пе-

реходит в разностное уравнение

где

а,— с, у, + у^ = , з = 1,..., N — 1,

а, = ^1,^-1/2 > 0, = ^1,^+1/2 > 0,

с,- = а, + + 6к

к01,,+1/2 — к01,,-1/2 + к( к00,,+1/2 + к00,,-1/2

= (их,,+1/2 + Мх,,-1/2] к + к [ (кхЯ + Пх),+1/2 — (кхЯ + Пх), —1/2

Я,+1/2 Я,-1/2

-3к2

Я,

+

,+1/2

Я,-1/2

2.2. Аппроксимация краевых условий

:17)

:1в)

Получим теперь разностное уравнение вида (17) в граничной точке ж = 0. Для этого рассмотрим элементарную ячейку ш0 = [0, к/2], примыкающую к левой границе области. На этой ячейке аналог интегрального соотношения (16) запишется как

/1/2

(^Ух) (к/2) — (к^) (0) — бМк/2М0) + 6М0М0) — 12 / &00У«х =

//2

//2

Я

= 2 + (кхЯ + Пх) (к/2) — (кхЯ + Пх) (0) — 6 Я«х

19)

Заметим, что краевое условие (10) для узла ж0 примет вид

(&1Ух)(0) — 6^01(0)^(0) = (кхЯ + Пх)(0),

и подчёркнутые члены в уравнении (19) сокращаются. Тогда, используя вышеприведённые аппроксимации, разностное уравнение для у в узле з = 0 запишем в виде

—с, У, + У,+1 = ,

(20)

где

= к1,,+1/2 > 0, с, = + 6к

к01 ,,+1/2 + кк00,,+1/2 2 Я,+1/2

= их,,+1/2к + к (кхЯ + Пх),+1/2 — 3к

Я,+1/2

(21)

Аналогично получается разностное уравнение в самом правом сеточном узле с номером з = N:

а У,-1 — с, У, = ,

(22)

в котором

к1,,-1/2 > 0, с, = а, + 6к

— к01,,-1 /2 + кк00,,-

00,,-1/2

= «х,,-1/2к — к (кхЯ + Пх),-1/2 — 3к

2 Я,-1/2 Я,-1/2

(23)

Таким образом, получена система из ^ + 1) линейных уравнений (17), (20), (22) относительно ^ + 1) неизвестных у,, з = 0,... , N. Для поиска решения этой системы используется метод прогонки.

а

,

2.3. О корректности метода прогонки

Рассмотрим коэффициенты системы уравнений (17), (20), (22). В силу (15) будут выполняться неравенства

^01,1/2 > 0, ко1,м-1/2 < 0. Поскольку всегда к00,3±1/2 > 0, то из формул (21), (23) следует

со > Ьо > 0, ем > ам > 0, (24)

что с учётом неравенств в (18) свидетельствует о положительности всех коэффициентов

а3 и Ь3.

Рассмотрим теперь условие диагонального преобладания [27, 28] во внутренних узлах:

|С3 |

а3 + Ь3 + 6К

к01,3±1/2 - к01,3—1/2 + К( к00,3±1/2 + к00,3-1/2

> 3 + Ь3. (25)

Если дно ровное, то к01 = 0, поэтому е3 > а3 + Ь3 > 0, т.е. имеет место строгое диагональное преобладание. В случае неровного дна для выполнения неравенства (25) достаточно потребовать, чтобы выражение в квадратной скобке было неотрицательным:

к01,3±1/2 - к01,3—1/2 + ^к00,3±1/2 + к00,3—1/2) > 0. (26)

к01,3±1/2 - к01,3—1/2 + К ( к00,3±1/2 + к00,3 — 1/2) = (Л , 3—2\ и2-3+2-Х

Для плоского откоса К(х) = ах + К0 (а = 0) левая часть этого неравенства принимает вид

И3—1/2 + И3+1/2 (4 + а2)И— 1/2И3±1/2

\ (и и Л , и/л , ^ И3—1/2 - И3 — 1/2И3±1/2 + И3±1/2 1

Х\ а (И3—1/2 - И3±1/2) + К(1 + а2) 3 ------3-!- \ .

[ И3 — 1/2И3±1/2 J

Знак данного выражения определяется знаком выражения в фигурных скобках:

а (И_,/2 - И±,/2) + + а2) (И—1/2 - 1 + ИИ2) >

\И 3+1/2 И3 —1/2/

—1

( К К \

> а (-Кцх,] + а(х3 - + К-0 - а(х3 + - Ки + К(1 + а2) = к (1 - а^),

где

= П 3+1/2 - П 3—1/2

Пх3 = к .

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

аПх, 3 < 1,

которое можно считать дискретным аналогом полученного ранее условия (14) и которое будет всегда выполняться для пологих волн, в частности, в состоянии покоя жидкости с невозмущённой свободной границей.

Вместе с тем выписанное неравенство (26) выполняется не для любого профиля дна. Рассмотрим простейший пример.

Пример 1 (о нарушении условия диагонального преобладания). Пусть дно задаётся функцией у = —Ь(ж), где

, , ( аж + Ьо, 0 < ж < жо, ( ,

Ь(ж) = ^ 0 (27)

[1, ж0 < ж < Я 47

1 — Ьо

Здесь 0 <ж0 < Ь, 0 < Ь0 < 1, а = - > 0, т.е. дно представляет собой плоский

жо

откос, сопрягающийся в точке ж0 с горизонтальной плоскостью.

Предположим, что точка сопряжения ж0 совпадает с узлом сетки, имеющим номер Тогда

к01,,0+1/2 = 0, к01,,0-1/2 = . 2 7/2-,

4 + а Н .70 —1/2

, _ 1 1 _ 1 + а2 1

к00,.0+1/2 = 4 Нз , к00,.о —1/2 = 4 + а2 Н

4 Я .0+1/2 4 + а Н.0 —1/2

и условие (26) примет вид

а 1 Ь /1 1 1 + а2 1 \

+ Ч4 Я^Т + 4 + а Я^ТУ

+ ^^4^- > 0,

4 + а2 Я.0—1/2 ^ Н0+1/2 4 + а2 Я.0—1/2

где Ь в отличие от величин Ь.+1/2, которыми обозначается глубина бассейна, — шаг равномерной сетки.

Нам достаточно показать, что неравенство (28) нарушается хотя бы в одном частном случае. В качестве последнего рассмотрим состояние покоя жидкости, когда П.+1/2 = 0, или, по-другому, Я.+1/2 = Ь.+1/2. При этом предположении получим

Я.0+1/2 = 1 Я.0 —1/2 = Ь.0 —1/2 = а ( ж0 — - ) + П0 = 1 — а - < 1,

и неравенство (28) запишется как

11

4 + а2 (1 — аЬ/2)

, и 1 + а2

—а + п-

1 — аП/2_

+ П > 0.

Для малых значений шагов П данное неравенство можно переписать следующим образом:

1 + аП + 0(П2) г о, П

4 + а2

или г

п

— а + П(1 + а2) + 0(П2)

+ 4 > 0

11

- +

4 4 + а2

+ 0(п2) >

а

4 + а2

При малых П неравенство выполняться не будет (правая часть положительна, а левая при П ^ 0 стремится к нулю).

Причина нарушения условия (26) состоит в том, что функция П(ж), заданная формулой (27), не является гладкой: в точке ж0 первая производная этой функции имеет разрыв первого рода. Пусть теперь функция П(ж) является достаточно гладкой. Тогда условие (26) можно переписать в виде

ПххЯ(4 — ПХ) + 2(4 + ПХ)(1 — Мх) + О(п) > 0.

2

Пренебрегая членом порядка О(П) и считая жидкость покоящейся, получим ограничение на форму дна:

4 + Ь2

Ь > -2_ х__(29)

Ьхх > 2Ь(ж)(4 — ЬХ). ( )

В силу неравенства

2 4 + *Х > 2

Ь(ж)(4 — ЬХ) Ьтах'

где Ьтах > 0 — максимальное значений функции Ь(ж), приходим к выводу о том, что для выполнения (29) достаточно потребовать, чтобы

2

Ьхх > — Ь-. (30)

Ьтах

Выполнение условия (30) достигается путём сглаживания функции у = — Ь(ж) в окрестности точек, в которых происходит потеря гладкости этой функции.

Пример 2 (о выполнении условия (30) для гладких функций). Сгладим кусочно-линейную функцию (27) в окрестности точки ж0 с помощью параболы:

(аж + Ь0, 0 < ж < жг, —а(ж — ж*)2 + 1, жг < ж < ж*, 1, ж* < ж < Ь,

где

т 1 — Ьо — аж*/2 а

жо < ж < 2жо < Ь, жг =-г--, а = —-г > 0,

а/2 2(ж* — жг)

ж*, жг — точки сопряжения параболы соответственно с горизонтальным участком дна и с плоским откосом, 0 < жг < ж0.

В рассматриваемом примере Ьтах = 1, поэтому для выполнения неравенства (30) на участке сглаживания необходимо выполнение неравенства а < 1 , или

ч , 4(1 — Ь0)

а < 2(ж — жг) = 4ж*--,

а

т.е. 4ж* > а + 4ж0. Таким образом, в данном примере достаточно выбрать

а

ж* > жо + 4. 2.4. Алгоритм решения полной задачи

Поскольку, как уже отмечалось, гиперболическая задача (12), (2) отличается от НЛ-модели только правой частью, для её решения будет использована хорошо зарекомендовавшая себя в этом классе задач явная схема типа предиктор-корректор второго порядка аппроксимации [29]. В настоящей работе мы не будем останавливаться на её подробном описании, а только отметим, что на шаге "предиктор" аппроксимируется

/* а }}

*+1/2, а на шаге корректор — вектор искомых

функций в целых узлах и..

Принцип численного решения расщеплённой системы на каждом временном шаге следующий. Сначала решается задача (9), (10) для у , где Я™, и™ берутся с предыдущего шага по времени или, если шаг первый, из начальных данных (3). Затем полученные

значения используются на шаге "предиктор" при вычислении f*+1/2 и и*. После этого снова решается эллиптическая задача, в которой вместо H, u используются ранее посчитанные H*,u*. Далее проводится заключительный этап — "корректор", на котором определяются искомые функции H™+1,u™+1.

3. Тестовые задачи

Вначале рассмотрим результаты решения тестовых задач о движении волн над неподвижным дном.

3.1. Распространение уединённой волны

Нетрудно показать, что в случае ровного дна (h = 1) при начальных данных

/ ч _a__/ ч UoVojx) . л

40<x) = cosh2(a(x - x0)) • U0(x) = ТГПф) (31)

система нелинейно-дисперсионных уравнений (1), записанная для безразмерных переменных, имеет аналитическое решение в виде солитона (рис. 1):

V(x,t) = ^^~, u(x,t) = ^HM-, (32)

cosh2 X H <x,t)

где a > 0 — амплитуда начальной волны (31), вершина которой располагается в точке

x = xo, U0 = л/1 + a, a = , X = a(x — x0 — U01). Отметим, что (32) является

2Uo

решением задачи Коши с начальным условием (31), но не точным решением краевой задачи с граничным условием (2). Поэтому сравнение численного решения с (32) имеет смысл только вдали от границ области.

Рис. 1. Профили свободной поверхности: 1 — в начальный момент времени, 2, 3 — при Ь = 200 с; 2 — аналитическое решение (32), 3 — численное решение при N = 400

В тестовой задаче был взят бассейн длиной Ь = 10000 м и глубиной Ьо = 100 м. Вершина волны (31) высотой а =10 м располагалась в точке жо = 8500 м (см. рис. 1). Волна движется влево до момента времени £ = 200 с. Из рисунка следует, что численное решение для N = 400 хорошо соответствует аналитическому. Численные эксперименты показали, что для протяжённых областей приближённое решение при измельчении сетки сходится к точному решению (32).

3.2. Накат уединённой волны на вертикальную стенку

Пусть в центре области с параметрами Ь = 3000 ми Ьо = 100 м в начальный момент времени располагается волна вида (31), которая при £ > 0 движется с постоянной скоростью влево, совершает накат на левую вертикальную стенку, отражается от нее и затем движется вправо. На рис. 2 изображены профили волны в конечный момент времени £ = 96 с для амплитуд а =10 и 50 м. Заметим, что при откате волны амплитудой 10 м форма профиля сохраняется, а с увеличением амплитуды картина течения усложняет-

Рис. 2. Профили свободной поверхности после отката: а = 10 м (а), 50 м (б)

Рис. 3. Зависимость относительного максимального заплеска на стенку от начальной амплитуды уединённой волны: 1 — НЛД-модель настоящей работы; 2 — модель мелкой воды; экспериментальные данные [24] (3) и [25] (4); 5 — НЛД-модель [30]

ся: головная волна немного понижается, за ней появляется дисперсионный "хвост" из волн малой амплитуды. Такой эффект при использовании НЛ-модели не наблюдается.

На рис. 3 для разных амплитуд а приведено сравнение полученного максимального заплеска Я на стенку с экспериментальными данными [24, 25], расчётами по НЛД-модели [30] и НЛ-модели. Видно, что при малых амплитудах (а/Ьо < 0.3) НЛД- и НЛ-модели дают достаточно хорошее соответствие с экспериментом, но при больших амплитудах НЛ-модель существенно занижает величину максимального заплеска.

3.3. Накат уединённой волны на плоский откос

В этой задаче волна вида (32) накатывается на плоский откос с углом а (рис. 4, а). Откос заканчивается вертикальной стенкой на глубине Ьге/ = 10 м, в то время как максимальная глубина составляет Ьо = 100 м. Рассмотрим картины отката волны амплитудой а =10 м в момент времени £ = 320 с при разных углах а (рис. 4, б). Из последнего рисунка следует, что при малых углах откоса картина отката усложняется: за головной волной появляется большое число волн дисперсионного "хвоста", причём амплитуды этих волн увеличиваются по сравнению со случаем горизонтального дна.

Рис. 4. Профили свободной поверхности: а — в начальный момент времени; б — при Ь = 320 с после отражения волны от вертикальной стенки; а = 2° (1), 5° (2), 10° (3)

Рис. 5. Зависимость относительного максимального заплеска от угла склона при амплитуде начального возвышения а = 1 м (а), 5 м (б), 10 м (в). Сплошные линии — результаты расчётов по НЛД-модели, штриховые — по модели мелкой воды

а

Рассмотрим теперь максимальные заплески на стенку в зависимости от угла а при разных начальных амплитудах волны а = 1, 5 и 10 м (рис. 5). При малой амплитуде максимальные накаты, полученные при помощи НЛ-модели и НЛД-модели, практически совпадают при всех рассмотренных значениях а. Для средней амплитуды наблюдаются небольшие различия полученных результатов, приблизительно равные при всех углах. Картина меняется в случае а = 10 м: при а > 5° НЛ-модель немного завышает значение заплеска относительно НЛД-модели, а при а < 5° начинает занижать, причём различие увеличивается с уменьшением угла и для очень пологого откоса достигает 15 м. Таким образом, можно сделать вывод, что при малых амплитудах в настоящей задаче различия между обеими моделями незначительны, а при амплитудах порядка 10 м и малых углах а важно учитывать нелинейно-дисперсионные свойства явления. Отметим также, что максимальный заплеск увеличивается с уменьшением угла склона.

4. Генерация волн подводным оползнем

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

Выделим в уравнении, описывающем дно, подвижную и неподвижную части:

г = -к(х,Ь) = кы(х) + hsl (х,Ь),

где отвечает за неподвижное дно акватории, а hsl — за оползень.

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

г = Ьг(х) = - ^X - + ^, (33)

где — глубина в точке х = 0, ^ — максимальная глубина, х = £ — точка вершины параболы. Акватория ограничена вертикальными стенками в точках х = 0 и х = 2£. Данная батиметрия соответствует условиям (15). Выберем следующие значения параметров:

= -10 м, ^ = -100 м, £ = 500 м. Форма оползня в начальный момент времени задаётся формулой

T

2n(x - x0)

cos '

|x - X0|< b/2,

z = (x)=< 2

0, |x - x0| >b/2,

где T — толщина, b — длина оползня, а x0 — начальное положение его центра.

Важными параметрами в законе движения являются также Cw — коэффициент присоединённой массы, Cd — коэффициент сопротивления воды, Cfr = tan в* — коэффициент трения скольжения, где в* — угол трения, y > 1 — отношение плотности оползня

к плотности воды. Влияние этих параметров на образование поверхностных волн в рамках модели мелкой воды исследовано в [13], здесь же, если не оговорено иное, берутся следующие их значения:

Сш = 1, Сл = 1, в* = 5°, 7 = 2, Т =10 м.

Для определения модельной задачи остаётся выбрать два параметра: длину Ь и начальное заглубление г0 оползня. Рассмотрим случай Ь = 100 м, г0 = -30 м. При таких параметрах оползень сдвигается с откоса вниз и в момент времени £ ^ 70 с останавливается почти в центре области. Отметим, что в некоторых случаях, например при уменьшении коэффициента трения, он может несколько раз пересечь центр водоёма, двигаясь то влево, то вправо, и лишь потом остановится.

Рассмотрим профиль образовавшейся волны в момент времени £ = 20 си марео-граммы, установленные на границах области за время £ = 100 с (рис. 6). На графиках результаты полученные по НЛД-модели и НЛ-модели, сравниваются с расчётами модели потенциальных течений (МПТ) [22], которую можно считать эталонной в данном классе задач [5]. Из рис. 6, а видно, что результаты НЛД-модели и МПТ почти совпадают правее оползня, который находится под впадиной, и левее него воспроизводят образование второстепенных волн. НЛ-модель немного точнее описывает впадину над оползнем, неплохо согласуется с МПТ правее оползня, но не описывает образование волн за ним. По мареограмме около левой границы можно сделать вывод, что первую волну понижения почти одинаково описывают все модели, но следующее за ней повышение значительно отличает НЛ-модель от остальных, которые показывают хорошую согласованность между собой. На мареограмме вблизи правой границы так же хорошо согласуются НЛД-модель и МПТ на всём временном промежутке, а НЛ-модель сильно упрощает картину, при этом значительно завышая амплитуду первой волны.

Проанализируем изменение результатов моделей при уменьшении длины оползня до Ь = 50 м (рис. 7,1). Из рисунка следует, что профиль в момент времени £ = 20 с у НЛД-модели сильно отличается от остальных: впадина над оползнем значительно больше, за оползнем появились аномально большие волны. Тем не менее это не сильно повлияло на заплески на стенку, где можно наблюдать удовлетворительное соответствие

Рис. 6. Профиль свободной поверхности при Ь = 20 с (а) и мареограммы на левой (б) и правой (в) границах области: 1 — НЛД-модель; 2 — модель потенциальных течений; 3 — НЛ-мо-дель. Длина оползня Ь = 100 м, начальное заглубление г0 = -30 м

Рис. 7. Профиль свободной поверхности при Ь = 20 с (а) и мареограммы на левой (б) и правой (в) границах области: 1 — результаты НЛД-модели; 2 — МПТ; 3 — НЛ-модели. Начальное заглубление г0 = -30 м, длина оползня Ь = 50 м (I), г0 = -60, Ь = 100 м (II)

с МПТ, хотя, несомненно, худшее, чем при Ь = 100 м. Похожая картина получается при увеличении начального заглубления г0 до 60 м (рис. 7,11).

Из данных примеров можно сделать вывод, что для адекватного описания профиля волны, образованной движением подводного оползня, при помощи предложенной НЛД-модели необходимо, чтобы оползень был длинный и при этом находился на малой глубине. Подобное ограничение для НЛД-модели Лью — Лайнетта получено в [21].

Заключение

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

На тестовой задаче о распространении уединённой волны показано хорошее соответствие численного решения аналитическому. В задаче о накате волны на вертикальную стенку при сравнении полученного решения с экспериментальными данными и результатами по НЛ-модели показана важность учёта нелинейно-дисперсионных свойств при больших начальных амплитудах. Значительные различия между НЛД- и НЛ-моделями обнаружены и при накате волны большой амплитуды на пологий плоский откос.

Проведена серия расчётов генерации волн подводным оползнем в модельной ограниченной акватории. При сравнении результатов, полученных по НЛД- и НЛ-моделям, с расчётами по модели потенциальных течений показано, что при оползнях достаточной длины и малых заглублениях НЛД-модель хорошо согласуется с МПТ, в отличие от НЛ-модели детально описывая картину течения. Однако в случаях оползней малой длины и больших заглублений предложенная НЛД-модель не позволяет адекватно описать процесс образования волн над оползнем, давая сильное расхождение как с МПТ, так и с НЛ-моделью.

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

[1] Watts P., Grilli S.T. Underwater landslide shape, motion, deformation, and tsunami generation // Proc. of the 13th Intern. Offshore and Polar Eng. Conf. Honolulu, Hawaii, 2003. Vol. 3. P. 364-371.

[2] Tappin D. R., Watts P., Grilli S. T. The Papua New Guinea tsunami of 17 July 1998: Anatomy of a catastrophic event // Natural Hazards and Earth System Sci. 2008. Vol. 8. P. 243-266.

[3] Шокин Ю.И., Чубаров Л.Б. О подходах к численному моделированию оползневого механизма генерации волн цунами // Вычисл. технологии. 2006. Т. 11. Спец. выпуск. Часть 2. С. 100-111.

[4] Федотова З.И., Чубаров Л.Б., Шокин Ю.И. Моделирование поверхностных волн, порожденных оползнями // Там же. 2004. Т. 9, № 6. С. 89-96.

[5] Елецкий С.В., Майоров Ю.Б., Максимов В.В. и др. Моделирование генерации поверхностных волн перемещением фрагмента дна по береговому склону // Там же. 2004. Т. 9. Спец. выпуск. Часть 2. С. 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] Pelinovsky E., Poplavsky A. Simplified model of tsunami generation by submarine landslides // J. Phys. Chem. Earth. 1996. Vol. 21, No. 12. P. 13-17.

[9] Tinti S., Bortolucci E., Vannini C. A block-based theoretical model suited to gravitational sliding // Natural Hazards. 1997. Vol. 16. P. 1-28.

[10] Heinrich P., Schindele F., Guibourg S., Ihmle P. Modeling of the february 1996 peruvian tsunami // Geophys. Res. Lett. 1998. Vol. 25, No. 14. P. 2687-2690.

[11] Savage S., Hutter K. The motion of a finite mass of granular material down a rough incline // J. Fluid Mech. 1989. Vol. 199. P. 177-215.

[12] 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. and Math. Modelling. 2011. Vol. 26, No. 1. P. 17-38.

[13] ХАкимзянов Г.С., Шокина Н.Ю. Численное моделирование поверхностных волн, возникающих при движении подводного оползня по неровному дну // Вычисл. технологии. 2010. Т. 15, № 1. С. 105-119.

[14] Watts P., Grilli S.T., Kirby J.T. et al. Landslide tsunami case studies using a Boussinesq model and a fully nonlinear tsunami generation model // Natural Hazards and Earth System Sci. 2003. Vol. 3, No. 5. P. 391-402.

[15] Дорфман А.А., Яговдик Г.И. Уравнения приближённой нелинейно-дисперсионной теории длинных гравитационных волн, возбуждаемых перемещениями дна и распространяющихся в бассейне переменной глубины // Численные методы механики сплошной среды: Сб. науч. тр. / ВЦ, ИТПМ СО АН СССР. 1977. Т. 8, № 1. С. 36-48.

[16] Shokin Yu.I., Fedotova Z.I., Khakimzyanov G.S. et al. Modelling surfaces waves generated by a moving landslide with allowance for vertical flow structure // Rus. J. of Numer. Anal. Math. Modelling. 2007. Vol. 22, No. 1. P. 63-85.

[17] 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.

[18] Железняк М.И., Пелиновский Е.Н. Физико-математические модели наката цунами на берег // Накат цунами на берег: Сб. науч. тр. / ИПФ АН СССР, 1985. С. 8-33.

[19] Вольцингер Н.Е., Клеванный К.А., Пелиновский Е.Н. Длинноволновая динамика прибрежной зоны. Л.: Гидрометеоиздат, 1989. 272 с.

[20] Алешков Ю.З. Течения и волны в океане. СПб.: Изд-во С.-Петербургского гос. ун-та, 1996. 226 с.

[21] 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.

[22] ХАкимзянов Г.С., Шокин Ю.И., Барахнин В.Б., Шокина Н.Ю. Численное моделирование течений жидкости с поверхностными волнами. Новосибирск: Изд-во СО РАН, 2001.

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

[24] Загрядская Н.Н., Иванова С.В., Нуднер Л.С., Шошин А.Н. Воздействие длинных волн на вертикальную преграду // Изв. ВНИИГ им. Б.Е. Веденеева: Сб. науч. тр. 1980. Вып. 138. С. 94-101.

[25] Манойлин С.В. Некоторые экспериментально-теоретические методы определения воздействия волн цунами на гидротехнические сооружения и акватории морских портов. Красноярск: Изд-во СО АН СССР, 1989. 45 с.

[26] Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1977.

[27] Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М.: Наука, 1978.

[28] ХАкимзянов Г. С., Черный С. Г. Методы вычислений. Ч. 2. Численные методы решения краевых задач для обыкновеных дифференциальных уравнений. Новосибирск: НГУ, 2005.

[29] Шокин Ю. И., Хлкимзянов Г. С. Схема предиктор-корректор, сохраняющая гидравлический скачок // Вычисл. технологии. 2006. Т. 11. Спец. выпуск. Часть 2. С. 92-99.

[30] Железняк М.И. Воздействие длинных волн на сплошные вертикальные преграды // Накат цунами на берег: Сб. науч. тр. / ИПФ АН СССР, 1985. С. 122-139.

[31] Диденкулова И.И., Пелиновский Е.Н. Цунамиподобные явления в российских внутренних водоёмах // Фундамент. и прикл. гидрофизика. 2009. № 3(5). С. 52-96.

Поступила в 'редакцию 22 мая 2012 г., с доработки — 30 июля 2012 г.