CC BY
30
Вычислительные технологии
Том 20, № 5, 2015
О численных методах решения задач о накате волн на берег. I. Сравнительный анализ численных алгоритмов для одномерных задач
Ю.И. Шокин1, А. Д. Рычков1, Г. С. Хлкимзянов1'2, Л. Б. Чубаров1'2'*
1Институт вычислительных технологий СО РАН, Новосибирск, Россия 2 Новосибирский государственный университет, Россия *Контактный e-mail: chubarovster@gmail.com
Проведен сравнительный анализ метода крупных частиц, метода "годуновско-го" типа, метода сглаженных частиц, комбинированного метода TVD + SPH, а также метода адаптивных сеток для численного моделирования наката длинных волн на берег в рамках теории мелкой воды. Значительное внимание уделено проверке качества алгоритмов и математической модели путем их сопоставления с аналитическими решениями и результатами лабораторных экспериментов. Показаны определенные преимущества метода TVD + SPH для решения задач наката волн на берег по сравнению с другими методами. Продемонстрированы возможности успешного применения этого метода для моделирования как слабо нелинейных (необрушающихся), так и существенно нелинейных (обрушающихся) волн. В последнем случае качество результата в значительной степени зависит от удачного подбора параметров, определяющих величину донного трения.
Ключевые слова: поверхностные волны, мелкая вода, накат на берег, численное моделирование.
Введение
Для моделирования наката волн на побережья с реальными очертаниями береговых линий и рельефами прилегающих участков акваторий и суши в рамках теории мелкой воды применяются различные численные методы. Их условно можно разделить на три основных класса: бессеточные методы сглаженных частиц (SPH-методы) [1-3], методы конечных разностей [4-6] и конечных элементов, а также комбинированные методы типа методов частиц в ячейках [7-9].
Обзор некоторых методов, позволяющих вполне адекватно описывать основные характеристики наката длинных поверхностных волн на берег в одномерной постановке, приведен в монографии [10]. Все упомянутые выше алгоритмы характеризуются значительным разбросом по точности, требуемым ресурсам, алгоритмической сложности. Очень важным фактором для рассматриваемого класса задач является способность того или иного численного метода определять такие специфические характеристики, как границы зон затопления, вторичные колебания линии уреза и т. п. Выбор численного метода, наиболее подходящего для решения прикладных задач о накате длинных
© ИВТ СО РАН, 2015
поверхностных волн на берег, требует тщательного тестирования рассматриваемой совокупности вычислительных схем.
Первая часть настоящей работы посвящена сравнительному анализу алгоритмов на примере одномерной задачи о накате одиночных волн на плоский откос. Численные и аналитические результаты решения этой классической тестовой задачи широко представлены в литературе. В статье приведена постановка одномерной задачи, дается краткое описание рассматриваемых численных алгоритмов и излагаются результаты некоторых тестовых расчетов. Вторая часть будет посвящена решению с помощью комбинированного метода ТУВ + БРИ модельной задачи о заплеске на плоский откос волн, распространявшихся до этого над участком дна постоянной глубины. Будут представлены результаты анализа особенностей заплесков волн различных конфигураций и амплитуд на берега различной крутизны. Особое внимание будет уделено исследованию процессов заплеска так называемых N-волн, имеющих самое непосредственное отношение к проявлениям у берега волн цунами, порожденных подводными землетрясениями, очаги которых расположены на небольшом расстоянии от берега. В третьей части на модельных примерах предполагается рассмотреть различные методики определения ха-рактертистик взаимодействия волн с берегом, применяемые для решении прикладных задач проблемы цунами.
1. Постановка задачи
В рамках модели мелкой воды рассматривается одномерная задача о накате волны на плоский откос, сопряженный с дном постоянной глубины к0. Используется декартова система координат Охх с вертикальной осью Ох, направленной вверх, и с координатной линией г = 0, совпадающей с невозмущенной свободной поверхностью слоя идеальной несжимаемой жидкости, ограниченного снизу неподвижным дном г = -Н(х), а сверху — подвижной свободной границей г = гц(х,Ь), где Ь — время. Система уравнений мелкой воды в этом случае записывается в виде законов сохранения:
искомыми величинами являются полная глубина слоя жидкости Н(х,Ь) = г/(х,1) + К(х) > 0 и и(х,Ь) — осредненная по вертикали от дна до свободной поверхности скорость; д — ускорение силы тяжести; х = Ьх — правая неподвижная граница области решения, имеющей слева подвижную границу х0(Ь), отделяющую воду от суши (точка уреза). Таким образом, в отличие от задач для уравнений мелкой воды в областях с фиксированными неподвижными границами, в задачах наката волн на берег одна из границ — линия уреза — является подвижной, заранее неизвестной и подлежащей определению. Поэтому в рассматриваемых задачах количество искомых функций увеличивается: наряду с поиском Н и и требуется искать дополнительную неизвестную х0(Ь).
Уравнения (1) дополняются краевыми
+ Еж = С, х0(Ь) <х<Ьх, Ь> 0,
(1)
где
н(хо(г),г) = 0, и(ьх,г) = 0, г > о,
и начальными
г](х, 0) = щ(х), и(х, 0) = и0(х), х00 ^х ^ Ьх (3)
условиями, где хоо = хо(0) — известная координата точки уреза в начальный момент времени, при этом предполагается, что волна движется справа налево.
В рассматриваемых ниже тестовых задачах рельеф дна и прилегающей суши задается функцией
7 / ч { - х tg/5 при 0 ^х ^ х3, г = -Н(х) = < , т (4)
I —м0 при хт ^ х ^ ьх,
где ß > 0 — угол наклона откоса; z0 > 0 — высота суши в точке х = 0; xs (z0 + h0) ctg ß — абсцисса основания склона. Значение z0 подбиралось таким, чтобы максимальный вертикальный заплеск волны на склон (ордината точки уреза на склоне) не превышал z0. Таким образом, при всех t > 0 будет выполняться неравенство x0(t) > 0.
2. Вычислительные алгоритмы
Для численного решения системы уравнений (1) в области х Е [0, Ьх] строилась равномерная разностная сетка хг = г Ах (г = 0, ... , N) с шагом Ах = Ьх/М, на которой система (1) решалась с помощью той или иной разностной схемы только в тех узлах сетки, в которых для значений полной глубины на п-м шаге по времени выполнялось неравенство Н™ > е, где е > 0 — заданная пороговая величина. В остальных узлах сетки все параметры течения полагались равными нулю. Далее, на каждом новом (п + 1)-м временном слое производился аналогичный анализ вычисленных значений полной глубины Н™"+1; если в некотором узле Xf она оказывалась меньше заданного порогового значения, то значения всех параметров течения левее этого узла (хг < х^) также полагались равными нулю. Такой подход известен в литературе как "метод улавливания скачка" [4] (в данном случае точки уреза х0 (¿)) в узел разностной сетки. В самой точке уреза расчет параметров течения на (п + 1)-м временном слое проводился с использованием в разностных схемах односторонних разностей по пространственной переменной, что обеспечивало возможность перемещения этой точки во времени, т. е. приближенное значение положения точки уреза на (п + 1)-м временном слое полагалось равным Xf.
2.1. Метод крупных частиц
Реализация метода крупных частиц проводится в два этапа. На первом (эйлеровом) этапе в системе (1) отбрасываются все конвективные члены и она записывается в виде
да ' 0
Ж + D = 0' D
дх
)
а на втором (лагранжевом) этапе решаются уравнения, полученные из (1), в которых сохранены только конвективные члены:
д q д F f
д д х
(Hu \ ^ Hu2 )
Для решения уравнений (5), (6) строится явная разностная схема первого порядка точности на разнесенной сетке, что обеспечивает, как показали расчеты, отсутствие
осцилляций численного решения. В основных узлах определяется величина Н, в сдвинутых вычисляется поток (Ни). Разностная схема первого этапа записывается следующим образом:
нг = н?, (Нй)г+1/2 = (Ниу+1/2 - дН?+1/2 ^ [(н - Ю+ - (н - ь)*г
где целый индекс относится к узлам основной сетки, дробный — к сдвинутым. На втором этапе для уравнений (6) используется противопоточная разностная схема:
^ = Н - A
Qi+1/2 - Чг-1/2
!
HiUi+i/2, если щ+1/2 > 0,
Qi+1/2 ]
Hi+iui+i/2, если Ui+1/2 < 0,
йг+1/2 = 2(Ни)г+1/2/(Нг+1 + Щ). Аналогичным образом записываются выражения для вычисления значений потоков
п+1 :
г+1/2:
( Ни)п+1
1 At г ~ ~ 1 { иг(Ни)г+1/2, если и > 0,
( Ни)^/2 = (Ни)г+1/2 - — Qi+1 - Qi , Qi+1 = I „
У Ui(HU)i+3/2, если и < 0,
Ui = (Ui+1/2 + Ui+3/2 )/2. Приближенное определение координаты точки уреза на (п + 1)-м временном слое
п+1 и
x0 осуществляется с помощью метода улавливания точки уреза в узел разностной сетки.
К достоинствам метода относятся его устойчивость к неточности задания рельефа дна, точное выполнение условия "lake at rest" (общепризнанный вычислительный тест "озеро в покое", проверяющий сохранение состояния покоя при задании нулевых начальных значений r/(x, 0) = 0, u(x, 0) = 0 во все последующие моменты времени при любом рельефе дна), отсутствие нефизических осцилляций численного решения, простота программирования и нетребовательность к вычислительным ресурсам.
Основным недостатком метода является его первый порядок точности. Однако при решении двумерных задач этот недостаток частично компенсируется тем, что приходится использовать достаточно подробную прямоугольную разностную сетку в силу необходимости детальной аппроксимации реальной береговой линии, а также рельефов суши и дна акватории.
2.2. Алгоритм типа схемы Годунова с TVD-свойством (HLL-метод)
Явную схему второго порядка точности, приведенную в [4], можно записать для уравнений мелкой воды (1) в следующем виде:
'П+1 _ П
q =q Ax
F(q? + 0.5Ax Q?) - F(q? - 0.5Ax Q?)
qn+1 + q n+1/2 qi + qi
qi =
, (-ЛП+1/2
+ tG, ,
2
qn+1 - qn T (Fn+1/2 Fra+1/^W ^Gn+1/2 qi = qi - AxlFi+1/2 - *-1/2) + TGi ■
Здесь
ОН = шт
/ап Ш0ё | —
Я? Я? - Я?- -
V Ах
Ах
О
п+-/2
(
0
дн1
+ -/2 Иг+-/2 — <Ч--/2
АХ
- 1
потоки F™+-:-///22 определяются на гранях ячеек из приближенного решения задачи Римана с помощью метода Хартена — Лакса — ван Лира (ИЬЬ):
рга+-/2 г+-/2
¥1
хк¥ь - хь¥к + хьхк(я« - яь)
Xя
Хь > 0,
Хь < 0 < X
к
Ел, Хк < 0,
где ¥ь = ¥(яь), ¥к = Е(яд), ЯЬ = яГ-/2, ЯК = Ят/2, Хь = и* - с*, Хд = и* + с*,
* = (Ни)Ь + (Ни)К * и Нь + Ни , с
Ньл/^Н1 + Н кл/дНя Нь + Нк
Достоинствами метода являются второй порядок точности, простота программирования и нетребовательность к вычислительным ресурсам. Выбор выражения вида (8) для аппроксимации правой части О представляет определенную проблему, если необходимо добиться точного выполнения условия "озеро в покое" в случае произвольной конфигурации дна.
2.3. Метод 8ГИ
Для уравнений мелкой воды в этом методе под "частицей" понимается столбик воды высотой Н, отсчитываемой от дна, и шириной т. Каждая г-я частица имеет массу Шг = Н^т^ и скорость щ. В рамках БРИ-подхода [1] система уравнений, описывающая движение -й частицы, взаимодействующей с окружающими ее соседними частицами и с дном, записывается в следующем виде:
а
= т
У^т (и - ,
зеР
¿НщШг (И
—т.
■Е-
зеР
9 Н , 9 Н2
2
((Нт а
2
0 (Хг
)
+ ^ + Щ) - дНг(1ь, - И)
АгШгз,
и
где суммирование проводится по всему ансамблю Р соседних частиц; А^Ш^ = А^Ш(Хг -х^; Ь) — производная по пространственной переменной от ядра сглаживающей функции для г-й частицы; Ь — длина сглаживания; П^ — необходимая для устойчивости счета искусственная вязкость, вычисляемая по формуле
П
Пг:,Ь СгуигуХгу и х < 0 | | ч ^г^гз
1 ^гз \
0,
ЩзХч > 0-
3
Здесь сгу = (Сг + ^)/2, Сг = л/дЩ,
и
г]
иг , Хг^ Хг , Нг^ (Нг + Ну)/2
аш е [0, 5].
В качестве сглаживающей функции обычно используется кубический сплайн Мона-гана [2], и тогда ядро принимает вид
Ж (з;Ь)
3з2 3з3
1 - т + 0 <"< 1 (2-^ 1 <,< 2,
2
3Ь
4 0,
| 1 ь
:ю)
2,
Длину сглаживания можно выбирать по-разному. Так, в [1] для одномерной модели мелкой воды рекомендуется задавать длину сглаживания постоянной, а в работе [11] она становится переменной и вычисляется по формуле Ьг^ = (тг + Wj)/2.
Достоинствами метода являются второй порядок точности и устойчивость к погрешностям задания батиметрии. Недостатки метода — приближенное выполнение условия "озеро в покое" в случае произвольной конфигурации дна, большие потребности в вычислительных ресурсах, необходимость подбора дополнительных параметров алгоритма (длина сглаживания, искусственная вязкость).
2.4. Метод ХУБ + 8ГИ
В основу алгоритма, который далее будет называться ТУВ + БРИ, положен метод, предложенный в работе [9] и представляющий собой оригинальную комбинацию метода БРИ и конечно-разностной схемы, обладающей ТУВ-свойствами. В узлах сетки, являющихся центрами "эйлеровых" ячеек, определены все искомые функции Н? и (Ни)1?, а также подвижные "лагранжевы жидкие" частицы, имеющие объемы V™ = Н?Ах, импульсы Р? = (Ни)?Ах и положения в пространстве ж?.
Алгоритм включает в себя два основных этапа и аналогично методу крупных частиц использует расщепление по физическим процессам. Первый этап состоит в определении характеристик частиц У^1, Р'?+1 и их положений х?+1 на следующем шаге по времени с помощью процедуры предиктор-корректор. Здесь используется стандартный БРИ-подход так, что движение каждой г-й частицы происходит только под действием гидродинамических сил, причем взаимодействующими считаются только две соседние частицы:
¿Нг о ^(Ни)г
сИ
сИ
п и ^
ах
¿хг (Ни)г
Н,
цг = Нг - к.
В этом случае длина сглаживания равна Ах.
Выражение для аналога градиента сил давления в соответствии с методом БРИ имеет вид
¿Г]
¿х
з=г+1
г+1 Н Н
Е (1^г - ^ 1А) - (|*° - х% Ах)
=г 1
з=г-1
I
где в качестве сглаживающего ядра Ш также используется выражение (10).
На втором этапе вычисляются потоки массы и импульса через границы ячеек на промежуточном временном слое (п +1/2) по формуле (9) путем решения задачи Римана
в сочетании с ТУВ-подходом. Здесь в системе уравнений (1) отбрасывается градиент гидродинамических сил
5 я + д¥ _ ( Н N _ _ (Ни \
т + дь _0, я _\ни)' Р ^Ни2 )
и применяется стандартная процедура конечно-разностной аппроксимации, используемая при построении схем годуновского типа, как было указано выше:
Я _ Чг - д^*г+1/2 - Р г-1/2,1 .
Здесь Я7+1 — значение вектора я, полученное после первого этапа. При приближенном решении задачи Римана (9) консервативные переменные я слева и справа от границ ячеек вычисляются на полуцелом слое по времени яГ+1/2 _ (Ч7+1 + Я?)/2, что обеспечивает второй порядок точности по пространству и времени.
Рассчитываются новые значения импульсов и объемов частиц, обусловленные их потоками через границы ячеек, и затем происходит возвращение "лагранжевых" частиц в исходное положение, т.е. в центры "эйлеровых" ячеек. Для определения положения точки уреза х0(Ь) также применялся упомянутый выше метод проектирования.
Достоинствами метода являются второй порядок точности, достаточная устойчивость, точное выполнение условия "озеро в покое" и нетребовательность к вычислительным ресурсам. В качестве недостатка следует отметить приближенный, однако более точный, чем в методе крупных частиц, расчет положения линии уреза.
2.5. Использование аналитических решений уравнений мелкой воды при постановке разностных краевых условий на подвижной линии уреза
В работе [6] предложен достаточно точный численный метод, использующий неравномерную подвижную разностную сетку и аналитическое решение задачи для расчета положения точки уреза. К сожалению, построить аналог этого метода для решения пространственных задач пока не удалось, однако его вполне возможно использовать в качестве "точного" для сравнительного анализа других методов при численном решении одномерных задач.
В рассматриваемом здесь методе используется явная схема предиктор-корректор [12] второго порядка аппроксимации, на шаге предиктор которой вычисляется вектор потоков
™ + Р™ г /1
г,+1/2 _ - 5( 7 № р - £С)) т/2, (11)
а на шаге корректор определяются полная глубина Н*г+1 и импульс (Ни)™+1:
( 7^+1 ( 7о)» (Р - (х*яГ) - (* - (х^яГ) ,
( 7 я)г - ( 7я)г + V_Л+1/2 V_Л-1/2 _ ^* (12)
т К 1
Обозначения взяты из работы [12], в частности, 7 — это величина, пропорциональная шагу неравномерной сетки, _ (х™+1 - х^ N, (х4)" _ (х^1 - х^ /т — скорость движения г-го узла подвижной сетки |х™} (г _ 0,..., N). В отличие от рассмотренных выше алгоритмов, здесь расчет ведется только в области, занятой водой. При этом на
каждом временном слое £ = ^ крайний левый узел х^ сетки совпадает с подвижной точкой уреза х0(¿"О-
При расчете величин на (п + 1)-м слое по схеме (11), (12) необходимо знать полную глубину и скорость в точке уреза, а также ее положение. Полная глубина определяется непосредственно из краевого условия (2): Н$+1 = 0. Для нахождения нового положения точки уреза и скорости ее перемещения используются результаты [6] аналитического исследования решений уравнений мелкой воды (1) в окрестности этой точки. В [6] решения задачи (1)-(3) построены в малой по времени и пространству такой окрестности для всех возможных режимов взаимодействия волны с берегом: 1) накат необрушающейся волны; 2) накат необрушающейся волны с вырождением (в момент времени Ь = касательные в точке уреза к свободной границе и рельефу дна совпадают); 3) накат с обрушением.
Способ определения положения х$+1 = х (%п+1) и скорости и$+1 = щ (^+1) точки уреза на новом временном слое зависит от режима взаимодействия. Для первого режима эти величины определяются в виде локально сходящихся рядов
Т2 т3 Т4 т2 т3
Хо+1 = Хо +ип0т + х02— +х03— + Х04-; + 0(т5), ип0+ = и% + Х02Т + Х03— +Х04— + 0(т4).
2 о 24 2 о
Разностные формулы для коэффициентов хо к (к = 2, 3,4) получены в [13].
Для второго и третьего режимов взаимодействия волны с берегом закон движения точки уреза определяется из решения системы двух обыкновенных дифференциальных уравнений
Х0 = и0 (г), щ = дЫ (х0 (г)) (13)
с соответствующими начальными условиями [6]. В случае плоского откоса (4) система (13) имеет точное аналитическое решение
т2
Xq+1 = Xq + ти* + —д tg р, Uq+1 = u* + т gtgp,
где и* = для второго режима и и* = и™ — 2Л/дН'? — для третьего.
Достоинствами метода являются теоретически обоснованное точное выполнение условия "озеро в покое" и отсутствие нефизических осцилляций численного решения. Применение метода проектирования подвижной границы в узел разностной сетки здесь не требуется, поскольку эта граница всегда попадает в левый узел сетки.
Недостаток метода — необходимость подбора двух дополнительных параметров алгоритма [14], отвечающих за выделение на каждом шаге по времени конкретного из трех возможных режимов взаимодействия волны с берегом.
3. Некоторые результаты расчетов
Как справедливо отмечено в сборнике трудов семинара исполнителей Национальной программы смягчения ущерба от цунами (National Tsunami Hazard Mitigation Program — NTHMP) [15], анализ численных моделей представляет собой непрерывный процесс, и даже проверенные численные модели должны подвергаться дополнительному тестированию по мере получения новых знаний, уточнения лабораторных материалов, разработки и совершенствования математических моделей и численных методов. Следует
формулировать новые тестовые задачи, учитывающие совершенствование представлений о механизмах генерации волн цунами и воздействия этих волн на берега и прибрежные сооружения. Необходимым является процесс регулярной проверки на этих обновленных задачах математических моделей и численных алгоритмов (используемых и вновь разрабатывамых) для решения теоретических и прикладных задач проблемы цунами.
Содержание настоящего раздела основано на результатах решения известных тестовых задач с помощью рассмотренных выше алгоритмов. Анализ полученных результатов позволяет обосновать выбор одного из этих алгоритмов (ТУВ + БРИ) как наиболее эффективного для обсуждаемого класса задач и продемонстрировать его возможности на ряде известных модельных примеров.
3.1. Сравнение алгоритмов
Первая серия расчетов проводилась для сравнения различных подходов к моделированию наката волны на берег. Цель численных экспериментов заключалась также в анализе на конкретном примере основных характеристик используемых алгоритмов и выборе наиболее перспективного для моделирования взаимодействия волн с берегом в реальных акваториях. Начальные данные (3) при моделировании наката уединенной волны на откос (4) задавались такими же, как в работе [16]:
где А — амплитуда начального возмущения; х,ш — начальная координата максимального смещения свободной поверхности. В расчетах полагалось
Необходимым этапом исследования свойств численных алгоритмов является анализ сходимости результатов на последовательности измельчаемых сеток. Такой анализ, выполненный для метода ТУВ + БРИ, демонстрируется на рис. 1, а, из которого следует, что с достаточной степенью доверия можно относиться к результатам численного решения уже на сетке, содержащей N _ 2000 узлов. Величина заплеска на склон при этом определялась как значение уровня воды в точке уреза гц(х0(Ь), ¿). Для сходимости других рассматриваемых в статье методов требуется, как правило, более детальная сетка. Так, например, для метода крупных частиц сходимость обеспечивается лишь при N _ 8000. Таким образом, проявляется еще одно достоинство метода ТУВ + БРИ, состоящее в довольно высокой скорости сходимости.
Сравнение различных методов моделирования динамики точки уреза при накате уединенной волны на плоский откос показывает (рис. 1, б), что значения продолжи-тельностей первого заплеска, рассчитанные каждым из рассмотренных алгоритмов, хорошо согласуются, в то время как на значения высот наката влияют диссипативные свойства вычислительных схем. Поэтому величина заплеска, определенная с помощью метода БРИ, оказывается наименьшей, в то время как наибольшая величина, практически не отличающаяся от "точной" [6], вычисляется годуновскими методами и методом
А _ 0.2 м, х,ш _ 25.6 м, 3 _ 15°, К0 _ 1 м, Ьх _ 45.6 м.
ТУВ + БРИ. Промежуточное положение занимает результат, полученный с использованием метода крупных частиц. При воспроизведении стадии отката по-прежнему близки (практически идентичны) результаты, полученные с помощью методов Годунова и ТУО + БРИ. Последующий за откатом процесс в "точном" решении представляется более сложной волновой картиной, включающей низкоамплитудные и короткопериодные колебания точки уреза.
Заметим, что сравнение методов крупных частиц и ТУВ + БРИ представляет для авторов особый интерес, поскольку первый из этих методов ранее был адаптирован к особенностям моделирования наката волн цунами, проверен на решении ряда тестовых задач и использован для определения характеристик заплесков на участки восточного побережья Японии и дальневосточного побережья России [7]. Оценивая результаты, полученные с помощью этих двух методов при решении рассматриваемой модельной
а б
Рис. 1. Динамика точки уреза, определенная с помощью различных методов: а — ТУБ + ЯГИ на сетке с числом узлов N = 250 (1), 500 (2), 1000 (3), 2000 (4), 4000 (5); б — ЯРИ (1), метода крупных частиц (2), ТУБ + ЯРИ (3), метода, представленного в работе [6], (4) и метода Годунова (5)
а б в
Рис. 2. Динамика точки уреза при многократном накате уединенной волны на откос с углом наклона £ = 2.5° (а), 5° (б), 10° (в). Метод ТУБ + ЯРИ
задачи, можно сказать, что точность метода крупных частиц оказалась существенно ниже и его сходимость обеспечивается лишь при N _ 8000. Поэтому предпочтение для дальнейших расчетов отдано методу ТУВ + БРИ, с помощью которого были получены все последующие результаты.
Моделирование многократных накатов волн на откос и их отражений от вертикальной стенки (рис. 2), замыкающей расчетную область справа, указывает на возникновение гидравлических скачков, движущихся в прямом и обратном направлениях. Сравнение вариантов, рассчитанных для разных углов наклона откоса, демонстрирует уменьшение временных интервалов между основными "пиками" заплесков для более крутого склона, связанное с тем, что такие волны проходят вдоль склона меньшее расстояние, а кроме того, большие их амплитуды вносят свой вклад в увеличение скорости распространения. Наличие вторичных "высокочастотных" пиков малой амплитуды, расположенных между основными колебаниями точки уреза, связано с заплесками волн, порождаемых в ходе чрезвычайно сложных множественных волновых взаимодействий и отражений как от склона, так и от стенки, сопровождаемых перераспределением энергии колебательного процесса.
3.2. Моделирование обрушающейся волны
Волны, описываемые в рамках модели мелкой воды, даже в случае дна постоянной глубины могут обрушаться в процессе своего движения. Как было показано выше, при многократном накате на откос возникает система гидравлических скачков, движущихся как в сторону откоса, так и в противоположном направлении. В таких случаях численные алгоритмы, как правило, генерируют нефизичные осцилляции, поэтому важное значение приобретает тестирование применяемых вычислительных методов на разрывных решениях. В настоящем разделе будет показано, как метод ТУВ + БРИ воспроизводит процесс появления скачка, его трансформацию вплоть до момента градиентной катастрофы и дальнейшее поведение разрывного решения при накате волны на противоположные вертикальные стенки и откате от них.
В тестовой задаче численное решение сравнивается с точным решением задачи Ко-ши [17]. Предполагается, что расчетная область включает участок постоянной глубины К(х) = К0 и при Ь _ 0 в интервале 0 < х < Ьх заданы начальные условия г](х, 0) _ щ(х), и(х, 0) _ и0(х). Точное решение задачи Коши для уравнений мелкой воды можно найти с помощью метода характеристик. В зависимости от начальных данных этот способ построения решения можно применять при всех Ь > 0 либо только до момента градиентной катастрофы
Рассмотрим волновые процессы, порожденные начальными возвышениями свободной поверхности в виде одиночных волн длиной Л:
В этом случае характеристики, уходящие с отрезка [хш — Л/2, х,ш] оси абсцисс, образуют сходящийся пучок, а с отрезка [хад, хт + Л/2] — расходящийся. Таким образом, в решении будут присутствовать волна разрежения и идущая перед ней волна сжатия, приводящая к градиентной катастрофе. При 0 < Ь < длина волны остается постоянной и равной длине Л начальной волны (15). Амплитуда также остается неизменной
|х — хш| < Л/2, |х — хт| > Л/2.
(15)
Г), м
0.15
0.10
0.05
0.00
20
-и, м/с 0.00
Ь -0.20
▲ ▲ ▲ 4
--3
• • • 2
--1
- -0.40
-0.60
30
х, м
0.00-
и С
Рис. 3. Расчет обрушающейся волны: а — сопоставление результатов расчетов по методу ТУБ + БРИ (маркеры) с аналитическим решением в виде волны Римана (сплошные линии) в момент времени £ = 5 с для возвышения свободной поверхности (1, 2) и скорости (3, 4); б — колебания свободной поверхности на левой стенке (1), правой стенке (2) и в точке х = хт первоначального положения вершины волны (15)
а
на промежутке 0 <Ь < однако вершина движется быстрее, чем точки, расположенные спереди и сзади, в результате чего профиль волны, движущейся влево со скоростью — \[дК~0, будет меняться со временем, укручиваясь на переднем фронте и выполаживаясь на заднем. Согласованная с возвышением (15) начальная скорость и0(х) вычисляется
по формуле _ _
и0(х) _ 2л/дК0 — 2л/д [щ(х) + К0]. (16)
На рис. 3, а показаны точное решение в момент времени _ 5 с и численное, полученное методом ТУВ + БРИ для начального возвышения с параметрами
А _ 0.2 м, х,ш _ 30 м, Л _ 10 м, К0 _ 1 м, Ьх _ 40 м.
Для этих значений градиентная катастрофа наступает при ¿* ^ 5.57 с, поэтому в момент времени _ 5 с передний фронт точного решения уже близок к вертикальному.
Из рис. 3, а видно, что метод ТУВ + БРИ очень хорошо воспроизводит аналитическое решение с большими градиентами, при этом осцилляции в численном решении не появляются. Динамика процесса многократных заплесков на стенки канала с дном постоянной глубины, показанная на рис. 3, б, подтверждает описанный выше сценарий эволюции волн возмущения, когда амплитуда быстро формирующихся ударных профилей убывает, а их длина возрастает (см. также [10]).
3.3. Проверка математических моделей и алгоритмов
Особого внимания заслуживают вопросы валидации используемых математических моделей и верификации реализующих эти модели алгоритмов на известных аналитических решениях и экспериментальных данных. Некоторые результаты такого рода исследований изложены в работах [7, 14]. Ниже представлены результаты оценки точности метода ТУВ + БРИ с помощью принятых в международном сообществе тестовых
задач. Наиболее подробное изложение постановок этих задач и результатов их решения с помощью хорошо зарекомендовавших себя алгоритмов и программ, созданных и используемых в различных коллективах северо-американских исследователей цунами в рамках Национальной программы снижения ущерба от цунами (National Tsunami Hazard Mitigation Program — NTHMP), приведено в трудах специального семинара, проведенного с 31 марта по 2 апреля 2011 г. в Техасском A&M университете. На этом семинаре, по существу, были приняты некоторые стандарты, удовлетворение которым стало обязательным для всех средств математического моделирования цунами, разработка и поддержка которых финансируется за счет NTHMP. Материалами для сравнения являются аналитические решения, результаты лабораторных экспериментов и полевых измерений.
Первоначально набор тестовых задач был сформирован несколько ранее. Его истоки связывают с проходившими в 1995 г. (Friday Harbor, Washington) и в 2004 г. (Catalina Island, California) семинарами по моделям наката длинных волн. Почти завершенный вид этот набор приобрел в отчете [18] и (в сокращенном виде) в опубликованной через год статье тех же авторов [19]. Необходимые для тестирования материалы представлены в сетевом пространстве в электронной форме [20]. Входные данные и результаты экспериментов размещены в Интернете (http://depts.washington.edu/clawpack/links/nthmp-benchmarks/).
Непосредственно с тематикой настоящей статьи связаны две тестовые задачи (BP1 и BP4). Первая из них предусматривает моделирование наката необрушающейся волны малой амплитуды A/h0 = 0.019 на плоский откос с углом наклона [ = arc cot (19.85) и сравнение полученных результатов (в виде профилей свободной поверхности в определенные моменты времени и мареограмм, рассчитанных в заданных точках расчетной области) с аналитическими решениями, опубликованными в статье [16]. Условиями тестовой задачи BP1 предполагается, что начало системы координат (x = 0) совпадает с начальным положением точки уреза, поэтому начало откоса имеет координату xs = h0 cot [. Кроме того, начальный профиль волны и ее начальная скорость несколько отличаются от (14):
r](x, 0) = A sech^-fX ^ , u(x, 0) = V(x, 0).
Величина L/h0 = arccosh (л/20^ /j интерпретируется как половина длины уединенной волны и определяет расстояние от кромки откоса (с координатой xs) до гребня волны (с координатой xw = xs + L).
Вторая тестовая задача (BP4) состоит в моделировании наката на склон с тем же углом наклона необрушающейся (A/h0 = 0.0185) и обрушающейся (A/h0 = 0.3) волн с последующим сопоставлением полученных результатов с материалами, полученными в экспериментальном лотке Калифорнийского технологического института (профили свободной поверхности в определенные моменты времени и максимальные значения вертикальных заплесков) для различных (более 40 реализаций) комбинаций значений величин A, h0. Условия этой задачи совпадают с изложенными выше для задачи BP1.
Результаты решения первой из упомянутых тестовых задач (рис. 4) позволяют говорить об успешной проверке алгоритма TVD + SPH, с достаточной точностью воспроизводящего аналитическое решение на сетке с числом узлов N = 4000.
Задача BP4 позволила оценить также используемую математическую модель путем сравнения результатов математического моделирования, выполненного с помощью
проверенного численного алгоритма, с данными лабораторных экспериментов, проводившихся в лотке длиной 31.73 м, глубиной 60.96 см и шириной 39.97 см, дно и боковые стенки которого выполнены из окрашенной нержавеющей стали. На рис. 5 представлены (в безразмерной форме) результаты математического и лабораторного моделирования вертикального заплеска К/Ы0 в зависимости от амплитуды начального возвышения А/ко. Графики показывают, что при малых амплитудах значения, полученные
T]/h0
0.04 0.03 0.02 0.01 0.00 -0.01
A o o o 4 - 3
M • • • 2 ----¡
Л
' * é « ' V > y Щ
é i # * Á г Ъ \ \ \
\
20
40
60
80 tí to
8 12 16 x/h0
Рис. 4. Аналитическое решение тестовой задачи ВР1 и ее численное решение, полученное методом ТУБ + ЯРИ: а — мареограммы, рассчитанные численно (маркеры) и аналитически (линии) в точках склона с координатами х/Но = 9.95 (кривые 1 и 2) и х/Но = 0.25 (3 и 4); б — профили свободной поверхности, полученные численно (маркеры) и аналитически (линии) в моменты времени ¿/¿0 =40 (1 и 4), 55 (2 и 5), 70 (3 и 6); ¿0 = \fhoJg
10.00 -=.
R/ha
1.00 -
0.10 ~А
0.01
"i 1 i 0.6 Alhо
Рис. 5. Аналитическое решение тестовой задачи BP4, ее численное решение, полученное методом TVD + SPH, и соответствующие экспериментальные данные: значения максимального вертикального заплеска, определенные численно без трения (маркеры 3) и с трением (п = 0.08, маркеры 4), аналитически (линия I); экспериментальные данные — маркеры 2. Вертикальная линия 5 соответствует экспериментальному критерию обрушения A/ho = 0.045 [15]
a
в ходе численного моделирования, хорошо воспроизводят данные физического моделирования и неплохо соответствуют максимальным значениям величины R/h0, определяемым по аналитической формуле
R/ho = 2.831 ^cot р (A/ho)5/4, (17)
предложенной впервые в [21] и применимой в пределах
(0.288 tanр)2 < A/ho < 0.478 (tanр)10/9. (18)
Однако с ростом амплитуды падающей на откос волны измеренные в экспериментах максимальные значения вертикального заплеска начинают все более отклоняться от теоретической кривой. Для того чтобы воспроизвести экспериментальные данные, уравнения движения модели (1) модифицируются за счет учета трения потока о дно и боковые стенки довольно узкого лабораторного канала. При этом предполагается, что трение квадратично зависит от скорости течения, а соответствующий коэффициент Шези определяется по формуле Маннинга. Такой диссипативной добавки с коэффициентом шероховатости п = 0.08 оказывается достаточно для хорошего согласования расчетных и экспериментальных данных.
Заметим, что близость к теоретической кривой результатов, полученных без учета трения, естественно, сохраняется значительно дольше, тем самым подтверждая сделанные при решении первой тестовой задачи BP1 выводы об успешной проверке алгоритма. Однако и эти результаты начинают расходиться после того, как в численном решении возникает эффект "обрушения" волны, выражающийся в предельном "укручении" переднего фронта и увеличении длины тыльной части волны.
Изложенные выше соображения подтверждаются графиками (рис. 6), иллюстрирующими динамику максимальных смещений свободной поверхности, полученных в экспериментах и в ходе расчетов с учетом трения и без него. Как следует из графиков на рис. 6, а, амплитуда необрушающейся волны (А = 0.0185 м) монотонно возрастает при движении по склону, в фазе заплеска движущаяся точка уреза достигает вертикального максимума, после чего волна переходит в фазу отката. Введение в модель трения (п = 0.08) практически не сказывается на фазе наката, но существенно замедляет откат волны. В случае обрушающейся волны (А = 0.3 м) ситуация значительно
0 20 40 60 t/t0 0 20 40 60 t/t0
Рис. 6. Результаты расчетов по методу ТУБ + ЯГИ для необрушающейся волны (а) с учетом (п = 0.08, линия 3) и без учета (линия 2) трения, а также для обрушающейся (б) с учетом (п = 0.0075 (линия 3), 0.05 (4), 0.08 (5)) и без учета (линия 2) трения; 1 — экспериментальные данные
меняется (рис. 6, б), что подтверждается и экспериментальными данными. Через некоторое время после начала движения волны начинает развиваться процесс ее обрушения, выражающийся, как уже было отмечено, в "укручении" переднего фронта волны и одновременном снижении амплитуды. Судя по графикам, расчетные и экспериментальные результаты несколько различаются в скоростях моделируемых волновых процессов так, что "обрушение" у расчетной волны начинается гораздо раньше, однако завершается оно примерно в то же время, что и у волны экспериментальной. Учет трения здесь практически не сказывается вплоть до завершения снижения максимального уровня свободной поверхности ("разрушения" волны) и начала его роста, но позволяет с большей точностью воспроизвести рост вертикальной отметки положения точки подвижного уреза на начальной стадии наката. Результаты расчетов, проведенных с различными значениями коэффициента шероховатости п (рис. 6, б), демонстрируют влияние определяемой им величины донного трения. Как видно, наибольшей близости к экспериментальным данным удается добиться при п = 0.08.
Заключение
Представленные в статье результаты сравнительного анализа некоторых методов моделирования наката длинных поверхностных волн на берег в рамках теории мелкой воды показали, что для решения задач этого класса наиболее перспективен метод TVD + SPH. В качестве тестовых задач рассматривались хорошо известные модельные примеры наката волн на плоский откос, либо обладающие аналитическими решениями, либо сопровождаемые достоверными результатами лабораторных экспериментов или результатами других авторов, полученными с помощью альтернативных численных алгоритмов.
Метод TVD + SPH, показавший в ходе тестовых испытаний возможность своего успешного применения для моделировании как слабо нелинейных (необрушающихся), так и существенно нелинейных (обрушающихся) волн, может послужить основой для разработки алгоритмов решения двумерных (в плане) задач наката длинных волн, в том числе — волн цунами, на берега со сложными конфигурациями береговых линий и с близкими к реальным рельефами дна прилегающих акваторий и суши.
Благодарности. Исследование выполнено при поддержке Российского научного фонда (проект № 14-17-00219).
Список литературы / References
[1] Leffe, M. De, Touze, D. Le, Alessandrini, B. SPH modeling of shallow-water coastal flows // Journal of Hydraulic Research. 2010. Vol. 48, Extra Issue. P. 118-125.
[2] Monaghan, J.J. Simulating free surface flows with SPH // Journal of Computational Physics. 1994. Vol. 110, No. 2. P. 399-406.
[3] Cho, Y.-S., Briggs, M.J., Kanoglu, U., Synolakis, C.E., Liu, P.L.-F. Runup of solitary waves on a circular island // Journal of Fluid Mechanics. 1995. Vol. 302. P. 259-285.
[4] Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001. 608 с. Kulikovskii, A.G., Pogorelov, N.V., Semenov, A.Yu. Mathematical aspects of numerical solution of hyperbolic systems: Monographs and Surveys in Pure and Applied Mathematics. Vol. 118. USA: Chapman & Hall/CRC, Boca Raton, FL, 2001. 560 p.
[5] Anastaious, K., Chan, C. Solution of the 2D shallow water equations using the finite volume method on unstructured triangular meshes // International Journal for Numerical Methods in Fluids. 1977. Vol. 24. P. 1225-1245.
[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] Шокин Ю.И., Бейзель С.А., Рычков А.Д., Чубаров Л.Б. Численное моделирование наката волн цунами на побережье с использованием метода крупных частиц // Мат. моделирование. 2015. Т. 27, № 1. С. 99-112.
Shokin, Yu.I., Beisel, S.A., Rychkov, A.D., Chubarov, L.B. Numerical simulation of the tsunami runup on the coast using the method of large particles // Mathematical Models and Computer Simulations. 2015. Vol. 7, No. 4. P. 339-348.
[8] Бейзель С.А., Рычков А.Д., Чубаров Л.Б. Сравнительный анализ некоторых алгоритмов моделирования наката волн цунами на берег //V Всерос. конф. с участием зарубежных ученых "Задачи со свободными границами: теория, эксперимент и приложения". Бийск, 29 июня — 4 июля, 2014. Тезисы докладов. Бийск: Изд-во Алтайского гос. техн. ун-та им. И.И. Ползунова, 2014. С. 16.
Beisel, S.A., Rychkov, A.D., Chubarov, L.B. Comparative analysis of some algorithms for the simulation of tsunami waves runup // V All-Russian conference with foreign participants "Free Boundary Problems: Theory, Experiment and Applications". Biisk, June 29 —July 4, 2014. Abstracts. Biisk: Publishing house of the Altai State Technical University named by I.I. Polzunov, 2014. P. 16.
[9] Храпов С.С., Хоперсков А.В., Кузьмин Н.М., Писарев А.В., Кобелев И.А. Численная схема для моделирования динамики поверхностных вод на основе комбинированного SPH-TVD подхода // Вычисл. методы и программирование. 2011. Т. 12. С. 282-297. Hrapov, S.S., Khoperskov, A.V., Kuzmin, N.M., Pisarev, A.V., Kobelev, I.A.
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.)
[10] Родин А.А., Пелиновский Е.Н. Динамика длинных волн в прибрежной зоне моря с учетом эффектов обрушения. Н. Новгород: Нижегород. гос. техн. ун-т им. Р.Е. Алексеева, 2014. 93 с.
Rodin, A.A., Pelinovsky, E.N. The dynamics of long waves in the coastal zone, taking into account the effects of the breaking. N. Novgorod: Nizhegorod. State Tehn. Univ. named by R.E. Alekseev, 2014. 93 p. (in Russ.)
[11] Паршиков А.Н. Применение решения задачи Римана в методе частиц // Журн. вычисл. математики и мат. физики. 1999. Т. 39, № 7. С. 1173-1182.
Parshikov, A.N. Application of a solution to the Riemann problem in the SPH method // Computational Mathematics and Mathematical Physics. 1999. V. 39, No 7. P. 1173-1182. (in Russ.)
[12] Хакимзянов Г.С., Шокина Н.Ю. Метод адаптивных сеток для одномерных уравнений мелкой воды // Вычисл. технологии. 2013. Т. 18, № 3. С. 54-79.
Khakimzyanov, G.S., Shokina, N.Yu. Adaptive grid method for one-dimensional shallow water equations // Computational Technologies. 2013. Vol. 18, No. 3. P. 54-79. (in Russ.)
[13] Шокина Н.Ю., Хакимзянов Г.С. Усовершенствованный метод адаптивных сеток для одномерных уравнений мелкой воды // Вычисл. технологии. 2015. Т. 20, № 4. С. 83-106.
Shokina, N.Yu., Khakimzyanov, G.S. An improved adaptive grid method for one-dimensional shallow water equations // Computational Technologies. 2015. Vol. 20, No. 4. P. 83-106. (in Russ.)
[14] Бейзель С.А., Шокина Н.Ю., Хакимзянов Г.С., Чубаров Л.Б., Ковырки-на О.А., Остапенко В.В. О некоторых численных алгоритмах расчета наката волн цунами в рамках модели мелкой воды. I // Вычисл. технологии. 2014. Т. 19, № 1. С. 40-62.
Beizel, S.A., Shokina, N.Yu., Khakimzyanov, G.S., Chubarov, L.B., Kovyrki-na, O.A., Ostapenko, V.V. On some numerical algorithms for computation of tsunami runup in the framework of shallow water model. I // Computational Technologies. 2014. Vol. 19, No. 1. P. 40-62. (in Russ.)
[15] National Tsunami Hazard Mitigation Program: Proceedings and Results of the 2011 NTHMP Model Benchmarking Workshop. Boulder: U.S. Department of Commerce, 2012. 436 p.
[16] Synolakis, C.E. The runup of solitary waves // Journal of Fluid Mechanics. 1987. Vol. 185, No. 6. P. 523-545.
[17] Вольцингер Н.Е., Клеванный К.А., Пелиновский Е.Н. Длинноволновая динамика прибрежной зоны. Л.: Гидрометеоиздат, 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.)
[18] Synolakis, C.E., Bernard, E.N., Titov, V.V., Kanoglu, U., Gonzalez, F.I. Standards, criteria, and procedures for NOAA evaluation of tsunami numerical models: NOAA technical memorandum OAR PMEL-135 / Pacific Marine Environmental Laboratory. Seattle, WA, May 2007. 60 p. Available at: http://nctr.pmel.noaa.gov/benchmark/SP_3053.pdf
[19] Synolakis, C.E., Bernard, E.N., Titov, V.V., Kanoglu, U., Gonzalez, F.I. Validation and verification of tsunami numerical models // Pure and Applied Geophysics. 2008. Vol. 165. P. 2197-2228.
[20] Gonzalez, F.I., LeVeque, R.J., Chamberlain, P., Hirai, B., Varkovitzky, J., George, D.L. (USGS) Validation of the GeoClaw Model: NTHMP MMS Tsunami Inundation Model Validation Workshop / GeoClaw Tsunami Modeling Group, University of Washington. 84 p. Available at: http://depts.washington.edu/clawpack/links/nthmp-benchmarks/geoclaw-results.pdf
[21] Synolakis, C.E. The Runup of Long Waves: Ph.D. Thesis. Pasadena, California Institute of Technology, 1986. 228 p. 91125.
Поступила в 'редакцию 17 июня 2015 г.
On numerical methods for solving run-up problems. I. Comparative analysis of numerical algorithms for one-dimensional 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 2 Novosibirsk State University, Novosibirsk, 630090, Russia
* Corresponding author: Chubarov, Leonid B., e-mail: chubarovster@gmail.com
The article presents the results of comparative analysis for some algorithms aimed at numerical modeling of the run up of long waves in the framework of shallow water theory.
© ICT SB RAS, 2015
Authors consider the problem of the waves run-up on a flat slope both including and excluding the forces of bottom friction. Considerable attention is paid to the verification and validation for algorithms of the mathematical model by comparing the calculated results with the known analytical solutions, calculations of other authors, as well as the results of laboratory experiments.
Numerical methods considered in this paper include the method of large particles, a Godunov type method, the method of smoothed particles, a method of smoothed particle hydrodynamics and finite difference scheme with TVD properties (TVD + SPH) and the method with a moving mesh using an analytical solution for the calculation of the flow parameters on the shoreline.
Some advantages of the TVD + SPH method for solving problems of the rolling waves on a beach are shown. Examples of a successful application of this method to simulate weakly nonlinear and essentially nonlinear (the breaking) waves are presented. In the latter case, the quality of the result is largely dependent on the successful selection of parameters defining the bottom friction.
The authors propose a continuation of the research presented in this paper for the simulation of more complex wave propagation phenomena and in extending the effective TVD + SPH method for the case of two spatial (horizontal) dimensions.
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 17 June 2015
