Численное моделирование наката волн цунами на берег произвольного профиля Текст научной статьи по специальности «Физика»

УДК 004.942

А. Г. Марчук \ П. С. Мошкалев

2

1 Новосибирский государственный университет ул. Пирогова, 2, Новосибирск, 630090, Россия

2 Институт вычислительной математики и математической геофизики СО РАН пр. Акад. Лаврентьева, 6, Новосибирск, 630090, Россия

2

E-mail: [email protected]; [email protected]

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ НАКАТА ВОЛН ЦУНАМИ НА БЕРЕГ ПРОИЗВОЛЬНОГО ПРОФИЛЯ

Предлагается новый метод для численного моделирования процесса наката длинных волн на берег. Для описания процесса распространения волн до точки уреза используются нелинейные уравнения мелкой воды. После этого во время процесса наката волны на берег используется специальный алгоритм для определения положения подвижной точки уреза и оценки параметров потока. Данный алгоритм основан на законах сохранения массы и энергии. С помощью этого алгоритма были проведены серии вычислений в одномерной прибрежной области океана. Был определен профиль берега, который дает максимальную высоту заплеска для волны с фиксированными начальными параметрами. В работе представлены результаты численного моделирования процесса наката волны цунами на берег с реальным профилем, расположенным в преф. Акита (Япония).

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

Одним из наиболее важных вопросов в прогнозе цунами является оценка высот наката волны цунами в различных точках побережья и затем определение возможных зон затопления в результате атаки такой волны. Численные методы моделирования распространения волн цунами в глубоком океане и в прибрежной зоне уже достаточно давно применяются широким кругом исследователей этого явления. Некоторые из них для нахождения заплеска (максимальной высоты над невозмущенным уровнем океана, достигнутой накатившейся волной) и возможных зон затопления применяют упрощенные соотношения, связывающие высоту заплеска с высотой волны у берега [1], в то время как эта высота существенно зависит от берегового профиля выше точки уреза. Самый простой подход, применяющийся при численном моделировании, это когда реальный берег заменяют вертикальной стенкой, на которой реализуется полное отражение волны. Другой подход состоит в ступенчатом представлении берегового профиля.

Постановка задачи

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

где и - скорость воды; п - смещение водной поверхности; Н - значение глубины, а g - гравитационная постоянная. Пусть в одномерной прибрежной области океана глубина зависит только от расстояния до береговой линии. Если фронт плоской волны параллелен берегу, то задача является одномерной и описывается уравнениями (1). От левой границы одномерной

Марчук А. Г., Мошкалев П. С. Численное моделирование наката волн цунами на берег произвольного профиля // Вестн. Новосиб. гос. ун-та. Серия: Информационные технологии. 2014. Т. 12, вып. 2. С. 55-63.

ut + иих + g Цх =

П + (и(П + H)) x = 0,

(1)

ISSN 1818-7900. Вестник НГУ. Серия: Информационные технологии. 2014. Том 12, выпуск 2 © А. Г. Марчук, П. С. Мошкалев, 2014

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

г

+

Рис. 1. Постановка задачи о накате волны на берег

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

и (х, 0) = п( х,0) = 0. (2)

Исходная волна генерируется принудительным движением воды на левой границе расчетной области:

) = 0,5п0 (1 + - п / 2)),

г е (0,2п / Ь) (3)

и (г) =п(г)4ё / Н.

После того, как волна полностью сформируется (это произойдет, когда значение г станет равным 2пЬ), ставятся граничные условия свободного выхода волн из расчетной области:

^ = 0, ^ = 0. (4)

Эх Эх

Численный алгоритм

Система уравнений (1) с начальными и граничными условиями (2)-(4) решается численно методом конечных разностей. В этом случае все переменные определены только в узлах расчетной сетки. Первоначально правая граница расчетной области находится в точке уреза, которой соответствует узел ' = М. Значения глубины и сеточного рельефа берега задаются массивом чисел Н. При этом в начальный момент количество расчетных узлов сетки равно М. Значения скорости и возвышения свободной поверхности на левой границе получаются исходя из граничных условий (3), а при г > 2п/Ь - из (4). Для вычислений во всех внутренних точках области (' = 2, ..., М - 1) используется явная разностная схема с центральными разностями

