О ЧИСЛЕННОМ МОДЕЛИРОВАНИИ ПОВЕРХНОСТНЫХ ВОЛН В БАССЕЙНЕ С ПОДВИЖНЫМИ НЕПРОНИЦАЕМЫМИ ГРАНИЦАМИ Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 24, № 4, 2019

О численном моделировании поверхностных волн в бассейне с подвижными непроницаемыми границами*

А. А. ПАллгинА, Г. С. Хлкимзянов^

Институт вычислительных технологий СО РАН, Новосибирск, Россия tkhak@ict.nsc.ru

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

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

Библиографическая ссылка: Палагина А.А., Хакимзянов Г.С. О численном моделировании поверхностных волн в бассейне с подвижными непроницаемыми границами // Вычислительные технологии. 2019. Т. 24, № 4. С. 70-107. БОТ: 10.25743/1СТ.2019.24.4.006.

Введение

Для численного решения задач о течениях жидкости со свободной поверхностью имеется большой набор алгоритмов, основанных на применении как сеточных, так и бессеточных методов. Среди этих задач наиболее трудными являются те, в которых подвижна не только свободная поверхность, но и некоторые другие фрагменты границы области, занятой жидкостью. Например, это задачи о поверхностных волнах, вызванных подводными [1, 2] и береговыми [3] оползнями, задачи исследования эффективности подвижных волнозащитных стенок и подвижных рабочих элементов волновых энергетических станций [4, 5], взаимодействия волн с полупогруженными [6] или полностью погруженными [7] в воду подвижными объектами, задачи генерации волн волнопродук-торами в лабораторных лотках и бассейнах [8, 9].

*Title translation and abstract in English can be found on page 107.

© ИВТ СО РАН, 2019.

Цели настоящей работы — анализ свойств и оценка возможностей вычислительного алгоритма, предназначенного для изучения поверхностных волн в ограниченном бассейне, часть непроницаемой границы которого может быть подвижной. За основу взят конечно-разностный алгоритм из работы [10], обобщенный здесь на задачи о потенциальных течениях жидкости при наличии подвижных непроницаемых границ. Закон движения этих частей границы считается заданным. Алгоритм опирается на использование криволинейных сеток, адаптирующихся ко всем частям границы, подвижным и неподвижным.

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

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

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

1.1. Постановка задачи в декартовой системе координат

Рассматривается плоскопараллельное потенциальное течение идеальной несжимаемой жидкости со свободной границей в ограниченном бассейне конечной глубины. Для описания течения используется декартова система координат с осью Ох, лежащей на невозмущенной свободной поверхности, и осью Ог, направленной вертикально вверх. Одно-связная область П(£), занятая жидкостью, ограничена слева криволинейной и подвижной в общем случае боковой стенкой Г^(£), а справа — неподвижной вертикальной стенкой Гг, положение которых задано уравнениями х = и х = Ь соответственно, где Ь — время. Снизу бассейн ограничен подвижным дном Г^(£), заданным функцией г = -к(х,Ь). Функции х = э(г^) и г = -к(х,Ь) считаются известными, а функция г = гц(х,Ь), описывающая свободную границу Г^(¿), является искомой.

В модели потенциальных течений жидкости со свободной границей искомые величины (потенциал скорости р и функция г/) удовлетворяют уравнению Лапласа

Ар = 0, х е ОД, (1)

кинематическому

и динамическому

Vt + Щх - V = 0, x G rf (t)

I I2

u

Pt + L2L + 9V = 0, x G rf (t)

условиям на подвижной свободной границе, условиям непротекания на подвижных ле-

вой стенке

и дне

St + VSZ - u = 0, x G г (t)

ht + uhx + v = 0, x G r^(i),

(5)

а также краевому условию на правой стенке Гг. Если правая граница является непроницаемой для жидкости, то на ней ставится условие непротекания

u(L,z,t) = 0, (L,z) G Г, t > 0.

(6)

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

гц(х, 0) = щ(х), u (x, 0) = uo(x),

(7)

по которому однозначно определяются [10] начальные значения для потенциала.

В приведенных выше формулах использованы следующие обозначения: x = (x,z); g = const — ускорение свободного падения; u — вектор скорости с компонентами u, v в направлении осей Ox, Oz соответственно; u = Vp = (<рх, pZ), |u|2 = p2x + p2y.

Подчеркнем, что в условиях (4), (5) во все моменты времени известно положение подвижных дна и левой стенки (т.е. функции z = -h(x,t) и x = s(z,t) заданы). Этим приведенная здесь постановка задачи отличается от постановки одной частной задачи из работы [11], в которой во все моменты времени была известна скорость точек дна, подвижного вследствие перемещения подводного оползня по плоскому откосу.

Далее будут использоваться безразмерные переменные

х

X

ho''

ho

z ~ h

ho' ho'

и

s

ho' v

и =

Vgho'

V =

Vgho'

ч>

JL L = —

ho' ho'

hoVghôJ

где к0 — некоторая характерная глубина, символом " обозначены безразмерные величины (во многих последующих формулах этот символ опущен).

1.2. Постановка задачи в подвижной системе координат

Область течения 0,(1) меняется со временем, поэтому при проведении расчетов используются подвижные сетки. Чтобы построить конечно-разностную схему на подвижной криволинейной сетке, вначале делается переход к новой постановке задачи — постановке в подвижной криволинейной системе координат, в которой все участки границы Г области 0(1) лежат на координатных линиях первого или второго семейства. Пусть преобразование координат

ж = г = г^1^2^) (9)

устанавливает в каждый момент времени Ь взаимно-однозначное непрерывно дифференцируемое соответствие между исходной (физической) областью 0(Ь) = 0(Ь) и Г и вычислительной областью в пространстве переменных д1, д2 — единичным квадратом С} = [0,1] х [0,1]. Будем предполагать, что боковые стенки и Гг области 0(1) являются образами при отображении (9) левой ^ и правой 7г сторон квадрата С}, нижняя сторона 'уь отображается на дно Г^), а верхняя ^^ — на свободную границу Гf(¿). Границу квадрата С} будем обозначать через 7 = 71 и 7г и 7ъ и ^^. Пусть Ц = С} \ 7.

В новой системе координат уравнение Лапласа (1), кинематическое (2) и динамическое (3) условия записываются с учетом формул обезразмеривания (8) в виде [10]:

где

W (*"W + ^© + if W + **© = 0 q = ^ ё (10)

Vt + v1r]qi — v = 0, q ё 7f, (11)

lui2

fit -uxt — vzt + — + T) = 0, q ё "Tf, (12)

22 2 1 h 1 = —J, k\2 = K21 = ——, K22 = -J, (13)

OO О О / \

9n = xqi + zqi, gi2 = 921 = xqixq2 + zqi zq2, g22 = xq2 + zq2, (14)

J = xqi zq2 — xq2 Zqi — якобиан преобразования (9), J > 0, первая контравариантная компонента v1 скорости определяется выражением

1 dq1 dq1 dq1 dq1 —xtzq2 + ztxq2 zq2 xq2

v = ~dt = ~ât +u^x+ v-ш =-J-+ uT — vir, (15)

при этом декартовы компоненты скорости вычисляются через потенциал по формулам

figl Zg2 — fiq2 Zgl — figl xg2 + fig2 xgl

u = ——-—-—-——, v =--——-——. (16)

J J

Замечание 1. В формулах (10)-(16) мы использовали для обозначения функций, зависящих от новых координат q1, q2, те же обозначения, что и в декартовых координатах x, z. Строго говоря, это разные функции, и следовало бы учитывать этот факт в обозначениях, например, добавлением черты сверху. Связь между значениями этих функций определяется заданными формулами преобразования координат. Например,

fi(q1, q2, t) = fi(x(q1, q2, t), z(q1, q2, t), t). (17)

Формула (17) означает, что в момент времени £ потенциал скорости в точке я принимает такое же значение, как в точке х, в которую точка я переходит при отображении (9) вычислительной области Q на физическую область 0,(1). Для упрощения обозначений черта над функциями в новых координатах будет опускаться, если это не приводит к недоразумениям.

Ранее [11] краевые условия для потенциала на подвижных частях границы ставились исходя из предположения, что известен вектор скорости каждой точки границы. Однако в задачах с криволинейными подвижными фрагментами границы известна, как правило, не скорость, а положение этих участков границы в каждый момент времени. Выведем для этого случая краевое условие для р. Вначале покажем, что из условия непротекания (5) следует равенство нулю второй контравариантной компоненты скорости в точках я е %, которая вычисляется по формуле

2 _ ¿д2 дд2 дд2 дд2 Хггд1 — хд1 хд1

(И дЬ дх дг 3 3 3 Используя формулу г = — Н(х,Ь) и равенства

г(д1, 0,г) = — Н (х(д1, 0,г),г) = —Цд1 ,г), (д1, 0) е ъ,

получаем выражения

= —Н = —Н — НхХг, V = —Н„1 = —кхх„1,

'18)

19)

(20)

подставляя которые в (18), приходим к равенству V2 = (Н + икх + ь) хд1, я е Отсюда с учетом (5) получаем, что

0.

(21)

1ь

Подставляя в (18) выражения (16) для декартовых компонент скорости и учитывая обозначения (13), (14), а также выражения (20) и (21), приходим к следующему краевому условию на подвижном дне для потенциала скорости:

к +к -од1 од2 я2=о

(х&д 1 — Хд1 Н г) = иь(д1,г).

д2=0

(22)

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

к 9* + к 9* од1 од2

д1=0

д1=0

при этом

X

и

(0,д2,1) = з ^(0,д2,1),г) = з(д2,1), (0,д2) е ц

= 0.

и

Условие непротекания (6) эквивалентно равенству

V

0

1г

и записывается в новых координатах в следующем виде:

к 9* + к 9* М1ТГТ +

од1 од2

(23)

(24)

(25)

(26) (27)

Я1 = 1

2

V

2. Численный алгоритм решения задачи в области с подвижными фрагментами границы

После перехода к новой системе координат необходимо построить конечно-разностную схему на равномерной прямоугольной сетке ((ь, покрывающей вычислительную область. Пусть ка (а = 1, 2) — шаги этой сетки в направлении осей Ода, ка = 1/Ма, Я5 — ее узлы, j = (]г, ]2) — мультииндекс, ]а = 0,... , Ма.

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

Пусть при Ь = ИП криволинейная сетка хП построена и на ней вычислены значения сеточных функций и . Вычисление решения на (п + 1)-м слое по времени, т. е. в момент времени = ^ + тп, состоит из нескольких этапов.

Этап 1. Новый шаг по времени. Сначала вычисляется шаг по времени

Тп кз

тш

П=0,...,Ы-1-1

|™П — гп I

\ХП + 1,М2 х31 ,М2\

(28)

где кзап — коэффициент запаса, 0 < кзап < 1. В пользу такого способа определения шага по времени указывают результаты теоретического анализа устойчивости линеаризованной задачи (см. ниже). Кроме того, опыт решения некоторых нелинейных задач с подвижными границами также подтверждает правомерность применения формулы (28) для достижения устойчивости вычислений. Нетрудно показать, что в случае сетки с вертикальными координатными линиями поверхностная волна проходит за время (28) путь, не превосходящий минимального расстояния между соседними координатными линиями второго семейства.

Этап 2. Новые значения потенциала на свободной границе. Далее вычис-

П + 1 и

