CC BY
45
Вычислительные технологии
Том 18, № 5, 2013
Моделирование поверхностных волн, генерируемых
>к
подводным оползнем в водохранилище*
О. И. Гусев, Н. Ю. ШокинА, В. А. Кутергин, Г. С. Хлкимзянов Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected], [email protected], [email protected]
Исследуется влияние дисперсии на картину поверхностных волн, возникающих при сходе подводного оползня в ограниченном водоёме. Для описания поверхностных волн используются как полная нелинейно-дисперсионная модель мелкой воды, так и новые приближённые модели типа модели Буссинеска для волн, генерируемых медленно спадающими либо малой высоты оползнями. Применяется единый для всех моделей подход к построению численных алгоритмов, основанный на аппроксимации расширенной системы уравнений, состоящей из системы уравнений гиперболического типа и уравнения эллиптического типа для негидростатической составляющей давления. Выполнено сопоставление численных результатов, полученных в рамках бездисперсионной модели мелкой воды, нелинейно-дисперсионных моделей и модели потенциальных течений со свободной границей.
Ключевые слова: подводный оползень, неровное дно, поверхностные волны, уравнения мелкой воды, нелинейно-дисперсионные уравнения, численный алгоритм.
Введение
Сход подводного оползня с крутого берегового склона водохранилища ГЭС может привести к возникновению поверхностных волн, опасных для береговых сооружений и плотины водохранилища. Подводные оползни генерируют более короткие волны [1], чем цунамигенные землетрясения, поэтому для описания таких волн необходимо учитывать дисперсию волн, отражающую в определённой степени неоднородность процесса в вертикальном направлении, т. е. использовать нелинейно-дисперсионные уравнения (НЛД-уравнения) [2, 3], способные воспроизводить дисперсию. Эти уравнения содержат смешанные производные высокого порядка, в связи с чем проблема конструирования эффективных численных алгоритмов для решения НЛД-уравнений не является тривиальной и ей в настоящее время уделяется пристальное внимание специалистов по вычислительной гидродинамике [4 - 8].
Стандартный подход к решению НЛД-уравнений заключается в вычислении полной глубины из уравнения неразрывности [9, 10] и решении системы обыкновенных дифференциальных уравнений относительно двух вспомогательных величин [4], которые используются затем в качестве правых частей для эллиптических уравнений относительно
* Работа выполнена при финансовой поддержке РФФИ (грант № 12-01-00721) и Программы интеграционных исследований СО РАН (проекты № 42 и № 117).
компонент скорости. В настоящей работе для решения системы НЛД-уравнений предлагается заменить её расширенной системой уравнений, состоящей из системы уравнений гиперболического типа, аналогичной системе уравнений мелкой воды первого гидродинамического приближения и отличающейся от последней лишь правой частью, и уравнения эллиптического типа для осреднённой по глубине дисперсионной составляющей давления. Этот подход к конструированию численных алгоритмов применяется здесь как для полных НЛД-уравнений, так и для новых приближённых моделей типа модели Буссинеска [11], описывающих поверхностные волны, генерируемые медленно сходящими либо малой высоты оползнями.
Представленная работа является продолжением исследований [12-14], посвящён-ных изучению влияния дисперсии на картину генерируемых оползнем поверхностных волн в прибрежной акватории морей. В отличие от указанных работ в нашем случае рассматриваются оползни в ограниченном водоёме. Исследовано влияние параметров, определяющих геометрию водоёма и движение квазидеформируемого оползня [15] по неровному склону, на величины максимальных заплесков на берег водохранилища. Выполнено сравнение с численными результатами, полученными по бездисперсионной модели мелкой воды [16] и модели потенциальных течений [17], а также с имеющимися экспериментальными данными [18-20] измерений высот волн, возникающих при движении твёрдых тел по плоскому подводному откосу.
1. Нелинейно-дисперсионные уравнения мелкой воды
Пусть слой несжимаемой жидкости ограничен снизу подвижным дном, заданным функцией г = -К(х,у,Ь), а сверху — свободной границей, описываемой функцией г = п(х,у,Ь), где Ь — время, х, у, г — координаты точки в декартовой системе координат Охуг, ось Ог которой направлена вертикально вверх, а плоскость Оху совпадает с невозмущённой свободной поверхностью.
В приближённых моделях мелкой воды искомыми величинами обычно являются полная глубина слоя жидкости Н = п + К > 0 и вектор скорости и = (и, у), связанный каким-либо образом с вектором скорости трёхмерного течения. Если в качестве и использовать осреднённую по глубине горизонтальную составляющую вектора скорости трёхмерного течения, то уравнение неразрывности приближённой модели принимает следующий вид:
Н + V ■ (Ни) = 0, (1)
где V = (д/дх, д/ду), V ■ и = их + уу.
В отличие от модели мелкой воды первого приближения в НЛД-модели имеет место [2] квадратичная зависимость давления р от вертикальной координаты
р = д
Н - (г + К)
(Н - (г + К)) Я2 + ( Н2 - М
Яг
где д — ускорение свободного падения,
Я2 = Б2К, Яг = Б (V ■ и) - (V ■ и)2
Б = д/дЬ + и - V — оператор полной производной. Если через Р обозначить среднее по толщине слоя давление (2)
п
1 Г , а- /Я2 Я
р = я/Р* = V "(."Я + ЯЯ
-н
то уравнение движения НЛД-модели примет вид [3]
щ + (и ■ У)и + Х1Нр) = д-К, Н
при этом
1
д = Нр
г=-Н 9 "I Н Н1 + Н2] • (6)
В отличие от классических уравнений мелкой воды в уравнения движения (5) входят смешанные производные по времени и пространственным переменным от компонент вектора скорости, что осложняет построение численного алгоритма. Используя известный в теории дифференциальных уравнений приём понижения порядка производных путём добавления к исходным неизвестным новых зависимых переменных, избавимся от этих смешанных производных высокого порядка, введя дополнительную зависимую переменную — дисперсионную составляющую
И3 и 2
V = -3- Я1 + — Д2 (7)
давления р, проинтегрированного по толщине слоя, и добавив к исходным уравнениям (1), (5) уравнение относительно новой зависимой переменной V. Для получения такого уравнения обозначим дисперсионную составляющую давления на дне через ф
Н2
ф = — Я1 + -Я 2. (8)
Тогда уравнение движения (5) можно записать в виде
„ -V - ф-К Би = -9— + -, (9)
где Би = (Бп,Бу). В этом случае функция ф выражается через новую переменную V и исходные зависимые переменные Н и и. В самом деле, из определений (7) и (8) следует
ф = ^ + НЯ2 • (10)
У2Н 4 2 ^ '
Если подставить сюда правую часть выражения для Я2
Я2 = Би --К + и ■ ((и --)-К) + В (11)
и вектор ускорения Би заменить правой частью уравнения (9), то после приведения подобных получим формулу
ф = 1 (И + НЯ + -V (12)
где
Я = -9-п --К + и ■ ((и --)-К) + В, В = Ны + 2 (и --К) , г = 4+ |-К|2 •
Уравнение для р выводится аналогично формуле (12) для функции ф. Из выражений (7), (8) следует
Н 3 Н
Р = 12 Яг + Нф, (13)
при этом
Яг = V- (Би) - 2 (V ■ и)2 + 2Vu XV«. (14)
В последней формуле использована бинарная операция для двумерных векторов, обозначенная так же, как операция векторного произведения трёхмерных векторов. Результатом действия этой операции на векторы Vu, Vu будет скаляр, определяемый по формуле Vu X V« = их У у иу «х.
Подставляя в (13) выражения (12), (14) и используя в (14) формулу (9), получим следующее уравнение для новой зависимой переменной р:
^ /V? (^ ^к) Vh\ / 2 г - 3 ^ / VII \\ ^
-нг ;-6Чн—+пНКг)) = Р (15)
где
Р = V- (gVп + —1 - Iя + 2 (V ■ и)2 - 2Vu XV«. г Нг
Итак, расширенная система НЛД-уравнений (КЬВ-модель) состоит из уравнения неразрывности (1), уравнения движения
и4 + (и ^)и + gVH = gVh + ^ - фVk (16)
Н
и уравнения (15) относительно дисперсионной составляющей р давления, проинтегрированного по глубине, при этом функция ф вычисляется по формуле (12).
Алгоритм расчёта на произвольном шаге по времени с номером п аналогичен описанному в [17] и состоит из двух шагов. На шаге предиктор сначала решается конечно-разностное уравнение, аппроксимирующее (15), и определяется функция рп. Коэффициенты уравнения вычисляются по известным значениям Нп и ип с п-го слоя по времени. Затем решается гиперболическая система уравнений (1), (16), при этом в правой части уравнения движения используются уже известные значения рп, фп, Нп и кп. Далее вновь численно решается уравнение (15) с использованием величин Н* и и*, вычисленных на предикторе. Найденные на этом этапе значения р* и ф* используются на шаге корректор для определения окончательных значений Нп+1 и ип+1 путём численного решения системы уравнений гиперболического типа, состоящей из уравнения неразрывности (1) и уравнения движения (16), записанного с дивергентной формой левой части:
•Н_2 2
где ии — тензорное произведение вектора и на себя. В случае ровного дна это уравнение принимает консервативный вид.
Преимуществом ЫЬЮ-модели является то, что система уравнений (1), (15), (16) допускает при численном решении расщепление на эллиптическую и гиперболическую части, что позволяет использовать для каждой из частей хорошо изученные алгоритмы [21]. Уравнение (15) является равномерно эллиптическим, и если коэффициент при р положителен (например, в случае горизонтального дна), то для нахождения численного решения этого уравнения с переменными коэффициентами можно построить
(Ни)4 + V ■ (Нии) + gV— = gHVh + V? - фVk, (17)
разностные схемы с положительно определёнными операторами. Левые части уравнений движения (16), (17) совпадают с левыми частями системы уравнений мелкой воды первого гидродинамического приближения, поэтому для численного решения гиперболической части (1), (16), (17) используется схема предиктор-корректор [22-24], хорошо зарекомендовавшая себя при исследовании волновых процессов в рамках бездисперсионной модели мелкой воды [25].
Приведём некоторые тривиальные решения выписанной системы уравнений, которые могут оказаться полезными при отладке численных алгоритмов. В случае горизонтального неподвижного дна
h(x,y,t) = h0 = const (18)
система уравнений (1), (15), (16) имеет решение "равномерный поток":
u = u0 = const, n = 0, H = h0, p = 0, (19)
для неровного неподвижного дна
h(x,y,t) = ho(x,y) (20)
существует решение "покой":
u = 0, n = 0, H = h0(x, y), p = 0, (21)
а для подвижного дна
h(x,y,t) = ho(x, y) + tb (22)
имеет место "вертикальный сдвиг":
u = 0, n = tb, H = ho(x, y), p = 0. (23)
Для одномерных течений, когда все зависимые переменные зависят только от одной горизонтальной координаты x и вектор скорости имеет только одну ненулевую компоненту u, система уравнений NLD-модели записывается в виде
Ht + (Hu), = 0, (24)
(Hu)t + ^Hu2 + g^ = gHhx + px - ^hx, (25)
( P. \ _ 6 f JL r - 3 ^ hx
4 - + Hri 1 = F (26)
где
1 "h,\ 6" 2 "(H + H" + p*4 , F = (gn, + "r) - ^ + 2uX,
x
Л = -д^х^х + + в, В = + г = 4 +
В случае горизонтального дна (18) уравнения (24)-(26) принимают наиболее простой вид
Я + (Яи)х = 0, (Яи)4 + (я«2 + д^) = Рх (- Н3Р = дПхх + 2«Х (27)
и имеют решение "уединённая волна":
п(х,£) = аоСовЬ 2( v/3aog (х - х0 - и0Ьь) ), Н(х,£) = к0 + п(х,Ь), (28)
V 2коЦо У
и(х,Ь) = ^0Н^, Р(х,Ь) = gH2(х,2) - К0 - ко^ои(х,Ь), (29)
где и0 — скорость перемещения волны, а0 — её амплитуда, х = х0 —
положение вершины волны при £ = 0.
2. Слабо дисперсионные уравнения мелкой воды
При выводе НЛД-уравнений (1), (5) предположение о малости амплитуды волн не использовалось. Приведём теперь приближённые НЛД-уравнения для волн малой амплитуды, допускающие в отличие от классических уравнений Буссинеска [26, 27] в качестве своего следствия закон баланса полной энергии [11]. В слабо дисперсионной модели типа модели Буссинеска (ЫЬВБ-модели) мы имеем [11] систему уравнений (1), (5) с модифицированными функциями (4) и (6)
Р = £ - (К Б ^ ■ и) + К Б2к) , д = g - (2 Б ^ ■ и) + Б2к) . (30)
Далее поступаем так же, как в случае полной НЛД-модели: записываем уравнение движения в виде (16) и конструируем уравнение для зависимой переменной р. В рассматриваемом случае формулы (7), (8) принимают вид
р = —Б и) + —Б2к, ф = НБ (hV■ и) + НБ2к, (31)
3 2 2
что приводит к изменению выражений (10), (13):
ф = 3р + НБ2к, р = —Б (hV■ и) + кф, (32)
^ 2к 4 ' у 12 1 ; 2^' 1 '
и небольшой коррекции формулы (12):
ф = ^6^ + НЯ + Vh^ . (33)
С учётом этого выражения, формулы (9), а также равенства
Б ^ ■ и) = (Бк) (V ■ и) + к (V ■ (Би) - (V ■ и)2 + 2Vu X V«) (34)
получаем уравнение вида (15) для дополнительной зависимой переменной р:
/ ^ _ (Ур- ук> Ук \ - бр /2 г-3 + / Ук\\ = р, (35)
Н Нг У ^Нк2 г + ДНкгу
где
^ „ { „ ЯVh\ 6Я Ч2 ч Бк ^ „
р = V ■ ( gVn +-) - + (V ■ и)2 - (V ■ и) — - 2Vu х V«.
Таким образом, структура уравнений полной ЫЬЮ-модели сохраняется и в случае слабо дисперсионной ЫЬЮБ-модели, для которой имеем систему уравнений гиперболического типа (1), (17) и равномерно эллиптическое уравнение (35). Поэтому описанный выше алгоритм численного решения полных ЫЬЮ-уравнений используется и для решения задач в рамках ЫЬЮБ-модели.
Следует отметить, что в одномерном случае гиперболическая часть уравнений слабо дисперсионной модели полностью совпадает с (24), (25), а уравнение для дисперсионной составляющей давления характеризуется более низкой, чем в (26), степенью нелинейности функции Н в коэффициентах и в правой части:
4 (Нг). - (Нр + (Нйг)*) = * (36)
при этом
^ ( ЯйЛ 6Д 2 Бй , 1/6. тт„ , , * = дп* +--- -¡— + и* - их—г~, ¥ = - ~т + НД + .Ж .
V Г / * йг П Г \ й /
Для горизонтального дна (18) изменится только последнее из уравнений (27):
(НО*- Нй!. = дп** +и*. (37)
3. Слабо дисперсионные уравнения мелкой воды в случае малых деформаций дна
Уравнения слабо дисперсионных течений над слабо деформируемым дном получаются [11] из уравнений ЫЬЮБ-модели при дополнительном предположении о слабой изменчивости дна во времени: й(я,у, ¿) = йо(я, у) + ей^я, у,£), где е — малое число [27]. В этой модели уравнение движения (5) принимает вид
и, + (и ■ У)и + = - йо) + дУйо,
Н
где
..N , й0гл2ъ \ _ (1
Р = ^ ^ У Б (йоУ ■ и) + 2 Б2^ , q = д - - Б (йоУ ■ и) + Б2^ . (38)
Из (38) следует, что формулы для функций . и ¥ в случае слабо деформируемого дна получаются из (31) простой заменой й на йо:
. = НгБ (^ ■ и) + НйоБ2йо, ¥ = |Б (йоУ ■ и) + НБ2йо. (39)
Подставляя в (34) функцию йо и используя уравнение движения, переписанное в виде
и, + (и -У)и + д*, = -Н¥Уйо, (40)
получим следующее выражение для ¥:
¥ + НД + V. -УйЛ , (41)
г \йо /
где
Я = ^^ ■ Vko + и ■ ((и ■ V)Vko), г = 4 + .
Таким образом, в ЫЬЮВ-модели уравнение неразрывности (1) остаётся прежним, уравнение движения принимает вид
Н2
(Ни)4 + V ■ (Нии) + gV— = gHVh + ^ - фVho, (42)
а уравнение для р преобразуется к следующему виду:
- ÍУ£^Уko) - 6- (Нк ^^+^\\ (43)
где
Р = V ■ (^П + ^^ - -6^ + (V ■ и)2 - (V ■ и) Б^ - 2Vu XV«. (44) \ г у к0г к0
Отметим, что в формуле (44) выражение Бк0 можно заменить на и ■ Vk0, поскольку функция к0(х, у) от времени не зависит, однако мы сохранили прежнее обозначение Бк0, чтобы подчеркнуть единообразие форм записи уравнений всех рассматриваемых здесь НЛД-моделей и возможность использования для них единого подхода к построению численных алгоритмов.
Для одномерных течений ЫЬЮВ-модель включает в себя уравнение неразрывности (24), модифицированное уравнение движения (25)
(Ни)4 + ^Ни2 + = gHkx + рх - фкю,х
и уравнение для дополнительной переменной р
4 т - 6.—+\=р.
VHr/ x r \Hhory,
где
I-, . . Д^о,Л 6R , 2 Dho , 1/6^
F = ( gnx H--— ) T--+ Ux - Ux^-, W = " "Г" + HR + ^xho,x
ho r ho r \ ho
x
d u 27 л i2 7 dho d h
R = —g^xho,x + u ho,xx, r = 4 + hox, ho,x = —;—, ho,xx =
o,xx 2 djx d^b
Отметим, что в случае неподвижного дна, в частности, неподвижного горизонтального дна, NLDB- и NLDD-модели совпадают.
Главными особенностями двух используемых нами слабо дисперсионных моделей являются наличие у них адекватных физике законов изменения энергии [11] и дивергентная форма записи левых частей всех уравнений. Таким образом, упрощённые модели сохраняют важные физические свойства полной NLD-модели, которая в свою очередь наследует эти свойства у трёхмерной модели гидродинамики. Это обстоятельство выгодно отличает используемые в настоящей работе слабо дисперсионные модели от других хорошо известных моделей типа модели Буссинеска (например, [26-28]).
4. Оползень в модельном водоёме с дном параболической формы
Вначале тестирование разработанных алгоритмов выполнялось на задаче о накате волн на берег водохранилища, форма дна которого не зависела от горизонтальной координаты у и задавалась дугой параболы (рис. 1)
г = к0 (х)
х
(кш - к) ( ^ - 1) + к, 0 < х < Ь,
(45)
где кад — глубина в точке х = 0, к — максимальная глубина водоёма. В расчёте использовались следующие значения параметров, определяющих форму дна: кад = -10 м, к = -100 м, £ = 500 м, Ь = 1000 м. С боков (в точках х = 0 и х = Ь) водоём был ограничен вертикальными непроницаемыми стенками с глубиной около них, равной кад.
Оползень представлялся квазидеформируемым телом, движущимся под действием сил тяжести, плавучести, трения о дно и сопротивления воды [16]. Параметрами, определяющими движение первоначально покоящегося оползня, являлись его высота Т и длина Ь (вдоль оси Ох), относительная плотность 7 (отношение плотности материала оползня к плотности воды), коэффициент присоединённой массы , угол трения 0*, коэффициент гидродинамического сопротивления С и начальное заглубление г0. В вычислительных экспериментах один или два из этих параметров варьировались, а другие оставались постоянными и принимали значения
Т =10 м, Ь = 200 м, 7 = 2, = 1, 0* = 5°, С = 1, г0 = -45 м. В расчётах использовалась следующая форма оползня:
2п(х - х0)
Мх)
Т 2 0,
1 + сое
Ь
, |х - х0| < Ь/2, |х - х01 > Ь/2,
(46)
(47)
200 400 600 800 1000
X, м
0 200 400 600 800 1000
X, М
2
а
Рис. 1. Движение оползня по неровному дну параболической формы: графики функции х(Ь) (а) и зависимости локального числа Фруда Ег от координаты х (б) при значениях параметров (46) (1), при Ь = 50 м (2), при Т = 1 м, в* = 10° (3)
где х0 — абсцисса вершины оползня при £ г0 = -45 м получаем
хо
е 1 -
¿0 - Н Н0 - Н
0. Соответственно г0 = Н0 (х0). При
109.13 м.
Функция Н в правой части уравнения движения всех моделей мелкой воды задавалась формулой
г = -Н(х,£) = Н0(х) + Н8[(х + х0 - х(£)), (48)
где х(£) — абсцисса вершины движущегося оползня, х(0) = х0. График функции х = х(£) изображен на рис. 1, а (линия 1). Видно, что при значениях параметров (46) оползень, начав движение из исходного положения, обозначенного на профиле дна кружочком 0, и пройдя точку х = е наибольшей глубины, поднимается по инерции на небольшую высоту противоположного склона и останавливается в положении, отмеченном кружочком 1. При этом движение оползня происходит с докритическими скоростями (см. линию 1 на рис. 1, б, изображающую зависимость локального числа Фруда от координаты х Ег(х) = х/^Но(х)).
Для выяснения влияния дисперсии на картину генерируемых волн были выполнены расчёты на основе модели плоскопараллельных потенциальных течений со свободной границей [29, 30], ЫЬЮ-модели и бездисперсионной модели мелкой воды [24]. В отсутствие экспериментальных данных в качестве "эталонных" использованы результаты, полученные по модели потенциальных течений, в которой учитываются вертикальные перемещения воды и нет ограничений на длину волны.
Из рис. 2, а видно, что для значений параметров (46) модель потенциальных течений и ЫЬЮ-модель дают очень близкие результаты. При уменьшении длины оползня падают скорость его движения, дальность распространения (см. линии 2 на рис. 1), генерируются более короткие поверхностные волны и различие результатов становится значительным (см. рис. 2, б). Таким образом, длина оползня существенно влияет
б
2.0 ть м 1.0
0.0
-1.0
-2.0
-3.0
/ / я
/ \ 1 \ / \ / А' / \ .у чЧ.
!/ \/ г ж
I / V 1 \ Д I 1 -----2 ..............3
1.0 Т],М
0.5 0.0 -0.5 -1.0
•'А \ ч \ ч У/\ Л /' 1 м /ч \ |\ /'\ ■. ! \ 1 Ч\ / Х ----\ / \ V__ /...—-с—"¿X_____V--'
' И \ / т ^
\ Ч / 1 : Ч / '■ Ч • 1
- V 1 1 А 1 / н 1 *,; п 1
*' \ / , , ,,,,,,,, I 1 I I
20
40
60
80 100
20
40 60
80 ^ 100 С
а
Рис. 2. Мареограммы, полученные при Ь = 200 м (а) и Ь = 50 м (б) в рамках КЬБ-модели (1), модели потенциальных течений жидкости (2) и бездисперсионной модели мелкой воды (3) (правый берег)
на адекватность численных результатов на основе КЬЮ-модели. Эта модель получена из уравнений Эйлера в предположении о том [3], что можно пренебречь всеми членами, имеющими порядок малости 0(^4), где ^ — параметр дисперсии, равный отношению характерных значений вертикального и горизонтального масштабов данного процесса. В работе [19] полагалось, что ^ = |г0|/А, где А — длина волны. Поскольку при использовании модели потенциальных течений и КЬЮ-модели возникают системы диспергирующих волн с разными длинами, то определить характерное значение А затруднительно. С другой стороны, общим для всех рассматриваемых моделей является формирование над движущимся оползнем сопровождающей его отрицательной волны ("впадины"), длина А которой сопоставима с протяжённостью оползня в горизонтальном направлении, т. е. с Ь. Если параметр дисперсии мерить по длине этой волны, то в соответствии с рис. 2 можно сделать вывод о том, что при = 0.002563 по сравнению с = 0.6561 КЬЮ-модель даёт лучшие результаты. В отличие от КЬЮ-модели по бездисперсионной модели мелкой воды получается большая скорость распространения волн, которые по истечению некоторого времени формируют одиночную волну, многократно меняющую направление своего распространения при отражениях от противоположных берегов водоёма. Таким образом, сравнение с результатами расчётов по "эталонной" модели показывает, что учёт дисперсии волн приводит к более точному воспроизведению волновой картины, возникающей при движении оползня. Отметим также, что несмотря на отличие волновых картин, воспроизводимых различными математическими моделями, максимальные значения заплесков на берег оказались для них сопоставимыми в довольно широких диапазонах изменения определяющих параметров.
На рис. 3 дано сравнение профилей свободной границы, полученных на основе N1/0- и КЬВБ-моделей. Последняя выведена [11] в предположении, что а = 0(и в NLD-модели можно пренебречь членами порядка 0(а^2), где а — параметр нелинейности. Графики на рис. 3, а соответствуют стандартным значениям параметров (46), на рис. 3, б — удвоенной высоте оползня. Отметим, что для оползня (47) закон его движения при изменении высоты Т не меняется [31], поэтому при разных Т траектории
а
б
Г|, м
-3
-1
-2
2
0
-8
0
200 400 600 800 1000
X, М
0 200 400 600 800 1000
X, м
Рис. 3. Профили свободной границы в момент времени Ь = 60 с, полученные при Т = 10 м (а) и Т = 20 м (б) в рамках КЬБ-модели (1) и КЬББ-модели (2)
а б
Рис. 4. Мареограммы, полученные при Т = 10 м, 9* = 5° (а) и Т = 1 м, 9* = 10° (б) в рамках КЬБ-модели (1) и КЬББ-модели (2) (правый берег)
движения оползня и графики локального числа Фруда изображаются одними и теми же линиями 1 на рис. 1. Однако амплитуда генерируемых волн существенно зависит от параметра Т, возрастая при его увеличении [16]. Из рис. 3, а видно, что для значений параметров (46) ЫЬЮБ-модель даёт в рассматриваемой задаче практически не отличимые от полной ЫЬЮ-модели результаты. При увеличении высоты оползня амплитуда генерируемых волн растёт и различие в рассчитанных по этим моделям профилях свободной границы увеличивается (см. рис. 3, б). Таким образом, для поверхностных волн умеренной амплитуды использование слабо дисперсионной ЫЬЮБ-модели может быть вполне оправданным, поскольку разработанный для этой модели численный алгоритм имеет большую устойчивость в силу понижения степени нелинейности в некоторых членах уравнения (35) по сравнению с уравнением (15) для ЫЬВ-модели.
Результаты на основе ЫЬЮВ-модели для параметров (46) значительно отличаются от полученных по ЫЬЮ-модели (рис. 4, а). Хорошее соответствие наблюдается только при условии слабой деформации дна, что видно по мареограммам на рис. 4, б, полученным при Т =1 ми угле трения 9* = 10°. Для этих параметров оползень имеет малую высоту, движется медленнее (см. линии 3 на рис. 1), проходит небольшое расстояние по склону (точка остановки изображена кружочком 3 на рис. 1, а) и генерирует волны небольшой амплитуды.
Аналогичные выводы относительно областей применимости трёх рассмотренных НЛД-моделей получены при решении задачи исследования величин заплесков на плотину ГЭС. В этом случае дно модельного водоёма имело параболическую форму в месте начального расположения оползня и было горизонтальным в окрестности плотины. На данной задаче вновь подтверждено: хотя классическая модель мелкой воды даёт и иную форму приходящей к плотине волны, предсказываемые ею вертикальные за-плески неплохо согласуются с результатами расчётов по более точным моделям.
Кроме того, проведенные вычислительные эксперименты подтвердили выводы работы [14] о том, что криволинейность склона оказывает существенное влияние на движение оползня и, как следствие, на величины максимальных заплесков волн на берег.
5. Сравнение с экспериментальными данными для плоского откоса
В лабораторных исследованиях, описанных в [18], подводный оползень имитировался движением по плоскому откосу полностью погружённого в воду твёрдого тела, имеющего плотность 1900 кг/м3. В продольном сечении форма тела была полуэллиптической с большой и малой полуосями эллипса равными 25 и 5 см соответственно. Плоский откос сопрягался справа с горизонтальным участком дна.
Сопоставление расчётных и экспериментальных данных было выполнено для разных углов наклона в плоского откоса. Здесь будут описаны некоторые результаты только одного эксперимента, проведённого для случая, когда откос ограничен слева вертикальной непроницаемой стенкой, т. е. рельеф дна описывается формулой
, / ч í hw — x tan в при 0 ^ x ^ £,
ho (x) = < , t , , \ (49)
[ hg при £ ^ x ^ L,
где hw — глубина около левой стенки (в точке x = 0), £ = (hw — hg)/tan в — абсцисса основания склона; при этом использовались следующие значения параметров: hg = —0.9 м, hw = —0.043 м, в = 15°. В этом случае горизонтальный участок начинается в точке £ = 3.198 м. В расчётах правая граница области решения представляла собой непроницаемую стенку и устанавливалась на удалении от начала координат L = 9.7 м, достаточном для того, чтобы отражённые от нее волны не успевали до конца расчёта достичь наклонного участка.
В экспериментах модель оползня, начав движение из точки с абсциссой x0 = 0.40 ми пройдя по откосу расстояние 2.2 м, мгновенно останавливалась в момент времени 2.7 с. В расчётах использовалась иная форма оползня, нежели в экспериментах [18]. Это связано с известным фактом [32], что все НЛД-модели весьма чувствительны к гладкости функции (48), задающей форму подвижного дна, а для полуэллипса функция (48) имеет неограниченные первые производные по x в точках пересечения его поверхности с неподвижным фрагментом дна (49). Поэтому численные эксперименты проводились для оползня (47) с гладкой поверхностью и с теми же высотой и площадью продольного сечения, что в лабораторном эксперименте. Такой же приём применялся, например, в [19] при сравнении численных результатов с экспериментальными данными, причём в указанной работе полуэллипс заменялся нефизичным бесконечно длинным "оползнем" с поверхностью, описываемой бесконечно дифференцируемой функцией. Это позволяло избежать появления высокочастотных колебаний, порождаемых разрывами производных функции (48). Кроме того, в наших численных экспериментах использовалось сглаживание кусочно-линейного дна (49) путём замены его в небольшой окрестности точки x = £ на дугу параболы.
В [18] отмечается, что характерные размеры бассейна допускали пренебрежение трением "оползня" о поверхность склона и сопротивлением воды. Тем не менее в численных расчётах эти эффекты учитывались с целью, чтобы траектории движения оползня по возможности лучше совпадали с экспериментальными. С учётом вышесказанного при численном моделировании использовались следующие значения параметров, определяющих движение оползня согласно уравнению из работы [16]:
T = 0.05 м, b = 0.785 м, y =1.9, Cw = 0.6, в, = 1°, Cd = 1, z0 = —0.15 м. (50)
Отметим, что для плоского откоса уравнение движения оползня [16] принимает более простой вид, чем в случае криволинейного склона, и записывается в виде обыкновенного дифференциального уравнения второго порядка, приведённого в [19].
В экспериментах [18] параметры волн фиксировались волномерами, установленными в четырёх точках. На рис. 5 приведены графики изменения поверхности воды в двух из этих четырёх точек. Если не принимать во внимание фазовый сдвиг, то видно, что полная ЫЬЮ-модель удовлетворительно описывает большее число волн, проходящих через точки замера, чем бездисперсионная модель мелкой воды. Вместе с тем при уходе оползня на глубину точность расчётов по ЫЬЮ-модели снижается.
аб
Рис. 5. Мареограммы в точках х = 0.9 м (а) и х = 1.9 м (б), полученные в рамках КЬБ-моде-ли (1), модели потенциальных течений жидкости (2), модели мелкой воды (3) и в лабораторных экспериментах [18] (4)
аб
с
Рис. 6. Мареограммы в точке х = 1.775 м: а — КЬБ-модель (1), модель потенциальных течений жидкости (2), модель мелкой воды (3), лабораторный эксперимент [19] (4); б — КЬБ-модель (1), КЬББ-модель (2), КЬББ-модель (3)
На рис. 6, а демонстрируется сравнение результатов расчётов с экспериментальными данными работы [19], в которой аналогично [18] исследовались закономерности волнообразования при движении твёрдой модели оползня по плоскому откосу (49) с тем же углом наклона в = 15°, но с иными значениями двух других параметров: Н = -1.05 м, = 0 м. Модель оползня также имела полуэллиптическую форму с полуосями 50 см и 5.2 см. В расчётах использовалась гладкая форма оползня (47), а вместо значений параметров (50) брались следующие:
Т = 0.052 м, Ь =1.57 м, 7 = 1.806, С,ш = 1, в* = 0°, С =1.0, ¿0 = -0.309 м. (51)
Видно, что волна повышения и следующая за ней волна понижения описываются в рамках КЬЮ-модели лучше, чем при использовании бездисперсионных уравнений мелкой воды. Для значений параметров (51) слабо дисперсионная КЬЮБ-модель также даёт удовлетворительные результаты (см. рис. 6, б) в отличие от упрощённой КЬВВ-модели, которая для рассматриваемых значений параметров даёт результаты, уступающие по точности даже полученным по бездисперсионной модели.
Заключение
Для численного моделирования поверхностных волн, возникающих при движении оползня по подводному склону ограниченного водоёма, использованы полная и приближённые нелинейно-дисперсионные модели. Разработан единый для всех рассмотренных моделей подход к построению численных алгоритмов, основанный на аппроксимации расширенной системы уравнений, состоящей из системы уравнений гиперболического типа и уравнения эллиптического типа для негидростатической составляющей давления. Показано, что учёт дисперсии волн существенно влияет на волновую картину. Показано также, что классическая модель мелкой воды даёт адекватные результаты по величинам заплесков на берег. Путём сопоставления численных результатов, полученных в рамках бездисперсионной модели мелкой воды, нелинейно-дисперсионных моделей
б
а
Рис. 7. Генерация поверхностных волн (а) оползнем пространственной формы, скользящим по плоскому откосу (б)
и модели потенциальных течений со свободной границей, обозначены области применимости как полной, так и приближённых НЛД-моделей. Выполнено сопоставление с экспериментальными данными для случая, когда оползень скользит по плоскому откосу и имеет неизменную форму поверхности в одном из горизонтальных направлений. Ввиду ограниченности объёма настоящей публикации результаты численных исследований для оползней пространственно неоднородной формы (рис. 7) и сравнений с экспериментальными данными работы [20] будут изложены в другой статье.
Список литературы
[1] Dütykh D., Días F. Energy of tsunami waves generated by bottom motion // Proc. R. Soc. A. 2009. Vol. 465. P. 725-744.
[2] Fedotova Z.I., Khakímzyanov G.S. Shallow water equations on a movable bottom // Rus. J. Numer. Anal. Math. Modelling. 2009. Vol. 24, No. 1. P. 31-42.
[3] Федотова З.И., Хлкимзянов Г.С. Анализ условий вывода НЛД-уравнений // Вычисл. технологии. 2012. Т. 17, № 5. С. 94-108.
[4] Dütykh D., Mítsotakís D., Beysel S., Shokína N. Dispersive waves generated by an underwater landslide // Numerical Methods for Hyperbolic Equations: Theory and Applications. An Intern. Conf. to Honour Prof. E.F. Toro / Eds. E. Vazques-Cendón, A. Hidalgo, P. García-Navarro. CRC Press, 2013. P. 245-250.
[5] Tappín D.R., Watts P., Gríllí S.T. The Papua New Guinea tsunami of 17 July 1998: Anatomy of a catastrophic event // Nat. Hazards Earth Syst. Sci. 2008. Vol. 8. P. 243-266.
[6] Gríllí S., Ioüalalen M., Asavanant J. et al. Source constraints and model simulation of the december 26, 2004 Indian Ocean tsunami //J. Waterway, Port, Coastal, Ocean Eng. 2007. Vol. 133, No. 6. P. 414-428.
[7] Environmental hazards. The Fluid Dynamics and Geophysics of Extreme Events / Eds. H.K. Moffatt and E. Shuckburgh. Singapore: World Sci., 2010. P. 273-300.
[8] Dütykh D., Katsaoünis T., Mítsotakís D. Finite volume schemes for dispersive wave propagation and runup //J. Comput. Phys. 2011. Vol. 230. P. 3035-3061.
[9] 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.
[10] Chübarov L.B., Fedotova Z.I., Shküropatskii 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.
[11] Dütykh D., Fedotova Z.I., Khakímzyanov G.S. On the equation of energy of some approximate models of wave hydrodynamics // Ibid. 2013. Vol. 28 (в печати).
[12] Chübarov L.B., Eletskij S.V., Fedotova Z.I., Khakímzyanov G.S. Simulation of surface waves generation by an underwater landslide // Ibid. 2005. Vol. 20, No. 5. P. 425-437.
[13] Shokin Yu.I., Fedotova Z.I., Khakímzyanov G.S. et al. Modelling surfaces waves of generated by a moving landslide with allowance for vertical flow structure // Ibid. 2007. Vol. 22, No. 1. P. 63-85.
[14] Beisel S.A., Chübarov L.B., Khakímzyanov G.S. Simulation of surface waves generated by an underwater landslide moving over an uneven slope // Ibid. 2011. Vol. 26, No. 1. P. 17-38.
[15] Beisel S.A., Chubarov L.B., Dutykh D. et al. Simulation of surface waves generated by an underwater landslide in a bounded reservoir // Ibid. 2012. Vol. 27, No. 6. P. 539-558.
[16] Хлкимзянов Г.С., Шокинл Н.Ю. Оценка высот волн, вызванных подводным оползнем в ограниченном водоёме // ПМТФ. 2012. Т. 53, № 5. С. 67-78.
[17] Численное моделирование течений жидкости с поверхностными волнами / Г.С. Хакимзянов, Ю.И. Шокин, В.Б. Барахнин, Н.Ю. Шокина. Новосибирск: Изд-во СО РАН, 2001.
[18] Елецкий С.В., Майоров Ю.Б., Максимов В.В. и др. Моделирование генерации поверхностных волн перемещением фрагмента дна по береговому склону // Совместный вып. журналов "Вычисл. технологии" и "Вестник КазНУ им. аль-Фараби". 2004. Т. 9, ч. II. С. 194-206.
[19] Grilli S.T., Watts P. Tsunami generation by submarine mass failure. I: Modeling, experimental validation, and sensitivity analyses //J. Waterway, Port, Coastal, Ocean Eng. 2005. Vol. 131, No 6. P. 283-297.
[20] Enet F., Grilli S.T. Experimental study of tsunami generation by three-dimensional rigid underwater landslides // Ibid. Vol. 133, No 6. P. 442-454.
[21] Гусев О.И. Об алгоритме расчёта поверхностных волн в рамках нелинейно-дисперсионной модели на подвижном дне // Вычисл. технологии. 2012. Т. 17, № 5. С. 46-64.
[22] ХАкимзянов Г.С., ШокинА Н.Ю. Некоторые замечания о схемах, сохраняющих монотонность численного решения // Там же. 2012. Т. 17, № 2. С. 78-98.
[23] Соммер А.Ф., ШокинА Н.Ю. О некоторых проблемах конструирования разностных схем на двумерных подвижных сетках // Там же. 2012. Т. 17, № 4. С. 88-108.
[24] ХАкимзянов Г.С., ШокинА Н.Ю. Метод адаптивных сеток для одномерных уравнений мелкой воды // Там же. 2013. Т. 18, № 3. С. 54-79.
[25] ХАкимзянов Г.С., ШокинА Н.Ю. Определение зоны затопления при накате длинных волн на берег // Тр. шестого Совещания российско-казахстанской рабочей группы по вычислительным и информационным технологиям. Алматы: Казак университета, 2009. С. 305-314.
[26] Peregrine D.H. Long waves on a beach // J. Fluid Mech. 1967. Vol. 27, pt. 4. P. 815-827.
[27] Дорфман А.А., Яговдик Г.И. Уравнения приближённой нелинейно-дисперсионной теории длинных гравитационных волн, возбуждаемых перемещениями дна и распространяющихся в бассейне переменной глубины // Численные методы механики сплошной среды. 1977. Т. 8, № 1. C. 36-48.
[28] Dias F., Milewski P. On the fully-nonlinear shallow-water generalized Serre equations // Phys. Lett. A. 2010. Vol. 374(8). P. 1049-1053.
[29] ХАжоян М.Г., ХАкимзянов Г.С. Численное моделирование взаимодействия поверхностных волн с подводными препятствиями // Вычисл. технологии. 2003. Т. 8, № 4. С. 108-123.
[30] ХАжоян М.Г. Численное моделирование поверхностных волн над подвижным дном // Там же. 2007. Т. 12, № 4. С. 96-105.
[31] ХАкимзянов Г.С., ШокинА Н.Ю. Численное моделирование поверхностных волн, возникающих при движении подводного оползня по неровному дну // Там же. 2010. Т. 15, № 1. С. 105-119.
[32] Федотова З.И., Чубаров Л.Б., Шокин Ю.И. Моделирование поверхностных волн, порожденных оползнями // Там же. 2004. Т. 9, № 6. С. 89-96.
Поступила в 'редакцию 26 августа 2013 г.