+и-ик-иик+ё ЧЬ -Пи = 0,

Аг ' 2 Ах 2 Ах

п"+1 + („ +л;) < - + и-1 Н'+1 +П"+1 - Н'-1 -пП-1 - (5)

Аг 2Ах 2 Ах

пП+1 - 2пП +П"-1

Ах2

Здесь к - это коэффициент искусственной вязкости, которая требуется для устойчивого счета в областях с большими градиентами функции возвышения водной поверхности. На очередном временном шаге сначала из первого уравнения (5) определяются значения скорости ия+1 во всех внутренних расчетных узлах. Затем из второго уравнения (5) можно найти значения возвышения свободной поверхности П+1 также во всех внутренних узлах области.

Пока возвышение свободной поверхности в точке уреза будет равно нулю или близко к нулю (< 0,01 м), там реализуются условия свободного выхода волн из области (4). Как только значение возвышения в этом расчетном узле превысит пороговое значение, включается специальный численный алгоритм для вычисления параметров течения при накате волны на сухой берег. Введем новую расчетную точку В, которая показывает текущее положение точки уреза на каждом шаге по времени. Пусть в некоторый момент времени пА(, когда волновой фронт достиг первоначального уреза (т. е. началось затопление берега), точки Е и А показывают положение водной поверхности в расчетных узлах сетки с номерами М - 1 и М соответственно (рис. 2). Величина возвышения свободной поверхности и скорость водного потока в расчетном узле с номером М - 1 вычисляются по разностной схеме (5). Для нахождения параметров течения в расчетном узле с номером М используется та же разностная схема, но центральные разности заменены левыми

UM+1 UM + un UM UM —1 + g nM nM—1 _ 0

~t~ и i.j : г t; ■ — и,

At Ax Ax

nn+1 —nn un+1 _ un+1 TT +Пп _ H —Пп

' Im ' IM I t U I ^n \ UM UM —1 I „,n+1 M 'IM nM —1 'IM —1 _ n -7"—+(HM +nM)-:-+ UM -:-_0

At Ax Ax

Рис. 2. Определение положения подвижной точки уреза

Для нахождения нового положения подвижной точки уреза D используется закон сохранения массы следующим образом. Найдем объем воды, прошедшей через сечение в узловой точке с номером M за один шаг по времени, который вычисляется как произведение скорости течения в этом узле на среднее арифметическое толщины водного слоя в моменты времени nAt и (n + 1)At. Здесь имеется в виду объем на единицу расстояния вдоль берега (в рассматриваемом сечении это площадь). Сложив его с уже имеющимся там объемом (площадь треугольника ABM) и затем разделив полученный суммарный объем на толщину водного потока ПМ+1 + HM, мы получим расстояние вдоль оси OX от точки C до новой точки уреза D. Горизонтальная проекция скорости точки уреза находится делением горизонтального расстояния, пройденного точкой уреза за последний

временной шаг, на величину этого шага Аг. Если после очередного шага по времени расстояние между проекциями точек С и Б становится больше пространственного шага Ах разностной схемы, то к набору неподвижных узлов расчетной сетки добавляется еще один узел с номером М + 1. При этом параметры течения в новой узловой точке находятся интерполяцией. Соответственно, при откате волны количество расчетных узлов будет уменьшаться и может стать меньше первоначального их числа, так как при откате оголяется дно ниже первоначального уреза. Но основная цель работы - нахождение максимального заплеска при накате длинной волны.

Результаты численных экспериментов

Используя данный алгоритм, мы провели серии экспериментов по накату волн цунами на берега различного профиля. Основной целью этих вычислений было тестирование нового алгоритма с учетом результатов других методов, а также нахождение максимального запле-ска волн цунами на берег в разных ситуациях.

Первая серия вычислений производилась на берег с постоянным уклоном по следующим параметрам (табл. 1).

Таблица 1

Значения параметров вычислительного эксперимента

Параметр Значение

Шаг по пространственным переменным, Ах 10 м

Начальное количество узлов в расчетной области, М0 400

Шаг по времени, Аг 0,2 с

Параметр длины волны в выражении (3), Ь 0,2

Начальная амплитуда волны, п0 2 м

Волна генерировалась на расстоянии 4 километра от начальной береговой линии. Профиль донной поверхности определялся следующим образом: от правого граничного узла сетки (начальной точки уреза) уклон донной поверхности был постоянным и равнялся tg(a) = 0,1 до уровня 100 метров. В дальнейшем он не менялся вплоть до левого узла сетки. В табл. 2 показана высота заплеска при накате волны цунами на берег с различным уклоном.

Таблица 2

Вертикальный заплеск волны при различных береговых уклонах

Угол /ё(а) 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

Максимальный заплеск 19,5 19,3 18,9 18,5 18,1 17,8 17,5 17,3 17,1 16,9

При распространении от левой границы до точки первоначального уреза амплитуда волны увеличивалась с 2 до приблизительно 8 метров. Из табл. 2 видно, что при увеличении угла наклона берега величина максимального заплеска уменьшается, что согласуется с результатами, полученными ранее как численными методами [2], так и аналитически [3-4]. Далее показан профиль водной поверхности в момент достижения волной своего максимального за-плеска (рис. 3).

Рис. 3. Профиль волны при уклоне ?ё(а) = 0,1

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

Н(г) = т-гё(а) + С-8т(пг/г0), 0 < г < г0.

В данном случае параметр г изменяется от начального узла сетки направо, а г0 - длина генерируемой волны (в данном случае она равна 250 м). Все остальные параметры остались такими же, как и в предыдущем случае. В этих расчетах максимальная высота волны в точке первоначального уреза составила 7,6 м. В табл. 3 показаны величина вертикального заплеска в зависимости от изменения положительного параметра С (вогнутая береговая поверхность) , а также береговой профиль и профиль свободной поверхности при параметре С = 4 (рис. 4).

Таблица 3

Заплеск волны в зависимости от параметра кривизны С

Параметр кривизны С 0,0 1,0 2,0 3,0 4,0 5,0 6,0 7,0

Высота заплеска 23,067 23,086 23,132 23,73 23,86 23,54 23,68 23,63

Рис. 4. Профиль волны при параметре С = 4

Далее показаны результаты вычислений при отрицательном значении параметра С (выпуклая береговая поверхность) (табл. 4), а также береговой профиль и профиль свободной поверхности при параметре С = -6 (рис. 5).

Таблица 4

Заплеск волны в зависимости от отрицательного параметра кривизны С

Параметр кривизны С -1,0 -2,0 -3,0 -4,0 -5,0 -6,0

Высота заплеска 22,58 22,82 22,40 22,37 21,88 20,86

Рис. 5. Профиль волны при параметре С = -6

Как видно из табл. 3 и 4 максимальный вертикальный заплеск волны цунами был достигнут при параметре С = +4.

Параметры течения во время наката и отката волны

Предлагаемый метод дает возможность вычислять и отслеживать параметры течения воды на протяжении всего процесса (накат волны цунами на сухой берег и затем ее откат с осушением затопленного во время фазы наката берега). В частности, скорость потока воды над затопленным берегом является одним из важнейших параметров этого процесса. Она определяет силу воздействия водного потока на береговые сооружения. Потенциальная сила потока на единицу ширины берега пропорциональна толщине водного слоя и квадрату скорости потока h XV2. В табл. 5 приводятся значения параметров водного потока в береговых расчетных узлах, которые были затоплены в процессе наката волны. В первой колонке приведены параметры исходной волны: длина волны, определяемая параметром Ь (см. формулу (3)), и начальная амплитуда волны Т]о. Во второй колонке представлены оцениваемые параметры волнового потока. В последующих колонках содержатся максимальные значения указанных параметров в береговых расчетных точках.