ляются значения потенциала ^ в сеточных узлах ^^ Е , лежащих на верхней стороне квадрата ( . Для этого используется разностная аппроксимация динамического условия (12):

~п+1

п

^ ^ + Ч

+ ТГП =0, j = (л,^),

31

1,

мл - 1,

(29)

где

х

У

ХП - хп-1

Тп-1

(30)

тп-1 — шаг по времени между временными слоями 1п-1 и ¿П, а иП и ьП вычисляются на

основе конечно-разностных аналогов формул (16):

иП

(

Ъд1 -

Рд2 Хд1

3

)П, *=(

Уд1Хд2 + фд2Хд1

3

)П

(31)

при этом производные по переменной 1 (вдоль верхней стороны квадрата ( ) аппроксимируются с использованием центральных разностей, а по 2 — односторонних разностей второго порядка.

2

В верхних левом и правом углах квадрата Q вместо (29) используются другие формулы. Если учесть свойство (25), то в узле я0 щ получается следующая аппроксимация динамического условия (12):

Ю.+ - ю'П 1 , ч

--- - 1 + ^) + ^п = 0, j = (О, N). (32)

Такая же формула применяется для вычисления Ю.+1 в правом верхнем углу квадрата Ящ N, если правая стенка является непроницаемой (условие (26)). В случае свободного прохода волн через эту границу можно использовать, например, неотражающее условие из [12] при j = (N^N2):

{

п+1 , ю- + гп (у -у/КЕи?) , и? > 0

ю ' ю?, и? < О, ()

где Кц, — глубина в месте расположения правой границы х = Ь.

Этап 3. Новые значения потенциала под свободной границей. После вычисления потенциала во всех узлах множества решается смешанная краевая задача для уравнения (10) с краевыми условиями (22) на 75, (23) — на 71, с условием непротекания (27) или неотражающим условием на 7Г и условием Дирихле на 7^ с использова-

и и П+1

нием в разностной задаче значений юу , вычисленных в узлах Яу Е на предыдущем этапе. Разностные уравнения для потенциала скорости получены интегроинтерполяци-онным методом [10] и имеют вид

\й=0 /у

У Я Е Qн \ 7М, (34)

где юП+1 — значение сеточной функции <р в узле, имеющем локальный номер к. Локальная нумерация узлов шаблона, соответствующего Яу, введена здесь для сокращения записи. Так, согласно табл. 1 вместо глобального номера ( 1, 2) используется локальный номер к = О, вместо (]1 — 1,]2) — к = 1, вместо (]1,]2 — 1) — к = 2 и т. д.

Для внутренних узлов Яу Е Qh разностные уравнения (34) имеют девятиточечный шаблон, состоящий из узлов с локальными номерами к = О,..., 8. Коэффициенты « ( к = 1,... , 8) этих уравнений приведены в первой строке табл. 2, в которой использованы следующие обозначения:

е = ^ — ^ с = ^ + 4К2к22 + \ак12, (35)

А,В,С иВ — центры Яп-ф^-^ Яп+ф^^ Яп+1/2,Ь+1/2 и Яп-1/2,п+1/2 ячеек с общей вершиной Яу о а = = 1, = о в = —1.

Таблица 1. Соответствие между глобальными номерами j узлов сетки и локальными номерами к узлов шаблона

j к j к j к

(3ъ32) 0 (31 + 1,32) 3 (31 + 1,32 — 1) 6

(л — 1,32) 1 (31,32 + 1) 4 (31 + 1,32 + 1) 7

(3ъ32 — 1) 2 (Л — 1,32 — 1) 5 (31 — 1,32 + 1) 8

Таблица 2. Коэффициенты ак уравнения (34) для внутренних и граничных узлов

Тип узла qj а1 (Х2 а4 а5 а6 а7 аз

Чп,12 Е Qh Ы + -НА + £Б ) Св + Сс -(& + ) (а (в (с Ся

qQ,j2 Е 71,ь, 0 < 32 < N-22 0 Чб Св + Сс 0 (в (с 0

qjuQ Е 7ь,ь, 0 <31 < 0 Сс -(& + Со) 0 0 (с Ся

*qN1,h Е 7г,Ь, 0 < 32 < N2 ^ + Со -и 0 а 0 0 Ся

qQ,Q 0 0 0 0 Сс 0

0 0 -Со 0 0 0 Сд

Для граничных узлов, не принадлежащих 7f,h, шаблон уравнений (34) является шеститочечным, за исключением угловых узлов qQ 0 и qWl 0, в которых он четырехточечный. Формально разностные уравнения в граничных узлах тоже можно записать в виде (34), занулив коэффициенты ак для тех узлов девятиточечного шаблона, которые не входят в шаблоны граничных узлов (табл. 2). Тогда для всех узлов qj Е ((h \7ffr коэффициент а0 в уравнении (34) будет вычисляться по одной и той же формуле

(36)

к=5

Отметим, что разностные уравнения в узлах правой границы выписываются только в случае задания на этой границе условия непротекания. Для этого случая формулы для вычисления коэффициентов уравнения (34) приведены в строках табл. 2, отмеченных символом "*". Если на правой границе задаются неотражающие условия, то, согласно [12], вместо (34) можно использовать формулу (33) для вычисления в граничных узлах qNlJ2 Е 7ф. Ниже мы приведем другую формулу, дающую меньшее отражение волн от правой границы.

Неоднородность некоторых уравнений (34) связана с подвижностью границ бассейна. В рассматриваемых здесь задачах подвижными могут быть левая стенка бассейна и фрагмент его дна. Будем считать, что вблизи боковых стенок дно является неподвижным. Тогда функции г-3 из (34) вычисляются по формуле

Г Пк2, j = (0,л), 1 < 32 <м2 - 1,

к

, j = (0,0),

2

j = (Л, 0), 0 < л <N1, 0 в остальных узлах,

(37)

при этом в формулах (22) для щ и (23) для щ производные по касательным направлениям к участкам 7ь и 71 границы 7 аппроксимируются с использованием центральных разностей, а производные по времени — аналогов формулы (30).

Этап 4- Новое положение свободной границы. Свободная граница на (п+1)-м слое по времени определяется путем решения уравнения (11), которое аппроксимируется явной противопоточной схемой

^п1 - у! , + у?! - у!- , (у1 -I* % Ъ+х - ч! ?, _ 0

тга ' 2 • }ц +2 • }ц ' 0

Л _!,..., N1 - 1, j _(л,М2), (38)

где

(V --^-, (39)

О т ' V >

3

компоненты вектора скорости и вычисляются по формулам (31), но с использованием

новых значении потенциала фа+1.

Формула (38) используется для вычисления Г]г^+1 только для "внутренних" узлов (0 <

л

]1 < N1), лежащих на верхней стороне квадрата Q. В силу равенства (25) уравнение (11) в угловом узле ^ принимает более простой вид ^ - V _ 0 и соответственно более простой вид при _ 0 будет иметь аппроксимирующее его разностное уравнение

^-^ _ ^. (40)

Если правая стенка является непроницаемой для жидкости, то на ней выполняется условие (26), поэтому при ]1 _ N1 также будем использовать разностное уравнение (40). В качестве неотражающего условия на правой границе наряду с формулой для потенциала (33) применялось следующее "мягкое" [13] краевое условие:

С1 _ ЛПм+-1. (41)

Этап 5. Сетка на новом временном слое. Все описанные выше вычисления выполнялись на сетке х?, соответствующей п-му слою по времени Ь _ ¿га. На данном этапе вычислительного алгоритма строится новая сетка 1. Для этого можно воспользоваться методом эквираспределения [5], однако этот метод, как и все другие дифференциальные методы, требует значительных временных ресурсов, поскольку координаты

п-\-1 и

узлов х^ приходится искать, решая итерационным методом систему нелинейных уравнений. В настоящей работе мы ограничимся более простыми и быстрыми алгебраическими методами, расчетные формулы для которых будут даны ниже при рассмотрении конкретных задач. Здесь укажем лишь на то, что при построении сетки область 0,(1) делилась на две непересекающиеся части с общей неподвижной вертикальной границей, лежащей на прямой х _ х*: область о1 (¿), ограниченную слева подвижной стенкой (см. разд. 3.2), и область 0Г (1) с подвижным дном. Сверху эти области ограничены подвижной свободной поверхностью. Остальные части границ являются неподвижными, в частности, в области дно неподвижно. Вначале сетка строится на границах областей 0,1 и 0Г, а затем — внутри них с помощью простых интерполяционных методов. На этом же этапе вычислений с использованием равенств (19) и (24) определяются новые положения левой границы и дна:

¿т 1_ н^х1,1) , 32 _ 0,...,N2, (42)

Щ+1 _ к Ч , л _0,...,^. (43)

Этап 6. Пересчет. После построения новой сетки необходимо повторить вычисления на этапах 2-4 с целью согласования величин фа+1 и г]п+1 с их областью определения — сеткой х?+1. Для этого вновь решаются уравнения (29), (32), (33), в которых

декартовы компоненты скорости (31) вычисляются по-прежнему на сетке х^ и с использованием значений потенциала с п-го слоя по времени, но теперь для скорости узлов (30) используется другая формула:

х?+1 - х?

= Х-1 - (44)

Тп

При повторном решении системы линейных уравнений (34) в их правых частях (37) учитывается новое положение левой стенки (42) и дна (43), а в формулах для коэффи-

и п + 1

циентов ак используются координаты узлов новой сетки х^ .

Еще одно изменение касается компонент скорости и, V, когда они используются для вычисления контравариантной компоненты (39) при пересчете нового положения свободной границы: декартовы компоненты скорости и?, г!1 вычисляются теперь по формулам (31) с верхним индексом (п +1) как для х, так и для (р.

3. Свойства численного алгоритма и результаты расчетов

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

3.1. Накат уединенной волны на вертикальную стенку

Наиболее распространенные прибрежные сооружения содержат в качестве одного из своих составных элементов вертикальную непроницаемую стенку. Поэтому уже много лет ведется изучение особенностей взаимодействия поверхностных волн с такой преградой. Первые работы этого направления, посвященные лабораторным [14-17], численным [18-20] и аналитическим [21, 22] исследованиям взаимодействия одиночной волны с вертикальной стенкой, дополнились позднее статьями, в которых более глубокий анализ особенностей процесса наката на вертикальную стенку и взаимодействия двух уединенных волн, в частности встречного столкновения двух одинаковых солитонов, был выполнен с использованием современных измерительных приборов [23-25], новых математических моделей и вычислительных алгоритмов [26-31]. Оказалось, что, несмотря на простую форму преграды в виде вертикальной стенки, процесс взаимодействия с ней уединенной волны большой относительной амплитуды имеет сложный характер: гребень волны движется с ускорением перед и после столкновения со стенкой, при этом амплитуда волны после соударения уменьшается, отражение волны от стенки приводит к искажению ее профиля и фазовому сдвигу, потере энергии уединенной волны за счет появления после соударения вторичных волн малой амплитуды ("дисперсионного хвоста") [21], движущихся за основной отраженной волной [26, 32].

Исследования в этом направлении продолжаются с неослабевающей интенсивностью до настоящего времени. Так, в одной из недавних работ [33] используется спектральный метод высокого порядка и отмечается, что максимальный заплеск группы волн на стенку может превышать их начальную амплитуду в шесть раз. В 2017 г. численно [34], а в 2019 г. экспериментально [25] вновь подтверждено, что при накате уединенной волны большой амплитуды давление и сила воздействия на стенку имеют два пика в моменты времени, не совпадающие с моментом максимального вертикального заплеска.

Тенденцией последнего времени, особенно заметной в крупных научных центрах, является переход к экспериментальным исследованиям в больших гидроволновых лотках [35] и бассейнах [36] с созданием условий, близких к натурным. Еще в 1980 г. в работе [15] говорилось о том, что "... интерес будет представлять также проведение новых крупномасштабных исследований воздействия длинных волн на вертикальную преграду". Такие эксперименты стали в настоящее время реальностью. Например, в работе [37] воздействие на вертикальную стенку высоких волн различной формы, в том числе и обрушающихся, изучалось в лотке длиной более 200 м, глубиной 7 ми шириной 5 м. В работах [38, 39] также представлены результаты крупномасштабных экспериментов по анализу взаимодействия регулярных и нерегулярных волн с вертикальной стенкой, в частности, в очередной раз подтвержден вывод работы [15] о том, что волны большой амплитуды вызывают двойное пиковое давление на стенку.

В настоящем подразделе излагаются результаты исследования свойств описанного выше численного алгоритма на задаче о накате уединенной волны, распространяющейся над горизонтальным дном z = -h0 = const, на вертикальную непроницаемую стенку. Эта задача не относится к классу задач с подвижными непроницаемыми границами, однако она позволяет более полно, чем в общем случае, изучить свойства численного алгоритма и выбрать для дальнейшего использования наиболее рациональные варианты реализации отдельных его частей. Кроме того, здесь сделана попытка объяснения причины появления двух пиков в хронограмме давления при накате на стенку волн большой амплитуды.

3.1.1. Согласованные начальные данные для модели потенциальных течений и моделей мелкой воды

При использовании модели потенциальных течений в задачах с уединенной волной, распространяющейся над горизонтальным дном, необходимо задать начальные данные (7) так, чтобы при t > 0 волна двигалась как солитон: с постоянной скоростью и без изменения своей формы. Однако для модели потенциальных течений не известны точные решения для уединенной волны в виде конечных формул, содержащих лишь элементарные функции, поэтому уединенную волну задают приближенно с той или иной погрешностью [40, 41]. При стремлении задать начальную уединенную волну с высокой точностью можно использовать результаты работ [42, 43].

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

В качестве примера рассмотрим здесь лишь две модели цепочки: представленную выше двумерную модель потенциальных течений и одномерную нелинейно-дисперсионную (НЛД-) модель, в которой роль скорости u(x, t) играет осредненная по толщине слоя жидкости горизонтальная составляющая вектора u(x, t). Уравнения НЛД-модели в случае горизонтального дна [44, 45]

1 ( H з \

Ht + (Hu)x = 0, rit + uux + gï]x = H , (45)

где H = ho + ï], Ri = uxt + uuxx — u2, дополняются начальными условиями

r](x, 0) = rjo(x), u(x, 0) = uio(x). (46)

Для того чтобы иметь возможность сопоставления численных результатов, полученных в рамках разных моделей, необходимо задавать для них одинаковые начальные условия. Но, как видно из формул (7), (46), начальные данные различаются (возвышение и вектор скорости в (7), возвышение и скорость в (46)), поэтому можно говорить лишь о желательности некоторого согласования начальных данных для рассматриваемых одномерной и двумерной моделей.

Поясним, что мы понимаем под согласованностью начальных данных для моделей разной пространственной размерности и как строятся эти согласованные начальные данные. Пусть начальные данные (46) для одномерной НЛД-модели заданы. Будем говорить, что начальные данные (7), (46) согласованы, если:

1) для модели потенциальных течений функция r]o(x) в (7) такая же, как в (46);

2) поле вектора скорости uo(x) в (7) является потенциальным (безвихревым);

3) после осреднения горизонтальная компонента скорости u0(x) совпадает с u0(x):

