CC BY
68
Вычислительные технологии
Том 19, № 1, 2014
О некоторых численных алгоритмах расчёта наката
• • Tjk
волн цунами в рамках модели мелкой воды. I
С. А. Бейзель1, Н.Ю. ШокинА1, Г. С. Хлкимзянов1, Л. Б. Чубаров1
О. А. КовыркинА2, В. В. Остапенко2 1 Институт вычислительных технологий СО РАН, Новосибирск, Россия 2Институт гидродинамики им. М.А. Лаврентьева СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected]
Представлен метод численного моделирования наката волн цунами на побережье, основанный на использовании модели мелкой воды сразу в двух приближениях: одномерном и двумерном. Вначале по двумерной модели с отражающим краевым условием на берегу рассчитывается распространение волны от источника к побережью. Параметры течения на некоторой изобате, полученные из этого расчёта, используются затем в качестве краевых условий для одномерного моделирования наката вдоль различных сечений, проведённых от этой изобаты до выбранной изолинии на суше. Описана процедура восстановления границы затопления суши по результатам решения одномерных задач. Дано сравнение трёх алгоритмов расчёта движения точки уреза, и приведены некоторые результаты моделирования Японского цунами 2011 года.
Ключевые слова: поверхностные волны, накат волн на берег, максимальный заплеск, линия уреза, плоский откос, неровный склон, уравнения мелкой воды, аналитические решения, численные методы.
Введение
Одной из не решённых до настоящего времени проблем численного моделирования является адекватное описание поведения волн цунами в зоне заплеска и надёжное определение границ затопления суши. Для получения удовлетворительного решения необходимы достаточно подробные сведения о рельефе рассматриваемого участка суши и батиметрии прилегающей акватории (с разрешением не ниже 50 м [1] для объектов, содержащих мелкомасштабные особенности, и 500 м для плавно изменяющихся высот и глубин). Однако доступные данные имеются, как правило, на неприемлемо редких сетках, что является одной из главных причин понижения точности результатов. Важное значение имеет также выбор модели источника начальной волны [2]. При решении задач наката имеются и чисто математические трудности. Так, с приближением волны к берегу линия уреза начинает смещаться, поэтому решение задачи о колебаниях свободной границы приходится искать в области с подвижной границей. Кроме того,
* Работа выполнена при поддержке Президентской программы "Ведущие научные школы РФ" (грант № НШ-6293.2012.9), Программы базовых фундаментальных исследований СО РАН (проект № 1У.38.2.1), Программы интеграционных исследований СО РАН (проекты № 37б и 117а), Программы Президиума РАН (проекты № 4.8 и 4.10), РФФИ (гранты № 12-05-00894-а и 12-01-00721-а), Совета по грантам Президента РФ для поддержки молодых российских учёных (грант МК-3477.2013.1).
вблизи береговой линии начинает расти число Фруда, достигая бесконечного значения на самой подвижной линии уреза, что создает большие проблемы при численном моделировании, поскольку в процессе взаимодействия волн с берегом происходят многократные переходы как от докритического режима к сверхкритическому, так и в обратном направлении — с образованием гидравлических скачков. Более того, на линии уреза собственные значения матрицы Якоби системы уравнений мелкой воды, которая чаще других используется при решении задач наката волн на реальное побережье, становятся кратными, а собственные векторы — линейно зависимыми, т. е. система теряет свойство сильной гиперболичности, присущее ей всюду, где полная глубина положительна.
Еще одна проблема, возникающая при моделировании выхода волн цунами на берег, связана с тем, что классическая система базисных законов сохранения теории мелкой воды, состоящая из законов сохранения массы и полного импульса [3-5], правильно передавая параметры прерывных волн, распространяющихся по жидкости конечной глубины [3], не допускает распространения прерывных волн по сухому дну. Точные решения, описывающие в рамках этой системы процесс течения воды по сухому дну, являются непрерывными волнами понижения (например, непрерывная волна понижения, возникающая при разрушении плотины над горизонтальным дном). Учёт донного трения, которое представляет собой распределённый источниковый член, не дающий вклада в условия Гюгонио на фронте прерывной волны, не может изменить эту ситуацию [6, 7]. В то же время лабораторные эксперименты [8-12] и результаты натурных наблюдений показали, что эти непрерывные решения существенно завышают скорость распространения передней кромки волны и заметно искажают профиль её поверхности. В опытах фронт волны, распространяющейся по сухому руслу, является существенно более крутым и на нём происходят обрушения, характерные для прерывных волн.
Перечисленные обстоятельства указывают на то, что процесс нелинейного взаимодействия длинных волн с реальной береговой линией является весьма трудным для численного моделирования, хотя подходов к решению проблемы оценки зон затопления предложено немало [13, 14]. Наибольшее распространение для расчёта взаимодействия волн с берегом получили конечно-разностные методы сквозного счёта, в которых нестационарная область течения с подвижной линией уреза вкладывается в область простой формы, например в прямоугольник, содержащий акваторию, линию уреза и прилегающую к берегу часть суши. В этой области с неподвижной границей решаются те же уравнения волновой гидродинамики, что и в исходной области течения, при этом либо суша покрывается тонкой плёнкой воды, удерживаемой силой трения, либо в точках суши скорость и полная глубина жидкости полагаются равными нулю. Для приближённого определения положения подвижной линии уреза проводится анализ вычисленных значений полной глубины воды на каждом шаге по времени [15], и линия уреза определяется, например, как граница области, в узлах которой полная глубина равна нулю [16] или не превышает заранее заданной малой величины [17, 18]. Этот же подход использовался для решения задач наводнения-осушения на основе методов конечных элементов и конечных объёмов [14, 19-25]. Недостатком такого подхода является довольно грубое определение положения подвижной линии уреза, а достоинством — возможность решения задач наката для береговой линии сложной формы.
Другой подход подразумевает счёт с выделением линии уреза, т. е. в области с частично подвижной границей [7, 26-30]. В этом методе область, занятая жидкостью, покрывается подвижной сеткой, одна из граничных координатных линий которой совпадает с линией уреза. Общий недостаток указанных работ — отсутствие подробного
описания постановки разностных краевых условий на подвижной линии уреза и их обоснования. Чаще всего ограничиваются простейшими вычислительными процедурами типа метода "экстраполяции на границу" [31]. В работе [32] выполнено аналитическое исследование решений одномерных нелинейных уравнений мелкой воды в окрестности границы вода — суша и разработаны новые аппроксимации краевых условий в подвижной точке уреза, существенно использующие полученный аналитически закон движения этой точки в процессе наката и отката волн на криволинейный склон. Данный подход трудно применить при решении плановых задач наката-отката, поскольку по мере продвижения волны на сушу даже плавная начальная линия уреза может принять чрезвычайно сложную форму, а первоначально односвязная область, занятая водой, может при затоплении части суши стать многосвязной за счёт образования "островов" и "озер".
В работе [33] предложен метод, позволяющий в рамках первого приближения теории мелкой воды моделировать процесс распространения прерывных волн по сухому дну. В основе этого метода лежит модифицированный закон сохранения полного импульса, в котором учитываются сосредоточенные потери импульса, связанные с образованием локальных турбулентно-вихревых структур на фронтах прерывных волн. Количественная оценка таких потерь получена путём вывода уравнений мелкой воды из уравнений Навье — Стокса [34] с учётом вязкости, влияние которой в областях турбулентного течения, описываемых прерывными волнами, резко возрастает. Возникающий при этом эвристический параметр подбирается путём согласования с результатами лабораторных экспериментов [11, 12]. В [18, 35] данный метод был применён для численного моделирования процесса распространения прерывных волн по сухому дну (включая их накат на наклонный берег) в одномерной и двумерной постановках.
В настоящее время численное моделирование на основе какого-либо одного метода часто заменяется комбинированными вычислительными алгоритмами. Так, в работе [36] реализована методика моделирования наката цунами на побережье реальной акватории, позволяющая повысить степень достоверности результатов моделирования волновой картины вблизи реального побережья за счёт предварительного исследования особенностей волновой картины вблизи берега для последовательности упрощённых модельных акваторий, наследующих базовые характеристики изучаемого фрагмента прибрежной зоны, и сопоставления на каждом этапе результатов, полученных с помощью различных вычислительных методов. На таких модельных задачах выявляются основные характеристики исследуемого явления, в дальнейшем претерпевающие детализацию и уточнение по мере приближения к реальным условиям, что может облегчить интерпретацию сложной картины наката на реальное побережье. В последнее время для более детального изучения процесса наката на выделенный малый участок побережья применяется также технология вложенных областей и вложенных расчётных сеток [36-38].
В работе [39] предлагается распространение волны цунами рассчитывать по плановой модели до некоторой изобаты вблизи берега и полученные на этой изобате параметры волны использовать затем для расчёта величин максимальных заплесков на берег с применением аналитических формул, полученных для плоского откоса. Отметим, что для плоского откоса аналитическая теория наката является весьма развитой. Аналитические решения системы одномерных нелинейных уравнений мелкой воды, описывающие накат и откат необрушающихся стоячих волн на плоский откос, найдены в работе [40]. Предложенный в ней подход в дальнейшем был развит в статьях [41, 42], в которых учтено влияние начальной формы волны на поведение траектории точки уре-
за и её скорости. В [43] получена аналитическая формула для вычисления максимума вертикального наката уединённой волны на относительно крутой плоский откос, сопряжённый с участком горизонтального дна. Установлено [44-47], что для относительно крутых плоских откосов максимальные значения высоты наката необрушающихся волн и глубины отката будут одинаковыми для линейной и нелинейной теорий мелкой воды, что позволяет рассчитывать максимальные заплески или откаты с помощью простых инженерных формул линейной теории. Подробно исследовано влияние формы набегающей волны на максимальные значения высоты и скорости наката на плоский откос, глубины и скорости отката [48-50]. Компактное описание подходов к получению аналитических решений задачи наката-отката дано в [31] ив обзорной статье [51].
Придерживаясь технологии [36] поэтапного решения сложных задач, в настоящей работе рассматривается решение задачи наката в упрощённой одномерной постановке, предполагающей зависимость батиметрии акватории и рельефа суши только от одной из горизонтальных координат. Суть методики в том, что распространение волны от источника к берегу рассчитывается по плановой модели мелкой воды с учётом пространственного характера распределения батиметрии, но на линии невозмущённого положения уреза устанавливается вертикальная стенка, отражающая набегающие волны. При этом выполняется такая модификация распределения глубин в прибрежной зоне, что глубины, меньшие некоторого заданного порогового значения hwaii, заменяются данным значением. Во время расчёта ведётся непрерывная запись параметров волн на изобате hin, расположенной на удалении от берега, достаточном для минимизации влияния отражённых волн на форму главной части падающей волны. Далее эти данные берутся в качестве краевых условий для одномерных расчётов вдоль нескольких сечений, проведённых от изобаты hin в сторону берега. Таким образом, в отличие от работы [39] "одномерный" накат определяется не по аналитическим формулам, а численно. В статье приводятся результаты расчётов наката с использованием трёх численных методов: сквозного счёта и с выделением линии разрыва на основе классической модели мелкой воды и метода сквозного счёта [18, 35], основанного на модифицированной модели мелкой воды.
1. Об алгоритмах расчёта движения линии уреза
Рассматриваемая методика расчёта наката волн цунами на выбранные участки побережья предполагает проведение предварительного 20-моделирования распространения волн в глобальной области, включающей сейсмический очаг и простирающейся почти до береговой линии. На этом этапе расчёт наката на берег не производится, поскольку вдоль береговой линии на заданной глубине hwaii устанавливается "вертикальная стенка", на которой принимается условие отражения. Начальное возвышение z = По(x,y) свободной поверхности воды предполагается равным вертикальному остаточному смещению поверхности дна, рассчитываемому в сейсмическом очаге по методике [2]. Расчёт процесса распространения волн производится с использованием разработанного ранее программного комплекса MGC [52], основанного на аппроксимации гидродинамической модели нелинейных уравнений мелкой воды.
Предполагается, что береговая линия не изрезана мелкомасштабными заливами и не экранируется от непосредственного воздействия очага цунами, а рельефы прилегающих акватории и суши не содержат мелкомасштабных особенностей. Затем в горизонтальной проекции этого участка проводятся сечения, соединяющие hcoast изолинию рельефа
суши с hin изобатой акватории (рис. 1, а, на котором чёрная кривая соответствует береговой линии, белая на суше — hcoast изолинии, белая в акватории — hin изобате). Обычно значение hcoast полагалось равным 20 м, а при аномально малых или больших значениях наката концы сечений на суше выбирались соответственно на 10- или 30-метровых изолиниях рельефа. Предполагается, что каждое такое сечение может пересекаться с изолиниями рельефа суши в затопляемой зоне только в одной точке. В более ограничительной форме математическая формализация этого предположения может выглядеть следующим образом. Пусть z = hbt ) — функция, описывающая рельеф дна и прилегающей суши, A1Б1 — какое-то из проведённых сечений (см. рис. 1, б) с конечными точками xAl на суше и xBl в море. Если отрезок A1Б1 задать параметрическим уравнением x(p) = (1 — p)xAl + pxBl (0 ^ p ^ 1, x = (x,y)), то рельеф суши и батиметрия дна в рассматриваемом сечении описываются функцией
z = hi (p) = hbt (x(p),y (p)).
Будем предполагать, что на суше (при h1 (p) > 0) функция z = h1 (p) является монотонно убывающей, что гарантирует существование единственного корня p* уравнения z* = hi (p), где 0 < z* < hcoast.
В ходе 2Б-расчётов определяется динамика возвышения свободной поверхности ("полные" мареограммы) в мареографных точках, установленных на "морских" окончаниях сечений (на глубине hin, м). Рассчитанные по 2Б-модели мареограммы используются затем в качестве краевых условий для одномерных расчётов наката волн на берег вдоль выбранных сечений, при этом начальные данные определяются из распределения начального 2Б-смещения свободной поверхности над соответствующим сечением, а из "полной" мареограммы используется только её главная часть — первая прошедшая к берегу волна, после попадания которой в расчётную одномерную область краевое условие на "морском" окончании сечения заменяется на неотражающее условие.
Рис. 1. Сечения для одномерных расчётов наката волн цунами (а) и граница заливания х*^) между двумя сечениями (б)
На основе таких одномерных расчётов вычисляются значения вертикального наката (высоты) и горизонтального заплеска (дальности) вдоль выбранных сечений. Далее с учётом цифрового рельефа суши определяется положение границы зоны затопления. Для приближённого определения этой границы поступаем следующим образом. Пусть и — значения максимальных вертикальных заплесков в двух соседних сечениях ЛхБх и Л2Б2 (см. рис. 1, б). Соединим концевые точки этих сечений отрезками прямых хА(д) = (1 — д)хА1 + дхАз и хв(д) = (1 — д)хв1 + дхв2 и для каждого значения параметра д (0 ^ д ^ 1) проведём вспомогательное сечение
х(р; д) = (1 — р)хА(д)+ Рхв (д). (1)
Полагаем, что в этих сечениях величина максимального вертикального заплеска линейно зависит от параметра д, т.е. г*(д) = (1 — д)г2 + д^. В силу сделанного выше предположения уравнение г*(д) = кы (х(р; д)) имеет единственный корень р2(д), что позволяет, используя уравнение (1), найти точку
х2(д) = (1 — р*(д))х^(д) + р2(д)хв(д), д е [0, 1], (2)
лежащую на границе зоны затопления и отвечающую значению параметра д, при этом точки х*(0) и х*(1) соответствуют вычисленным значениям максимальных заплесков г2 и в рассматриваемых сечениях ЛхБх и Л2Б2. Формула (2) задает в параметрической форме фрагмент границы зоны затопления между любыми двумя соседними сечениями. Сшивая эти фрагменты, получаем искомую границу для всего рассматриваемого участка.
Отметим, что при равенстве заплесков г2 = г2 кривая (2) совпадёт с дугой соответствующей изолинии рельефа суши. Если же рельеф суши между соседними сечениями является плоским откосом, то линия (2) будет отрезком прямой.
Далее остановимся на описании алгоритмов, используемых для расчёта наката в одномерных сечениях.
1.1. Метод сквозного счёта
Проведём через одно из сечений вертикальную плоскость и введём в ней локальную декартову систему координат Охг так, чтобы уравнение свободной поверхности покоящейся жидкости имело вид г = 0. Пусть г = — к(х) — известная функция, задающая в рассматриваемом сечении рельеф дна и прилегающей суши. Система уравнений мелкой воды, описывающая течение несжимаемой жидкости со свободной границей, имеет следующий вид:
и + fx = о, 0 <жо(г) <х<ь, £>0, (3)
где
" =( Ни), £ (и)^ НИ +^Я2/0 , С(Х, и)Ч А
£ — время, и(х,£) — осреднённая по вертикали от дна до свободной поверхности скорость, Н = п + к — полная глубина, п(х,£) — отклонение свободной поверхности от невозмущённого уровня, д — ускорение свободного падения, х0(£) — искомая координата подвижной точки уреза. Предполагается, что эта точка является левой подвижной границей области решения, а справа область ограничена неподвижной границей х = Ь.
Уравнения (3) дополняются краевым
Н(хо(£),£) = 0, £ > 0, (4)
и начальными
п(х, 0) = п0(х), и(х, 0) = и0(х), х00 ^ х ^ Ь (5)
условиями, где х00 = х0(0) — известная координата точки уреза в начальный момент времени. О постановке краевых условий на правой границе х = Ь будет сказано ниже при обсуждении результатов конкретных вычислительных экспериментов. Очевидно, что начальное условие для отклонения свободной поверхности от невозмущённого уровня п можно заменить условием для полной глубины
Н(х, 0) = Н0(х) = п0(х) + к(х), х00 ^ х ^ Ь. (6)
При использовании рассматриваемого в настоящем разделе метода сквозного счёта линия уреза не выделялась и уравнения (3) решались во всей области х е (0, Ь) с помощью схемы предиктор-корректор [53] на равномерной сетке Xj = ^'Дх с шагом Ах = Ь/Ж (^ = 0,... , N). На каждом шаге по времени производился анализ вычисленных значений полной глубины Н™+1 на новом временном слое, и если в некотором узле Xj она оказывалась меньше заданного порогового значения ен > 0, то полная глубина в этом узле полагалась равной нулю. Скорость течения и™+1 вычислялась по разностным формулам схемы [53] только в тех узлах, в которых выполнялось неравенство Н™+1 > ен, в других же узлах она, как и полная глубина, занулялась.
Отметим, что для применения метода сквозного счёта необходимо модифицировать и начальные данные задачи, расширив область их определения на весь отрезок [0, Ь] так, что при постановке разностной задачи полагаем Н° = и® = 0 во всех узлах, лежащих левее начальной точки уреза х00, а также в узлах, в которых выполняется неравенство Н0^) < ен .В остальных узлах начальные данные вычисляются с использованием заданных значений начальной скорости и начальной полной глубины: Н° = Н0^-), и0 = и0^-). Очевидно, что при этом в окрестности начальной точки уреза будет возникать нефизичное течение. Так, если в начальный момент времени жидкость покоилась (и0(х) = 0), а свободная граница была невозмущённой (п0(х) = 0), то даже при малом искусственном "осушении дна" справа от начальной точки уреза в последующие моменты времени вместо состояния покоя возникнет течение с колебаниями уровня в окрестности точки х00. При уменьшении значения ен амплитуда колебаний уменьшается, но возрастает возможность потери устойчивости вычислений [16].
Путём значительного усложнения алгоритма обработки точки уреза можно исключить появление указанного численного эффекта [14] или существенно ослабить его, выбирая пороговое значение ен на каждом шаге по времени в зависимости от профилей свободной границы и дна в окрестности точки уреза, однако при этом затрудняется перенос таких алгоритмов на двумерный случай. По этой причине метод сквозного счёта, используемый в настоящей работе, реализован в простейшем классическом варианте с выбранным заранее пороговым значением ен. В ряде работ проводились численные эксперименты с целью определения оптимального значения ен для конкретных задач (см., например, [54]). В рассматриваемых ниже задачах наката использовано значение, равное одному проценту от безразмерной амплитуды (отнесённой к глубине) набегающей волны.
1.2. Расчёт движения линии уреза с использованием подвижных сеток
В расчётах с выделением подвижной точки уреза применялась схема предиктор-корректор [53] на адаптивной сетке {х0} = 0,..., N), левый узел х0 которой совпадал на каждом временном слое £ = £0 с точкой уреза х0(£0). Разностные краевые условия в этой точке использовались для вычисления величин х^+1, Н0+1 и и0+1 на (п + 1)-м слое по времени. Полная глубина определялась непосредственно из краевого условия (4):
Н га+1 Н0
0.
(7)
Для определения нового положения х0+1 точки уреза и скорости её перемещения и0+1 использовались результаты [32] аналитического исследования решений уравнений мелкой воды (3) в окрестности подвижной линии уреза. В [32] решения задачи (3)-(6) построены в малой по времени и пространству окрестности точки уреза для всех возможных режимов взаимодействия волны с берегом: 1) накат необрушающейся волны при условиях Нх(х0(£0), £0) = 0, Нх(х0(£0), £0) < то; 2) накат необрушающейся волны при совпадении в момент времени £ = £0 касательных в точке уреза к свободной границе и рельефу дна (Нх(х0(£0),£0) = 0); 3) "накат с обрушением" (Нх(х0(£0),£0) = то). Для первого режима компоненты вектора решения построены в виде локально сходящихся рядов. Используя их частичные суммы, можно получить следующие формулы:
0+1
т
хо + тио - у дп
п
Х,0,
и
о+1
и
тдпю,0 + 0 9 (2ио,0Нж,0 + и0^ож,0)
где т — шаг по времени, п0о, иП0, НП0, Л.хх,0 — односторонние разностные производные в узле х0 подвижной сетки {х°}, например
о
Но Но А ;
где А = х0 — х0.
Для второго режима закон движения точки уреза определяется из системы двух обыкновенных дифференциальных уравнений [32], в качестве приближённого решения которой можно взять следующее:
о+1
х0 + тип + у 9^;
о
х,0,
и
о+1
т
и0 + тд^о,0 + 2 дио^ож,0.
:ю)
Для третьего случая полагается [32]
о+1
х0 + ти* + — 9^;
о
х,0,
и
0+1
0 т 0
и* + тд^0,0 + у ди*^0ж,0,
11)
при этом
и* = и0 — 2
12)
В численных расчётах будем выделять конкретный режим взаимодействия волны с берегом в соответствии со значениями разностной производной (9). Если она удовлетворяет условию
т ^ |Н°01 ^ М,
13)
2
0
2
0
где 0 < т ^ М, т, М — заданные числа, то считаем, что реализуется первый режим, а при выполнении условия
|ЯП.о|<т — (14)
второй. В качестве критерия возникновения третьего случая будем использовать неравенство | |
|НП,о|>М. (15)
При удачном выборе значений параметров т и М формулы (8), (10), (11) могут дать точные значения для положения точки уреза и её скорости. Покажем это на примере движения волны понижения по сухому горизонтальному дну. Предполагается, что справа от плотины, расположенной в точке с координатой ж^, полная глубина слоя покоящейся жидкости равна Нг, а слева воды нет. После разрушения плотины возникает волна понижения
0
Н (ж,*)
при 0 ^ ж ^ ж0(*), + ж ^ при Жо(*) ^ ж ^ жг (*),
Нг
0
и (ж, *)
2
- Д л/^НТ -0
'ж ж<^
'16)
при жг(*) ^ ж ^ Ь,
при 0 ^ ж < ж0(*), при ж0(*) ^ ж ^ жг (*), при жг(*) ^ ж ^ Ь,
'17)
18)
ограниченная слева подвижной точкой контакта с сухим дном
ж0(*) = ж^ - 2^х/уНТ, * > 0,
а справа — звуковой характеристикой
жг(*) = ж^ + 1л/д~Н~г, * > 0,
при этом в подобласти ж0(*) < ж < ж^ течение является сверхкритическим, а в "точке уреза" ж = ж0(*) число Фруда равно бесконечности:
Иш
|и(ж,*)|
11ш 2 (_ 11=00,
^ дН (ж,*) V 5
где 5 = ж — ж0(*) > 0.
Пусть в момент времени * = скорость и положение "точки уреза" определены по точным формулам (17), (18):
и
—2л/дН
г, ж0
ж^—2гЛ/дж
а пороговое значение т выбрано так, что удовлетворяется неравенство (14), т.е.
1 .
Л < т.
9д*2
'19)
(20)
Тогда в численном алгоритме новая скорость и новое положение "точки уреза" будут определяться по формулам (10)
n+1 - - -^x/gHr, ХП+1 = хП + тиП = Xd - 2(tn + т) VgHT.
un = U = — 2 л/ , Жп = X,
0
n
Такие же значения дают формулы точного решения в момент времени £ = £га + т.
Аналогичный результат справедлив для формул (11), описывающих приближённое решение сразу после распада разрыва в точке уреза. Если в рассматриваемой задаче о разрушении плотины формулы (11), (12) применить для расчёта скорости и положения "точки уреза" на первом шаге по времени, то получим
м1 = -2л/дЯг, жО = ж^ + тм* = ж^ - 2тЛ/дЯТ,
что совпадает с точными значениями м(ж0(т),т), ж0(т), вычисленными по формулам (17), (18) в момент времени £ = т, при этом условием применимости формул (11) будет выполнение неравенства (15), т.е.
Н > М. (21)
А v у
Основываясь на вышеизложенном и формулах (20), (21), в расчётах по методу адаптивных сеток были использованы следующие выражения для пороговых значений:
т = стА, М = СМ, (22)
где А = жП — жП — первый шаг неравномерной сетки {жп}, ст, СМ — некоторые постоянные.
1.3. Применение модифицированной модели мелкой воды для расчёта наката прерывных волн на наклонный берег
Дифференциальным следствием классической системы базисных законов сохранения теории мелкой воды (3), состоящей из законов сохранения массы и полного импульса, является уравнение
u
ut + ( у + gH
описывающее изменение локального импульса каждой частицы жидкости вдоль её линии тока. Поскольку система (3) не допускает распространения прерывных волн по сухому дну, которые наблюдаются в лабораторных экспериментах [8-12], в работах [33, 34] было предложено для такого моделирования наряду с законом сохранения массы
Ht + (Hu)x = 0 использовать модифицированное уравнение импульса
gH 2'
(23)
(Hu)t + Hu2 +
2
+ Y
ut + (т + *H) x
g(H + y )hx
(24)
представляющее собой линейную комбинацию законов сохранения полного и локального импульсов. В уравнении (24) y = = const, где h0 — некоторая характерная для
x
x
данного течения постоянная глубина, 8 — малый безразмерный параметр, выбираемый путём согласования с результатами лабораторных экспериментов. Вводя обозначения
д = Ни, Н = Н + 7, Я = Н + Я = Н + 27, перепишем уравнение (24) в виде закона сохранения модифицированного импульса
д* + (#«2 + д™) = дН (25)
В работе [18] построена полунеявная условно устойчивая разностная схема первого порядка, аппроксимирующая модифицированную систему законов сохранения теории мелкой воды (23), (25). В этой схеме подавление нефизических осцилляций на фронтах прерывных волн достигается путём введения специальной искусственной вязкости в модифицированное уравнение импульса (25). Моделирование волны, распространяющейся по сухому берегу, осуществляется методом сквозного счёта [18], аналогичным описанному в разделе 1.1.
2. Верификация алгоритмов
В данном разделе приводятся результаты тестовых испытаний описанных выше трёх методик расчёта наката. Рассматриваются задача наката на плоский откос, для которой имеется возможность сравнения с точным аналитическим решением, и задача наката на неровный откос, заданный аналитической формулой и являющийся модельным для некоторых типов реального дна.
2.1. Накат на плоский откос
Рассматривается задача о набегании уединённой волны на плоский откос, сопрягающийся с горизонтальным дном глубины Л,0 = 1. Рельеф дна и прилегающей суши описывается функцией
,, , , ¿0 — жtgв при 0 ^ ж ^ Жs, — Л(ж) = < s (26)
1 —1 при ж. ^ ж ^ Ь,
где в > 0 — угол наклона откоса, ¿0 > 0 — высота "суши" в точке ж = 0, ж. = (г0+1)с^ в — абсцисса основания склона. Значение ¿0 подбиралось таким, чтобы его не превышал максимальный вертикальный заплеск уединённой волны.
В качестве начальных данных (5), (6) задавались следующие:
П0(ж) = й0сЬ-2 ^Л//3а0д (ж — ж,шЛ , Н0(ж) = Л,(ж) + П0(ж), «0(ж) = — V П°(ж) , (27) V 2и ) Н0 (ж)
где а0 — амплитуда подходящей к склону начальной волны, ж = жад — абсцисса её вершины, ж е [0, Ь]. Отметим, что для "полных" нелинейно-дисперсионных уравнений [55] задача Коши с начальными данными (27) имеет решение в виде уединённой волны, распространяющейся над горизонтальным дном влево с постоянной скоростью
2
и = \/#(1 + а0) без изменения своей формы. Для нелинейных уравнений мелкой воды (3) профиль волны, движущейся над горизонтальным дном, будет меняться со временем, укручиваясь на переднем фронте и выполаживаясь на заднем. Наличие откоса (26) приводит к дополнительной трансформации формы волны при её распространении и накате.
Приведённые ниже результаты расчётов получены при значениях
жад = ж, + 20, Ь = ж, + 40, (28)
что гарантировало расположение основной части начальной волны над горизонтальным участком дна и независимость до некоторого момента времени картины наката волны на откос от краевых условий на правой границе. В силу достаточной удалённости правой границы от начального положения уреза, в точке ж = Ь допустима постановка условий непротекания. Для вычисления возвышения и скорости и*^1 на правой границе (в узле сетки ж* = Ь) использовались также модифицированные неотражающие условия из работы [56], полученные на основе дисперсионного соотношения линеаризованной системы уравнений мелкой воды:
пП+1 - пП + - 4пП-1 + пП-2 пП - 2пП-1 + пП-2 =0 т +С 2Дж 2 (Дж)2 0
га+1 га+1 пга+1 .юга+1
-1 - д„ -1 = 0, (29)
гга+1И*-1 П* П* -1
Дж Дж
где с = и* + \fgHNj, Дж — шаг равномерной сетки. Эта формула обобщается очевидным образом и на случай подвижных неравномерных сеток.
В начальный момент времени точка уреза жоо определялась как корень уравнения П0(ж) + Л,(ж) = 0. В силу пренебрежимой малости величины п0(ж00), координату начальной точки уреза можно положить равной значению ж00 = в.
Первая серия расчётов проводилась для плоского откоса с углом наклона в = 2.88° [43] и г0 = 0.5. В этих экспериментах определялось влияние амплитуды а0/Л,0 уединённой волны (27) на величину максимального вертикального заплеска Я/Л,0. На рис. 2, а маркерами изображены точки в плоскости переменных (в, а0/Л,0), для которых выполнялись эксперименты (в первой серии в = 2.88°, а0/Л,0 = 0.005 ^ 0.05), а на рис. 2, б в виде различных маркеров представлены результаты этих численных экспериментов, демонстрирующие зависимость максимального вертикального заплеска Я/Л.0 на плоский откос от параметра 5, совпадающего с правой частью точного решения [43]
5/4
£=2-831^(£) ■ (30>
В системе координат (5, Я/Л.0) графиком аналитической зависимости (30) является прямая. При оценке отклонения маркеров от этой прямой надо иметь в виду, что формула (30) получена в [43] при условии
0° > (0.288 в)2 (31)
и в предположении, что волны, взаимодействующие с плоским откосом, являются необ-рушающимися. Такое возможно, если при накате выполнено условие [43]
О0 < 0.82^в)10/9, (32)
Я-0
а при откате [29] —
а0 < 0.479^ 0)10/9. (33)
Если необходимо определять только границу зоны затопления, то представляет интерес лишь процесс наката, поэтому в качестве ограничения для начальных данных в этом случае можно использовать формулу (32). Однако для более полной оценки последствий цунами необходимо знать скорости потока на затопляемой части суши как при накате, так и при откате волны, а также полное время нахождения под водой тех или иных участков побережья, т. е. надо иметь представление в целом о взаимодействии волны с берегом. В этом случае необходимо учитывать более жёсткое ограничение (33). Далее будем считать, что условиями применимости формулы (30) являются ограничения (31) и (33). На рис. 2, а приведены границы области, в которой можно использовать формулу (30), и видно, что часть численных экспериментов проводилась для параметров, лежащих за пределами этой области. Соответствующие таким экспериментам результаты изображены на рис. 2, б светлыми маркерами, а результаты расчётов для параметров 0, а0/Л,0, взятых из области применимости аналитического решения (30), представлены тёмными маркерами. Расчёты показали, что с увеличением амплитуды набегающей волны заплеск возрастает, при этом с использованием метода адаптивных сеток получаются большие значения максимального вертикального заплеска, чем в случаях применения методов сквозного счёта для классических (3) и модифицированных (23), (25) уравнений мелкой воды.
Другая серия экспериментов, выполненная при постоянном значении амплитуды начальной волны а0/Л,0 = 0.01 и разных значениях угла наклона откоса (0 = 1° ^ 10°), показала уменьшение заплесков при возрастании угла 0, при этом методы сквозного
а б
Рис. 2. Накат уединённой волны на плоский откос. а — нижняя (1) и верхняя (2) границы области, в которой справедлива формула (30), определяемые соотношениями (31), (33); б — значения максимальных заплесков, полученные по формуле (30) (сплошная линия), в экспериментах [43, 47] (+), с помощью численных методов сквозного счёта (а, а), адаптивных сеток (•, о) и схемы сквозного счёта для модифицированных уравнений мелкой воды (■, □)
счёта демонстрируют меньшую зависимость максимального заплеска от параметра 0, чем метод адаптивных сеток.
Эксперименты для высоких волн (а0/Л0 = 0.1 ^ 0.256) и крутых склонов (0 = 10° ^ 15°, г0 = 1.3) подтвердили, что заплески, полученные с использованием метода адаптивных сеток, всегда превышают значения заплесков, вычисленные двумя другими методами.
Сравнение численных результатов с аналитическим решением в той области параметров, для которых справедлива формула (30), показывает, что ни один из рассмотренных методов не имеет явного преимущества и все численные результаты неплохо согласуются с точным решением. Более того, если, например, в методе адаптивных сеток для каждой амплитуды и каждого угла наклона откоса подбирать оптимальные значения границ т и М, то можно получать значения заплесков, практически совпадающие с аналитическими. Такой подбор здесь не проводился, и во всех расчётах принимались одни и те же значения т = 10-3 и М = 10. Аналогичное замечание относится к методам сквозного счёта, при использовании которых величины максимальных заплесков сильно зависят соответственно от значений £я и 7. Эти параметры также не подбирались и во всех расчётах были постоянными: £я = 0.01а0/Л0, 7 = 0.02.
Анализ результатов расчётов (см. светлые маркеры на рис. 2, б) для углов наклона откоса и амплитуд набегающей волны, не принадлежащих области (31), (33), показывает, что метод адаптивных сеток завышает значения вертикального заплеска по сравнению с методом сквозного счёта, а схема сквозного счёта для модифицированных уравнений мелкой воды — занижает. С другой стороны, результаты, полученные с использованием модифицированных уравнений мелкой воды, в среднем наиболее близки к экспериментальным данным. Этот факт не является неожиданным, поскольку в уравнении для модифицированного импульса (25) в косвенной форме (через эмпирический коэффициент 7) учитывается влияние на течение турбулентной вязкости в поверхностном слое на фронте размазанной прерывной волны. Отметим, что уже в работе [29] подчёркивалось, что численные методы, основанные на классических уравнениях волновой гидродинамики для идеальной жидкости, всегда дают превышение максимальных заплесков над экспериментальными данными, и предпринималась попытка связать этот эффект с тем, что в используемых математических моделях вязкость жидкости не учитывалась.
2.2. Накат на криволинейный откос
Взаимодействие волн с криволинейным береговым склоном имеет более сложный характер, чем с плоским откосом, и аналитическими методами изучено недостаточно полно, поэтому для детального моделирования процессов наката-отката в этом случае необходимо использование методов численного моделирования [57]. Здесь будет рассмотрена задача о набегании уединённой волны на криволинейный откос, сопрягающийся с горизонтальным дном глубины Л0 = 1:
0 + ( —---^ ) Л [с(х — £)] при 0 ^ х ^ Хз
—Л(х) = { 2 ^ 2 7 1 4 ^ 1 (34)
— 1 при хз ^ х ^ Ь,
где ¿г > 0 — высота суши в левой бесконечно удалённой точке, ¿о — высота суши в точке х = 0 (0 < х0 < ¿г),
tg 0о . 1 , 1 + ¿г
" С = т"1п
¿г - (¿о - 1)/2' 2с ¿г - ¿о'
х* = 2С, 0о — максимальный угол наклона неровного склона, достигаемый в точке перегиба С- Согласно уравнению (34), в начальный момент времени точка уреза имеет абсциссу
* 1 1 Л ^ 1 - ¿о
Хоо = С - 7Г 1п 1 +
2с \ ¿г
при этом 0 < хоо < С- В расчётах использовались значения параметров ¿г = 0.15, ¿о = 0.14.
На рис. 3 изображены графики функции (34) для двух значений угла 0о, а в табл. 1 приведены геометрические характеристики криволинейных профилей дна, при этом 0оо означает угол наклона откоса в начальной точке уреза хоо, а в3 — то же, в точке сопряжения х* криволинейного откоса с горизонтальным участком. Видно, что при 0о < 11° локальные углы наклона профиля в окрестности начальной точки уреза не превышают 5° и стыковка участков дна в точке сопряжения (0* < 0.38°) достаточно гладкая.
Результаты расчёта наката на откос (34) сравнивались с данными, полученными для плоского откоса, аппроксимирующего криволинейный и проходящего (см. штри-
г \ ' 0 . 4
-1.0
0 10
20
30 х/Нп
Рис. 3. Профили дна и суши (сплошные линии) в задаче о накате уединённой волны на неровный откос (34) при ео = 5° (1) и 9° (2)
Ь
z
0
Таблица 1. Зависимость формы дна от заданного в точке перегиба £ угла наклона во
до Хоо воо £ х§ в* хт Ь в
3 15.70 1.35 26.26 52.51 0.10 72.51 92.51 1.56
5 9.41 2.26 15.73 31.46 0.17 51.46 71.46 2.60
7 6.70 3.16 11.21 22.41 0.24 42.41 62.41 3.64
9 5.20 4.08 8.69 17.38 0.31 37.38 57.38 4.69
11 4.23 5.00 7.08 14.16 0.38 34.16 54.16 5.75
ховые линии на рис. 3) через точку (хз, —1) сопряжения криволинейного участка дна с горизонтальным и начальную точку уреза (х00, 0):
х — х00
при 0 ^ х ^ хз
— Л(х) = ^ х5 — х00 (35)
— 1 при хз ^ х ^ Ь.
Значения угла наклона плоского откоса 0 приведены в табл. 1.
Начальное возмущение задавалось в виде уединённой волны (27), при этом абсциссы вершины волны и правой границы определялись по формулам (28). В табл. 2 приведены значения величины Я/а0, характеризующей степень усиления высоты заплеска по отношению к амплитуде набегающей волны. Видно, что как на криволинейном, так и на плоском склоне значения коэффициента усиления уменьшаются с ростом угла наклона 00 и возрастают при увеличении амплитуды а0. Для пологого криволинейного откоса (34) заплески оказываются несколько большими, чем для аппроксимирующего его плоского откоса. Хотя при малых 00 отличия по высотам незначительны, тем не менее дальности заплеска могут различаться в разы, что связано с тем, что в отличие от плоского откоса в случае криволинейного (34) волна продвигается по склону, выпо-лаживающемуся по мере удаления от берега, и поэтому проходит по суше боольшее расстояние (вдоль оси Ох), чем при плоском откосе. Сильно различаются также картины взаимодействия волны с берегом. На рис. 4 показана динамика свободной поверхности при накате уединённой волны на криволинейный откос и плоский склон. Видно, что в том и в другом случае скорость волны с приближением к берегу замедляется. После наката волны на плоский откос она скатывается с него с образованием основной отражённой волны. Для криволинейного же склона отражённая волна возникает в момент времени, когда на суше ещё продолжается накат, и во время пребывания воды на пологом криволинейном склоне появляются несколько отражённых волн (см. рис. 4, а), чего не наблюдалось для плоских откосов.
Интересная особенность взаимодействия уединённой волны с пологими откосами — колебательный характер процесса наката-отката (см. рис. 4) с затухающей амплитудой колебаний. Для криволинейного откоса наблюдаются колебания большей амплитуды, боольшего периода и медленнее затухающие, чем для плоского. С увеличением крутизны откосов колебательный характер движения точки уреза практически исчезает и в процессе взаимодействия волны с берегом образуется только одна отражённая волна. Однако на криволинейном склоне колебания точки уреза наблюдаются для боольших значений углов 00, чем на аппроксимирующем его плоском откосе.
Таблица 2. Зависимость максимального вертикального заплеска от наклона дна и амплитуды набегающей волны. Метод адаптивных сеток
00 К/а0
неровный склон (34) линейный склон (35)
а0/Н0 = 0.01 а0/Н0 = 0.02 а0/Н0 = 0.01 а0/Н0 = 0.02
3 5.62 6.34 5.58 6.30
5 4.40 5.15 4.37 5.11
7 3.71 4.36 3.74 4.42
9 3.26 3.90 3.35 3.94
11 2.95 3.52 3.07 3.60
Рис. 4. График поверхности г = п(х, £) при накате уединённой волны на криволинейный откос (а), на плоский откос (б); ао/Ло = 0.02, = 5°, метод адаптивных сеток
Таблица 3. Зависимость максимального вертикального заплеска от наклона дна и амплитуды набегающей волны. Метод сквозного счёта для классических (3) и модифицированных (23), (25) уравнений мелкой воды
Я/ао
0о неровный склон (34) линейный склон (35)
ао/Ло = 0.01 ао/Ло = 0.02 ао/Ло = 0.01 ао /Ло = 0.02
Метод еквозного счёта
3 4.83 5.02 4.63 4.79
5 4.08 4.49 4.11 4.27
7 3.62 4.00 3.70 3.86
9 3.28 3.54 3.38 3.60
11 3.02 3.36 3.23 3.46
Схема для модифицированных уравнений мелкой воды (23), (25)
3 4.90 4.50 4.97 4.80
5 4.18 4.32 4.20 4.55
7 3.59 4.14 3.63 4.18
9 3.16 3.82 3.28 3.82
11 2.88 3.53 3.02 3.51
В табл. 3 приведены результаты расчётов с использованием методов сквозного счёта для классических и модифицированных уравнений мелкой воды. Видно, что в данных случаях для пологих откосов величины заплесков получаются заметно меньшими, чем при использовании метода адаптивных сеток. Однако для более крутых склонов результаты вполне сопоставимы. Интересной особенностью является и то, что методы сквозного счёта дают меньшую зависимость высоты заплеска от крутизны склона, чем метод адаптивных сеток.
3. Результаты моделирования наката на участок реального побережья
В настоящем разделе приводятся некоторые результаты моделирования наката на берег волн цунами, вызванных сильнейшим землетрясением (М = 9.0) 11 марта 2011 г., эпицентр которого находился в точке с координатами 38.9° N, 141.7° Е. Высота заплес-ков на некоторых участках восточного побережья Японии достигала 37 м. Отмечалась дальность заливания до 10 км.
Для моделирования генерации цунами и распространения волн от области источника к берегу была выбрана область, простирающаяся с запада на восток от 135°42'Е до 149°00'Е и с юга на север от 33°00^ до 44°30'N. Двумерное моделирование выполнялось в подобласти, отделённой от берега вертикальной стенкой, установленной на глубине Л.№ап = 30 м. Использовалась равномерная сетка с шагом 15 географических секунд и размерностью 3193 х 2761 узлов.
Для расчёта характеристик наката волны цунами и определения области затопленного волной побережья авторами были выбраны (см. рис. 1, а) две расположенные недалеко друг от друга подобласти ("южная" и "северная") с достаточно гладкой береговой линией. В каждой из них построено по семь сечений: 51 ^ Б7 и N1 ^ N7 . Профили склонов вдоль сечений в своих подводных фрагментах в основном близки к линейным. Что касается рельефа суши, то можно выделить две группы сечений, в первой из которых профили суши почти линейны и не содержат особенностей, а во второй — в своих прибрежных участках содержат пологие фрагменты (рис. 5), что может приводить, как отмечалось в разделе 2.2, к резкому отличию динамики затопления.
Для каждого сечения задавалась входящая через его мористое окончание волна = рассчитанная с помощью 2В-моделирования. Для сечений Б6 и N4 ма-
реограммы входящей волны изображены на рис. 6 линиями 1. Для этих сечений амплитуды головной входящей волны равны примерно 2.3 и 2.6 м соответственно. На основе 2В-моделирования показано, что на вертикальной стенке, установленной на глубине hwa.il = 30 м, эти волны усилились почти в три раза и привели к максимальным заплес-кам в 6 и 7.2 м (см. мареограммы г(¿) = п(х№а11,^) на рис. 6, изображённые линиями 2). Линии 3 на рис. 6 представляют зависимость вертикального смещения точки уреза от времени г(¿) = п(х0(£),£), рассчитанную методом адаптивных сеток. Для сечений Б6 и N4 максимальный вертикальный заплеск волн на берег составил соответственно 12.64 и 14.7 м, т. е. вдвое превысил заплеск на вертикальную стенку и примерно в 5.5 раза — амплитуду входящей волны.
М
-80
-40
0
0
5000
10000
15000 м 20000
Рис. 5. Профили рельефов дна и прилегающей суши в сечениях Б6 (1) и N4 (2)
Хотя максимальные вертикальные заплески в сечениях Б6 и N4 различаются примерно на 10%, горизонтальные заплески (дальности заливания) различаются существенно — более чем в полтора раза — и составляют соответственно 2268 и 3439 м. Причина последнего состоит в том, что в сечении N4 имеется более протяжённый по сравнению с Б6 пологий участок прибрежной суши. Так, в сечении N4 суша имеет высоту 5 м на удалении более 2 км от точки уреза спокойной воды, а в Б6 — около 1 км. Соответственно и время нахождения суши под водой в процессе наката-отката
Рис. 6. Вертикальное (3) и горизонтальное (4) смещение точки уреза в сечениях Б6 (а) и N4 (б) при входе волны (1) через мористое окончание сечений (изобата 100 м); 2 — колебания свободной поверхности на вертикальной стенке, рассчитанные по двумерной модели
шщ 10.9 II
12.3*— 12.
9 Я 8.54 Г )
НО. 7 140.8
37.2
N
37.1 3736.9-
Зб.а-
Щ 1489^ I
4221*— 327226^
814*/ 402ч Г )
1-10.9 141 Е 141.1 140.7 140.В 1-10.9
аб
141 Е 141.1
Рис. 7. Зона затопленного берега. Числа на границе зоны затопления указывают максимальный вертикальный заплеск (а) и максимальную дальность затопления (б) в метрах; N — северная широта, Е — восточная долгота
первой волны составляет примерно 29 мин для сечения N4 и 16.6 мин для Б6. Горизонтальное смещение точки уреза показано на рис. 6 линиями 4, которые представляют собой графики функции /(¿) = хоо — хо(£), принимающей положительные значения при выходе волны на берег и отрицательные — при осушении дна. Соответствующие оси изображены справа от графиков.
Наличие весьма протяжённого и очень пологого участка в сечении N4 является, видимо, причиной того, что в фазе отката первой волны присутствует эффект повторного наката, после которого первая волна покидает сушу с осушением дна до глубины —10.7 м и отходом кромки воды на расстояние 1180 м от начальной точки уреза.
На рис. 7 приведена криволинейная граница зоны затопления, определённая по результатам расчётов во всех сечениях на основе методики, описанной в начале раздела 1. Серый цвет — зона акватории, белый — суши, не залитой водой, тёмный — зона затопления. Маркеры на границе зоны затопления указывают точки, до которых дошла волна вдоль выбранных сечений. Видно, что для рассматриваемого участка побережья высота заплесков в зависимости от рельефа суши колеблется в пределах от 8 до 14. 7 м, а дальность заливания — от 377 до 4221 м. При этом по высотам заплесков имеет место не только качественное соответствие численных результатов данным наблюдений, но и хорошее количественное согласование. Так, вблизи сечения Б6 есть точка, в которой наблюдалась высота заплеска в 11.2 м, а расчёт для этого сечения дает значение 12.64 м.
Заключение
В статье представлен комбинированный метод численного моделирования наката волн цунами на побережье, основанный на использовании двумерной модели мелкой воды для расчёта распространения волны от источника землетрясения к побережью и одномерного моделирования наката вдоль сечений, проведённых от некоторой прибрежной изобаты до выбранной изолинии рельефа на суше. При этом результаты двумерного моделирования используются в качестве входных данных для одномерных расчётов. Выполнено сравнение трёх численных методов для расчёта наката: сквозного счёта и с выделением линии разрыва на основе классической модели мелкой воды и конечно-разностного метода, основанного на модели мелкой воды с модифицированным уравнением импульса. Результаты тестовых расчётов показали, что метод адаптивных сеток можно использовать для определения границы затопления с "запасом", т. е. он дает границу, до которой волна может дойти. Метод для модифицированных уравнений мелкой воды позволяет дать оценку снизу для границы затопления, т. е. определить линию на суше, которую волна при накате пройдёт. И, наконец, рассмотренный вариант метода сквозного счёта может использоваться в качестве дополнительного контрольного механизма.
Изложенная методика применялась для расчёта зон затопления многих участков побережья и показала удовлетворительное согласование численных результатов с данными натурных наблюдений. При этом надо иметь в виду, что методика может использоваться только для участков с простой формой береговой линии и расположенных "напротив" источника. Ввиду ограниченности объёма настоящей публикации, другая часть работы, посвящённая двумерной методике расчёта наката, позволяющей выполнять моделирование в условиях изрезанной береговой линии, наличия разномасштабных заливов и бухт, косого или касательного подхода волны, будет изложена в отдельной статье.
Список литературы
[1] Petukhin A., Yoshida K., Miyakoshi K., Irikura K. Tsunami simulation for the 2011 Great Tohoku earthquake (Mw9.0), Japan, using seismic inversion source model and fully nonlinear tsunami model // 15th World Conference on Earthquake Engineering. Lisbon, Portugal, 2012. http://www.iitk.ac.in/nicee/wcee/article/WCEE2012_3644.pdf
[2] Гусяков В.К. Остаточные смещения на поверхности упругого полупространства // Условно-корректные задачи математической физики в интерпретации геофизических наблюдений. Новосибирск: ВЦ СО АН СССР, 1978. C. 23-51.
[3] Стокер Дж.Дж. Волны на воде. М.: Иностранная литература, 1959.
[4] Куликовский А.Г., Погорелов Н.В., Семёнов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001.
[5] Остапенко В.В. Гиперболические системы законов сохранения и их приложение к теории мелкой воды. Новосибирск: Изд-во НГУ, 2004.
[6] Whitham G.B. The effects of hydraulic resistance in the dambreak problem // Proceed. Royal Society. Ser. A. 1955. Vol. 227, No. 1170. P. 399-407.
[7] Судобичер В.Г., Шугрин С.М. Движение потока воды по сухому руслу // Изв. СО АН СССР. Сер. техн. наук. 1968. Т. 13, вып. 3. С. 116-122.
[8] Martin J.C., Moyce W.J. An experimental study of the collapse of liquid columns on a rigid horizontal plane // Phil. Trans. Roy. Soc. London. 1952. Vol. 244A. P. 312-324.
[9] Dressler R.F. Comparison of theories and experiments for the hydraulic dam-break wave // Intern. Assoc. Sci. Hydrology. 1954. Vol. 3, No. 38. P. 319-328.
[10] Stansby P.K., Chegini A., Barnes T.C. The initial stages of dam-break flow // J. Fluid Mech. 1998. Vol. 374. P. 407-424.
[11] Букреев В.И., Гусев А.В., Малышева А.А., Малышева И.А. Экспериментальная проверка газогидравлической аналогии на примере задачи о разрушении плотины // Изв. РАН. Механика жидкости и газа. 2004. № 5. С. 143-152.
[12] Букреев В.И., Гусев А.В. Начальная стадия генерации волн при разрушении плотины // Докл. АН. 2005. Т. 401, № 5. С. 1-4.
[13] ВольцингЕр Н.Е., Клеванный К.А., ПЕлиновский Е.Н. Длинноволновая динамика прибрежной зоны. Л.: Гидрометеоиздат, 1989.
[14] Medeiros S.C., Hagen S.C. Review of wetting and drying algorithms for numerical tidal flow models // Intern. J. Numer. Meth. Fluids. 2013. Vol. 71. P. 473-487.
[15] Hibberd S., Peregrine D.H. Surf and runup on a beach: A uniform bore //J. Fluid Mech. 1979. Vol. 95, pt 2. P. 323-345.
[16] Федотова З.И. Обоснование численного метода для моделирования наката длинных волн на берег // Вычисл. технологии. 2002. Т. 7, № 5. С. 58-76.
[17] Kobayashi N., DeSilva G.S., Watson K.D. Wave transformation and swash oscillation on gentle and steep slopes //J. Geophys. Res. 1989. Vol. 94, No. C1. P. 951-966.
[18] Борисова Н.М., Остапенко В.В. О численном моделировании процесса распространения прерывных волн по сухому руслу // Журнал вычисл. математики и матем. физики. 2006. Т. 46, № 7. С. 1322-1344.
[19] Bradford S.F., Sanders B.F. Finite-volume model for shallow-water flooding of arbitrary topography //J. Hydraul. Eng. 2002. Vol. 128. P. 289-298.
[20] Delis A.I., Kazolea M., Kampanis N.A. A robust high-resolution finite volume scheme for the simulation of long waves over complex domains // Intern. J. Numer. Meth. Fluids. 2008. Vol. 56. P. 419-452.
[21] Bokhove O. Flooding and drying in discontinuous Galerkin finite-element discretizations of shallow-water equations. Part 1: One dimension //J. Sci. Comput. 2005. Vol. 22-23. P. 47-82.
[22] Westerink J.J., Feyen J.C., Atkinson J.H. et al. A basin to channel-scale unstructured grid hurricane storm surge model applied to Southern Louisiana // Monthly Weather Rev. 2008. Vol. 136, No. 3. P. 833-864.
[23] Ern A., Piperno S., Djadel K. A well-balanced Runge —Kutta discontinuous Galerkin method for the shallow-water equations with flooding and drying // Intern. J. Numer. Meth. Fluids. 2008. Vol. 58. P. 1-25.
[24] Bunya S., Kubatko E.J., Westerink J.J., Dawson C. A wetting and drying treatment for the Runge-Kutta discontinuous Galerkin solution to the shallow water equations // Comput. Meth. Appl. Mech. Eng. 2009. Vol. 198, No. 17-20. P. 1548-1562.
[25] KarnA T., de Brye B., Gourgue O. et al. A fully implicit wetting-drying method for DG-FEM shallow water models, with an application to the Scheldt Estuary // Ibid. 2011. Vol. 200, No. 5-8. P. 509-524.
[26] Лятхер В.М., МилитЕЕВ А.Н. Расчёт наката длинных гравитационных волн на откос // Океанология. 1974. Т. 14, № 1. С. 37-42.
[27] Yeh G.-T., Chou F.-K. Moving boundary numerical surge model //J. Waterway, Port, Coastal, Ocean Division. 1979. Vol. 105. P. 247-263.
[28] Lynch D.R., Gray W.G. Finite element simulation of flow in deforming regions //J. Comput. Phys. 1980. Vol. 36. P. 135-153.
[29] Pedersen G., Gjevik B. Run-up of solitary waves // J. Fluid Mech. 1983. Vol. 135. P. 283-299.
[30] Liang Q., Borthwick A.G.L. Adaptive quadtree simulation of shallow water flows with wet-dry fronts over complex topography // Comput. Fluids. 2009. Vol. 38. P. 221-234.
[31] Madsen P.A., Fuhrman D.R. Run-up of tsunamis and long waves in terms of surf-similarity // Coastal Eng. 2008. Vol. 55. P. 209-223.
[32] Bautin S.P., Deryabin S.L., Sommer A.F. et al. Use of analytic solutions in the statement of difference boundary conditions on a movable shoreline // Rus. J. Numer. Anal. Math. Model. 2011. Vol. 26, No. 4. P. 353-377.
[33] Борисова Н.М., Гусев А.В., Остапенко В.В. О распространении прерывных волн по сухому руслу // Изв. РАН. Механика жидкости и газа. 2006. № 4. С. 135-148.
[34] Остапенко В.В. Модифицированные уравнения теории мелкой воды, допускающие распространение прерывных волн по сухому руслу // ПМТФ. 2007. Т. 48, № 6. С. 22-43.
[35] Борисова Н. М. О моделировании процесса набегания прерывной волны на наклонный берег // СибЖВМ. 2007. Т. 10, № 1. С. 43-60.
[36] Gusyakov V.K., Fedotova Z.I., Khakimzyanov G.S. et al. Some approaches to local modelling of tsunami wave runup on a coast // Rus. J. Numer. Anal. Math. Model. 2008. Vol. 23, No. 6. P. 551-565.
[37] Flather R.A. Existing operational oceanography // Coastal Eng. 2000. Vol. 41. P. 13-40.
[38] TsuNAMis-Assessing the Hazard for the UK and Irish Coasts. London: Department for Environment, Food and Rural Affairs, 2006. http://www.defra.gov.uk
[39] Choi B.H., Kaistrenko V., Kim K.O. et al. Rapid forecasting of tsunami runup heights from 2-D numerical simulations // Natural Hazards Earth Syst. Sci. 2011. Vol. 11. P. 707-714.
[40] Carrier G.F., Greenspan H.P. Water waves of finite amplitude on a sloping beach // J. Fluid Mech. 1958. Vol. 4, No. 1. P. 97-109.
[41] Carrier G.F., Wu T.T., Yeh H. Tsunami run-up and draw-down on a plane beach // Ibid. 2003. Vol. 475. P. 79-99.
[42] Kanoglu U. Nonlinear evolution and runup-rundown of long waves over a sloping beach // Ibid. 2004. Vol. 513. P. 363-372.
[43] Synolakis C.E. The runup of solitary waves // Ibid. 1987. Vol. 185. P. 523-545.
[44] Мазова Р.Х., ПЕлиновский Е.Н. Линейная теория наката волн цунами на берег // Изв. АН СССР. Физика атмосферы и океана. 1982. Т. 18, № 2. С. 166-171.
[45] Кайстренко В.М., ПЕлиновский Е.Н., Симонов К.В. Накат и трансформация волн цунами на мелководье // Метеорология и гидрология. 1985. № 10. C. 68-75.
[46] Pelinovsky E.N., Mazova R.Kh. Exact analytical solutions of nonlinear problems of tsunami wave run-up on slopes with different profiles // Natural Hazards. 1992. Vol. 6, No. 3. P. 227-249.
[47] Synolakis C.E. Tsunami runup on steep slopes: How good linear theory really is // Ibid. 1991. Vol. 4, No. 2-3. P. 221-234.
[48] Диденкулова И.И., Заибо Н., Куркин А.А. и др. Накат нелинейно деформированных волн на берег // Докл. АН. 2006. Т. 410, № 5. C. 676-678.
[49] Диденкулова И.И., Куркин А.А., ПЕлиновский Е.Н. Накат одиночных волн различной формы на берег // Изв. РАН. Физика атмосферы и океана. 2007. Т. 43, № 3. С. 419-425.
[50] Диденкулова И.И., ПЕлиновский Е.Н. Накат длинных волн на берег: влияние формы подходящей волны // Океанология. 2008. Т. 48, № 1. С. 5-10.
[51] Synolakis C.E., Bernard E.N., Titov V.V. et al. Validation and verification of tsunami numerical models // Pure and Appl. Geoph. 2008. Vol. 165. P. 2197-2228.
[52] Shokin Yu.I., Babailov V.V., Beisel S.A. et al. Mathematical modeling in application to regional tsunami warning systems operations // Computational Sci. and High Performance Computing III. Notes on Numerical Fluid Mechanics and Multidisciplinary Design. 2007. Vol. 101. P. 52-69.
[53] ХАкимзянов Г.С., ШокинА Н.Ю. Метод адаптивных сеток для одномерных уравнений мелкой воды // Вычисл. технологии. 2013. Т. 18, № 3. С. 54-79.
[54] Khan A.A. Modeling flow over an initially dry bed // J. of Hydraulic Res. 2000. Vol. 38, No. 5. P. 383-388.
[55] Гусев О.И., ШокинА Н.Ю., Кутергин В.А., ХАкимзянов Г.С. Моделирование поверхностных волн, генерируемых подводным оползнем в водохранилище // Вычисл. технологии. 2013. Т. 18, № 5. С. 74-90.
[56] ХАкимзянов Г.С., Шокин Ю.И., Барахнин В.Б., ШокинА Н.Ю. Численное моделирование течений жидкости с поверхностными волнами. Новосибирск: Изд-во СО РАН, 2001.
[57] Yeh H., Liu P., Synolakis C.E. Long-Wave Runup Models. Singapore: World Sci. Publ., 1996.
Поступила в редакцию 4 декабря 2013 г.