В этой серии расчетов наклон берега tё(a) не менялся и был равен 0,1. Анализ приведенных в табл. 5 данных позволяет заключить, что скорость водного потока при накате волны достигает своего максимума в заключительной фазе процесса наката примерно в последней четверти от общего количества затапливаемых расчетных узлов (считая от первоначальной точки уреза). В частности, максимальная скорость водного потока во время наката цунами высотой (у берега) 6 м была зафиксирована в 6-м расчетном узле (от уреза) из всех 20, до которых дошла волна. В этом же численном эксперименте скорость водного потока при откате волны во всех точках была больше, чем при накате. При этом максимум этой скорости наблюдался вблизи точки первоначального уреза и был равен 24 м/с. Интегральные параметры hXV и hXV2 также достигали своего максимума вблизи точки первоначального уреза из-за достаточно большой толщины водного потока в этом месте.

Таблица 5

Параметры водного потока во время наката и отката волны

Параметры потока Зарегистрированные значения 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

Ъ = 0,05 % = 4 ^шах 4,67 4,50 4,26 3,93 3,60 3,44 3,02 2,80 2,25 1,79 1,54 0,09

^шш -15,3 -13,9 -13,1 -12,0 -11,2 -10,1 -9,23 -8,15 -6,97 -5,54 -3,66 -0,52