т/о(х)

—— uo(x, z) dz = uo(x), Ho(x) = ho + щ(x). (47)

Ho(x) J -h0

Для построения согласованных начальных данных будем использовать формулы [46], позволяющие реконструировать (восстановить) с определенной точностью компоненты вектора скорости uo(x) двумерной задачи по начальным данным для одномерной (45):

, , _ , . f Ho (x)2 (z + ho)2 \ _ц /ч / ч / ï\—//\ uo(x, z) = uo(x) + I —6---2-)Uo(x), vo(x, z) = —(z + ho)uo(x). (48)

Очевидно, что при такой реконструкции требование (47) удовлетворяется. Кроме того, поле вектора скорости uo(x) будет потенциальным:

dvo, . duo n te(x) — (x) "

Таким образом, начальное условие (7) с компонентами (48) вектора скорости согласовано с (46).

Приведем два примера согласованных начальных данных. Пусть начальные функции в условии (46) заданы в виде

rqo(x) = aosech2(X ), (49)

lo(x) = , (50)

Ho(x)

где Х = к(х — хо), ао — амплитуда начальной волны, х = хо — положение ее вершины,

и0 = у/д(ао + ко), к = ^

ко У 4(ао + ко)'

Тогда НЛД-уравнения (45) имеют точное решение г](х, Ь) = щ(х — и0¿), и(х, Ь) = ио(х — и0Ь), описывающее уединенную волну, движущуюся с постоянной скоростью и0 над горизонтальным дном. Поскольку для функций (49), (50) справедливы равенства

и0(х) = ^Щх)Ло(х), и'0>(х) = (Но(х)'П^(х) — 2 (71о(х))2) , (51)

г}0(х) = —2кт]о(х) 1апЬ(Х), ^(х) = 2к2щ(х)(2 — , (52)

о о ао

то, в соответствии с выражениями (48), согласованные начальные данные получаются при задании компонент скорости в (7) по формулам

ио(х, х) = и0(х)

1 + 3(г + ко)2\ Но(х) (2ао — 3щ(х)) + 4 (щ(х) — ао) 'Цо(х)

+ и 4 Н2(х) )

4 4 Н2(х) ) ко(ао + ко)

ь0(х, г) = у/3а0д тр(тК(^ + ко)ЬапЪ.(Х). (53)

1 2( х)

Интересно отметить, что после элементарных преобразований формулы (53) для вычисления компонент начальной скорости совпадают с применявшимися ранее в работе [10], в которой они выведены другим способом с использованием результатов [47]. Как показали численные расчеты в рамках модели потенциальных течений, формулы (53) дают при Ь > 0 уединенную волну, движущуюся с постоянной скоростью, однако при ао/ко > 0.5 за основной волной возникает движущийся вместе с ней цуг волн малой амплитуды.

Во втором примере начальные данные (46) имеют конечный носитель и описывают одиночную волну длины А > 0:

^х) = ( !(1 + С°8(Х)), |х — х 1 ^ А/2, ио(х) = Щрх). (54)

[ 0, |х — х01 >А/2, Н(х)

Здесь использованы те же обозначения, что и в формулах (49), (50), за исключением одного: к = 2ж /А. При |х — х| < А/2 справедливы формулы (51), а вместо (52) теперь следует использовать выражения г]'0(х) = —0.5 као вт(Х), г]'(х) = 0.5 к2 (ао — 2г]о(х)). Таким образом, согласованность условий (7) с (46) будет иметь место, если компоненты начальной скорости в модели потенциальных течений вычисляются по формулам

щ(х, г) = щ(х) + ио-

ко (I (г + ко)2\ Н(х) (ао — 2щ(х))+4('Цо(х) — а0)щ(х)

(I (г + ко)2 \ 1з н2 (х) )

2

Но(х)\3 Н2(х) ) (А/ж)

Уо(х, г) = ио а°, , ^ (г + ко)ып(Х). (55)

но(х) А

Замечание 2. Если начальные данные (7) для модели потенциальных течений заданы по формулам (55) и первой из формул (54), то очевидно, что функции щ(х) и ьо(х, г) являются непрерывными. Легко проверить, что горизонтальная компонента скорости ио(х, г) имеет разрыв первого рода на линиях х = хо ± А/2, пренебрежимо малый в случае длинных волн.

3.1.2. Устойчивость разностной схемы

Для исследования устойчивости разностной схемы, представленной в разд. 2, будем использовать линеаризованный вариант задачи (1)-(3), (5) в бесконечной в горизонтальном направлении полосе, ограниченной снизу горизонтальным дном h0 = const, а сверху прямой z = 0, найти решение p(x,z, t) уравнения Лапласа

A p = 0 - то < x < -h0 < z < 0 (56)

и функцию z = Tj(x, t), описывающую свободную границу, такие чтобы выполнялись линеаризованные кинематическое и динамическое условия

Vt -v\ z=Q = 0, pt + 9V\Z=0 = 0, (57)

где v = pz, и условие непротекания на дне:

Pzl=-h0 = 0. (58)

Решением задачи (56)-(58) являются периодические волны

r](x, t) = а0 sm(ut — kx), (59)

p(x, z, t) =-0 cosh (k(z + h0)) cosiut — kx), (60)

ш cosh (kh0)

где a0 — амплитуда волны, к = 2ц/А — волновое число, А — длина волны, ш — ее частота,

ш(к) = ±y/g!i~^ta:nh(kho). (61)

Далее в этом подразделе будем использовать в (61) знак "+", что соответствует волне, движущейся вправо со скоростью

ш(к) tanh( kho) г-j- . ,

vp^ = — = ^^Т", ^ = (62)

При применении к линейной задаче и ее дискретизации только по времени схема, рассмотренная в разд. 2, может быть записана в виде (56), (58) и следующих аппроксимациях условий (57):

^(,, 0) -У"(х, 0) +^{х) = 01 гГЧ*) Г(х) „^ 0) = 0 (63)

Т т

Устойчивость этой полудискретной схемы будем исследовать на классе функций вида

рп(х, г) = рпФоегкх С08к{ф + Иа)], г]п(х) = рпЕаегкх. (64)

Эти функции являются решением задачи (56), (58), (63) при условии, что множитель перехода р удовлетворяет уравнению р2 — 2(1 — т2ш2/2)р + 1 = 0. Тогда из необходимого условия устойчивости (|р| < 1) получаем следующее ограничение на шаг по времени:

2

г < -ггт. (65)

шахш(к) к

Рис. 1. График функции ^(5), заданной формулой (67) и определяющей верхнюю границу в условии устойчивости (66)

Далее будем предполагать, что используется равномерная сетка, шаг которой в горизонтальном направлении обозначим через Ах. Самые короткие волны (64), которые можно воспроизвести на сетке с шагом Ах, имеют длину Л = 2Ах (к = ж/Ах) и выражаются формулой г]п(^Ах) = (-1)^рпЕ0. Следовательно,

таахш(к) = тах ш(к)=ш(-Ж\ = I , Х ^т0).

к К ' о<к<-к/Ах \АхУ Ах у ко \Ах)

Таким образом, условие устойчивости (65) можно записать в виде

со т/Ах < ^ (8), (66)

где введена функция

4

F ($) = -\1-т^, (67)

w Xi ж tanh(^8) v У

зависящая от параметра

5 = h0/Ах. (68)

Этот параметр характеризует степень разрешимости сетки в горизонтальном направлении по отношению к характерной глубине h0. Поведение функции F(8) (рис. 1) свидетельствует о том, что на более мелких сетках условие устойчивости (66) становится менее ограничительным.

Огрубленное условие устойчивости (66) имеет вид неравенства

^ < inf F(8) = lim F(8) = -, (69)

Ах s>о г^о ж

которое для безразмерных переменных (8) принимает вид

-

Т < -АХ. (70)

ж

Аналог (28) условия (70) использовался при решении нелинейных задач на неравномерных подвижных сетках.

3.1.3. Неотражающие условия

Для решения (59), (60) линейной задачи (56)-(58) выполняются равенство

<Рг + = 0, х = Ь, — ко < г < 0, (71)

а также кинематическое и динамическое условия (57) в точке х = (Ь, 0) и условие непротекания (58) при х = (Ь, — к0), где х = Ь — произвольное вертикальное сечение области.

Выпишем неотражающие условия на выходной границе х = Ь для нелинейной разностной задачи, требуя, чтобы аналоги перечисленных равенств и условий выполнялись и в нелинейном случае. Таким образом, будем считать, что разностное динамическое условие (29) справедливо вплоть до выходной границы. На вертикальной границе выполняется равенство XI = 0, поэтому условие (29) упрощается:

'П+1 п

