CC BY
28
Вычислительные технологии
Том 20, № 4, 2015
Усовершенствованный метод адаптивных сеток для одномерных уравнений мелкой воды
Н.Ю. ШокинА*, Г. С. Хлкимзянов
Институт вычислительных технологий СО РАН, Новосибирск, Россия *п1па.8Ьок1па@1сЬ.пБс.ги
Представлен модифицированный алгоритм реализации явной схемы предиктор-корректор, основанный на новом способе вычисления правой части уравнений мелкой воды на неровном дне. Предложен новый способ выбора схемных параметров, гарантирующий выполнение ТУБ-свойства для усовершенствованной схемы предиктор-корректор. Описан способ получения консервативных схем на подвижных сетках, основанный на выборе соответствующих значений схемных параметров в схеме предиктор-корректор. На основе известных аналитических решений уравнений мелкой воды в окрестности границы вода — суша разработаны уточненные разностные краевые условия в подвижной точке уреза, аппроксимирующие аналитические решения с большей точностью, чем использованные ранее. Приведены результаты расчетов с помощью усовершенствованного метода адаптивных сеток.
Ключевые слова: нелинейные уравнения мелкой воды, конечно-разностная схема, адаптивная сетка, поверхностные волны, накат на берег.
Введение
Среди многочисленных задач о течениях жидкости со свободной поверхностью важное место занимают задачи моделирования процессов распространения длинных поверхностных волн в больших акваториях и взаимодействия этих волн с берегом. В настоящее время для решения таких задач интенсивно используются методы численного моделирования на основе иерархии вложенных друг в друга математических моделей [1] и оптимальных вычислительных алгоритмов. Нелинейные модели мелкой воды занимают в этой иерархии центральное место. Поэтому разработка и исследование численных методов для решения практически важных задач в рамках моделей мелкой воды является актуальной проблемой современной вычислительной гидродинамики [2].
Данная работа является продолжением исследований, начатых в [3], и посвящена некоторым проблемам численного моделирования в рамках бездисперсионной модели мелкой воды (КЬБШ-модели). Для неровного дна КЬБШ-система является неоднородной, поэтому многие конечно-разностные схемы, хорошо зарекомендовавшие себя в случае горизонтального дна, теряют важное свойство сохранения простейших аналитических решений КЬБШ-модели [4]. В настоящей работе показано, что для схемы предиктор-корректор указанная проблема решается за счет особой аппроксимации правой части уравнения импульса.
© ИВТ СО РАН, 2015
Известно, что решение уравнений КЬБШ-модели может стать разрывным даже при задании гладких начальных данных, что во многих случаях приводит к появлению ос-цилляций численного решения [5 - 7]. Здесь предложен новый способ выбора схемных параметров, основанный на исследовании дифференциального приближения разностной схемы, гарантирующий выполнение для нее ТУВ-свойства.
В работе [3] важное свойство консервативности схем на подвижных сетках для нелинейного скалярного уравнения достигалось с помощью соответствующего выбора схемного параметра. В данной работе такой же прием применен для автоматического построения консервативных схем, аппроксимирующих уравнения мелкой воды на подвижных сетках.
Методы расчета процессов наводнения-осушения, необходимые для определения зон затопления при выходе волны цунами на берег, могут быть разделены на две основные группы [8]. К первой группе относятся методы сквозного счета с прямым применением существующих алгоритмов на фиксированных пространственных сетках, где в каждый момент времени необходимо вычислять положение линии уреза. Ко второй группе принадлежат методы на подвижных сетках, адаптирующихся так, чтобы узлы на границе численной области всегда совпадали с подвижной границей. Для расчета течений жидкости в задачах наводнения-осушения используются методы конечных объемов [9, 10], конечных элементов [11, 12, 13], конечных разностей [14, 15], метод частиц [16] и другие.
В настоящей работе для расчета зон затопления используются подвижные адаптивные сетки и положение границы вода — суша определяется с использованием в окрестности линии уреза точных аналитических решений уравнений мелкой воды на временном промежутке, равном шагу по времени разностной схемы. Приведены уточненные формулы для определения положения и скорости движения точки уреза и результаты численного решения задач о накате волны на вертикальную стенку и плоский откос.
1. Новая процедура вычисления правой части
В работе [3] предложена и исследована конечно-разностная схема на подвижной сетке для решения одномерных уравнений мелкой воды
£ — время; и(х, — скорость; г/(х, — возвышение свободной поверхности относительно невозмущенного уровня; к(х) — глубина дна, отсчитываемая от невозмущенной свободной границы; Н = 'ц+к — полная глубина; д — ускорение свободного падения. Предполагается, что для системы уравнений (1) поставлена некоторая начально-краевая задача на отрезке х Е [0, I] и временном промежутке Ь Е [0, Т].
Координаты узлов подвижной сетки, покрывающей в момент времени Ь = ¿га отрезок [0, /], обозначены через х™ (] = 0,... , N). Эти узлы являются образами при некотором взаимно-однозначном отображении х = узлов ^ равномерной сетки с шагом
к = 1/Ы на отрезке д Е [0, 1]. В схеме предиктор-корректор [3] на шаге предиктор вычисляется вектор потоков в центрах ячеек:
щ + ^ = О,
(1)
где
который используется затем на шаге корректор для определения полной глубины и импульса по явным формулам:
(3"Г1 - V")? + - (*Ц)"),+1/2 - - (Г'Ц)");-1/2 = е*, (3)
г К 7'
где г — шаг по времени; Г" = ¥(""); Р"+1/2 = (Ьц)"+1/2'
Ц" + Ц" Ц" _Ц" _+
Ц" = Ц7 + Ц7+1 Ц" = Ц7+1 Ц7 . ™ = _. ™ = х',7 + х',7+1.
"•+1/2 = 2 ' 4,7+1/2 = К . = г . Х',^'+1/2 = 2 .
- X,- ^7 + 1 /2 + 1/2 / и
• + 1 ^ _ „," . Т" _ ^7+1/2^ ^7-1/2,
• + 1/2 = К = Ж<" + 1/2' ^ = 2 . ^+1/2 = \, £ (ЯК,)"+1/2
(у (як°! )"+1/0.
К"+1 - К" / ч
К" = •+1 7 . К" = у,7+1/2 = К . ' 7 = К(Х7 )
Матрицы Ь, Я, Л из уравнения (2) задаются следующими формулами:
Г" ( -Л"'+'/2 1 \ ™ = 3+!/» ( -1 М
Ь^1/2 = (г" )2 I Л" ' ' 7+1/2 = 2 I Л " Л"
(С7+1/2) V -Л1,7+1/2 1 / 2 \ -Л1,7+1/2 Л2,7+1/2 /
/ Л",7+1/2 0 \ V 0 Л",7+1/2 )
д" _ Д" __С Д"
^ 7+1/2 = Л7+1/2 Х',7+1/2С' Л7+1/2 ,
0 ¿,7+1/2
где £ — единичная матрица; А"7-+1/2 (& =1, 2) — собственные значения матрицы Л"+1/2, аппроксимирующей матрицу Якоби системы (1):
Л"+1/2 =( 2И"1 )=(ЯЛЬ)"+1/2 ' (4) \ "7 "7+1 + Ул7+1/2 2а7+1/2 /
^",7+1/2 = (м - с)"+1/2 ' ^",7+1/2 = (м + с)"+1/2 ' (5)
с"+1/2 = ^(«"+1/2)2 - М"М"+1 + ##"+1/2 — > 0 •
Для вычисления элементов диагональной матрицы
1 + ^/о 0
/ 1 + *}+1/2 0 \ V 0 1 + ^2+1/2 У
Т)"
®,+1/2 ' 0 1+^
в работе [3] предложена следующая формула для схемных параметров:
Г 0' I ^7^+1/2 I < | ^7^+1/2-^^ I и ^7?+1/2 ■ #7+1/2-, — 0'
^7+1/2 = ^ ^0,7+1/2(1 - Й+1/2^ 1 ^^+1/2 1 > 1 #7+1/2-^ 1 и ^^+1/2 ■ #7+1/2-^ — 0'
I ^0,7+1/2' ^?+1/2 ■ ^^+1/2-5^ < 0'
где к = 1' 2; £7+1/2 = #7+1/2-/#7+1/2; ^ = 5",7+1/2 = дД",7+1/2;
^1 \" . ^ _ 1 _1.
#7+1/2 = (|дД^ |(1 - С^ )"+1/2. С",7+1/2 = ) "+1/2. 00,7+1/2 = С"
й,7+1/2
Pkj+i/2 — компоненты вектора P™+1/2; X%j+1/2 = (\k — xt)™+1/2 —диагональные элементы матрицы Л™+1/2; ж = т/h. В расчетах использовалась модифицированная формула (6), полученная [3] в результате энтропийной коррекции схемы на основе метода дифференциального приближения и гарантирующая неотрицательность аппроксимационной вязкости в окрестности звуковой точки волны разрежения.
Величины C^ j+1/2 естественно называть локальными числами Куранта, соответствующими собственным значениям Х^ j+1/2. Аналогично случаю равномерных неподвижных сеток шаг по времени т выбирается из соображений устойчивости вычислений на подвижных сетках, т. е. так, чтобы выполнялось условие
max C£j+1/2 < 1.
к, J
В работе [3] правая часть уравнения (3) вычислялась по формуле
0
*
^ V 9Н*К,3
поэтому на шаге предиктор дополнительно к двум уравнениям системы (2) приходилось использовать еще одно уравнение, которое получалось в результате аппроксимации первого уравнения системы (1) и служило для вычисления полной глубины Н* в целых узлах. Для этого дополнительного уравнения необходимо было вычислять в целых узлах величины и б^ с использованием сложных формул, отличных от (5), (6). Такие дополнительные вычисления увеличивали время решения одномерных задач, и особенно существенно — двумерных.
Рассмотрим более экономичную процедуру вычисления правой части уравнения (3):
О* = ( д(як°0);+1/0 , (7)
где
ип+1/2 = П] + 1 + П]-1 + nj + 1 + nj-1 1П+1/2 = к] + 1 П]-1 + ,Ь] +1 ,Ь]-1
пз = 4 ; пя,3 = 4к .
Интересно, что несмотря на то, что теперь в (7) присутствуют величины с (п + 1)-го временного слоя, тем не менее схема (2), (3) по-прежнему остается явной. В модифицированном алгоритме вычислять величины Н* не нужно, а шаг корректор (3) выполняется отдельно для Н™+1 и (Ни)™+1, т. е. вначале из уравнения неразрывности находятся величины Н™+1, а затем по явной формуле определяется импульс (Ни)™+1 с использованием в правой части (7) уравнения движения вычисленного значения Н™+1.
Итак, в отличие от [3], теперь на шаге предиктор решаются два уравнения вместо трех, что приводит к заметному сокращению временных затрат. При этом все положительные свойства схемы предиктор-корректор, отмеченные в [3], сохраняются и в случае вычисления правой части по формуле (7). В частности, остается верным утверждение о сохранении состояния покоя.
Лемма 1. Пусть неравномерная сетка {хявляется неподвижной и на п-м слое по времени жидкость находится в состоянии покоя:
= 0, и[Ч = 0. (8)
Тогда при использовании схемы предиктор-корректор (2), (3), (7) состояние покоя сохранится на (п +1 )-м слое.
Доказательство. Из условий (8) и условия неподвижности сетки XI = 0 следует
(ЛР - сс)™+1/2 = 0
поэтому на шаге предиктор
^ 3+1/2
1 з + 13+1
Тогда из формулы (3) шага корректор с учетом (7), а также равенств Щ+1 = Щ и = З™ для неподвижной сетки следует что = 0, и™+1 = 0. Таким образом, состояние покоя сохраняется и на (п + 1)-м слое по времени. □
2. Новый прием переключения схем
Формула (6) регулирует переключение с одной схемы на другую в зависимости от локального поведения вектора Р = Сщ. Так, например, при выполнении первого из условий в (6) в обеих соседних точках (1^-1/2 и (1^+1/2 и для обоих значений к = 1 и к = 2 расчет предикторных величин (2) выполняется с единичной матрицей "Р = 8, т. е. схема (2), (3) будет в этом случае локально (в узле ) совпадать со схемой Лакса — Вендроффа на подвижной сетке. Для других соотношений между компонентами вектора Р в соседних узлах получаются либо известные, либо новые разностные схемы. В частности, при
0^+1/2 = ^0^+1/2, ^3-1/2 = ®\,з-1/2, ^ = 1, 2, (9)
схема (2), (3) в узле qj превращается в противопоточную схему первого порядка аппроксимации.
Переключение по формуле (6) фактически означает выбор той или иной схемы на основе анализа поведения разностной производной щ от вектора решения, т. е. производных Нд и (Ни)д от полной глубины и полного импульса. Поскольку Ня = г]я + кя, при таком способе переключения поведение функции к(х), описывающей дно, может оказывать сильное влияние на выбор схемы. Но на основе опыта многочисленных расчетов движения поверхностных волн над неровным дном авторы данной работы пришли к выводу от том, что переключатель схем должен анализировать не изменение полной глубины Нд, а изменение возвышения свободной границы — т.е. производную г]д, а также производную от скорости щ. При новом способе переключения все формулы сохраняют свой прежний вид, только в выражениях д^+1/2 вместо компонент вектора Р
теперь используются компоненты р)1/^ ^+1/2 нового вектора Р, не содержащие производную Нд .
Из определения вектора Р следует, что
Рп = (ги ^ = 1 ( -сщ + Нщ\п +/кЛп ( -1 \
3+1/2 = ( 9Ь+1/2 =(С-+1/2)Ч + Нщ)]+1/2 +и )Ш/А 1) .
Таким образом, вектор Р представлен в виде суммы двух векторов, первый из которых зависит от поведения решения (от производных г]д, щ), но не зависит от изменений глубины кя. Этот первый вектор и берется в качестве нового вектора Р.
3. Консервативные схемы на подвижных сетках для уравнений мелкой воды
Известно [6, 17], что на разрывных решениях неконсервативные схемы могут расходиться, поэтому консервативность является одним из ключевых требований, предъявляемых к разностным схемам. Для неоднородных систем уравнений вида (1) под консервативностью разностной схемы понимается выполнение разностных аналогов уравнений баланса, выражающих собой законы изменения так называемых консервативных переменных. Для системы уравнений мелкой воды (1) такими переменными являются полная глубина Н и полный импульс Ни. Уравнения баланса для этих переменных описывают закон сохранения массы и закон изменения полного импульса. В случае ровного дна закон изменения полного импульса превращается в закон сохранения.
Построение консервативных схем на подвижных сетках осложняется тем, что формулы для скоростей движения узлов и для вычисления площадей ячеек подвижной сетки должны быть строго согласованными. Особая аппроксимация указанных величин гарантирует выполнение разностного аналога геометрического закона сохранения [18, 19], являющегося необходимым условием согласованности формул. В работах [3, 18] предложена оригинальная процедура получения консервативных схем для скалярного закона сохранения, основанная на выборе схемного параметра , соответствующего конкретной схеме. Этот прием можно использовать и для систем уравнений.
Например, согласно [18], для получения противопоточной схемы на подвижной сетке для уравнений мелкой воды достаточно в консервативной схеме предиктор-корректор (2), (3) задать параметры по формулам (9). Приведем вид полученной таким
способом консервативной противопоточной схемы на подвижной сетке в эквивалентных между собой дивергентной и недивергентной формах записи.
Итак, если в схеме предиктор-корректор (2), (3) использовать параметры вк, определенные равенствами (9), то
К/
К/
К/
КЛ
РЛ = — 5, РЛ2 = — (Л+ - Л") , ЯРЛ2£ = — (Д+ - А~) , = —
и схема превращается в противопоточную схему, записанную в дивергентной форме:
( - ( /и)П _
1
+К
СП | СП
г.?+1 + Ь'
(Д+ - Д-)
.+1/2
К \
(иП+1 - иП) - 0г4и)П+1/2 + К (^С)П+1/2 )
С + 1 (Д+ -
2
2
.-1/2 ,П
К
иП - иП-0 - (®4и)П-1/2 + 2 (Я5£0)П-1/2
где правая часть вычисляется по формуле (7)
ъ = Ы / 1/|А1| 0 = Г ^ 0 1/|А2|
), 5 (ов 2), Л (олг)
- (10)
С**,
Ак
Ак ± I Ак
Д
±
Л±
П
2
П
при этом справедливы равенства
Л = Л+ + Л-, А+ + Д" = А = А -хг£.
Для сокращения записи индексы п и ] + 1/2 здесь опущены.
Левая часть консервативной противопоточной схемы (10) записана в дивергентной форме. Путем эквивалентных преобразований эту схему можно записать в недивергентной форме, аналогичной записи противопоточных схем для скалярных уравнений. Для этого используются равенство
СП
1 4,3+1/2
(Ли,
п
^ 3+1/2
справедливое для матрицы А, введенной по формуле (4), и тождество
1
где
(х4и)
= (х^и)? + 2 [(х*и?^+1/2 + (х*и9)Г^-1/2 г _ (х*и)"+1/2 - (х*и)"-1/2 .
ч,з
к
х
_ /у' ь
1,]+1/2 х 1,3-1/2
к
В результате разностное уравнение (10) принимает следующий вид:
( 3и)^1 - ( 3и)? + 1 т +2
Аид - [А+ -А-) щ - хи + ПвСС
1 + 2
Аид + {А+ -А-) щ - хи - ПБСС
3-1/2
- (х!ди)
3+1/2
С*.
+
Далее, с использованием равенства А = А+ + А и геометрического закона сохранения [18]
3-+1 = 3- + тх^ (11)
получается окончательный вид записи противопоточной схемы в недивергентной форме:
и^1 - и?
+
37+1
(л+и);_ 1/2 + (л- и,)
п
1) 3+1/2
С* - 1 ((П8ССТ+1/2 - (ПБСС)^ 1/2
:12)
Явная противопоточная схема (12) реализуется следующим образом. Сначала из первого уравнения системы (12) находится полная глубина Н™+х во всех узлах сетки, а затем, с использованием новых значений Нп+1 в правой части (7), вычисляется полный импульс ( Ни)™+1. Кроме того, что эта схема консервативная, еще одним ее преимуществом является то, что она сохраняет состояние покоя жидкости.
Лемма 2. На неподвижной сетке противопоточная схема (12) сохраняет состояние покоя (8).
Доказательство. Сохранение состояния покоя г]г^+1 = 0, и™+1 = 0 на (п + 1)-м слое по времени следует из формулы (12) с учетом условий (8), равенств х1 = 0 и Щ+1 = Щ для неподвижной сетки и формулы для правой части (7). □
71
П
1
1
4. Тестирование схем на точном решении уравнений мелкой воды
Для разностных уравнений, аппроксимирующих систему нелинейных уравнений мелкой воды, не удается строго обосновать некоторые свойства решений, присущие численным решениям линеаризованных систем или скалярных нелинейных уравнений. Поэтому особую важность приобретает исследование свойств схем на численных тестах. В проблемах, связанных с накатом волн на берег, одним из самых распространенных тестов является задача наката волн на плоский откос, сопрягающийся с горизонтальным дном постоянной глубины [20 - 22]. При численном решении этой задачи начальная волна располагается над горизонтальным дном и некоторое время движется над ним в сторону откоса, при выходе на который и начинается непосредственно процесс наката волны на берег. Для достоверной оценки зон затопления необходимы не только качественные алгоритмы для определения движения линии уреза, но и надежные численные методы для моделирования начального этапа движения волны, в том числе и над горизонтальным дном. Приведем некоторые результаты такого моделирования с использованием рассмотренных выше методов на адаптивных сетках.
Итак, в рассматриваемой тестовой задаче предполагается, что дно горизонтальное, h(ж) = Ко, и при t = 0 заданы начальные условия
?/(ж, 0) = ?/0(ж), и(ж, 0) = и0(ж). (13)
Численное решение сравнивается с точным решением задачи Коши (1), (13), представленным в [5]. В отличие от нелинейно-дисперсионных уравнений [23, 24], для которых задача с точным решением в виде уединенной волны, распространяющейся над ровным дном, входит в список обязательных тестов, точное решение [5] уравнений бездисперсионной модели мелкой воды редко используется для тестирования численных алгоритмов решения. Авторы данной работы полагают, что причиной непопулярности этого теста является недостаточно подробное описание процедуры вычисления точного решения в [5]. Поэтому ниже предложен другой способ получения формул точного решения. Запишем уравнения мелкой воды (1) в инвариантах Римана [7]
3г + s г + 3s _
П + Гж = 0, St + ^ = 0, (14)
где
г = и — 2с; s = и + 2с; с = л/ </Н; Я = т/ + h0. (15)
Решение уравнений (14) ищется в виде распространяющейся справа налево г-волны [25], т. е. в предположении, что s = s0 = const. Полагая, что на бесконечности жидкость покоится, получаем s0 = 2л/fih0 = 2с0. Тогда из формулы и + 2с = s0 следует, что скорость выражается через полную глубину:
и(ж, i) = 2С0 — 2^Я(ж, i) = 2С0 — 2^0 Мж i) + h0], (16)
зд(ж) = и (ж, 0) = 2 С0 — 2\J (7Н0 (ж) = 2 С0 — 2^ [г/0 (ж) + Л^]. (17)
Если ввести новую зависимую переменную
3 г + s 0
то уравнение для инварианта г преобразуется в уравнение Хопфа + = 0, для которого решение задачи Коши находится методом характеристик и задается неявно формулой [25]
р(ж, ¿) = ро(ж — (19)
где р0 — начальная функция,
( ) 3 г(ж, 0) + в о 3(Ио(ж) — У#Я°(ж))+2со 2 з ^^ (20) Ро(х) =-4-= —-4---= 2 Со — Яо0*0. (20)
Тогда из (19) следует, что для заданных ж, £ величина р(ж, ¿) является корнем р уравнения
р = 2со — 3^ (Ко + г/о(ж — р£)). (21)
Найденное р используется для определения возвышения свободной поверхности жидкости г/(ж,¿). С учетом равенств (15), (16), (18), получаем, что
ч(х, г) = ()' — Ко. (22)
Итак, точное решение уравнений (1) находится по следующему алгоритму. При £ = 0 задается начальное возвышение г/о(ж) в виде непрерывной функции. Согласованная с ним начальная скорость ио(ж) вычисляется по формуле (17). В момент времени £ в точке ж находится корень уравнения (21). По формуле (22) определяется возвышение г/(ж,¿), после чего с использованием (16) вычисляется скорость м(ж,¿).
С помощью описанного метода характеристик можно найти точное решение либо при всех £ > 0, либо только до момента градиентной катастрофы ¿*. Рассмотрим далее начальное возвышение в виде одиночной волны длины Л [20]:
1 , ! 2"^(ж жад) \ 1 + СОв (---)
Л
а
г?о(ж) = { 2
0, |ж — жад| > Л/2
|ж жад1 ^ Л/2, (23)
где жад — абсцисса вершины волны при £ = 0, а — ее амплитуда. При 0 < £ < длина волны, движущейся влево со скоростью — о, и ее амплитуда остаются постоянными и равными соответствующим значениям начальной волны (23). Профиль движущейся волны будет деформироваться со временем так, что ее передний фронт укручивается, а задний — выполаживается, поскольку характеристики ж/ = о( ж), уходящие с отрезка [жад — Л/2, жад], образуют сходящийся пучок, а с отрезка [жад, жад + Л/2] — расходящийся. Таким образом, в решении будет присутствовать волна разрежения и идущая перед ней волна сжатия, приводящая к градиентной катастрофе.
На рис. 1, а показан график точного решения у = г^(ж, ¿) при £ = 5 с (сплошная тонкая линия). Безразмерные (полученные делением на Ко) значения параметров, определяющих начальную волну, были следующими:
а = 0.2, жад = 30, Л =10, Ко = 1.
Для этих значений градиентная катастрофа наступает при ~ 5.57 с, поэтому в момент времени = 5 с передний фронт точного решения уже стал близким к вертикальному.
Отметим, что уравнение (21) имеет единственное решение р при Ь < ¿*, поскольку функция
¡(р;х, г) = р + д (ко + г/о(х - рЬ)) - 2со (24)
является монотонно возрастающей функцией переменной при всех х и г е [0,г*). На рис. 2 показан график функции /(р; х, ¿) при фиксированном значении Ь = 5 с их е [0, Ь], при этом Ь = 40 и по оси Ор откладывается безразмерная величина р = р/Со. Сплошная линия изображает линию уровня / = 0. Из рисунка видно, что уравнение = 0 действительно имеет при каждом фиксированном значении переменной х единственный корень . Для нахождения корня можно использовать какой-либо
0.20 П
0.15 Н 0.10 0.05 Н 0.00
9 12 15 18 х 21
10 20
30
40
Рис. 1. Движение одиночной волны (23) по горизонтальному дну: а — графики точного решения (1) в момент времени £ = 5 си численного, полученного по схеме предиктор-корректор (2) и противопоточной схеме (3); б — траектории узлов адаптивной сетки
-1.0 Р/со -1.2
Рис. 2. График функции (24) при £ = 5 с
а
г
5
4
3
2
1
0
0
X
метод решения нелинейных уравнений. В качестве начального интервала, содержащего корень, можно взять [р1, р2], где
На рис. 1, а также изображено численное решение, полученное с помощью двух схем на адаптивных сетках, покрывающих отрезок [0, Ь]. Адаптивные сетки строились методом эквираспределения [18] с использованием управляющей функции
с коэффициентами «о = 10, «1 = 10. Из рис. 1, б, на котором изображена динамика сетки в расчете по схеме предиктор-корректор, видно, что при = 0 узлы сгущаются в зоне начальной волны, а с ростом времени происходит непрерывное автоматическое перераспределение узлов с усилением сгущения в зоне переднего укручивающегося фронта волны. Такое поведение адаптивной сетки приводит к повышению точности расчета. Как видно из рис. 1, а, графики точного решения и численного, полученного с помощью схемы предиктор-корректор, визуально неразличимы, хотя количество узлов сетки было небольшим ( N = 100). Использование неподвижных равномерных сеток с тем же числом узлов ведет к существенному понижению точности расчетов, что согласуется с выводами работ [3, 18, 19], в которых рассматривалось решение скалярных уравнений на адаптивных сетках.
Что касается противопоточной схемы, то применение адаптивных сеток также приводит к повышению точности расчетов по сравнению с равномерными сетками, хотя и не в такой степени, как для схемы предиктор-корректор (см. рис. 1, а, на котором, в частности, различимо размазывание границ волны и падение ее амплитуды при использовании противопоточной схемы).
5. Метод адаптивных сеток в задачах наката волн на берег
В работе [22] метод адаптивных сеток применялся для расчета наката волн в одномерных сечениях, проведенных от некоторой изобаты до выбранной изолинии на суше. Для расчета движения точки уреза использовались аналитические решения уравнений мелкой воды [20], представленные в виде бесконечных рядов для одного из возможных случаев взаимодействия волны с берегом и в виде решения системы обыкновенных дифференциальных уравнений для двух других. При численном решении [22] использовались несколько первых членов из полученных рядов, а для решения систем ОДУ применялся метод Эйлера второго порядка точности. В настоящей работе предлагается алгоритм расчета, в котором используются большее число членов в частичных суммах рядов и более точные методы решения обыкновенных дифференциальных уравнений.
5.1. Уточненные формулы для расчета движения точки уреза
Рассмотрим задачу о накате волн на берег, в которой требуется искать неизвестную функцию ж = жо(£), описывающую закон движения точки уреза, и решение системы уравнений (1) при ж ^ жо(£). Далее предполагается, что начальные данные для уравнений (1) заданы не при = 0, а при = о, т. е. известны начальное положение точки уреза жо (£о) = жоо, начальные функции при ж ^ жоо
I /у-П _ Т)'П
п л , I п I , /7 + 1 Н
«>¿+1/2 = 1 + «о | ^"+1/2 | + «1--
(25)
м(ж, i0) = м0 (ж), Н (ж, i0) = Н0(ж), (26)
и необходимо найти ж = ж0(£), «(ж, i), Н(ж, i) при i > i0. Дополнительное краевое условие в точке уреза заключается в требовании, чтобы полная глубина в этой точке равнялась нулю во все моменты времени:
Н i) = 0. (27)
В зависимости от начальных данных (26) возможны три режима взаимодействия волны с берегом [20]. Первый случай — накат на берег необрушающейся волны:
Нж U=ico(i)= (28)
Покажем, что для скорости движения точки уреза U(¿) = м(ж0(£), ¿) выполнено равенство
U(¿)=ж0(*), ¿0. (29)
Здесь и далее точка сверху используется для обозначения производной по времени от функции, зависящей только от одной переменной ¿. Введем новую независимую переменную £ = ж — ж0(^) и запишем систему (1) в новых координатах (£, ¿):
Н + (м — ж0)Н + Hwg = 0, Mt + (м — ж0+ gHg = ghg, £ ^ 0, (30)
где
К(ж) = h (ж(£, i)) = h(f, i), H(ж, i) = H(ж(£, i), i) = #(f, i), м(ж, i) = м(ж(£, i), i) = , i).
В новых переменных (£, i) точка уреза неподвижна и имеет при всех i ^ t0 одну и ту же координату = 0, при этом
Я(0, i) = 0. (31)
Используя в первом уравнении системы (30) при £ = 0 тождество (31), а также его следствие Ht(0, ¿) = 0 и условие (28), получаем, что ж0 — м|^=0 = 0, что совпадает с равенством (29).
Закон движения точки уреза ищется в виде степенного ряда по времени
ж0(*) = ж00 + ж01(* — U) + ж02 (t — ^ + ж0з (t —*0) + ж04— ^ +0 ((i — ¿0)5) . (32)
2 6 24
Из равенства (29) следует, что
(i — i )2 (i — i )3 U(i) = ж0(*) = ж01 + ж02(i — ¿0) + ж03 ( 0 0) + ж04 ( г 0) +0 ((i — h)^ , (33)
2 6
т. е. ж01 = U(¿0) = м0(ж0(¿0)) = м0(ж00) = м00, где через м00 обозначена скорость жидкости в точке уреза в начальный момент времени.
Для нахождения остальных коэффициентов можно воспользоваться рекуррентными формулами, предложенными в [20]. Однако это неудобно при численной реализации алгоритма, поскольку один и тот же член ряда необходимо аппроксимировать с разным порядком в зависимости от того, сколько членов ряда берется в частичной сумме для использования в разностном граничном условии. Более предпочтительным является использование явных формул для вычисления коэффициентов ж0& (fc > 2) через начальные данные задачи. Получим эти формулы.
Для вычисления коэффициента х02 рассмотрим второе уравнение системы (30) при £ = 0 и Ь = ¿0: и(¿0) = -дг]'0(х00). С другой стороны, из (33) следует, что (У(¿0) = х02, поэтому
х02 = -д Г]0(х00). (34)
Здесь и далее штрихом обозначается производная по х от функций, зависящих только от пространственной координаты.
Для нахождения х03 используем уравнения, полученные в результате дифференцирования первого уравнения системы (30) по переменной £, а второго — по ¿:
+ (и -0С0,)Н.^ + 2Н?щ + Нщ£ = 0, йы + (щ -х0)щ + (и -000)^1^ + дН^ = дх0Ы1. (35)
Полагая в этих уравнениях сначала = 0, а затем — = 0, приходим к формуле
х0з = и(¿0) = 0 (2и0Н0 + Щ0Ы1)
■ (36)
х=хоо
Для вычисления коэффициента х04 будем использовать уравнения (30), (35), а также следующие уравнения, получающиеся из указанных путем дифференцирования по переменным или :
Щы + (и - хс0)Н^ + (и - х0)Н^ь + 2Н^щ + 2Н?и^ + Ньи^ + Ни?0 = 0,
иш + (иЙ - х0)щ + 2(щ - хс0)щь + (и - х0)ща + дН^л = # [х0Ы" + (хх0)2Ы"], (37) Щг + (щ)2 + (и - х0)щц + = дк''.
Полагая в этих уравнениях £ = 0 и учитывая равенства (29), (31), а также их следствия и(£) = х0(£), и(Ь) = х0(£), Н4(0, ¿) = 0, получаем систему уравнений в точке £ = 0:
Щи + 2Н# щ + 2Щи& = 0, Н0 + 2Н? щ = 0, и0 + (Щ )2 + = рЫ7,
[7 + дЩы = # + (х0)2Ы"].
Исключая из последнего уравнения величину с использованием первых трех уравнений и полагая Ь = ¿0, получаем выражение для коэффициента х04 через входные данные задачи:
х04 = и(¿0) = д - п0Ы' - 2Н0 (д^ + 3(и0)2)
■ (38)
х=хоо
Если в момент времени Ь = Ь0 условие (28) не выполняется, то возможны две ситуации: либо Нх(х00, ¿0) = 0, либо Нх(х00, ¿0) = то. Первый случай соответствует накату необрушающейся волны, для которой касательная в точке уреза к свободной границе совпадает в момент времени = 0 с касательной к поверхности дна, а второй — накату обрушающейся волны. В обоих случаях закон движения точки уреза определяется как решение системы обыкновенных дифференциальных уравнений [20]
х0 = и (г), и = дЫ (х0(г)) (39)
с начальными условиями
х0(Ь 0) = х00, и (¿0) = и*, (40)
при этом и* = и00 для волны с касанием и и* = и00 - 2^/дН00 для обрушающейся волны.
5.2. Алгоритм расчета
Для численного решения поставленной задачи о накате волны используется описанная выше схема предиктор-корректор на адаптивной сетке {ж™} (j = 0,... , N), левый узел Жд которой совпадает на каждом временном слое t = i™ с подвижной точкой уреза ж0(¿™). Предполагается, что в момент времени t = i™ скорость, полная глубина и положение точки уреза известны, и требуется определить эти величины на следующем временном слое Г+1 = Г + т.
Для того чтобы с помощью схемы предиктор-корректор определить значения возвышения свободной поверхности и скорости и™+1, необходимо сначала построить сетку {ж™+1}, а для этого надо знать новое положение точки уреза, т.е. крайнего левого узла Правая граница области считается неподвижной, поэтому ж*^+1 = L. В этом узле задается разностное неотражающее краевое условие [22], с помощью которого определяются значения г^1 и
В невырожденном случае (28) для вычисления координаты ж^1 и скорости и^1 используются разностные аналоги частичных сумм рядов (32) и (33):
г2 г3 г4 г2 г3
жд+1 = жд + идт + Ж02ТТ +Ж037Г +Ж04^Г, ид+1 = ид + Ж02Т + Ж03 7Т +Ж047Г. (41)
2 6 24 2 6
Разностные формулы для коэффициентов ж0& (fc = 2, 3, 4) получаются в результате аппроксимации производных из правых частей равенств (34), (36), (38) с помощью односторонних разностей соответствующего порядка в левом узле неравномерной сетки. Из первой формулы (41) видно, что для достижения порядка точности О(г5) в определении положения точки уреза ж^1 достаточно аппроксимировать коэффициент ж02 с погрешностью О (К3), ж03 — О (К2), а для ж04 можно ограничиться формулами первого порядка аппроксимации, при этом предполагается, что закон предельного перехода имеет вид ж = const. Поэтому в формуле для коэффициента
ж02 = ^Г (42)
будем использовать конечные разности третьего порядка
-11 /70? + 18rf{ - 9??™ + 2?$ ™ -11жд + 18жд - 9жд + 2жд
_ ж™ |1 Q^n П„» О».» га = -I0 — g --I2 ' -Ч3 жп '
''9,0 = пь, , ж9,0
3к ' 3к
Как и в [18], здесь предполагается, что узлы ж™ неравномерной подвижной сетки являются образами узлов qj Е [0, 1] равномерной неподвижной сетки с шагом к = 1/ Ж при некотором преобразовании координат ж = ж(д, ¿), которое в каждый момент времени £ взаимно однозначно отображает отрезок [0, 1] на [ж0(£), £]. Используя отображение ж = ж( , ) , можно выразить производные по переменной ж через производные по д, например, = щ/жд. На дискретном уровне это означает, что задачу аппроксимации производных на неравномерной сетке ж™ можно заменить более простой задачей их аппроксимации на равномерной сетке qj, что и использовано в формуле (42). В соответствии с этим аппроксимационная формула для коэффициента (36) будет иметь следующий вид:
1
ж°3 = #(ж» )2
2и™ Н™ + I К™ _ К™ |
2 "д,0Нд,0 + и0 I Кдд,0 Кд,0 ж™
V жд,0 /
(43)
Как было сказано выше, входящие в эту формулу производные достаточно аппроксимировать односторонними разностями второго порядка:
га = -3< + 4u" - u" un = -3Щ + 4En - E2Q n = -3X" + 4X" - x" — 2h , Eq'° — 2h , Xq'° — 2h , n 2h" - 5h" + 4h" - h" n 2x" - 5x" + 4x" - x"
кяя,о к2 , хяя,о к2 ■
Поскольку производные из формулы (38) необходимо вычислять лишь с первым порядком, аппроксимационную формулу для х04 проще выписать сразу на неравномерной сетке {хгп}:
x04 — g
(Uo) hxxx,0 Vx,0hxx,0 2 Ех,0 (gVxx,0 + 3(иж,о) )
(44)
где для первых разностных производных используются формулы вида
Е n_е n
En0) — ^-0, (45)
x о n n
ji x0
а для вторых и третьей — известные формулы первого порядка на неравномерной сетке [17].
Для плоского откоса
7/ ч Г -х tg6, при 0 ^х ^ х3, , л
г — —к(х) — < т (46)
[ —1 при х3 ^ х ^ Ь,
сопряженного с горизонтальным дном постоянной глубины к0, полученные формулы существенно упрощаются в силу того, что к' — tgв, к'' — к111 — 0. В формуле (46) в > 0 — угол наклона откоса, г0 > 0 — высота "суши" в точке х — 0, х3 — (г0 + 1)е1§ б1 — абсцисса основания склона. Значение г0 подбирается так, чтобы максимальный вертикальный заплеск волны на склон не превосходил этого значения.
Как было показано выше, для второго и третьего режимов взаимодействия волны с берегом закон движения точки уреза определяется из решения системы двух обыкновенных дифференциальных уравнений (39). В работе [20] для решения этой системы использовался метод Эйлера. В настоящей работе применяется метод Рунге — Кутты четвертого порядка точности с начальными условиями х0(Ьп) — хП, и(¿п) — и*, где и* — иП для второго режима и и* — иП — 2л^дН'П — для третьего.
В случае плоского откоса (46) система (39) имеет точное аналитическое решение
2
x'n+1 — x" + ти* + —д tg0, u°n+l — и* + тд tg0 (47)
и оба указанных метода (Эйлера и Рунге — Кутты) дают одинаковые результаты.
Как и в [22], в численных расчетах конкретный режим взаимодействия волны с берегом определяется в соответствии со значением разностной производной Н'П 0 в узле х'П подвижной сетки. В [22] эта производная вычислялась с помощью односторонней разности первого порядка (45), а в данной работе используется формула второго порядка
_ — 3НП + 4НП — нп
Нх,0 -
_ Q /yn I /I /yn _ /yn
о ^-Л-11 Jb2
Если эта разностная производная удовлетворяет условию т ^ |Я™0| ^ М, где т, М — заданные числа (0 < т ^ М), то считается, что реализуется первый режим, а при выполнении условия | Я™0 | < т — второй. В качестве критерия возникновения третьего случая используется неравенство | Я™ 0 | > М.
Полученные выше формулы для определения положения подвижной точки уреза ж£+1 и скорости ее движения зависят от режима взаимодействия волны с берегом, тогда как полная глубина в крайнем левом узле сетки равна нулю для всех режимов: ЯО^1 = 0.
Лемма 3. Если положение точки уреза и ее скорость определяются по формулам (41)-(44), то схема предиктор-корректор (2), (3), (7) сохраняет состояние покоя (8).
Доказательство. При выполнении условия покоя (8) реализуется первый режим взаимодействия волны с берегом (т ^ | Я™0 | ^ М). Легко проверить, что при использовании формул (41)-(44) точка уреза в момент времени ¿га+1 имеет нулевую скорость = 0) и прежнее положение (ж^+1 = ж^). Тогда из краевого условия Я01+1 = 0 получаем, что и возвышение свободной границы останется в этой точке неизменным: = = 0. Далее доказательство такое же, как в лемме 1. □
5.3. Некоторые результаты расчетов наката уединенной волны на плоский откос
Рассмотрим результаты численного решения задачи о накате на плоский откос уединенной волны, заданной в начальный момент времени = 0 по формулам:
?7о(ж) = а0 еЬ-2 Г ^^(ж - жад Л , Я0 (ж) = к(ж) + г/0(ж), щ (ж) = тттг, (48) V 2 и0 ) Я0 (ж)
где ао — амплитуда волны, ж = жад — абсцисса ее вершины, [70 = л/#(к0 + ао), ж € [0, Ь]. Плоский склон задавался по формуле (46) со значением = 1.
В начальный момент времени вершина волны (48) была на достаточном удалении от основания склона ж5 и от правой границы ж = Ь, что гарантировало расположение основной части начальной волны над горизонтальным участком дна и независимость до некоторого момента времени картины наката волны на откос от краевых условий на правой границе. Как показали расчеты, при длительном движении над горизонтальным дном волна может обрушиться. Поэтому если требуется в рамках модели мелкой воды изучить взаимодействие недеформированной уединенной волны с откосом, то волну нельзя устанавливать в начальный момент времени слишком далеко от откоса.
При малой амплитуде ао (например, ао = 0.01) начальная уединенная волна имеет большую эффективную длину, поэтому в численных экспериментах выбирались следующие значения абсцисс вершины начальной волны и правой границы вычислительной области:
= ж^ + 30, Ь = + 30.
Для более коротких начальных волн (например, при ао = 0.1) точка ж = жад располагалась ближе к откосу и длина области была меньше:
= ж^ + 20, Ь = + 20.
На рис. 3 показаны графики максимального усиления волны (максимального за-плеска на плоский откос, отнесенного к амплитуде ао начальной волны). Маркеры 1 соответствуют уточненным формулам для расчета движения точки уреза, полученным в настоящей работе. Маркеры 4 получены при использовании тех же формул, что и в работах [20, 22], т. е. при использовании в частичных суммах (41) только первых трех членов. Для малой начальной амплитуды ао полученные результаты практически одинаковые (рис. 3, а). Для начальной амплитуды, увеличенной в десять раз, различие становится заметным, особенно при малых углах наклона откоса (рис. 3, б). Расчеты с другими амплитудами набегающих волн подтвердили вывод о том, что при малых углах наклона откоса и больших значениях ао уточненные формулы дают результаты, отличающиеся от применения прежних формул.
Во многих работах, посвященных исследованию наката волн на берег, выполняются тесты на сравнение численных результатов со значениями, полученными по аналитической формуле [21]
1/4
(49)
К
— = 2.831
ао
ао
V
Но )
Однако к ней надо относиться с осторожностью, поскольку данная формула имеет очень узкую область применимости. Согласно работам [21, 26] условиями применимости являются следующие ограничения:
^ > (0.288 в)2 Но
(50)
^ < 0.479^ в)10/9. Но
(51)
2
3
^ 4
..... ..... ..... .....
I
7 6 5
- \ : \ - V -•- 1
- > - - ~[г 4 ' - / N1 N - ^ \ -------- - - 2 3 у— 4
/Г"! ^ ~ // 1 1 N ~ // ¡4 " и 1 ^
" " -
" ¿^-^-*-
'..... ..... ..... ..... .....
12 15
12 15
9, град
а
4
3
2
0
3
6
9
0
3
6
9
Рис. 3. Максимальное усиление волны в зависимости от угла наклона откоса: ао/Ьо = 0.01 (а), 0.1 (б). 1 (уточненные формулы для расчета движения точки уреза), 4 (формулы из [20, 22]) — численный расчет для плоского откоса; 2 — аналитическая зависимость (49); 3 — численный расчет для вертикальной стенки, установленной на глубине Ь„ац = ао
На рис. 3 показаны величины (линии 2) максимального усиления волны, полученные по формуле (49). Если считать, что нижняя граница области применимости (49) определяется выражением в правой части неравенства (50), то почти вся линия 2 на рис. 3, а (за исключением небольшой верхней части, соответствующей углам в окрестности (9=1°) лежит в области применимости.
Известно, что величина максимального усиления Д/а0 при накате уединенной волны на вертикальную стенку больше двух [27], а на плоский склон — больше, чем на вертикальную стенку [14]. Однако для углов, больших 11°, формула (49) дает значение меньше двух (рис. 3, а, линия 2), т.е. даже меньше, чем при накате на вертикальную стенку (маркеры 3). Это вызывает сомнение в физической состоятельности данной формулы при а0/к0 = 0.01 и б1 > 11°. Соответственно, аналитическая формула (49) применима только для малых углов (9, и тогда в ограничении (50) знак ^ означает, что выражение в правой части должно быть существенно меньше выражения слева. Таким образом, даже при такой малой начальной амплитуде набегающей волны можно делать сравнение с аналитической формулой лишь для небольшого диапазона углов , который невозможно определить с помощью неравенства (50).
Если такой диапазон существует, то графики на рис. 3, а показывают, что улучшенные формулы для расчета движения точки уреза (маркеры 1) дают численные значения заплесков, чуть более близкие к аналитическим, чем формулы из [20, 22] (маркеры 4). Для амплитуды а0/к0 = 0.1 аналитическая формула (49) не применима ввиду ограничений (50), (51), поэтому не имеет смысла использовать ее для сравнения с численными результатами. Чтобы подчеркнуть это, аналитическая зависимость (49) изображена на рис. 3, б штриховой линией 2.
Интересной особенностью поведения зависимости максимального усиления волны Д/а0 от угла в является немонотонный характер ее поведения при больших амплитудах (маркеры 1 на рис. 3, б). Причиной немонотонности является обрушение волны
П
0.4
0.2
0.0 40
80 60 £
п
0.08
0.06
0.04-
0.02-
0.00
2 1
- ; з
■ 1 . 1 . 1 1 1 и 11 1» 11
1 " 1 - 1 - 1 1 11 11 1 1 1 *
1 1 1 1 ! % ч _ / 1
60
80
100
120
140
Рис. 4. Накат уединенной волны (48) с амплитудой 00/^0 = 0.1 на плоский откос с углом наклона в = 1°: а — возвышение свободной поверхности г = г)(х, £); б — возвышение свободной поверхности в начальный момент времени £ = 0 (1), в момент обрушения волны £ =
12.64 с (2), в момент начала наката £ = 29.97 с (3)
а
х
на плоском откосе. На рис. 4, а показан процесс наката и последующего отката волны на плоский откос. Заметим, что после основного отката вблизи первоначальной точки уреза образуется волна отрицательной полярности со значительной амплитудой, генерирующая впоследствии целую серию накатов-откатов с затухающими амплитудами, что физически объяснимо, но не отмечалось в предыдущих публикациях по рассматриваемой теме. Из этого же рисунка видно, что при приближении начальной волны к точке уреза ее передний фронт укручивается.
Для более детального исследования процесса наката на этой стадии на рис. 4, б показаны профили свободной поверхности в начальный момент времени, в момент обрушения волны и в момент начала наката. К моменту обрушения амплитуда волны (линия 2) стала больше начальной (линия 1), а ее передний фронт стал почти вертикальным, при этом волна еще не дошла до первоначальной точки уреза (левый кружок на оси Ох), а по склону ее вершина прошла путь, чуть больше чем треть его длины, считая от начала склона (правый кружок на оси Ох). После момента обрушения амплитуда и скорость волны уменьшаются и к точке уреза фактически подходит волна со значительно меньшей высотой, чем до момента обрушения, а в данном случае даже меньшей, чем первоначальная амплитуда ао (линия 3). Это и является причиной того, что при малых углах наклона в заплеск на берег может оказаться меньшим, чем при больших углах. Отметим, что для малых начальных амплитуд ао такого не происходит. Здесь волна очень длинная, имеет пологие склоны и не успевает обрушиться даже при малых углах в, т. е. на длинном склоне.
Заметим также, что для рассматриваемой амплитуды ао/Но = 0.1 обрушение на плоском склоне происходит при угле его наклона в* ^ 6°. Формула (51) дает пороговое значение в* ^ 12°. Есть и другие теории, дающие иные критерии обрушения, согласно которым получаются иные пороговые значения в*. В рассматриваемом случае расхождение пороговых значений в численных и аналитических расчетах можно объяснить тем, что к склону фактически подходит уже деформированная волна с более крутым передним фронтом, чем у начальной волны (48).
Во многих работах при решении реальных задач вместо наката на берег моделируется накат на вертикальную стенку, установленную на береговой линии. Это значительно упрощает вычислительный алгоритм, но приводит к погрешностям при вычислении максимального усиления волны, а следовательно, и границ затопления суши.
В случае одномерной задачи для замены берега стенкой в начальной точке уреза хоо = z0ctgв ставится вертикальная стенка и на ней задается условие непротекания и(хоо,Ь) = 0, означающее, что подвижность точки уреза игнорируется. Дно (46) в окрестности стенки модифицируется так, что если глубина становится меньше некоторого заданного значения Н„ац, то полагается, что Н(х) = Н„ац. Иными словами, дно в окрестности стенки становится горизонтальным с глубиной Н„ац, т.е. дно "срезается" и выравнивается в окрестности точки хоо. Там, где глубина больше чем Н„ац, форма дна остается прежней.
Были выполнены многочисленные расчеты задачи наката волны на вертикальную стенку. Некоторые результаты представлены на рис. 3 (маркеры 3). Видно, что, как и в случае задачи наката на берег, зависимость максимального усиления волны К/ао от углов наклона откоса имеет монотонный характер при малых амплитудах и немонотонный — при больших. В последнем случае причина по-прежнему заключается в обрушении волны на откосе еще до подхода к вертикальной стенке. При малых амплитудах и больших углах откоса значения максимального усиления волны при накате на стенку
Рис. 5. Мареограммы, полученные аналитически (линии) и численно (маркеры) в точках, расположенных справа от начального положения точки уреза на расстояниях 0.25 (1, 3) и 9.95 (2, 4)
и берег практически неразличимы, а для больших амплитуд они различаются даже при больших углах наклона. Самое большое отличие наблюдается при 9 ^ 6° (см. рис. 3, б). Здесь накат на берег значительно больше, чем на вертикальную стенку.
К настоящему времени сложился большой набор тестовых задач [28], принятых в качестве обязательных тестов для проверки достоверности результатов, получаемых с помощью численных алгоритмов расчета наката волн на берег, а также для сравнения характеристик различных алгоритмов [29]. Один из таких тестов предусматривает сравнение численных результатов с точным аналитическим решением задачи о накате уединенной волны малой амплитуды ао/Но = 0.019 на пологий плоский откос (46) с углом наклона 9 = аге1ап(1/19.85) ~ 2.88°.
На рис. 5 приведены мареограммы, зафиксированные двумя виртуальными мареографами, установленными на расстояниях 0.25 и 9.95 справа от начального положения точки уреза хоо. Для первого мареографа (х = хоо + 0.25), расположенного вблизи начальной точки уреза, имеется промежуток времени, когда происходит осушение дна. На графике этот факт изображается горизонтальной линией, проведенной на высоте ^ ~ -0.012594. В точке расположения второго мареографа (х = хоо + 9.95) осушения не происходит. Из графиков видно, что численные результаты, полученные с помощью схемы предиктор-корректор, визуально неотличимы от аналитического решения в фазе наката и в большей части фазы отката. Отличие можно увидеть только во временном промежутке от Ь = 75 до Ь = 100, где Ь — безразмерное время, полученное делением размерного на величину ¿о = л/д/Но. Это различие можно объяснить тем, что аналитическая формула для точного решения, выведенная для необрушающихся волн, на этом временном промежутке становится неверной из-за обрушения волны в конце фазы отката. В самом деле, согласно выводам работы [26], обрушения волны при откате не происходит, если выполняется условие (51). Но в рассматриваемом здесь случае это условие нарушается, поскольку ао/Но = 0.019, а правая часть неравенства (51) принимает значение, равное 0.0173.
Заключение
Рассмотрен усовершенствованный метод адаптивных сеток для численного решения задач о распространении и накате на берег поверхностных волн, описываемых в рамках модели мелкой воды в одномерном приближении. Представлен модифицированный алгоритм реализации явной схемы предиктор-корректор, основанный на новом способе вычисления правой части уравнений мелкой воды на неровном дне, дающий при сохранении порядка аппроксимации экономию времени счета по сравнению с известной версией алгоритма [3] и гарантирующий сохранение схемой состояния покоя жидкости при переходе с одного временного слоя на другой. Предложен новый способ выбора схемных параметров, гарантирующий выполнение TVD-свойства для усовершенствованной схемы предиктор-корректор. Предложен способ получения различных консервативных схем на подвижных сетках, основанный на выборе соответствующих значений схемных параметров для схемы предиктор-корректор, представляющей собой каноническую форму двухслойных явных схем для уравнений мелкой воды. В качестве примера приведена консервативная противопоточная схема на подвижной сетке в дивергентной и недивергентной формах записи.
На основе известных аналитических решений уравнений мелкой воды в окрестности границы вода — суша разработаны уточненные разностные краевые условия в подвижной точке уреза, аппроксимирующие аналитические решения с большей точностью, чем использованные ранее [20]. Доказано, что при использовании полученных краевых условий разностная схема предиктор-корректор на адаптивной сетке сохраняет состояние "спокойной воды" во все моменты времени, если в начальный момент жидкость находилась в состоянии покоя и имела невозмущенную свободную границу. Это одно из преимуществ разработанных краевых условий по сравнению с известными методами сквозного счета, для которых сохранение состояния покоя в задачах наката является, как правило, проблематичным.
Полученные результаты найдут применение при решении двумерных задач в рамках классической модели мелкой воды, а также в алгоритмах решения нелинейно-дисперсионных уравнений [23, 24].
Благодарности. Исследование выполнено при поддержке Российского научного фонда (проект № 14-17-00219).
Список литературы / References
[1] Шокин Ю.И., Федотова З.И., Хакимзянов Г.С. Иерархия моделей гидродинамики длинных поверхностных волн // Докл. АН. 2015. Т. 462, № 2. С. 168-172.
Shokin, Yu.I., Fedotova, Z.I., Khakimzyanov, G.S. Hierarchy of nonlinear models of the hydrodynamics of long surface waves // Doklady Physics. 2015. Vol. 60, No. 5. P. 224-228.
[2] Toro, E., Garcia-Navarro, P. Godunov-type methods for free-surface shallow flows: A review // Journal of Hydraulic Research. 2007. Vol. 45, No. 6. P. 736-751.
[3] Хакимзянов Г.С., Шокина Н.Ю. Метод адаптивных сеток для одномерных уравнений мелкой воды // Вычисл. технологии. 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.)
[4] Vazquez-Cendon, M.E. Improved treatment of source terms in upwind schemes for the shallow water equations in channels with irregular geometry // Journal of Computational Physics. 1999. Vol. 148, No. 2. P. 497-526.
[5] Вольцингер Н.Е., Клеванный К.А., Пелиновский Е.Н. Длинноволновая динамика прибрежной зоны. Л.: Гидрометеоиздат, 1989. 272 c.
Voltsinger, N.E., Pelinovskii, E.N., Klevannyi, K.A. The Long-wave Dynamics of the Coastal Zone. Leningrad: Gidrometeoizdat, 1989. 272 p. (in Russ.)
[6] LeVeque, R.J. Numerical Methods for Conservation Laws. Berlin: Birkhauser Verlag, 1992. 214 p.
[7] Toro, E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics: a Practical Introduction. Berlin: Springer-Verlag, 2009. 724 p.
[8] Lynch, D.R., Gray, W.G. Finite element simulation of flow in deforming regions // Journal of Computational Physics. 1980. Vol. 36. P. 135-153.
[9] Hou, J., Liang, Q., Zhang, H., Hinkelmann, R. An efficient unstructured MUSCL scheme for solving the 2D shallow water equations // Environmental Modelling & Software. 2015. Vol. 66. P. 131-152.
[10] Hou, J., Liang, Q., Simons, F., Hinkelmann, R. A stable 2D unstructured shallow flow model for simulations of wetting and drying over rough terrains // Computers & Fluids. 2013. Vol. 82. P. 132-147.
[11] Duran, A., Marche, F. Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms // Computers & Fluids. 2014. Vol. 101. P. 88-104.
[12] Funke, S.W., Pain, C.C., Kramer, S.C., Piggott, M.D. A wetting and drying algorithm with a combined pressure/free-surface formulation for non-hydrostatic models // Advances in Water Resources. 2011. Vol. 34, No. 11. P. 1483-1495.
[13] Karna, T., de Brye, B., Gourgue, O., Lambrechts, J., Comblen, R., Legat, V., Deleersnijder, E. A fully implicit wetting-drying method for DG-FEM shallow water models, with an application to the Scheldt Estuary // Computer Methods in Applied Mechanics and Engineering. 2011. Vol. 200, No. 5-8. P. 509-524.
[14] Li, Y., Raichlen, F. Non-breaking and breaking solitary wave run-up // Journal of Fluid Mechanics. 2002. Vol. 456. P. 295-318.
[15] Flouri, E.T., Kalligeris, N., Alexandrakis, G., Kampanis, N.A., Synolakis, C.E.
Application of a finite difference computational model to the simulation of earthquake generated tsunamis // Applied Numerical Mathematics. 2013. Vol. 67. P. 111-125.
[16] Шокин Ю.И., Бейзель С.А., Рычков А.Д., Чубаров Л.Б. Численное моделирование наката волн цунами на побережье с использованием метода крупных частиц // Матем. моделирование. 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.
[17] Самарский А.А. Теория разностных схем. М.: Наука, 1989. 616 с.
Samarskii, A.A. The Theory of Difference Schemes. USA: Marcel Dekker, Inc., 2001. 788 p.
[18] Хакимзянов Г.С., Шокина Н.Ю. Некоторые замечания о схемах, сохраняющих монотонность численного решения // Вычисл. технологии. 2012. Т. 17, № 2. С. 78-98. Khakimzyanov, G.S., Shokina, N.Yu. Some notes on monotonicity preserving schemes // Computational Technologies. 2012. Vol. 17, No. 2. P. 78-98. (in Russ.)
[19] Shokina, N.Yu. To the problem of construction of difference schemes on movable grids // Russian Journal of Numerical Analysis and Mathematical Modelling. 2012. Vol. 27, No. 6. P. 603-626.
[20] 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.
[21] Synolakis, C.E. The runup of solitary waves // Journal of Fluid Mechanics. 1987. Vol. 185. P. 523-545.
[22] Бейзель С.А., Шокина Н.Ю., Хакимзянов Г.С., Чубаров Л.Б., Ковырки-на О.А., Остапенко В.В. О некоторых численных алгоритмах расчета наката волн цунами в рамках модели мелкой воды. 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.)
[23] Гусев О.И., Шокина Н.Ю., Кутергин В.А., Хакимзянов Г.С. Моделирование поверхностных волн, генерируемых подводным оползнем в водохранилище // Вычисл. технологии. 2013. Т. 18, № 5. С. 74-90.
Gusev, O.I., Shokina, N.Yu., Kutergin, V.A., Khakimzyanov, G.S. Numerical modelling of surface waves generated by underwater landslide in a reservoir // Computational Technologies. 2013. Vol. 18, No. 5. P. 74-90. (in Russ.)
[24] Шокин Ю.И., Бейзель С.А., Гусев О.И., Хакимзянов Г.С., Чубаров Л.Б., Шо-кина Н.Ю. Численное исследование дисперсионных волн, возникающих при движении подводного оползня // Вестник ЮУрГУ. 2014. Т. 7, № 1. С. 121-133.
Shokin, Yu.I., Beisel, S.A., Gusev, O.I., Khakimzyanov, G.S., Chubarov, L.B., Shokina, N.Yu. Numerical modelling of dispersive waves generated by landslide motion // Bulletin of the South Ural State University. 2014. Vol. 7, No. 1. P. 121-133. (in Russ.)
[25] Рождественский Б.Л., Яненко Н.Н. Системы квазилинейных уравнений и их приложения к газовой динамике. М.: Наука, 1978. 688 с.
Rozhdestvenskiy, B.L., Yanenko, N.N. Systems of Quasilinear Equations and Their Application to Gas Dynamics. Moscow: Nauka, 1978. 688 p. (in Russ.)
[26] Pedersen, G., Gjevik, B. Run-up of solitary waves // Journal of Fluid Mechanics. 1983. Vol. 135. P. 283-299.
[27] Cooker, M.J., Weidman, P.D., Bale, D.S. Reflection of a high-amplitude solitary wave at a vertical wall // Journal of Fluid Mechanics. 1997. Vol. 342. P. 141-158.
[28] 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.
[29] Horrillo, J., Grilli, S.T., Nicolsky, D., Roeber, V., Zhang, J. Performance benchmarking tsunami models for NTHMP's inundation mapping activities // Pure and Applied Geophysics. 2015. Vol. 170, No. 3-4. P. 1333-1359.
Поступила в редакцию 15 июля 2015 г.
106
H. K. W^OKHHa, r. C. XaKHM33HOB
An improved adaptive grid method for one-dimensional shallow water equations
Shokina, Nina Yu.*, Khakimzyanov, Gayaz S.
Institute of Computational Technologies SB RAS, Novosibirsk, 630090, Russia * Corresponding author: Shokina, Nina Yu., e-mail: [email protected]
An improved adaptive grid method is considered for the numerical solution of the problems on propagation and run-up of surface waves, described by the one-dimensional shallow water model. The modified algorithm for the realization of the explicit predictor-corrector scheme is presented, which is based on the new way of computation of the right-hand side of the shallow water equations. The algorithm provides savings in computational time in comparison with its earlier version while preserving the approximation order. Also the preservation of the state of rest is guaranteed in transition from one time level to a next one. A new method for choosing the scheme parameters on the basis of the analysis of the differential approximation is suggested that guarantees the satisfaction of the TVD-property for the improved predictor-corrector scheme. The presented method for construction of different conservative schemes on moving grids is based on an appropriate choice of the scheme parameters for the predictor-corrector scheme, which represents the canonical form of the two-layer explicit schemes for the shallow water equations. As an example, a conservative upwind scheme on moving grid is provided in the divergent and non-divergent forms. The properties of the upwind scheme and the predictor-corrector scheme on dynamically adaptive grids are demonstrated for the exact solution of the nonlinear shallow water equations.
Using the known analytical solutions of the shallow water equations in the vicinity of the water-land boundary the improved difference boundary conditions are obtained at the moving waterfront point. These boundary conditions approximate the analytical solutions with a higher accuracy than the conditions used in the earlier works. It is proved that if a fluid is at rest and has a non-perturbed free boundary at the initial time moment, then the difference predictor-corrector scheme on adaptive grid preserves the state of rest at all subsequent time moments when the newly obtained conditions are used. This is one of the advantages of the developed boundary conditions in comparison with the known shock-capturing methods, where the preservation of the state of rest is usually problematic for the run-up problems. The numerical experiments have shown that for the run-up problems the substitution of a slope by a vertical wall in the initial position of the waterfront point leads to the significant change of the wave amplification in the case of very smooth slopes even if a wall embedding is small.
It is expected that the obtained results will be used for solving two-dimensional problems in the framework of the classical model of shallow water, as well as in the algorithms for solution of nonlinear dispersive equations.
Keywords: nonlinear shallow water equations, finite-difference scheme, adaptive grid, surface waves, run-up.
Acknowledgements. The study was supported by the Russian Science Foundation (Project No. 14-17-00219).
Received 15 July 2015
© ICT SB RAS, 2015
