CC BY
32
Вычислительные технологии
Том 20, № 5, 2015
0 численных методах решения задач о накате волн на берег. II. Опыт решения модельных задач
Ю.И. Шокин1, А. Д. Рычков1, Г. С. Хлкимзянов1'2, Л. Б. Чубаров1'2'*
1 Институт вычислительных технологий СО РАН, Новосибирск, Россия
2 Новосибирский государственный университет, Россия *Контактный e-mail: chubarovster@gmail.com
Настоящая статья продолжает цикл работ по численному моделированию процессов заплеска поверхностных волн на берег. В ней изложены результаты решения ряда модельных задач, полученные в рамках теории мелкой воды для длинных волн различных конфигураций, амплитуд и углов наклона склонов. Продемонстрированы возможности применения вычислительного алгоритма, основанного на комбинированном методе TVD + SPH, для моделирования как слабо нелинейных (необрушающихся), так и существенно нелинейных (обрушающихся) волн. Полученные результаты сравниваются с результатами других авторов. Обсуждаются возможные причины сходства и различия сопоставляемых результатов.
Ключевые слова: поверхностные волны, мелкая вода, накат на берег, численное моделирование.
Введение
Одно из наиболее интересных и сложных явлений, изучаемых в рамках математических моделей волновой гидродинамики, — взаимодействие волн с берегом и прибрежными сооружениями. Именно к этому классу принадлежат задачи, связанные с заключительным, подчас носящим разрушительный характер этапом существования волн цунами — их заплеском на берег, проявляющимся в виде наката и отката волн. Чрезвычайно важным для решения прикладных задач прогнозирования ущерба от цунами и смягчения катастрофических последствий для населения прибрежных зон, а также производственных объектов является умение с необходимой точностью и заблаговременностью определять экстремальные величины заплесков, дальностей проникновения волн на берег и протяженностей зон осушения, глубин заливания и продолжительности пребывания участков суши под водой. Безусловно, эти характеристики существенно зависят от распределений глубин акваторий, прилегающих к защищаемым участкам суши, от рельефов этих участков, их насыщенности защитными сооружениями, растительностью, другими структурами естественного или рукотворного происхождения, наконец, от степени изрезанности береговой линии, направления подхода волны и влияния на нее других природных явлений — приливов, ветров, штормов и т. п.
Конечно же, задачи о накате волн на берег с реальным рельефом дна и суши должны рассматриваться в пространственной постановке, а в определенных условиях — с учетом эффектов обрушения волн, турбулентности, образования двухкомпонентной водо-
© ИВТ СО РАН, 2015
воздушной смеси и т. п., однако для рассмотрения физических основ процессов взаимодействия волн с сушей вполне допустимо использовать приближенные гидродинамические модели мелкой воды в одномерной постановке.
Создание надежных алгоритмов для решения задач такого класса с учетом упомянутых выше разнообразных сопутствующих факторов требует не только теоретического анализа математических моделей и применяемых для их реализации численных методов, но и проверки их на тестовых задачах. Одной из "канонических" задач такого типа, по мнению авторов работы [1], является задача определения параметров заплеска длинных волн на плоский откос с заданным углом наклона, накатывающихся на него со стороны участка с дном постоянной глубины. Принято считать, что такая задача с достаточной полнотой отражает основные физические аспекты рассматриваемого явления. Содержание настоящей статьи можно определить как попытку продемонстрировать возможности рассмотренного ранее [2] комбинированного метода TVD + SPH для решения упомянутой выше "канонической" задачи в условиях, воспроизводящих заплес-ки длинных волн различных конфигураций и амплитуд на склоны различной крутизны, выявить базовые особенности исследуемых явлений как для слабо нелинейных (необру-шающихся), так и для существенно нелинейных (обрушающихся) волн, а также оценить устойчивость (робастность) алгоритма в широком диапазоне изменения параметров такой задачи.
Статья состоит из двух частей, в первой из них кратко изложена постановка задачи, а во второй — представлены результаты численного решения двух типов задач. Первая задача связана с воспроизведением волновых процессов, экспериментально исследованных в Большом волновом канале Ганноверского университета [3], а вторая — с моделированием заплеска на откос так называемых Ж-волн. Полученные результаты решения первой задачи сопоставляются с результатами численного моделирования, выполненного группой исследователей под руководством Е.Н. Пелиновского [4]. Результаты численного решения второй задачи сравнивались с результатами аналитических исследований аналогичных процессов [1]. Несмотря на очевидную "модельность" рассмотренных N-волн, следует заметить, что они имеют самое непосредственное отношение к волнам цунами, порожденным подводными землетрясениями, очаги которых расположены на небольшом расстоянии от берега. В таких случаях протяженность трассы распространения волны может оказаться недостаточной для того, чтобы начальное возмущение, направленное в сторону берега, трансформировалось в волну повышения солитонного типа, и к берегу подходит N-волна зачастую с лидирующим фрагментом понижения (денивиляции) уровня. Полученные результаты численного моделирования таких волн качественно согласуются с выводами работы [1].
1. Постановка задачи и вычислительные алгоритмы
Как и в предыдущей статье авторов [2], в рамках модели мелкой воды рассматривается одномерная задача о накате волны на плоский откос, сопряженный с горизонтальным дном глубины h0. Используется декартова система координат Oxz с вертикальной осью Oz, направленной вверх, и с координатной линией z = 0, совпадающей с невозмущенной свободной поверхностью слоя идеальной несжимаемой жидкости, ограниченного снизу неподвижным дном z = — h (х), а сверху — подвижной свободной границей z = ^ (x,t), где t — время. Система уравнений мелкой воды в этом случае записывается в виде законов сохранения
Я, + ¥х = О, 1> 0, (1)
где
)' Г ={ни-2 2/2) , О = (9Н}гх) ,
д — ускорение силы тяжести. Искомыми величинами являются полная глубина слоя жидкости Н(х, = г](х, + 7(х) > 0 и и(х, ¿) — скорость, осредненная по вертикали от дна до свободной поверхности.
В статье [2] рассматривался канал, в котором плоский откос располагался слева от участка с дном постоянной глубины. Эта же схема используется и в настоящей работе. Однако в некоторых численных экспериментах для более удобного сопоставления полученных результатов с результатами других авторов применяется модифицированная форма расчетной области, в которой откос располагается справа. В таком случае область решения ограничивается слева непроницаемой стенкой, расположенной в точке х = 0, а справа — подвижной границей хо (¿), отделяющей воду от суши (точка уреза). Эта граница заранее неизвестна, поэтому количество искомых функций увеличивается.
Задача для уравнения (1) замыкается постановкой начальных и краевых условий, которые в предположении, что начальное положение х0 (0) точки уреза известно и волна движется слева направо, принимают следующий вид:
н (х0(г),г) = 0, и (0,г) = 0, 0; (2)
г] (х, 0) = г]0 (х), и (х, 0) = и0 (х) , 0 < х < х0 (0). (3)
Рельеф дна и прилегающей суши задается функцией
!
= 7 ( ) = I — 70 при 0 < х < х3
( ) ' —70 + (х — х3) tg[ при х8 <х < Ьх,
где [ > 0 — угол наклона откоса; х3 > 0 — заданная абсцисса основания склона (кромки откоса); Ьх = х3 + (г0 + 70) ctg[; г0 > 0 — высота суши в точке х = Ьх. Значение г0 подбирается так, чтобы максимальный вертикальный заплеск волны на склон (ордината точки уреза на склоне) не превышал х0. Таким образом, при всех Ь > 0 будет выполняться неравенство х0 (¿) < Ьх.
В исследованиях, результаты которых представлены в настоящей статье, авторы для моделирования наката волн на склон использовали два численных метода второго порядка точности. Первый из них основан на методе, предложенном в работе [5], и представляет собой оригинальную комбинацию метода сглаженных частиц (БРИ) и конечно-разностной схемы, обладающей ТУБ-свойствами. Этот метод (ниже он будет называться ТУБ + БРИ) авторами статьи был дополнен алгоритмом расчета положения подвижной точки уреза. Результаты сравнительного анализа методов моделирования наката длинных поверхностных волн на берег в рамках теории мелкой воды [2] показали, что наиболее перспективным для решения задач этого класса является именно метод ТУБ + БРИ.
Второй метод, условно называемый "точным" [6], основан на использовании подвижной разностной сетки и применении аналитического решения задачи в малой окрестности точки уреза для расчета ее положения и скорости движения.
2. Некоторые результаты расчетов
Приведем результаты решения двух модельных задач, первая из которых состоит в определении характеристик заплеска на откос одиночных волн повышения и понижения (положительной и отрицательной полярности) в широком диапазоне начальных амплитуд, а вторая — заплеска их комбинации, так называемых Ж-волн, в которых лидирующей может быть как волна понижения, так и волна повышения. Результаты решения каждой из рассмотренных задач сопоставляются с результатами других авторов. В первом случае это качественное и количественное сравнение с результатами группы профессора Е.Н. Пелиновского [4], а во втором — качественное сравнение с аналитическими результатами S. Tadepalli и C.E. Synolakis [1]. Для того чтобы облегчить сопоставление с результатами работы [4], используются отмеченное выше "правое" расположение склона и начальное распределение скоростей, обеспечивающее распространение волн по направлению к склону. Вторая задача решается в расчетной области с "левым" расположением склона.
2.1. Моделирование наката одиночных волн на плоский откос
Ниже излагаются результаты численного моделирования наката одиночных волн на склон и сопоставления их с результатами группы Е.Н. Пелиновского [4, 7, 8] включающими аналитические решения (для необрушающихся волн), материалы лабораторных и численных исследований. Эксперименты выполнялись в Большом волновом канале Ганноверского университета с использованием волнового лотка, обустроенного в этом канале и состоящего из участка постоянной глубины h0 = 3.5 м длиной xs = 250 м, сопряженного с откосом 1:6, который располагался у правой границы лотка [3]. Численные эксперименты авторами работы [4] выполнялись с помощью пакета программ CLAWPACK [9].
Начальные данные задавались соотношениями
ц (х, 0) = А cosh 2 (х — xw)/L , u (х, 0) = 2 \Jд (h0 + ц (х, 0)) — л/gh
(5)
Как показано в [10-12], они соответствуют точному решению (волна Римана) уравнений мелкой воды
Н (х,Ь) = Но [х - V (Н) Ь] , У (Н) = 3^дН - Н0 (х) = 1г0 + у (х, 0).
Авторы работы [4] интерпретируют величину Ь как половину длины начального возвышения, полагая в своих расчетах Ь =11 м. Экстремальные значения смещений свободной поверхности (амплитуды волн) в начальный момент находились над точкой с координатой х,ш = 50 м, для волн положительной полярности они изменялись в диапазоне от 0.05 до 3.5 м, а для отрицательной — от -0.05 до -3.49 м.
По сравнению с расчетами, выполненными авторами работы [4], при подготовке настоящей статьи дополнительно моделировались волны с меньшими начальными амплитудами — 0.001, 0.0025, 0.005, 0.01 м. Следует заметить, что пределы применимости теории мелкой воды, определяемые возможностью использования аналитической формулы Синолакиса [13], в условиях проведенных авторами расчетов ограничиваются интервалом изменения начальных амплитуд 0.008064 < А < 0.228976576 м. Однако результаты, полученные для амплитуд, выходящих за пределы указанного интервала, могут оказаться полезными для проверки новых алгоритмов моделирования заплес-ков в рамках модели мелкой воды. По крайней мере, эти результаты явно указывают
на то, что используемые численные методы могут работать даже за теоретическими границами применимости аппроксимируемой ими математической модели.
Выбирая форму представления результатов, авторы следуют образцам, предложенным в [4]. К сожалению, некоторые особенности постановки задач численного моделирования в работе [4] не позволяют выполнить непосредственное количественное сопоставление результатов, однако их качественный анализ вполне возможен.
Первая серия графиков (рис. 1) представляет профили свободной поверхности, рассчитанные для того же, что ив [4], набора амплитуд (рис. 2, 3, 5 и 6 из [4]). На каждом из графиков изображены волны в начальный момент времени (кривые 1), вблизи
б
2
0.5-
0.4-
0.3 ■
0.2-
0.1 •
3
7\ 4 Л 2 1
IX У \ ¿А У
—1
50
100
150 200
х, м
в
50
100
150 200
X, м
250
300
V
> \ У / /
) Л у4 \ к/
250
300
1.41.21
0.80.60.40.20-
3
1
/
1 \ А 2/ /
/ \ Л / /
/ У \ 4 / /
—^
I 1 1
50
100
150 200
X, м
250
300
Рис. 1. Профили свободной поверхности при накате на плоский откос волн положительной полярности с начальными амплитудами А = 0.1 м (а), 0.5 м (б), 1.5 м (в), 3.5 м (г) в начальный момент времени, в момент подхода волны к кромке склона, в момент почти максимального вертикального заплеска и, наконец, в стадии подхода отраженной откосом волны к противоположной от него границе расчетной области. Значения соответствующих моментов времени указаны в тексте
а
г
кромки склона (кривые 2, для амплитуд 0.1 и 0.5 м они соответствуют моменту времени £ = 30 с, для А = 1.5 и А = 3.5 м — £ = 20 с); в момент почти максимального вертикального заплеска (кривые 3, соответствуют для А = 0.1 и А = 0.5 м времени £ = 40 с, для А = 1.5 м — £ = 35 си для А = 3.5 м — £ = 30 с); наконец, в стадии подхода отраженной откосом волны к противоположной от него границе расчетной области (кривые 4, соответствуют для А = 0.1 и А = 0.5 м — £ = 70 с, для А =1.5 м — £ = 60 с, а для А = 3.5 м — £ = 50 с).
Надо заметить, что при малых амплитудах наблюдается не только качественное, но и количественное согласие с результатами из [4], в то время как с ростом начальных амплитуд различия в профилях свободной поверхности, сформировавшихся после взаимодействия волн со склоном, существенно возрастают.
Влияние нелинейности сказывается уже для волны самой малой амплитуды, передний фронт которой постепенно становится круче, но обрушения (градиентной катастрофы) здесь не наступает (рис. 1, а), амплитуда волны, проходящей над участком дна постоянной глубины, практически не снижается. После отражения от склона волна существенно изменяет свою форму, превращаясь, по сути, в ^-волну, содержащую фрагмент понижения уровня ниже нулевого.
С ростом начальных амплитуд волн воздействие нелинейных эффектов на их характеристики усиливается. Под их влиянием, проходя над участком канала постоянной глубины, передние фронты волн становятся существенно круче, в то время как их задние фронты выполаживаются так, что увеличиваются длины волн и снижаются их амплитуды. В результате к кромке откоса подходят и накатываются на откос уже "обрушенные" волны. Эти соображения подтверждаются также рис. 2, представляющим характеристики падающих и отраженных волн, рассчитанных над участком лотка постоянной глубины (рис. 2, а) и над кромкой откоса (рис. 2, б). Из рис. 2, а видно, что падающие и отраженные волны разделены во времени, причем падающие волны имеют вертикальный передний фронт, а отраженные от откоса волны — гладкую форму.
а
Рис. 2. Мареограммы, рассчитанные над участком виртуального лотка постоянной глубины в точке х = 150 м (а) и над кромкой откоса в точке х = 250 м (б) при заплеске на плоский откос волн с начальными амплитудами А = 0.1 м (I), 0.5 м (2), 1.5 м (3), 3.5 м (4)
Над кромкой откоса разделения волн во времени не наблюдается, и здесь падающая волна, не успев полностью пройти над точкой х = х8, встречает волну, отражающуюся от склона.
Сопоставление результатов, представленных на рис. 3, а, с результатами из [4] (рис. 8) демонстрирует не только количественные, но и существенные качественные различия. Так, рассчитанные с помощью метода ТУБ + ЯРЫ высоты заплесков оказываются значительно больше, в то время как продолжительность наката — существенно короче. Кроме того, в отличие от результатов работы [4], зависимости величин максимального (по абсолютной величине) отката от начальных амплитуд являются монотонными, а зависимости вертикальных заплесков для всех амплитуд — гладкими. Наконец, если в [4] с ростом начальной амплитуды абсолютные значения величин отката изменяются немонотонно и уменьшаются, устремляясь к нулю (что плохо согласуется с представлениями о физике моделируемых процессов) для самых больших начальных амплитуд, начиная с А = 2 м, то в представленных на рис. 3, а результатах эти значения пусть медленно, но монотонно увеличиваются. Имеют место и различия в скоростях, выражающиеся в том, что процессы наката и отката для одних и тех же условий в [4] происходят существенно медленнее.
Для возможного использования в качестве тестовых материалов некоторые данные, полученные при решении рассматриваемой задачи, представлены в виде таблицы.
Следующие результаты получены в ходе моделирования наката на откос волн отрицательной полярности (А < 0) — волн понижения уровня, начальные положения которых совпадали с начальными положениями волн положительной полярности. Совпадали также с точностью до зеркального отражения и формы начальных волн. Как и в случае волн повышения, результаты, полученные методом ТУБ + ЯРЫ, близки к результатам из [4] при малых амплитудах, но существенно различаются для | А| > 1 м. Уже при первичном анализе рассчитанных профилей легко заметить, что для таких амплитуд времена прихода экстремальных значений уровня в одни и те же точки не совпада-
Рис. 3. Динамика величин вертикального заплеска на склон волн положительной (а, А > 0) и отрицательной (б, А < 0) полярности с начальными амплитудами А = ±0.1, ±0.3, ±0.7, ±1.5, ±2.5 м (кривые 1-5 соответственно), А = 3.5 м (а) и А = -3.49 м (б) — (6)
Основные характеристики заплеска волн положительной полярности на плоский откос, рассчитанные двумя вычислительными алгоритмами
А, м Отношение максимального наката к начальной амплитуде Отношение максимального отката
к начальной амплитуде к максимальному накату
ТУБ + БРИ [4] ТУБ + БРИ [4] ТУБ + БРИ [4]
0.05 4.868 4.7354 1.104 1.0053 0.2268 0.2365
0.1 5.499 4.6779 1.024 1.0015 0.1862 0.2211
0.3 5.658 3.5827 0.682 0.7363 0.1205 0.2055
0.5 5.0096 2.9408 0.5778 0.5733 0.1153 0.1962
0.7 4.503 2.5674 0.4829 0.4672 0.1072 0.1820
1 3.9676 2.2246 0.396 0.3653 0.0998 0.1642
1.5 3.4581 1.9056 0.2983 0.1877 0.0863 0.0983
2 3.1606 1.7105 0.2547 0.0899 0.0806 0.0526
2.5 2.9644 1.5826 0.2144 0.0229 0.0723 0.0144
3 2.8293 1.4976 0.189 0.0065 0.0668 0.0044
3.5 2.7244 1.3807 0.1663 0.0009 0.0611 0.0007
ют, поэтому представленные на рис. 4 профили соответствуют отличным от выбранных в [4] временам (рис. 10-12, 14 из [4]). На каждом из графиков волны изображены в начальный момент времени (кривые 1), вблизи кромки склона (кривые 2, соответствуют моменту времени £ = 20 с); в момент почти наибольших вертикальных отката (кривые 3, соответствующие для А = -0.1 м моменту времени £ = 41 с, для А = -0.5 м — £ = 43 с, для А = —1 м — £ = 46.5 с, для А = -3.49 м — £ = 53.86 с) и наката (кривые 4, соответствующие для А = —0.1 м моменту времени £ = 44 с, для А = —0.5 м £ = 48 с, для А = —1 м £ = 51.2 с, дляА = —3.49 м £ = 58.8 с); наконец, в стадии подхода отраженной откосом волны к противоположной от него границе расчетной области (кривые 5, соответствующие для А = —0.1 м, —0.5 ми —1 м моменту времени £ = 75 с, а для А = —3.49 м — £ = 85 с).
Так, при А = 1.0 м в расположении профилей наката и отката отмечается некоторое опережение по сравнению с результатами работы [4], а при А = —3.49 м — отставание от них. Численный эксперимент с такой амплитудой волны отрицательной полярности трудно считать осмысленным из-за очевидной неприменимости в этом случае теории мелкой воды, хотя возможность проведения такого рода вычислений свидетельствует о высокой работоспособности использованного алгоритма. На рис. 4, г видны существенные отличия от близкого по содержанию рисунка из работы [4] для всех моментов времени, проявляющиеся в гораздо более сложных формах падающей и отраженной волн.
Общее представление о динамике уреза при накате на плоский откос волн отрицательной полярности дают результаты (рис. 3, б и рис. 15 из [4]), полученные методом ТУО + БРИ на сетке с числом узлов, равным 6000. Они, впрочем, оказались практически идентичными полученным на сетке с числом узлов 4000. Графики на этом рисунке демонстрируют монотонный рост величины отката, как первичного, так и вторичного, однако величина наката, возрастая вплоть до начальной амплитуды 2.5 м, затем уменьшается для амплитуд А = —3 и —3.499 м.
Далее с результатами работы [4], представленными в обобщенной форме, сопоставляются результаты, полученные авторами с помощью метода ТУБ + БРИ и алгоритма, основанного на определении положения точки уреза с использованием аналитических
Рис. 4. Профили свободной поверхности при накате на плоский откос волн отрицательной полярности с начальными амплитудами А = -0.1 м (а), -0.5 м (б), -1 м (в), -3.49 м (г) в начальный момент времени, в момент подхода волны к кромке склона, в моменты почти максимальных вертикальных отката и наката, наконец, в стадии подхода отраженной откосом волны к противоположной от него границе расчетной области. Значения соответствующих моментов времени указаны в тексте
а
в
2
решений уравнений мелкой воды при постановке разностных краевых условий на подвижной линии уреза [6]. Значения величин из работы [4] были получены частично из ее текста, а частично оцифровкой приведенных в ней графиков.
На рис. 5 кривые 1-3 иллюстрируют зависимость в стадии наката волн положительной полярности максимальных значений величин К/А от начальной амплитуды. Эти величины, рассчитанные с помощью алгоритмов метода ТУБ + ЯРЫ, существенно выше результатов из работы [4] практически для всех рассмотренных амплитуд, но близки к результатам, полученным с помощью методики работы [6], которые, однако, оказыва-
О 0.5 1 1.5 2 2.5 3 3.5
А, М
Рис. 5. Относительные характеристики заплеска волн положительной полярности в зависимости от их начальных амплитуд; максимальные значения величин Я/А в стадии наката (1-3) и в стадии отката (4-6) рассчитаны с помощью метода ТУБ + БРИ (1, 4), методов [6] (2, 5) и [4] (3, 6); серым цветом выделена область применимости аналитического решения [13]
ются несколько меньше при самых больших амплитудах. Надо отметить, что если для абсолютных величин наката Я имеет место монотонный рост, то соответствующие относительные величины Я/А сначала (вплоть до амплитуды 0.25 м) возрастают и только затем начинают монотонно снижаться, что может быть объяснено усилением воздействия нелинейных эффектов, которые приводят к росту крутизны переднего фронта, уменьшению высоты волн, подходящих к откосу, и к увеличению их длины. В результатах из [4] экстремум не обнаруживается из-за отсутствия данных для малых амплитуд. Существенно меньшие значения из работы [4] можно объяснить тем, что в используемом ее авторами компьютерном пакете CLAWPACK реализован неявный алгоритм, абсолютная устойчивость которого позволяет вести расчет с довольно большим шагом по времени, приводящим к существенному росту численной диссипации, что и проявляется в некотором снижении качества обсуждаемых результатов.
Рассматривая характеристики заплесков в стадии отката, можно отметить некоторое изменение ситуации. Здесь для малых амплитуд (почти до А =1 м) значения, полученные с помощью алгоритма ТУБ + БРИ, монотонно возрастают и практически совпадают с результатами из [4], которые затем начинают резко уменьшаться, устремляясь к нулю при самой большой амплитуде. Результаты, рассчитанные алгоритмом ТУБ + БРИ, напротив, продолжают устойчивый монотонный рост, оставаясь во всем диапазоне изменения амплитуд меньше результатов, полученных по методике из работы [6]. Анализ графиков, представляющих максимальные значения величин в стадии отката, прежде всего указывает на сохраняющееся аномальное стремление результатов из [4] к нулю для А > 1 м и на их близость к результатам, полученным с помощью алгоритма ТУБ + БРИ, вплоть до амплитуды А = 1 м. Поведение последних качественно близко к распределению, рассчитанному по методике [6], которое, однако, как и ранее, расположено значительно выше. Результаты, рассчитанные для малых начальных амплитуд, как и в стадии наката, демонстрируют наличие экстремума,
несколько смещенного в сторону самых малых значений А. В работе [4] необходимые для сравнения результаты, рассчитанные в обсуждаемом диапазоне изменения аргумента, к сожалению, отсутствуют.
Отношения максимальных значений |Д| в стадии отката к максимальным значениям Я, аналогичным величинам в стадии наката, уменьшаются с ростом начальной амплитуды подходящей к откосу волны. По-прежнему результаты из [4] устремляются к нулю, расчетные данные, полученные по методике [6], принимают наибольшие значения, а результаты, рассчитанные алгоритмом ТУБ + ЯРЫ, демонстрируют наличие экстремума в области очень малых аргументов.
Аналоги рассмотренных выше характеристик заплеска для волн отрицательной полярности показывают, что, как и в случае волн повышения, результаты, полученные с помощью алгоритма ТУБ + ЯРЫ, превышают результаты работы [4], что в большей степени проявляется в характеристиках наката, в то время как соответствующие характеристики отката оказываются очень близкими. Оба набора относительных величин для стадии наката (рис. 6) обладают максимумами в диапазоне малых начальных амплитуд, после достижения которого соответствующие распределения монотонно уменьшаются. Относительные характеристики стадии отката экстремумов не имеют и монотонно уменьшаются с ростом начальной амплитуды.
Нужно отметить некоторую особенность изменения величины Я, рассчитанной с помощью алгоритма ТУБ + ЯРЫ, состоящую в том, что вблизи начальной амплитуды А = 2.5 м монотонный рост максимального вертикального положения подвижной точки уреза сменяется столь же монотонным его понижением. Можно предположить, что таким образом проявляется специфика воспроизведения применяемым алгоритмом нелинейных эффектов, приводящих к раннему обрушению подходящих к откосу волн большой амплитуды с одновременным снижением высоты их положительных фрагментов. Этот вывод находит свое подтверждение как на профилях соответствующих свобод-
Рис. 6. Относительные характеристики заплеска волн отрицательной полярности в зависимости от их начальных амплитуд; максимальные значения величин К/А в стадии наката (I, 2) и в стадии отката (3, 4) рассчитаны с помощью метода ХУБ + SPH (I, 3) и метода [4] (2, 4); серым цветом выделена область применимости аналитического решения [13]
ных поверхностей (см. рис. 4), так и на графиках динамики подвижной точки уреза (см. рис. 3, б). Распределения отношений максимальных значений |Д| в фазе отката к максимальным значениям |Д| в фазе наката для обоих наборов данных качественно весьма схожи. Они сначала довольно резко уменьшаются и принимают минимальное значение для начальной амплитуды А = —1 м, а затем медленно возрастают. При этом распределения из работы [4] оказываются выше.
Выполненный анализ особенностей заплесков волн различной полярности (волн повышения и понижения) показал, что волны положительной полярности приводят к большим величинам вертикального заплеска К, однако скорости перемещения точки уреза при накате волн отрицательной полярности оказываются значительно больше и это приводит к более крутым задним фронтам волн, что качественно совпадает с выводами авторов работы [7]. Заметим также, что при накате волн положительной полярности вертикальные смещения точек уреза в фазе наката существенно превышают смещения в фазе осушения, в то время как в случаях волн отрицательной полярности эти величины оказываются гораздо ближе. Немонотонность форм свободной поверхности над точкой излома рельефа дна (см. рис. 2, б), отсутствующая в результатах, представленных в [7], является следствием взаимодействия подходящей и отраженной от склона волн.
Подводя итог сопоставления с результатами, изложенными в работе [4] и в определенной степени связанными с условиями физических экспериментов в ганноверском лотке, надо еще раз отметить, что даже использование известного компьютерного кода, созданного высокопрофессиональными специалистами, требует детального исследования реализуемого этим кодом алгоритма, понимания его особенностей и правильного задания определяющих параметров. В противном случае можно попасть в ситуацию, когда присущие алгоритму численные эффекты серьезно исказят значения вычисляемых величин, затруднят адекватную интерпретацию полученных результатов и приведут к далеким от реального физического смысла выводам.
Следует также отметить, что используемый в методе ТУБ + БРИ явный алгоритм обеспечивает необходимую точность расчета и подтверждает правильность выбора численных схем подобного класса для моделирования тонких эффектов, возникающих при заплеске волн даже на простые модельные склоны.
2.2. Моделирование наката на плоский откос N-волн различной конфигурации
В завершающей серии расчетов моделировался заплеск на склон Ж-волн, в начальный момент состоящих из расположенных рядом двух косинусоидальных фрагментов положительной и отрицательной полярности (возвышения и понижения или понижения и возвышения), которые задавались соотношениями
0,
Х^ + хп ^
0 —— ^с ^^ ^с^, ^с^^ ^с —— ^ ^,
гц(х, 0) = <
А ¥
А ¥
(
X
сое 2п
V
сое / 2п-
2
хк+х
т
х —
2
+ 1
+ 1
— % — ,
— % —
(6)
и(х, 0)
0, 0 < х < хп, хт < х < Lx,
¡—-— гц(х, 0)
9(ho + А) ——:—;—^т, Хп < х < хк,
ho + <q(x, 0)
т-IT v(x, 0) ^ ^
-v 9(ho - А)-—:—7—ttv, %к < х < Хт;
ho + rj(x, 0)
величина А в (6) и (7) равна 0.2 м.
Откос располагался у левой границы расчетной области, а волна распространялась справа налево. Рельеф дна и прилегающей суши задавался функцией
7 / ч Г — в при 0 < х < х3, /гЛ
г = — п (х) = < ю\
-h0 при xs < х < Lx,
где, как и в (4), ß > 0 — угол наклона откоса, а z0 > 0 — высота суши в точке х = 0, xs = (z0 + h0) ctg ß — абсцисса основания склона, h0 = 1 м, z0 = 2 м, хт — хк = хк — хп = 10 м, %п — %s = 15 м. Для выбранного значения z0 неравенство х0 (t) > 0 будет выполняться при всех t > 0.
Интерес к исследованию характеристик заплеска на берег такого типа волн объясняется, в частности, тем, что в задачах, связанных с изучением трансформации сейсмоген-ных цунами в прибрежной зоне, используются, как правило, модели очагов подводных землетрясений [14, 15], порождающих цунами и приводящих в широком диапазоне своих параметров к формированию смещений поверхности воды, имеющих в поперечных сечениях форму Ж-волн. В случаях, когда эти очаги расположены на небольшом расстоянии от берега, протяженность трассы распространения может оказаться недостаточной для того, чтобы начальное возмущение, унаследовавшее от сейсмического очага фрагмент понижения уровня, как правило, направленный в сторону берега, трансформировалось в волну повышения солитонного типа. Тщательное аналитическое исследование задач о заплеске N-волн на плоский откос, сопряженный с участком дна постоянной глубины, представлено в работе [1], авторы которой рассматривали отличающиеся от использованных в настоящей статье начальные формы (6) волн ("обобщенные", "равнобедренные" и "двойные" Ж-волны ), облегчающие получение точных аналитических решений.
Анализ графиков (рис. 7, а) показывает, что при накате на откос с углом наклона ß =15° различных типов N-волн, в том случае, когда лидирующим является фрагмент отрицательной полярности, величина вертикального заплеска оказывается больше, но амплитуды последующих колебаний точки уреза — меньше, чем у волны положительной полярности.
Аналогичные результаты для угла наклона склона ß = 5° (рис. 7, б) демонстрируют изменение общих характеристик заплеска волны для более пологого откоса. В случае, если лидирующей является волна отрицательной полярности, следующий за ней фрагмент возвышения успевает за счет своей более высокой скорости ее догнать и в результате взаимодействия фрагментов различных полярностей эффект осушения при первом заплеске на склон практически исчезает. Этот результат подтверждается и сравнением рассчитанных величин заплеска N-волн при четырех значениях угла наклона склона, когда лидирующими являются волны отрицательной полярности (рис. 8). Соответствующие графики показывают, что при уменьшении угла наклона разница во времени между приходом на склон волны отрицательной полярности и следующей за ней волны повышения сокращается.
-0.5 -|-------------- -0.25 -|--------------
0 10 20 30 40 50 60 70 0 10 20 30 40 50 60 70
г, с t,c
Рис. 7. Динамика величин вертикального заплеска Ж-волн с лидирующей волной положительной (сплошная линия) и отрицательной (штриховая линия) полярностей. /3 = 15° (а); 5° (б)
ф • • 1 1 / 2 3 4
• •
1(1 -----.
¡: 1
1:!
1 ! л •• • •
• ! • | \ г • ф
• • 1: : • / г ;
• • г *' у \ • • • * г / / \
• ? V/ у • 1 и _ 3 : / / л \ \
' и •• 1 * • • \ / \ N.
ч %
0 10 20 30 40 50 60 70
г, с
Рис. 8. Динамика величин вертикального заплеска Ж-волн с лидирующей волной отрицательной полярности. /3 = 2.5, 5, 10, 15° — кривые 1-4 соответственно
Рассмотренные ситуации, как уже было отмечено, сходны с поведением волн цунами при их взаимодействии с берегом, часто начинающемся с осушения прибрежной зоны — отката линии уреза, за которым следует ее стремительное возвращение — накат волны на берег. Возможность такого варианта развития катастрофического процесса определяется различными факторами: механизмом источника возмущения, его удаленностью от "защищаемого" участка побережья, ориентацией относительно береговой линии, ее изрезанностью, наконец, характеристиками рельефа дна вдоль трассы распространения волны и в прибрежной зоне, типом рельефа суши, примыкающей к месту выхода волны на берег. Эти обстоятельства стимулируют интерес к исследованию ^-волн. Тем не менее в большинстве работ, связанных с моделированием заплеска волн цунами
Рис. 9. Динамика величин вертикального заплеска на откос с углом наклона /3 = 15° волн различной конфигурации: одиночной волны положительной полярности (4), одиночной волны отрицательной полярности (2), Ж-волны с лидирующей волной повышения (3) и волной понижения (1)
на берег, в качестве начального возвышения рассматриваются одиночные волны соли-тонного типа положительной или отрицательной полярности.
Приведенные на рис. 9 графики демонстрируют различия величин вертикального заплеска для волн различной конфигурации, подходящих к склону с [3 = 15°. Это, во-первых, одиночные волны положительной и отрицательной полярности, представляющие собой соответствующие фрагменты начального возвышения (6), расположенные так, чтобы совпадали точки положения их передних фронтов, и, во-вторых, ^-волны двух типов, форма которых также задавалась соотношениями (6). Излагаемые здесь результаты дополняют результаты статьи [1], авторы которой в силу ограничений применяемых ими аналитических методов приводят только значения максимальных за-плесков, в то время как методика численного моделирования позволяет выявить более тонкие детали изучаемых процессов.
Из представленных на рис. 9 графиков следует, что одиночная волна повышения и ^-волна с лидирующим фрагментом положительной полярности приводят к практически совпадающим максимальным значениям вертикального заплеска в стадии наката на склон, но различаются в фазе отката, в то время как аналогичные характеристики для волны понижения и ^-волны с лидирующим фрагментом отрицательной полярности близки только в начальной фазе отката, когда экстремальное значение заплеска по своему абсолютному значению оказывается меньше у ^-волны. Однако ее вертикальный заплеск в стадии наката превышает аналогичную характеристику одиночной волны более чем в два раза. Наблюдаются различия и в последующих колебаниях точки уреза. Особенно заметны они в первом случае, когда величина второго вертикального заплеска в стадии наката для ^-волны с лидирующим фрагментом положительной полярности оказывается в несколько раз больше, чем для одиночной волны повышения.
Наиболее интересными являются результаты сравнения ^-волн разной полярности между собой, качественно совпадающие с выводами авторов [1] (естественно, для случая "обобщенных" ^-волн) о том, что экстремальные характеристики заплеска ^-волны с лидирующим фрагментом понижения оказываются больше в стадиях первичных наката и отката.
Заключение
1. Результаты, представленные в статье, продемонстрировали возможности комбинированного численного метода TVD + SPH в определении основных характеристик заплеска волн различной конфигурации при решения "канонической" задачи о накате длинных волн на плоский откос. Этот метод оказался работоспособным и эффективным не только в области малых амплитуд, теоретически допускающих использование приближения мелкой воды, но и за пределами ее теоретической применимости в области обрушающихся волн.
2. Анализ особенностей заплесков одиночных волн различной полярности (волн повышения и понижения) показал, что волны повышения приводят к большим величинам вертикального заплеска по сравнению с волнами понижения, и это качественно совпадает с выводами авторов работы [4]. Показано также, что взаимодействие падающей и отраженной от склона волн над точкой излома рельефа дна приводит к немонотонности форм свободной поверхности, влияющей на характеристики заплеска.
3. Моделирование заплеска на склон волн более сложной формы, так называемых N-волн, продемонстрировало, что величины вертикального заплеска таких волн с лидирующим фрагментом отрицательной полярности существенно превышают аналогичные величины для подобных волн с лидирующим фрагментом положительной полярности. Результаты сравнения заплесков различных типов N-волн, а также их сопоставления с близкими по форме одиночными волнами качественно согласуются с аналитическими решениями авторов статьи [1].
4. Перспективными направлениями исследований авторы считают создание и тщательный анализ версии алгоритма TVD + SPH для моделирования заплеска катастрофических волн различных конфигураций в пространственной постановке с учетом реальных свойств акваторий и прилегающих к ним побережий.
Благодарности. Исследование выполнено при поддержке Российского научного фонда (проект № 14-17-00219).
Список литературы / References
[1] Tadepalli, S., Synolakis, C.E. The run-up of N-waves on sloping beaches // Proceedings of Royal Society of London. Series A: Mathematical and Physical Sciences. 1994. Vol. 445, No. 1923. P. 99-112.
[2] Шокин Ю.И., Рычков А.Д., Хакимзянов Г.С., Чубаров Л.Б. О численных методах решения задач о накате волн на берег. I. Сравнительный анализ численных алгоритмов для одномерных задач // Вычисл. технологии. 2015. Т. 20, № 5. С. 214-232. Shokin, Yu.I., Rychkov, A.D., Khakimzyanov, G.S., Chubarov, L.B. On numerical methods for solving run-up problems. I. Comparative analysis of numerical algorithms for one-dimensional problems // Computat. Technologies. 2015. Vol. 20, No. 5. P. 214-232. (in Russ.)
[3] Denissenko, P., Didenkulova, I., Rodin, A., Listak, M., Pelinovsky, E. Experimental statistics of long wave runup on a plane beach // Journal of Coastal Research. 2013. Special Issue No. 65. P. 195-200.
[4] Rodin, A., Didenkulova, I., Pelinovsky, E. Numerical study for run-up of breaking waves of different polarities on a sloping beach // Extreme Ocean Waves. Eds: E. Pelinovsky, C. Kharif. Switzerland: Springer International Publishing, 2016. P. 155-172. DOI: 10.1007/9783-319-21575-4.
[5] Храпов С.С., Хоперсков А.В., Кузьмин Н.М., Писарев А.В., Кобелев И.А. Численная схема для моделирования динамики поверхностных вод на основе комбинированного SPH-TVD подхода // Вычисл. методы и программирование. 2011. Т. 12. С. 282-297. Hrapov, S.S., Khoperskov, A.V., Kuzmin, N.M., Pisarev, A.V., Kobelev, I.A.
Numerical scheme for simulating the dynamics of surface water on the basis of the combined SPH-TVD approach // Computational Methods and Programming Techniques. 2011. Vol. 12. P. 282-297. (in Russ.)
[6] Bautin, S.P., Deryabin, S.L., Sommer, A.F., Khakimzyanov, G.S., Shokina, N.Yu.
Use of analytic solutions in the statement of difference boundary conditions on a movable shoreline // Russian Journal of Numerical Analysis and Mathematical Modelling. 2011. Vol. 26, No. 4. P. 353-377.
[7] Диденкулова И.И., Пелиновский Е.Н., Диденкулов О.И. Накат длинных уединенных волн различной полярности на плоский откос // Изв. РАН. Физика атмосферы и океана. 2014. Т. 50, № 5. С. 604-611.
Didenkulova, I.I., Pelinovsky, E.N., Didenkulov, O.I. Run-up of long solitary waves of different polarities on a plane beach // Izvestiya, Atmospheric and Oceanic Physics. 2014. Vol. 50, No. 5. P. 532-538.
[8] Родин А.А. Численные расчеты наката обрушенных одиночных волн на плоский откос // Тр. Нижегородского гос. техн. ун-та им. Р.Е. Алексеева. 2013. № 1 (98). С. 36-43. Rodin, A.A. Numerical calculations of run-up of solitary breaking waves on a plane beach // Transactions of Nizhni Novgorod State Technical University named by R.E. Alekseev. 2013. № 1 (98). P. 36-43. (in Russ.)
[9] Clawpack-5 home page. Available at: http://www.clawpack.org/
[10] Вольцингер Н.Е., Клеванный К.А., Пелиновский Е.Н. Длинноволновая динамика прибрежной зоны. Л.: Гидрометеоиздат, 1989. 272 c.
Voltsinger, N.E., Klevannyy, K.A., Pelinovskiy, E.N. Long-wave dynamics of the coastal zone. Leningrad: Gidrometeoizdat, 1989. 272 p. (in Russ.)
[11] Пелиновский Е.Н., Родин А.А. Нелинейная деформация волны большой амплитуды на мелководье // Доклады Академии наук. 2011. Т. 438, № 3. С. 337-340. Pelinovsky, E.N., Rodin, A.A. Nonlinear deformation of a large-amplitude wave on shallow water // Doklady Physics. 2011. Vol. 56, No. 5. P. 305-308.
[12] Пелиновский Е.Н., Родин А.А. Трансформация сильно нелинейной волны в мелководном бассейне // Изв. РАН. Физика атмосферы и океана. 2012. Т. 48, № 3. С. 343-349. Pelinovsky, E.N., Rodin, A.A. Transformation of a strongly nonlinear wave in a shallow-water basin // Izvestiya, Atmospheric and Oceanic Physics. 2012. Vol. 48, No. 3. P. 383-390.
[13] Synolakis, C.E. The runup of long waves: Ph. D. Thesis. Pasadena: California Institute of Technology, 1986. 228 p. 91125.
[14] Okada, Y. Surface deformation due to shear and tensile faults in the half space // Bulletin of the Seismological Society of America. 1985. Vol. 75. P. 1135-1154.
[15] Гусяков В.К. Остаточные смещения на поверхности упругого полупространства. Условно-корректные задачи математической физики в интерпретации геофизических наблюдений. Новосибирск: ВЦ СО АН СССР, 1978. C. 23-51.
Gusiakov, V.K. Residual displacements on the surface of elastic half-space. Conditionally correct problems of mathematical physics in interpretation of geophysical observations. Novosibirsk: Computer Center, 1978. P. 23-51. (in Russ.)
Поступила в 'редакцию 31 августа 2015 г.
On numerical methods for solving run-up problems. II. Experience with model problems
Shokin, Yuriy I.1, Ryohkov, Alexandr D.1', Khakimzyanov, Gayaz S.1'2, Chubarov, Leonid B.1,2'*
institute of Computational Technologies SB RAS, Novosibirsk, 630090, Russia 2Novosibirsk State University, Novosibirsk, 630090, Russia
* Corresponding author: Chubarov, Leonid B., e-mail: chubarovster@gmail.com
In this article we consider numerical simulation of the run-up of long surface waves on sloping beaches. Using a method based on the combination of TVD and SPH, the authors use the shallow water approximation to solve the well known model problem of the run-up of a wave approaching from an area of constant depth towards a beach with a constant slope. The numerical method works well and is effective not only in the case of small amplitudes, but also for modelling of breaking waves that is outside of the theoretical limits of applicability of the shallow water theory. An analysis of the run-up of waves of various configurations and amplitudes on the slopes of different steepness was performed. Basic characteristics of the phenomenon under study are reproduced both for weakly nonlinear (nonbreaking), and for essentially nonlinear (breaking) waves. Particular attention is paid to the study of the processes of run-up of the so-called N-waves having a direct relation to tsunami waves generated by underwater earthquakes located at a small distance from the shore.
Qualitative and partial quantitative comparison with the results of numerical calculations of other authors [4] are presented. The differences in the results due to differences in numerical algorithms are highlighted. The results of modeling the W-waves and comparing them with the run-up of solitary waves are qualitatively consistent with the conclusions of [1].
Keywords: surface waves, shallow water, run-up, numerical simulation.
Acknowledgements. The study was supported by the Russian Science Foundation (Project No. 14-17-00219).
Received 31 August 2015
© ICT SB RAS, 2015