ю'-:' - — ш'п I иП,

--- — ^ + ^ + ^ = 0, j = ( Мъм2). (72)

Т~п 2

Вид уравнения (38) зависит от знака компоненты скорости V1 и при з1 = N1 условие (38) не имеет смысла, если V1 < 0. По этой причине в качестве неотражающего условия при V1 < 0 будем использовать модифицированное уравнение (38):

г,п+1 _ гП ,пп _ 'П

^-^ + (ь1)' • ПМ1 . ^ -1 —У, = 0 при (V1)' > 0,

Тп пп+1_ - к1 j = ( N1,М2). (73)

' Ш! ' п / п п __^

—1-1 — ^ = 0 при (V 1)п < 0,

п

Для длинных волн фазовая скорость , входящая в качестве коэффициента в уравнение (71), мало отличается от величины с0, аналогом которой в нелинейном случае является выражение л/дН. Поэтому для нелинейной задачи вместо (71) будем использовать следующий его аналог, записанный в подвижной системе координат для безразмерных переменных:

тп+1 - тп ,-

--— — <43 + \Щип = 0, j = (N1, з2), з2 = 1,..., N2 — 1. (74)

2

Ез 1 V "3 "3

п

Наконец, последним уравнением на выходной границе будет конечно-разностный аналог условия (58):

тШ+о = тШ^. (75)

Результаты решения модельных задач показывают, что применение на выходной границе условий (72)-(75) дает существенно меньшее отражение волн при их прохождении через эту границу, чем использованных ранее условий (33), (41). Иллюстрирующий пример приведен в разд. 3.3.

3.1.4. Начальное приближение для итерационного процесса

Система уравнений (34) относительно новых значений потенциала тп+1 решается итерационным методом. Для запуска итерационного процесса необходимо задать начальное итерационное приближение т>п+1'°, например, как в [10]

тп+10 = тп. (76)

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

ф.п+1,0 = ^п ,п - , (77)

Тп-1

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

3.1.5. Исследование свойств численного алгоритма на задаче о накате уединенной волны на вертикальную стенку

Для расчета взаимодействия поверхностных волн со стенкой выбрана область длиной Ь/к0 = 20 с горизонтальным дном к(х, Ь) = к0 = 1 м. В начальный момент времени по формулам (49), (53) задавалась уединенная волна, вершина которой имела абсциссу х0/к0 = 10. Сетка была равномерной в направлении оси Ох. Вдоль вертикальных координатных линий второго семейства узлы подвижной сетки также расставлялись равномерно от дна до свободной границы на каждом шаге по времени.

Была проведена большая серия вычислительных экспериментов с целью исследования степени влияния схемных и физических параметров на свойства разработанного численного алгоритма. В настоящем разделе представлены некоторые результаты расчетов при варьировании амплитуды а0 набегающей волны, размера сетки N1 х и числа из критерия остановки итерационного процесса расчета потенциала:

\\фп+1т -фп+1т-1\\ <£. (78)

Здесь используется равномерная норма, т — номер итерационного приближения.

Расчеты показали, что выбор шага тп по формуле (28) со значением кзап = 0.95 гарантирует устойчивость вычислений для значений амплитуд а0/к0 < 0.7.

Результаты расчетов на последовательностях измельчающихся сеток указывают на сходимость численного решения во всем диапазоне амплитуд набегающей волны. В качестве примера, подтверждающего это свойство, на рис. 2, а показаны хронограммы точки пересечения свободной границы со стенкой при накате на нее уединенной волны с относительной амплитудой а0/к0 = 0.4. Видно, что графики, полученные на сетках 400 х 20 и 800 х 40, почти неразличимы. Сходимость при измельчении сетки проверена также на профилях свободной границы в заданные моменты времени и на величинах максимальных заплесков на стенку. Для волн большой амплитуды скорость сходимости при измельчении сетки падает, поэтому для численного исследования этих волн необходимо использовать мелкие сетки. На таких сетках хорошо описывается особенность, впервые отмеченная в теоретической работе [21] и характерная для взаимодействия со стенкой волн с амплитудой а0/к0 > 0.4, — образование после их отражения от стенки "дисперсионного хвоста" из коротких волн малой амплитуды (рис. 2, а).

Графики на рис. 2, б изображают зависимость от Ь количества итераций т, необходимых на каждом шаге по времени для решения системы уравнений (34) методом последовательной верхней релаксации при использовании для остановки итерационного процесса критерия (78) се = 10-8. Видно, что даже на грубой сетке ( хМ2 = 200 х 5) начальное итерационное приближение (77) дает существенный выигрыш в количестве итераций по сравнению с приближением (76). Этот вывод остается справедливым для всех амплитуд и расчетных сеток любой разрешимости.

t, с t, с

Рис. 2. Накат уединенной волны с относительной амплитудой aо /ho = 0.4 на вертикальную стенку: a — графики колебаний уровня воды на стенке, полученные на сетках с числом узлов Ni х N2 = 100 х 5 (I), 200 х 10 (2), 400 х 20 (3), 800 х 40 (4); б — количество итераций га, необходимое на каждом шаге по времени для вычисления потенциала р при задании начального итерационного приближения по формулам (76) — I, (77) — 2

В ходе численных экспериментов были обнаружены и некоторые другие свойства рассматриваемого алгоритма. Первое из них очевидно: для заданных амплитуды а0 и сетки N1 х N меньшее значение е требует большего числа итераций. Кроме того, установлено, что для заданных значения е и размера сетки чем больше амплитуда а0, тем больше потребуется итераций для вычисления рп+1, а при фиксированном значении а0 число итераций возрастает на тех расчетных шагах, когда волна взаимодействует со стенкой. Показано, что для каждой сетки и каждой начальной амплитуды а0 можно подобрать значение параметра такое, что дальнейшее его уменьшение ведет лишь к увеличению числа итераций, фактически не повышая точность численного решения на заданной сетке, например точность определения максимального вертикального за-плеска на стенку. Так, для сетки N1 х N = 100 х 5 результаты расчетов с е = 10-8 и е = 10-5 практически не различаются, однако в последнем случае число итераций уменьшается в 5-6 раз по сравнению с данными, представленными на рис. 2, б.

Одной из главных характеристик взаимодействия поверхностных волн со стенкой является величина Я максимального вертикального заплеска. На рис. 3 изображены графики зависимости величины Я от амплитуды а0 набегающей уединенной волны, полученные в настоящей работе (сплошная линия), в экспериментах (темные маркеры) и расчетах (светлые маркеры) других авторов. Для относительной амплитуды а0/к0 < 0.55 результаты расчетов максимального заплеска на основе предложенного нами алгоритма практически совпадают с результатами расчетов [19, 26, 34], а с экспериментальными данными [16] они согласуются вплоть до а0/к0 < 0.65.

Замечание 3. В работе [21] получена формула для вычисления максимального заплеска на вертикальную стенку уединенной волны малой амплитуды а = а0 /к0:

Я = 2а + \а2 + 3а3. (79)

к0 2 4

Я/И0 2.4-

■ 1

• 2

▲ 3

▼ 4

□ 5 О 6 д 7 -8

V 9

0 0.1 0.2 0.3 0.4 0.5 0.6 а0/Н0

Рис. 3. Зависимость максимального вертикального заплеска К на вертикальную стенку от амплитуды ао набегающей волны: эксперименты: 1 — [15], 2 — [14], 3 — [16], 4 — [48]; расчеты: 5 — [19], 6 — [18], 7 — [26], 8 — настоящая работа, 9 — [34]

Позже выяснилось [5, 26, 49], что эта простая аналитическая формула дает хороший прогноз даже для волн большой амплитуды. Более того, для таких амплитуд отклонение величины К из (79) от некоторых известных экспериментальных данных получается меньшим, чем при использовании многих нелинейных математических моделей.

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

р = -р\<рг + ^ + дЛ Ух е ВД, (80)

+ Т" +

где р — плотность жидкости. На рис. 4 изображены хронограммы давления в точке х = (Ь, -к0), являющейся основанием вертикальной стенки. Видно, что при набегании

Рис. 4. Хронограммы давления в основании вертикальной стенки для значений амплитуды набегающей волны ао/Но = 0.1 (1), 0.2 (2), 0.3 (3), 0.4 (4), 0.5 (5), 0.6 (6)

Таблица 3. Зависимость максимального давления ртах в основании вертикальной стенки от амплитуды ао набегающей уединенной волны

ао/ho Pmax/Pgho

Настоящая работа [34]

0.1 1.186 1.189

0.2 1.334 1.351

0.3 1.461 1.482

0.4 1.578 1.606

0.5 1.708 1.737

0.6 1.848 1.868

уединенной волны с амплитудой а0 > 0.4 в хронограмме давления появляются два локальных максимума (пика давления), при этом чем больше амплитуда волны, тем больше расхождение максимальных значений пиков. Двухпиковость давления имеет место и в других точках стенки.

К сожалению, совсем мало работ, содержащих табличные данные для давления, которые можно было бы использовать для сопоставления результатов. В статье [34] такие данные имеются, причем вычислены они довольно точным полуаналитическим методом. В табл. 3 приведены максимальные значения давления в основании стенки, полученные в настоящей работе и в [34]. Можно сказать, что имеет место хорошее соответствие, поскольку расхождение результатов не превышает 2 %.

На рис. 5 изображены хронограммы возвышения свободной поверхности на стенке и давления в ее основании при набегании волны с относительной амплитудой a0/h0 = 0.6. Видно, что первый пик давления появляется перед максимальным заплеском на стенку, а второй — после него, причем второе пиковое давление меньше первого. Эти свойства хронограмм давления обнаружены как в лабораторных экспериментах [15, 16, 24, 25, 38], так и в численных [20, 26, 29, 34, 49]. Заметим, что в экспериментах указанные свойства хронограмм давления были впервые описаны в статье Н.Н. Загрядской [15], опубликованной в 1980 г., но не попавшей в списки цитированной литературы в более поздних публикациях зарубежных исследователей, считающих своим открытие свойства двухпиковости давления (например, в статье 2019 г. [25]: "This is first proof from the experimental observation"). Такая же ситуация, к сожалению, и с численными рас-

Рис. 5. Хронограммы возвышения д/Но свободной поверхности на стенке (1), давления р/рдИо в основании стенки (2), ускорения ач/д точки, в которой свободная поверхность пересекает вертикальную стенку (3), и давления р/рдк0, вычисленного по формуле (83) (4)

четами: мы не нашли ссылок в зарубежной печати (за исключением [49]) на одну из первых работ М.И. Железняка [20], в которой достаточно подробное изучение свойств хронограмм давления выполнено с помощью численных методов.

Интересно отметить, что причина появления двух пиков в хронограммах давления до сих пор не объяснена в полной мере ни экспериментаторами, ни вычислителями. В статьях [20, 25, 26, 29] выдвигается гипотеза о том, что двухпиковость давления может быть связана с немонотонным характером ускорения частиц жидкости на стенке в фазах наката и отката. Впервые на это указано в работе М.И. Железняка [20], однако связь давления с ускорением не приведена. Выведем такую связь в случае НЛД-модели мелкой воды (45). Для этой модели давление на дне вычисляется по формуле [20, 46]

2

Н

р\,=-Но = Р9Н - РНГК1. (81)

Из уравнения неразрывности (45) следует, что

Щ + иНх

их =--Н-' (82)

поэтому Кг = иХ1 + иихх + их (Щ + иНх) /Н = ((Них)г + и (Них)х) /Н. Вновь используя (82), получим К = — (Ни + щНх + иНхЛ + и (Нйх)х) /Н.