Ъ = 0,05 % = 5 ^шах 4,69 4,71 4,80 4,77 4,69 4,60 4,48 4,38 4,22 4,00 3,76 3,43 2,95 2,23 1,35

^шш -16,4 -15,5 -14,9 -13,9 -13,3 -12,4 -11,8 -10,9 -10,1 -9,21 -8,24 -7,14 -5,90 -1,36 -1,95

Ъ = 0,05 % = 6 I/ ^ шах 4,80 5,04 5,45 5,63 5,72 5,76 5,74 5,73 5,66 5,52 5,31 5,19 4,96 4,63 4,26 3,73 3,12 2,30

к ■ ^ шш -24,3 -22,3 -20,0 -19,0 -17,7 -16,5 -15,6 -14,5 -13,7 -12,8 -12,0 -11,1 -10,2 -9,31 -8,28 -7,25 -5,81 -1,13

Ъ = 0,03 % = 5 ^шах 2,11 2,09 2,62 2,71 2,54 2,20 2,17 2,06 1,53 1,32 1,04 0,25

^шт -14,9 -13,9 -13,2 -12,4 -11,4 -10,4 -9,31 -8,21 -6,99 -5,61 -3,84 -1,05

Ъ = 0,05 Ло = 5 (А'^Отах 31,8 28,9 26,1 23,1 20,6 17,9 15,2 12,9 10,2 8,12 6,15 4,40 2,91 1,56 0,47

-17,1 -12,3 -38,6 -34,3 -30,6 -26,8 -23,2 -19,8 -16,4 -13,2 -10,2 -7,31 -1,63 -2,28 -0,43

тах)2 144 127 115 101 88 76,8 63,4 52,0 41,8 31,6 24,5 17,5 11,6 6,8 3,0 0,5

А^тш)2 603 543 461 407 341 292 240 197 156 120 88,0 60,9 37,5 19,5 7,0

Ъ = 0,05 Ло = 6 (/? ^Отах 45,4 43,7 38,9 36,1 33,3 30,0 26,7 23,4 21,4 17,8 15,2 12,7 10,4 8,07 6,01 4,16 2,58 1,17

(/? ^Отах -60,2 -55,7 -51,7 —17,3 -13,6 -39,3 -35,7 -31,9 -28,2 -24,9 -21,3 18,2 -14,9 -12,0 -9,07 -6,39 -3,89 -1,73

А^тш)2 245 227 215 184 183 160 139 122 107 93,9 72,6 63,2 50,6 39,1 27,8 19,3 12,0 6,4

А^тш)2 877 788 702 627 552 490 423 371 314 267 221 179 144 108 80,0 53,5 32,7 16,1