Выписанная формула для Кг справедлива при отсутствии в потоке препятствий. Однако хорошо известно [50], что накат на вертикальную непроницаемую стенку можно рассматривать как процесс взаимодействия двух уединенных волн одинаковой амплитуды, распространяющихся в противоположном направлении. Таким образом, можно считать, что стенка, препятствующая течению, отсутствует. По этой причине формула (81) будет верной и в основании стенки х = (Ь, —ко), при этом

и1 х=(Ь,-Но) = 0' х=(Ь,-Но) = °.

Следовательно, Кг = — щ/Н и

р\х=(ь,-к0) = руН + рН^ = Р(ко + ч)(д + у). (83)

Из этой формулы видно, что давление зависит как от возвышения свободной границы г/, так и от ее ускорения а^ = щ, причем нелинейным образом. Поскольку давление (81) в НЛД-модели с большой точностью [49] воспроизводит давление в модели потенциальных течений, то можно предположить, что для последней характер зависимости (83) тоже имеет место. Косвенно это подтверждается графиком 4 на рис. 5. Хронограмма 4 изображает зависимость (83), в которую вместо возвышения г/ и ускорения а^ из НЛД-модели подставлены соответствующие величины, вычисленные в рамках модели потенциальных течений. Для этой модели ускорение а^ (Ь, Ь) имеет два локальных максимума (линия 3 на рис. 5), при этом, как и для давления (линия 2 на рис. 5), первый максимум наступает раньше максимального заплеска, а второй — позже. Отличием хронограммы ускорения от хронограммы давления является превышение второго максимума над первым. Еще одной особенностью хронограммы ускорения является наличие локального минимума примерно в тот же момент времени, когда происходит максимальный заплеск. Видно, что после достижения максимальной высоты столб жидкости начинает почти свободное падение вниз, т. е. с ускорением, немного меньшим ускорения свободного падения . Такое поведение свободной поверхности на стенке характерно только для волн большой относительной амплитуды [23, 25, 29].

3.2. Генерация поверхностных волн подвижной вертикальной стенкой

Рассмотрим теперь задачу о распространении и трансформации волн над неровным дном, состоящем из горизонтального участка на глубине К0 = 21.8 см и нескольких плоских откосов, последний из которых примыкает на глубине = 4.7 см к вертикальной непроницаемой стенке (рис. 6), расположенной на расстоянии Ь = 23.23 м от начала лотка. Координаты вершин ломаной г = -к(х), описывающей дно лотка, приведены в табл. 4. В лабораторных экспериментах [51] волна генерировалась в левой части лотка движением щита поршневого волнопродуктора. Рассматриваемая задача включена в набор обязательных тестов [52] для сравнения возможностей различных численных алгоритмов с использованием имеющихся экспериментальных хронограмм возвышения в точках установки волномеров. Последние расставлены в лотке так, как указано в табл. 5 (см. также рис. 6). Кроме того, известны измеренные в экспериментах значения максимальных заплесков на правую вертикальную стенку для случаев, когда относительная амплитуда основной волны, прошедшей через волномер В4, была равна 0.05 (вариант А), 0.3 (В), 0.7 (С).

Задача использовалась для верификации различных численных алгоритмов, например конечно-разностного метода решения НЛД-уравнений [53], конечно-разностного метода для бездисперсионных уравнений мелкой воды и БРИ-метода для уравнений

Рис. 6. Правая часть гидроволнового лотка с установленными в нем волномерами и фрагмент расчетной сетки Ы1 х А2 = 400 х 20 в момент времени £ = 15 с (вариант В)

Таблица 4. Координаты вершин ломаной, описывающей дно гидроволнового лотка

Вершина ломаной Х0 Х1 Х2 Х3 Х4

Абсцисса вершины х, м 0 15.04 19.4 22.33 23.23

Ордината вершины х, м -0.218 -0.218 -0.136 -0.116 -0.047

Таблица 5. Расстояния от волномеров В4—Вю до правой стенки лотка и глубины в местах расположения волномеров

Волномер В10 в9 в8 Вт в6 в5 В 4, вариант А В 4, вариант В В 4, вариант С

Расстояние, м 0.43 0.9 2.37 3.83 6.01 8.19 10.59 9.17 8.83

Глубина, см 8 11.6 12.6 13.6 17.7 21.8 21.8 21.8 21.8

Эйлера [54] и др., при этом во многих работах процесс генерации волн подвижной стенкой волнопродуктора не воспроизводился. Вместо этого сразу в начальный момент времени левее волномера В4 задавалась уединенная волна вида (49), (50), при этом ее амплитуда а0 подбиралась такой, чтобы в расчетах высота волны в момент ее прохода через волномер В4 как можно лучше соответствовала экспериментально измеренным значениям. При втором подходе на неподвижной левой границе, совпадающей с вертикальным сечением установки волномера В4 (см. рис. 6 и табл. 5), задавались условия, имитирующие вход через эту границу волны с формой, известной из экспериментов [51]. Недостаток этого подхода заключается в отсутствии экспериментальных данных для скорости в точке установки волномера В4, необходимых для численного моделирования.

Поскольку рассмотренный в разд. 2 алгоритм позволяет проводить расчеты в областях с подвижными границами, то имеется возможность полностью воспроизвести лабораторный эксперимент, а именно при численном моделировании создать волну с помощью подвижной стенки, движущейся по тем же законам, что и в [51]. В расчетах единственным отличием от лабораторного эксперимента было то, что движение стенки "численного волнопродуктора" всегда начиналось из точки х = 0 (рис. 7).

Подвижность левой вертикальной стенки расчетной области учитывалась за счет применения подвижных сеток, отслеживающих положение стенки на каждом шаге по времени. Наибольший сдвиг стенки 5 тах = шахз(£) равен 0.4 м (рис. 7, б, вариант С).

Исходя из этого выбиралось некоторое х*, 5тах < х* < Ь, и в подобласти (£) С расположенной между вертикальными границами х = х* их = Ь, строилась равномерная в горизонтальном направлении сетка х^ с абсциссами узлов

Ь х*

х* ¿2 =х* + (л -N1 )Дх, л = М* ,...,N1, 32 = 0,...,Ж2, Дх = Ь - , (84)

где N1 — заданное число разбиений отрезка [5п, х*], 1 < N1 < N, = п), координаты х^не меняются со временем.

Рис. 7. Траектория х = s(t) движения стенки волнопродуктора при численном моделировании: вариант A (I); B (2); C (3)

В левой подобласти С П(£), заключенной между стенкой х = и неподвижной границей х = х*, строилась неравномерная сетка, шаги которой в горизонтальном направлении изменяются по закону геометрической прогрессии, что приводит к следующей формуле для абсцисс узлов:

¿к* = ^ + ' , л = 0,...,N*1, 32 = 0,...,N2, (85)

(86)

1 - (Сга)л

„га _ „га I A VS 7

(Сnf1-1 - (Сга)

где (га — зависящий от времени корень нелинейного уравнения

1 - (C)Wl* _ ж* - sn

(£n)N*-1 - (£n)N* = Дх .

Легко проверить, что при выборе числа N**, удовлетворяющего условию

Г _

^tNI < N* <N1 - Ъ Х , (87)

Г 1 1 1 т»* _ о vy

и jy ■'max

существует единственный корень (га уравнения (86), при этом (га > 1, т. е. горизонтальные шаги сетки плавно возрастают при удалении от подвижной стенки. Кроме того, ячейки, примыкающие к сечению х = х* слева, имеют такой же шаг Ах в горизонтальном направлении, как и ячейки равномерной сетки (84) справа от этого сечения, что обеспечивает гладкую стыковку сеток, построенных в областях Qi(t) и (t).

Ординаты узлов хга вычислялись по одной и той же формуле во всей области Q(i):

пга + h"

^ = -hra + J2AZI, л = 0,...,NU J2 = 0,...,N2, =

п + hra

га _ 7 га I • А га • _ п лт ■ _ п лт А ~га _ ''Л Л

Таким образом, в подобласти П (¿) узлы сетки вначале движутся и в вертикальном и в горизонтальном направлении, а после момента остановки стенки волнопродуктора — только в вертикальном. В правой подобласти Пг (¿) они передвигаются исключительно в вертикальном направлении.

На рис. 6 изображен фрагмент сетки размером N1 х М2 = 400 х 20, построенной при х* = 2 м, N1* = 40. Для этих данных условие (87) выполняется, а из (84) следует, что Ах ~ 5.897 см. Наибольшее значение (п = 1.0217 знаменателя геометрической прогрессии (корня уравнения (86)) получается в варианте С после прекращения движения стенки волнопродуктора. В этом случае возле нее возникает минимальный в горизонтальном направлении шаг х™^ — х^2 = 2.553 см (]2 = 0,... , Ы2).

На рис. 8 изображены хронограммы возвышения свободной границы, измеренные двумя волномерами в лабораторных и численных экспериментах для варианта А движения стенки волнопродуктора. Видно, что расчеты на основе модели потенциальных течений достаточно хорошо воспроизводят экспериментальные данные, в частности уменьшение длины и увеличение амплитуды волны при ее движении в мелкой части лотка в сторону правой стенки, постепенное формирование "дисперсионного хвоста" по мере обратного движения волны над неровным дном в сторону больших глубин после отражения от правой стенки. Заметим, что в случае горизонтального дна "дисперсионный хвост" после отражения от стенки волны малой амплитуды а0/к0 = 0.05 не возникает. Таким образом, причиной образования "дисперсионного хвоста" в варианте А является неровность дна.

2, см

1.0

0.5

0.0

0

4

8

12

16

20

24 I, С

Рис. 8. Трансформация одиночной волны при ее движении над кусочно-линейным дном и отражении от вертикальной стенки: записи волномеров В4 (1, 3) и Вд (2, 4) в расчете (линии 1, 2) и эксперименте (маркеры 3, 4)

Анализ графиков на рис. 8 показывает, что расчет несколько завышает амплитуду распространяющихся волн, особенно после отражения от правой стенки. Это связано с тем, что в используемой математической модели не учитываются трение о дно и вязкость жидкости. Тем не менее рассчитанная величина максимального вертикального заплеска на стенку Д/К0 = 0.1353 находится в хорошем соответствии с экспериментальным значением К/К о = 0.1257. Отметим, что в случае горизонтального дна и а0/К0 = 0.05 значение коэффициента наката К/а0, характеризующего степень роста высоты заплеска по отношению к амплитуде набегающей волны, примерно равно двум (см. рис. 3 и формулу (79)). Для неровного дна в варианте А получается существенно большее значение К/а0 = 2.51, а в варианте В имеем значение К/а0 ~ 7, при этом волна близка к обрушению около правой стенки. В варианте С волна обрушается еще до подхода к правой стенке, вследствие чего получается меньшее значение коэффициента наката К/а0 ~ 1.8.

3.3. Взаимодействие одиночной волны с поверхностными волнами, генерируемыми подводным оползнем

В работе [55] показано, что японское мегацунами 2011 г. нельзя считать чисто сейсмическим и должен был сработать какой-то дополнительный механизм, внесший существенный вклад в образование таких аномально высоких волн, какие были зарегистрированы датчиками в открытом море и в нескольких пунктах побережья. В [55] предлагается и численно обосновывается гипотеза ("двойного источника" происхождения поверхностных волн) о том, что наряду с вертикальной подвижкой дна после землетрясения произошел мощный подводный оползень, усиливший волну цунами до катастрофической.

Мы тоже рассматриваем задачу о распространении цунами вблизи берега при воздействии на него оползневых волн, но в отличие от [55] здесь исследуется модельная ситуация, когда формы оползня и подводного склона задаются аналитическими формулами. На этой модельной задаче легче исследовать свойства численного алгоритма. Кроме того, численные эксперименты могут быть воспроизведены другими исследователями, поскольку не требуют знания данных из [55] о реальном рельефе дна и прилегающей суши, нуждаясь лишь в небольшом числе определяющих параметров.

Для расчетов выбрана северо-восточная часть акватории Черного моря вблизи места исторического оползня [56]. Реальный профиль дна в сечении по меридиану 37.5° в. д. (штриховая линия на рис. 9) был заменен модельным (сплошная линия), заданным с помощью гладкой монотонно убывающей функции, хорошо аппроксимирующей большую часть континентального склона:

/ \ - к — - к — г /

кы(х) =----1----tann[c(х — £)],

(89)

где к+ < 0, к- < 0 — предельные значения ординаты модельного рельефа соответственно в правой и левой бесконечно удаленных точках,

2 tan в0 к_-к

+

1 к — к+

£ = - 1П --;-,

2с к_ — кт

О0 — максимальный угол наклона модельного рельефа; £ — точка перегиба; кш — глубина в точке х = 0, в которой при численном моделировании устанавливалась вертикальная непроницаемая стенка. Все расчеты выполнялись при следующих значениях параметров:

к

+

-2100 м, к- = -25 м, О0 = 16°, к,, = -30 м.

(90)

Длина Ь расчетной области варьировалась от 40 до 75 км.

Модельный подводный оползень, первоначально располагавшийся на некоторой заданной глубине при Ь > 0 скользил вниз по неровному склону (89) до остановки. Нижняя граница г = — к(х, ¿) = кы(х) + кБ\(х, ¿) области П(£), занятой жидкостью, изменялась со временем в силу подвижности и деформируемости оползня, поверхность которого описывается функцией г = к3\(х, ¿) > 0. Начальная форма модельного оползня задавалась по формуле

ка\ (х, 0)

Т 2

1 + сое

^2тт(х — х0) ^

|х — х0| < Ь/2, |х — х0| > ь/2,

(91)

где х0, Т,Ь — абсцисса вершины оползня при £ = 0, его толщина и длина (вдоль оси Ох), при этом = кы(х°). При Ь > 0 форма оползня и его положение определялись путем решения уравнения движения [57], в которое в качестве параметров входят следующие

0

2, М

Рис. 9. Модельный рельеф дна, заданный аналитической формулой (89) (сплошная линия), реальный рельеф (штриховая) и невозмущенная поверхность воды (пунктирная)

величины: Cw, Cfr, C¿ — коэффициенты присоединенной массы, трения и гидродинамического сопротивления соответственно, 7 = ps/pw > 1, pw — плотность воды, ps — плотность оползневой массы, Cfr = tan в*, в* — угол трения.

Варьирование в некоторых диапазонах этих определяющих параметров, а также размеров оползня и его заглубления позволило установить характерные для модельного рельефа (89) особенности движения оползня и распространения генерируемых им волн. В качестве примера на рис. 10 изображены графики z = изменения свободной

поверхности при создании волн подводным оползнем, начинающим свое движение по склону с глубины z0 = -1500 м. Были взяты следующие значения параметров:

Т = 20 м, b = 6 км, Cw = 1, Cd = 1, 7 = 2, в* = 1o. (92)

Подводный оползень, движущийся по склону (89), генерирует сложную систему волн: к берегу распространяется волна понижения (волна отрицательной полярности), другая волна отрицательной полярности сопровождает оползень и движется вместе с ним в глубоководную часть акватории, длинная волна повышения движется в сторону моря, все более удаляясь от оползня. Отметим, что такие же волновые картины получаются и для склонов иной, чем (89), формы [2].

Кроме того, постепенно отставая от оползня, в сторону моря движется также цуг коротких волн — "дисперсионный хвост". Такой цуг волн хорошо описывается и при использовании для расчетов нелинейно-дисперсионной модели мелкой воды [2, 49], но бездисперсионная модель его не воспроизводит. Причиной образования "дисперсионного хвоста" не может быть лишь неровность подводного склона (89), поскольку возникновение этого цуга коротких волн наблюдается и в случае плоского откоса как в лабораторных экспериментах, так и в численных расчетах [58, 59]. Более того, этот цуг имеет место даже в случае движения тела по плоскому горизонтальному дну, что подтверждается недавно проведенными расчетами и лабораторными экспериментами [60]. Видимо, образование "дисперсионного хвоста" связано лишь с особенностями обтекания

a б

Рис. 10. Графики свободной поверхности в задаче о распространении волн, возникших при подводном оползне: а — применение новых неотражающих условий (72)—(75), б — неотражающих условий (33), (41) из [10]

контура оползня потоком воды. Известно, например, что при стационарном обтекании объектов, расположенных на дне, за ними возникает волнистая свободная граница — ондулярный бор [61-63].

На задаче о генерации волн подводным оползнем численно исследовались неотражающие условия, предложенные в подразд. 3.1.3. График поверхности г = г](х, ¿), изображенной на рис. 10, а, получен при использовании на выходной границе Ь = 40 км новых неотражающих условий (72)-(75), а на рис. 10, б — условий из [10]. В последнем случае заметно сильное отражение от границы. При использовании новых неотражающих условий волны как положительной, так и отрицательной полярности свободно проходят через правую границу х = Ь, при этом картина течения при х < 40 км визуально совпадает с той, которая получается при проведении расчета в более протяженной области с Ь = 75 км, что также подтверждает высокое качество неотражающих условий

(72)-(75).

На этой же задаче выявлена сильная восприимчивость алгоритма к такому свойству криволинейной сетки, как вытянутость ее ячеек. Например, если при Ь = 40 км использовать равномерную в направлении оси Ох сетку с числом ячеек N1 = 1000 и равномерную по г с N2 = 30, то в начальный момент времени, когда г/ = 0, отношение Я3(х) постоянного шага Ах = Ь/N1 к переменному шагу Дг(х) = к(х, 0)/N2 будет изменяться согласно данным (90) от 40 около левой границы до 4/7 возле правой. Таким образом, вблизи берега ячейки сетки являются сильно вытянутыми в горизонтальном направлении. Такая сетка приводит к сильному замедлению скорости сходимости итерационного процесса для вычисления потенциала скорости и понижению точности вычислений [5]. Численные эксперименты показали, что удовлетворительные результаты получаются только при условии 0.1 < Я3(х) < 10 (0 < х < Ь).

Для предотвращения сильной вытянутости ячеек сетки мы применяем следующий прием при ее конструировании. Ординаты узлов х-1 изменяются со временем и вычисляются при каждом п по формуле (88). Абсциссы узлов стационарны, при этом первый шаг сетки в горизонтальном направлении Дх1/2 полагается равным вертикальному в состоянии покоя, поэтому = Дх1/2 = А= (]2 = 0,... , N2). Для других ячеек сетки шаги в горизонтальном направлении плавно возрастают по закону геометрической прогрессии со знаменателем £ > 1. Следовательно, для вычисления абсцисс узлов получаются следующие формулы:

1 — СЛ

х31,32 = Дх1/2 1 — ^ , к =0,...^Ъ 32 = 0,...,N2,

где — корень уравнения,

1 — Ь

(93)

•ш •

1 — С Дх1/2

Это уравнение имеет единственный корень £ > 1 при условии N1 < N2 ■ Ь/к На рис. 11 изображен график отношения К3 шагов (х^1,о — х^1-1,о)/Дг<°1 Сл = 1,... , для покоящейся в начальный момент времени воды при указанных выше значениях N1, N и параметров (90). Для этого случая имеем = 1.0054, а максимальная вытяну-тость ячеек К3, имеющая место около правой границы, равна 3.08, что допустимо для используемого алгоритма с точки зрения его экономичности и точности вычислений.

Для того чтобы выяснить степень влияния волн, генерируемых подводным оползнем, на величины максимальных заплесков на берег одиночной волны, проведена боль-

2, М

0

-500 -1000 -1500 -2000

Я

8000

^ - ' '

** -1

/ Ч ' \ ' N | \ - ' ' ----

ч ч

16000

24000

32000 ^ „ 40000 х, м

Рис. 11. Модельный рельеф дна (1) и график отношения К3(х) горизонтального размера ячейки сетки к вертикальному (2)

шая серия экспериментов, в которых варьировались значения всех параметров, определяющих оползень и начальную одиночную волну. Последняя задавалась по формулам (55) и первой из формул (54), в которых полагалось = -кы(х0).

Некоторые результаты расчетов представлены на рис. 12. На рис. 12, а показана динамика свободной поверхности при накате на берег одиночной волны с параметрами а0 = 1 м, Л =10 км, х0 = 20 км. Видно, что неровности дна приводят к распаду приближающейся к берегу одиночной волны на две, в силу чего возникает двойной за-плеск на берег. На рис. 12, б показано изменение со временем свободной границы при сходе оползня, заданного параметрами (92), с глубины = -600 м. Видно, что качественная картина возникших волн идентична той, что показана на рис. 10, а. Наконец, на рис. 12, в при тех же значениях определяющих параметров изображена динамика свободной поверхности при одновременном движении оползневых волн и одиночной. Видно, что в результате взаимодействия одиночной волны с головной оползневой волной и с оползневым "дисперсионным хвостом" волновая картина стала более сложной, чем в первых двух случаях, и к берегу подходит пакет волн малой длины. В итоге на берегу наблюдается последовательность чередующихся накатов и откатов, при этом как и в случае единственного оползневого источника при "двойном источнике" первой к берегу подходит волна понижения.

На рис. 13 изображены графики максимальных значений вертикальных залесков при накате одиночной волны, оползневых волн и волн от "двойного источника". Эти графики получены при значениях параметров (92) и а0 = 1 м, Л =10 км. Варьируемыми параметрами была глубина с которой начинал свое движение оползень, и удаленность х0 вершины одиночной волны от берега. В этой серии расчетов получено, что с какой бы глубины ни стартовал оползень, он оказывает тем меньшее влияние на величину максимального заплеска, чем дальше от берега находилась первоначально одиночная волна. Связано это с тем, что дальняя одиночная волна накатывает на берег много позже момента взаимодействия оползневой волны, когда последняя фактически не оказывает уже никакого влияния на процесс наката одиночной волны. И наоборот, если одиночная волна находится первоначально в таком положении, что она подходит к берегу по волне понижения, сгенерированной оползневой волной, то заплеск на берег оказывается значительно большим, чем при накате только одиночной волны. Эта ситуация схожа с накатом на берег классической Ж-волны с лидирующей волной

0

понижения: именно в этом случае величина заплеска ^-волны значительно превосходит заплеск одиночной волны положительной полярности [64].

40000 0

Рис. 12. Графики свободной поверхности при накате на берег одиночной волны (а); волн, возникших при подводном оползне (б); волн от "двойного источника" (в)

Таблица 6. Значения максимальных вертикальных заплесков при накате на берег поверхностных волн от различных источников

х0, км

Накат 20 25 30 35 40

*0, м ж0, км оползневой Накат одиночной волны, м

волны, м 4.4 3.8 3.5 3.3 3.2

Накат волн от "двойного источника" ,м

-600 9.2 2 6.1 6 5 4 3.9

-800 10.1 1.5 5.8 5.5 4.7 3.9 3.8

-1000 11.1 1.2 5.2 5 4.4 3.8 3.7

-1200 12 0.9 4.9 4.6 4.1 3.7 3.6

Я, м 6

1

-1200

\1

- у/ """'

-.......Л..А..-.......^.......................

.........г......:........х.......:................:

-1000

-800 о „ -600

2 , М

Рис. 13. Максимальные заплески на берег К при накате волн различного происхождения: сплошные линии 1-5 — накат волн, возникших при совместном движении оползня с глубины г0 и одиночной волны с амплитудой ао = 1 ми начальным положением гребня в точке Хо = 20 км (1), 25 км (2), 30 км (3), 35 км (4), 40 км (5); штриховые прямые 1-5 — одиночной волны; сплошная линия 0 — волн, возникших при движении оползня

Некоторые результаты расчетов сведены в табл. 6. В третьей колонке приведены максимальные заплески оползневых волн, в четвертой строке справа — максимальные заплески одиночной волны. Приведены также максимальные заплески при одновременном движении оползневых волн и одиночной волны. Из таблицы видно, что значения максимальных заплесков при совместном движении больше, чем значения максимальных заплесков одиночной волны и волн, вызванных только подводным оползнем.

Заключение

В работе выполнено обобщение конечно-разностного алгоритма из [10] на случай движения поверхностных волн в ограниченном бассейне, у которого боковые стенки или дно могут перемещаться по заданным законам.

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

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

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

Разработанный алгоритм планируется использовать для моделирования процесса работы волновых энергетических станций нового типа [65], в которых энергия поверхностных волн передается погруженной в воду горизонтальной гибкой пластине.

Благодарности. Исследование выполнено при поддержке Программы Президиума РАН № 27 "Фундаментальные проблемы решения сложных практических задач на основе высокопроизводительных вычислений".

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

[1] Хакимзянов Г.С., ШЪкина Н.Ю. Оценка высот волн, вызванных подводным оползнем в ограниченном водоеме // Прикладная механика и техническая физика. 2012. Т. 53, № 5. С. 67-78.

Khakimzyanov, G.S., Shokina, N.Yu. Evaluation of the height of waves generated by an underwater landslide in a confined water reservoir // J. of Applied Mechanics and Technical Physics. 2012. Vol. 53, No. 5. P. 690-699.

[2] Гусев О.И., Шокина Н.Ю., Кутергин В.А., Хакимзянов Г.С. Моделирование поверхностных волн, генерируемых подводным оползнем в водохранилище // Вычисл. технологии. 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 // Comput. Technologies. 2013. Vol. 18, No. 5. P. 74-90. (In Russ.)

[3] Lynett, P.J., Liu, P.L.-F. A numerical study of the run-up generated by three-dimensional landslides // J. of Geophysical Research. 2005. Vol. 110. C03006.

[4] Коробкин А.А., Стуколов С.В., Стурова И.В. Движение вертикальной стенки, закрепленной на пружинах, под действием поверхностных волн // Прикладная механика и техническая физика. 2009. Т. 50, № 5. С. 132-142.

Korobkin, A.A., Stukolov, S.V., Sturova, I.V. Motion of a vertical wall fixed on springs under the action of surface waves // J. of Applied Mechanics and Technical Physics. 2009. Vol. 50, No. 5. P. 841-849.

[5] Khakimzyanov, G., Dutykh, D. Numerical modelling of surface water wave interaction with a moving wall // Communications in Computational Physics. 2018. Vol. 23, No. 5. P. 1289-1354.

[6] Kashiwagi, M. Non-linear simulations of wave-induced motions of a floating body by means of the mixed Eulerian-Lagrangian method // Proc. of the Institution of Mechanical Engineers. Pt C: J. of Mechanical Engineering Science. 2000. Vol. 214, No. 6. P. 841-855.

[7] Zhao, H., Freund, J.B., Moser, R.D. A fixed-mesh method for incompressible flow-structure systems with finite solid deformations // J. of Computational Physics. 2008. Vol. 227. P. 3114-3140.

[8] Нуднер И.С., Семенов К.К., Хакимзянов Г.С., Шокина Н.Ю. Исследование взаимодействия длинных морских волн с сооружениями, защищенными вертикальными экранами // Фундамент. и прикл. гидрофизика. 2017. Т. 10, № 4. С. 31-43.

Nudner, I.S., Semenov, K.K., Khakimzyanov, G.S., Shokina, N.Yu. Investigations of the long marine waves interaction with the structures protected by the vertical barriers // Fundamental and Applied Hydrophysics. 2017. Vol. 10, No. 4. P. 31-43. (In Russ.)

[9] Нуднер И.С., Семенов К.К., Лебедев В.В. и др. Численная модель гидроволновой лаборатории для исследования взаимодействия морских волн с гидротехническими сооружениями // Вычисл. технологии. 2019. Т. 24, № 1. С. 86-105.

Nudner, I.S., Semenov, K.K., Lebedev, V.V. et al. Numerical model of the hydrowave laboratory for studying the interaction of sea waves with hydrotechnical structures // Comput. Technologies. 2019. Vol. 24, No. 1. P. 86-105. (In Russ.)

[10] Численное моделирование течений жидкости с поверхностными волнами / Г.С. Хакимзянов, Ю.И. Шокин, В.Б. Барахнин, Н.Ю. Шокина. Новосибирск: Изд-во СО РАН, 2001. 394 с.

Numerical Simulation of Fluid Flows with Surface Waves / G.S. Khakimzyanov, Yu.I. Shokin, V.B. Barakhnin, N.Yu. Shokina. Novosibirsk: Izd-vo SO RAN, 2001. 394 p. (In Russ.)

[11] Хажоян М.Г. Численное моделирование поверхностных волн над подвижным дном // Вычисл. технологии. 2007. Т. 12, № 4. С. 96-105.

Khazhoyan, M.G. A numerical simulation of surface waves above moving bottom // Comput. Technologies. 2007. Vol. 12, No. 4. P. 96-105. (In Russ.)

[12] Хажоян М.Г., Хакимзянов Г.С. Численное моделирование взаимодействия поверхностных волн с подводными препятствиями // Вычисл. технологии. 2003. Т. 8, № 4. С.108-123.

Khazhoyan, M.G., Khakimzyanov, G.S. Numerical simulation of the surface waves and underwater obstacles interaction // Comput. Technologies. 2003. Vol. 8, No. 4. P. 108-123. (In Russ.)

[13] Волков К.Н. Разработка и реализация алгоритмов численного решения задач механики жидкости и газа // Вычисл. методы и программирование. 2007. Т. 8. С. 40-56. Volkov, K.N. Development and implementation of numerical algorithms for solving the problems of fluid dynamics // Numer. Methods and Programming. 2007. Vol. 8. P. 40-56. (In Russ.)

[14] Maxworthy, T. Experiments on collisions between solitary waves //J. Fluid Mech. 1976. Vol. 76, pt 1. P. 177-185.

[15] Загрядская Н.Н., Иванова С.В., Нуднер Л.С., Шошин А.И. Воздействие длинных волн на вертикальную преграду // Изв. ВНИИГ им. Б.Е. Веденеева. 1980. Т. 138. С. 94101.

Zagryadskaya, N.N., Ivanova, S.V., Nudner, L.S., Shoshin, A.I. Action of long waves on a vertical obstacle // Bulletin of VNIIG. 1980. Vol. 138. P. 94-101. (In Russ.)

[16] Давлетшин В.Х. Силовое воздействие одиночных волн на вертикальные сооружения // Совещание по цунами: Тез. докл. Горький: ИПФ АН СССР, 1984. С. 41-43. Davletshin, V.H. Force action of solitary waves on vertical structures // Soveshchanie po Tsunami: Tez. dokl. Gorky: IPF AN SSSR, 1984. P. 41-43. (In Russ.)

[17] Seabra-Santos, F.J., Temperville, A.M., Renouard, D.P. On the weak interaction of two solitary waves // Europ. J. Mechanics / B Fluids. 1989. Vol. 8, No. 2. P. 103-115.

[18] Chan, R.K.C., Street, R.L. A computer study of finite-amplitude water waves //J. Comput. Phys. 1970. Vol. 6, No. 1. P. 68-94.

[19] Fenton, J.D., Rienecker, M.M. A Fourier method for solving nonlinear water-wave problems: application to solitary-wave interactions //J. Fluid Mech. 1982. Vol. 118. P. 411-443.

[20] Железняк М.И. Воздействие длинных волн на сплошные вертикальные преграды // Накат цунами на берег: Сб. науч. тр. Горький: ИПФ АН СССР, 1985. С. 122-140. Zheleznyak, M.I. Influence of long waves on vertical obstacles // Soveshchanie po Tsunami: Tez. dokl. Gorky: IPF AN SSSR, 1985. P. 122-140. (In Russ.)

[21] Su, C.H. Mirie, R.M. On head-on collisions between two solitary waves //J. Fluid. Mech. 1980. Vol. 98, No. 3. P. 509-525.

[22] Byatt-Smith, J.G.B. The reflection of a solitary wave by a vertical wall //J. Fluid Mech. 1988. Vol. 197. P. 503-521.

[23] Chen, Y.Y., Kharif, C., Yang, J.H. et al. An experimental study of steep solitary wave reflection at a vertical wall // Europ. J. Mechanics / B Fluids. 2015. Vol. 49. P. 20-28.

[24] Umeyama, M. Experimental study of head-on and rear-end collisions of two unequal solitary waves // Ocean Engineering. 2017. Vol. 137. P. 174-192.

[25] Chen, Y.Y., Li, Y.-J., Hsu, H.C., Hwung, H.H. The pressure distribution beneath a solitary wave reflecting on a vertical wall // Europ. J. of Mechanics / B Fluids. 2019. Vol. 76. P. 66-72.

[26] Cooker, M.J., Weidman, P.D., Bale, D.S. Reflection of a high-amplitude solitary wave at a vertical wall // J. Fluid Mech. 1997. Vol. 342. P. 141-158.

[27] Madsen, P.A., Bingham, H.B., Liu, H. A new Boussinesq method for fully nonlinear waves from shallow to deep water // J. Fluid Mech. 2002. Vol. 462. P. 1-30.

[28] Craig, W., Guyenne, P., Hammack, J. et al. Solitary water wave interactions // Phys. Fluids. 2006. Vol. 18, No. 5. Paper 57106.

[29] Chambarel, J., Kharif, C., Touboul, J. Head-on collision of two solitary waves and residual falling jet formation // Nonlinear Proc. Geophys. 2009. Vol. 16. P. 111-122.

[30] Dutykh, D., Katsaounis, T., Mitsotakis, D. Finite volume schemes for dispersive wave propagation and runup //J. Comput. Phys. 2011. Vol. 230. P. 3035-3061.

[31] Touboul, J., Pelinovsky, E. Bottom pressure distribution under a solitonic wave reflecting on a vertical wall // Europ. J. Mechanics / B Fluids. 2014. Vol. 48. P. 13-18.

[32] Протопопов Б.Е. Численный анализ трансформации уединенной волны при отражении от вертикальной преграды // Изв. АН СССР. Механика жидкости и газа. 1990. № 5. С. 115-123.

Protopopov, B.E. Numerical analysis of the transformation of a solitary wave reflected by a vertical wall // Fluid Dynamics. 1990. Vol. 25, iss. 5. P. 746-753.

[33] Akrish, G., Rabinovitch, O., Agnon, Y. Extreme run-up events on a vertical wall due to nonlinear evolution of incident wave groups// J. Fluid Mech. 2016. Vol. 797. P. 644-664.

[34] Paprota, M., Staroszczyk, R., Sulisz, W. Eulerian and Lagrangian modelling of a solitary wave attack on a seawall // J. of Hydro-Environment Research. 2018. Vol. 19. P. 189-197.

[35] Riefolo, L., Contestabile, P., Vicinanza, D. Seiching induced by bichromatic and monochromatic wave conditions: experimental and numerical analysis in a large wave flume // J. of Marine Science and Engineering. 2018. Vol. 6, No. 2. Paper 68.

[36] Kihara, N., Niida, Y., Takabatake, D. et al. Large-scale experiments on tsunami-induced pressure on a vertical tide wall // Coastal Engineering. 2015. Vol. 99. P. 46-63.

[37] Hofland, B., Kaminski, M.L., Wolters, G. Large scale wave impacts on a vertical wall // Coastal Engineering Proc. 2010. No. 32.

[38] Frandsen, J.B., Berube, F. Large scale experimental wave impact on walls in the Quebec Coastal Physics Laboratory // Proc. of the ASME 2015 34th Intern. Conf. on Ocean, Offshore and Arctic Engineering. Newfoundland, Canada. 2015. Vol. 6. P. 1-10.

[39] Frandsen, J.B., Tremblay, O.G., Xharde, R. Preliminary investigations of wave impact on vertical walls with/without parapets and toe protection on deformable beach // Proc. of the 6th Intern. Conf. on the Application of Physical Modelling in Coastal and Port Engineering and Sci. (Coastlab16). Ottawa, Canada. 2016. P. 1-15.

[40] Laitone, E.V. The second approximation to cnoidal and solitary waves // J. Fluid Mech. 1960. Vol. 9, No. 3. P. 430-444.

[41] Tanaka, M. The stability of solitary waves // Phys. Fluids. 1986. Vol. 29, No. 3. P. 650-655.

[42] Clamond, D., Dutykh, D. Fast accurate computation of the fully nonlinear solitary surface gravity waves // Comput. & Fluids. 2013. Vol. 84. P. 35-38.

[43] Dutykh, D., Clamond, D. Efficient computation of steady solitary gravity waves // Wave Motion. 2014. Vol. 51, No. 1. P. 86-99.

[44] Железняк М.И., Пелиновский Е.Н. Физико-математические модели наката цунами на берег // Накат цунами на берег: Сб. науч. тр. Горький: ИПФ АН СССР, 1985. С. 8-33. Zheleznyak, M.I., Pelinovsky, E.N. Physico-mathematical models of the tsunami run upon a beach // Tsunami Climbing a Beach: Collection of Sci. Papers. Gorky: IPF AN SSSR, 1985. P. 8-33. (In Russ.)

[45] Федотова З.И., Хакимзянов Г.С., Гусев О.И. История развития и анализ численных методов решения нелинейно-дисперсионных уравнений гидродинамики. I. Одномерные модели // Вычисл. технологии. 2015. Т. 20, № 5. C. 120-156.

Fedotova, Z.I., Khakimzyanov, G.S., Gusev, O.I. History of the development and analysis of numerical methods for solving nonlinear dispersive equations of hydrodynamics. I. One-dimensional models problems // Comput. Technologies. 2015. Vol. 20, No. 5. P. 120-156. (In Russ.)

[46] Khakimzyanov, G., Dutykh, D., Fedotova, Z., Mitsotakis, D. Dispersive shallow water wave modelling. Pt I: Model derivation on a globally flat space // Communications in Comput. Physics. 2018. Vol. 23, No. 1. P. 1-29.

[47] Овсянников Л.В. Параметры кноидальных волн // Проблемы математики и механики: сб. науч. тр., посвящ. памяти акад. М.А. Лаврентьева. Новосибирск: Наука, 1983. С.150-166.

Ovsyannikov, L.V. Parameters of cnoidal waves // Problems of Mathematics and Mechanics. Novosibirsk: Nauka, 1983. P. 150-166. (In Russ.)

[48] Манойлин С.В. Некоторые экспериментально-теоретическое методы определения воздействия волн цунами на гидротехнические сооружения и акватории морских портов. Красноярск, 1989. 45 с. (Препр. / ВЦ CO АН СССР. № 5).

Manoylin, S.V. Some experimental and theoretical methods of estimation of tsunami wave action on hydro-technical structures and seaports. Krasnoyarsk, 1989. 45 p. (Prehrint / VZ SO AN SSSR. No. 5.) (In Russ.)

[49] Khakimzyanov, G., Dutykh, D., Gusev, O., Shokina, N. Dispersive shallow water wave modelling. Pt II: Numerical simulation on a globally flat space // Communications in Comput. Physics. 2018. Vol. 23, No. 1. P. 30-92.

[50] Стокер Дж.Дж. Волны на воде. Математическая теория и приложения. М.: Изд-во иностр. лит., 1959.

Stoker, J.J. Water waves: The mathematical theory with applications. New York: Interscience, 1957.

[51] Kanoglu, U., Synolakis, C.E. Long wave runup on piecewise linear topographies // J. of Fluid Mechanics. 1998. Vol. 374. P. 1-28.

[52] Synolakis, C.E., Bernard, E.N., Titov, V.V. et al. Validation and verification of tsunami numerical models // Pure and Applied Geophisics. 2008. Vol. 165. P. 2197-2228.

[53] Chubarov, L.B., Fedotova, Z.I., Shokin, Yu.I., Einarsson, B.G. Comparative analysis of nonlinear dispersive models of shallow water // Intern. J. of Comput. Fluid Dynamics. 2000. Vol. 14, No. 1. P. 55-73.

[54] Dao, M.H., Xu, H., Chan, E.S., Tkalich, P. Modelling of tsunami-like wave run-up, breaking and impact on a vertical wall by SPH method //J. Natural Hazards and Earth System Sciences. 2013. Vol. 13. P. 3457-3467.

[55] Tappin, D.R., Grilli, S.T., Harris, J.C. et al. Did a submarine landslide contribute to the 2011 Tohoku tsunami? // Marine Geology. 2014. Vol. 357. P. 344-361.

[56] Khakimzyanov, G.S., Gusev, O.I., Beisel, S.A. et al. Simulation of tsunami waves generated by submarine landslides in the Black Sea // Russ. J. of Numerical Analysis and Mathematical Modelling. 2015. Vol. 30, No. 4. P. 227-237.

[57] Beisel, S.A., Chubarov, L.B., Dutykh, D. et al. Simulation of surface waves generated by an underwater landslide in a bounded reservoir // Russ. J. of Numerical Analysis and Mathematical Modelling. 2012. Vol. 27, No. 6. P. 539-558.

[58] Елецкий С.В., Майоров Ю.Б., Максимов В.В. и др. Моделирование генерации поверхностных волн перемещением фрагмента дна по береговому склону // Вестн. КАЗНУ им. Аль-Фараби. Сер.: Математика, механика, информатика. 2004. Т. 9, спецвыпуск, ч. 2. С. 194-206.

Eletskij, S.V., Maiorov, Yu.B., Maximov, V.V. et al. Simulation of surface waves generation by a moving part of the bottom down the coastal slope // Vestn. KAZNU im. Al'-Farabi. Ser.: Matematika, Mekhanika, Inform. 2004. Vol. 9, special issue, pt 2. P. 194-206. (In Russ.)

[59] Shokin, Yu.I., Fedotova, Z.I., Khakimzyanov, G.S. et al. Modelling surfaces waves of generated by a moving landslide with allowance for vertical flow structure // Russ. J. of Numerical Analysis and Mathematical Modelling. 2007. Vol. 22, No. 1. P. 63-85.

[60] Whittaker, C.N., Nokes, R.I., Lo, H.-Y. et al. Physical and numerical modelling of tsunami generation by a moving obstacle at the bottom boundary // Environmental Fluid Mechanics. 2017. Vol. 17, No. 5. P. 929-958.

[61] Хажоян М.Г., Хакимзянов Г.С. Численное моделирование обтекания ступеньки потоком идеальной несжимаемой жидкости // Прикл. механика и техн. физика. 2006. Т. 47, № 6. С. 17-22.

Khazhoyan, M.G., Khakimzyanov, G.S. Numerical modeling of ideal incompressible fluid flow over a step //J. of Appl. Mechanics and Techn. Physics. 2006. Vol. 47, No. 6. P. 785-789.

[62] Букреев В.И. Ондулярный прыжок при обтекании открытым потоком порога в канале // Прикл. механика и техн. физика. 2001. Т. 42, № 4. С. 40-47.

Bukreev, V.I. Undular jump in open-channel flow over a sill // J. of Applied Mechanics and Technical Physics. 2001. Vol. 42, No. 4. P. 596-602.

[63] Букреев В.И., Гусев А.В. Волны за ступенькой в открытом канале // Прикл. механика и техн. физика. 2003. Т. 44, № 1. С. 62-70.

Bukreev, V.I., Gusev, A.V. Waves behind a step in an open channel // J. of Applied Mechanics and Technical Physics. 2003. Vol. 44, No. 1. P. 52-58.

[64] Shokin, Yu.I., Rychkov, A.D., Khakimzyanov, G.S., Chubarov, L.B. A combined computational algorithm for solving the problem of long surface waves runup on the shore // Russ. J. of Numerical Analysis and Mathematical Modelling. 2016. Vol. 31, No. 4. P. 217-227.

[65] Carmigniani, R.A., Benoit, M., Violeau, D., Gharib, M. Resonance wave pumping with surface waves // J. of Fluid Mechanics. 2017. Vol. 811. P. 1-36.

Поступила в 'редакцию 4 апреля 2019 г., с доработки — 17 мая 2019 г.

Numerical simulation of surface waves in a basin with moving impermeable boundaries

Palagina, Anna A., Khakimzyanov, Gayaz S.*

Institute of Computational Technologies SB RAS, Novosibirsk, 630090, Russia * Corresponding author: Khakimzyanov, Gayaz S. e-mail: khak@ict.nsc.ru

In numerical simulation of fluid flows with a free surface, the most difficult problems are those in which not only the free boundary is mobile, but also some other parts of the boundary of the region occupied by the liquid. For example, these are the problems of surface waves caused by underwater and coastal landslides, the generation of waves by wavemakers in laboratory flumes and tanks, the problem of waves run-up on the shore, the interaction of waves with moving wave protection walls and problems addressed moving semi-submerged or fully immersed objects.

The purpose of this paper is to analyze the properties and evaluate the capabilities of the computational algorithm based on the finite-difference approximation of the equations of potential fluid flows with a free boundary and designed to study surface waves in a confined basin, part of the impermeable boundary of which can be mobile. The algorithm relies on the use of curvilinear grids that adapt to all parts of the boundary, both moving and stationary.

New initial data are proposed for the problem of motion of a solitary or a single wave, consistent with the initial data for the shallow water equations of the second long-wave approximation. New non-reflecting conditions have been developed that allow waves to be emitted across the boundary of a flow region with very little reflection. A new initial approximation is proposed for the iterative process of calculating the potential of the velocity vector. By using this approach we significantly reduce the number of iterations at each time step. The original stability condition of the linearized difference scheme is derived. The reasons for the appearance of two peaks in the chronograms of pressure when the long waves of large amplitude roll onto a vertical wall are discussed.

The capabilities of the numerical algorithm are demonstrated on the problem of generating waves by a moving wall travelling in the initial part of the flume. The results

© ICT SB RAS, 2019

of the calculations well reproduce the experimental data, in particular, the decrease in the length and increase in the amplitude of the wave when it moves towards the shallow part of the flume, as well as the formation of a "dispersion tail" as the waves reverse motion after reflection from the vertical wall installed at the end of the flume.

The developed algorithm was used to study the process of generation for surface waves by an underwater landslide moving along an uneven bottom, and the interaction of these landslide waves with a single surface wave moving towards the shore. It is shown that the surface waves caused by an underwater landslide significantly affect the process of rolling of a single wave on the shore and it could increase its maximum splash.

Keywords: surface waves, movable boundary, underwater landslide, potential flows, finite-difference method, movable grid, stability, non-reflective condition, results of calculations.

Cite: Palagina, A.A., Khakimzyanov, G.S. Numerical simulation of surface waves in a basin with moving impermeable boundaries // Computational Technologies. 2019. Vol. 24, No. 4. P. 70-107. (In Russ.) DOI: 10.25743/ICT.2019.24.4.006.

Acknowledgements. The study was supported by the Program of Presidium of RAS No. 27 "Fundamental problems of solving sophisticated practical problems using supercomputers".

Received on April 4, 2019 Received in revised form on May 17, 2019