Моделирование наката волны на берег с реальным профилем

Проведен также вычислительный эксперимент по накату цунами на реальный профиль берега. Береговой профиль (рис. 6) был взят из Отчета по воздействию цунами 26 мая 1983 г. на побережье преф. Акита 1. В данном месте (неподалеку от д. Минехама) был зафиксирован накат волны с вертикальным заплеском, равным 14,08 м. Согласно результатам моделирования по описанному здесь методу для достижения такого значения волна с периодом 65 с должна иметь высоту 8 м перед накатом на берег. Волновой профиль при накате такой волны показан на рис. 6-7.

Ш №

Рис. 6. Вертикальный профиль берега вблизи населенного пункта Минехама (преф. Акита, Япония)

Рис. 7. Моделирование цунами в Японском море 26.05.1983 на побережье преф. Акита. Показан профиль волны цунами в момент наибольшего заплеска

1 Akita Prefecture Office Report. 1984. ¡Survey Report of Tsunami of May 26 1983 Along the Coasts of Akita Prefecture. Japan: Akita Prefecture Office, 1984. (in Jap.)

Заключение

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

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

1. Titov V. V., Gonzalez F. Implementation and Testing of the Method of Splitting Tsunami (MOST) // Technical Memorandum ERL PMEL-112, National Oceanic and Atmospheric Administration. Washington DC, 1997.

2. Марчук А. Г. Метод расчета наката длинных гравитационных волн на наклонный берег // Эволюция цунами от очага до выхода на берег. М.: Радио и связь, 1982. C. 64-68.

3. Пелиновский Е. Н. Нелинейная динамика волн цунами. Горький: ИПФ АН СССР, 1982. 226 с.

4. Пелиновский Е. Н. Накат волн цунами на берег. Горький: Ин-т прикладной физики АН СССР, 1985.

Материал поступил в редколлегию 26.03.2014

A. G. Marchuk, P. S. Moshkalov

NUMERIC SIMULATION OF THE TSUNAMI RUNUP PROCESS ON THE SHORE

WITH ARBITRARY PROFILE

In the paper the new method for numerical simulation of the long wave run-up process is proposed. Nonlinear shallow water equations are used for description of the wave propagation up to the water edge point. Then the special algorithm used for estimation of the flow parameters and the location of the moving water edge. It is based on the energy and mass conservation laws. Several series of one-dimensional computations were carried out. The shore profile, that gives the maximum run-up height for the fixed initial wave parameters have been found. Results of modeling of the tsunami run-up on the real shore in the Akita prefecture (Japan) are presented.

Keywords: shallow water equations, tsunami, long wave runup process, difference scheme.

References

1. Titov V. V., Gonzalez F. Implementation and Testing of the Method of Splitting Tsunami (MOST). Technical Memorandum ERL PMEL-112, National Oceanic and Atmospheric Administration, Washington DC, 1997.

2. Marchuk A. G. Metod rascheta nakata dlinnyh gravitacionnyh voln na naklonnyj bereg. [The method for calculation of the long wave run-up on a sloping beach]. Evolution of tsunami from a source till run-up on a shore. Moscow, Radio i Svyaz, 1982, p. 64-68. (In Russ.)

3. Pelinovsky E. N. Nelinejnaya dinamika voln cunami [The nonlinear dynamics of tsunami waves]. Gorky, Institute of Applied Physics of the USSR Academy of Sciences, 1982. (In Russ.)

4. Pelinovsky E. N. Nakat voln cunami na bereg [Tsunami waves run-up to the shore]. Gorky, Institute of Applied Physics of the USSR Academy of Sciences, 1982, 226 p. (In Russ.)