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

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

Том 10, № 2, 2005

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДВИЖЕНИЯ УЕДИНЕННОЙ ВОЛНЫ НАД ПОДВОДНЫМ

ПРЕПЯТСТВИЕМ

К. Е. Афанасьев, Е.Н. Березин Кемеровский государственный университет, Россия e-mail: [email protected], [email protected]

This paper is an extension of the authors' work on interaction of a solitary wave with a solid body partially submerged into fluid. This paper provides an analysis of the kinematical and dynamical characteristics of the wave during its passage over an underwater obstacle located on the horizontal bottom. Variation of the dynamic load during the interaction of surface waves with solid boundaries is considered. Fully nonlinear formulation of the problem is solved using the boundary element method.

В последнее время изучению взаимодействия уединенных волн с подводными препятствиями уделяется большое внимание. По этой тематике имеется много работ, в которых приводятся результаты как натурных экспериментов [1, 2], так и теоретических расчетов, выполненных на основе аналитических, полуаналитических и численных методов [3, 4, 5, 6, 7].

При исследовании взаимодействия поверхностных волн с подводными препятствиями необходимо учитывать множество характеристик, таких как стратификация жидкости, трение на границах, возможность отрыва слоев жидкости с поверхности тела, возможность обрушения поверхностных волн и т.д. Исследование указанных явлений требуют прямого численного моделирования ввиду существенной нелинейности процесса. Однако, если канал имеет большие размеры, а слой жидкости не является тонким, то можно пренебречь влиянием вязкости [8] и осуществлять моделирование на основе моделей идеальной жидкости. Так, например, в работах [9, 3] в рамках модели потенциальных течений достаточно детально изучены различные режимы обтекания препятствий типа полуцилиндра, одиночного выступа и ступеньки [10, 11, 7]. Для задачи о взаимодействии уединенной волны с частично и полностью погруженным в жидкость препятствием получено хорошее соответствие численных результатов с экспериментальными [5, 6].

Изучение взаимодействия жидкости с препятствием имеет важное значение при проектировании гидротехнических сооружений. По строительным нормам и правилам, действующим в Российской Федерации [12], распространяющимся на речные и морские гидротехнические сооружения, основным параметром при проектировании гидротехнического сооружения служит нагрузка, создаваемая волнением поверхности воды.

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

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

В настоящей работе, являющейся продолжением работы [17], представлены результаты расчетов волнового движения жидкости при движении уединенной волны над расположенным на дне препятствием прямоугольного сечения в виде одиночного выступа. Описывается алгоритм численного решения и исследуется влияние ширины и высоты препятствия, а также амплитуды волны на основные характеристики возникающего течения жидкости. Для численного исследования используется метод граничных элементов — один из эффективных методов решения задач гидродинамики идеальной жидкости [13, 18].

Методика решения нестационарных задач гидродинамики идеальной жидкости со свободными границами методом граничных элементов изложена в работе [13]. Для описания траекторий частиц (точек) свободной границы применяется метод Лагранжа. При этом интеграл Коши — Лагранжа и кинематическое условие на свободных границах жидкости записываются в виде обыкновенных дифференциальных уравнений первого порядка. Задача является нестационарной, и для ее решения применяется метод Эйлера с автоматическим выбором шага по времени. Начальная форма уединенной волны и распределение потенциала на ней получены из численного решения нелинейной стационарной задачи [19].

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

В расчетной области Д (рис. 1), ограниченной поверхностями Сх и Г!, Г2, решается уравнение Лапласа

Аф = 0, х(х,у) е Д. (1)

Граница Сх является свободной поверхностью жидкости, Г! — твердая граница бассейна, Г2 — твердая граница препятствия. На твердых границах выставляется условие непротекания

дП = 0, х(х,у) е Гь Г2. (2)

На свободной границе выполняются кинематическое и динамическое условия

пх

— = Уф, х(х,у) е Сх; (3)

Пф - 2 |УФ12 + у = 0, х(х,у) е Сх. (4)

В начальный момент времени £ = 0 граница Сх соответствует уединенной волне, на которой распределение потенциала задается из решения нелинейной стационарной задачи

Рис. 1. Схема расчетной области.

[19]. Краевая задача (1)-(4) для потенциала скоростей записана в безразмерном виде, где в качестве характерных величин выбираются ускорение свободного падения д и глубина бассейна Н. Все геометрические размеры и обозначения показаны на рис. 1.

2. Алгоритм движения по времени

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

Пусть в некоторый момент времени ¿к заданы положение свободной границы С и распределение потенциала фк на ней. Необходимо решить уравнение Лапласа (1) с условием фк на С и условием (2) на Г^ Г2. Новое положение свободной границы и распределение потенциала на ней для момента времени £к + т можно вычислить, используя условия (3) и (4), дискретный аналог которых расписывается по схеме Эйлера следующим образом:

хк+1 = хк + (уф)к т,

фк+1 = фк + (0, 5|(Уф)к|2 - ук)т, к = 0,1, 2..., (5)

где хк, фк — значения функций на к-м шаге по времени. В результате на каждом шаге по времени ¿к решается смешанная краевая задача для уравнения Лапласа (1) с граничными условиями (2) и (5). Описанный алгоритм представляет собой явную схему Эйлера первого порядка точности, которая, как известно из теории численных методов, является условно устойчивой и требует наличия некоторых дополнительных условий на выбор шага по времени. В данной работе используется хорошо зарекомендовавший себя алгоритм выбора шага по времени, описанный в работах [13, 17].

3. Численное моделирование

На каждом шаге по времени методом граничных элементов [18, 20] решается граничное интегральное уравнение (записанное по границе Г = С1 и Г1 и Г2), к которому сводится уравнение Лапласа внутри области Д. Применяя третью формулу Грина, в плоском случае можно записать

в (е) ф (£) + / ф(х) ^Г(х) = / дщхус (^ х) ^(х), (6)

г г

где Г(х) — граница области Д; ф(х) — гармоническая функция; п(х) — внешняя по отношению к области Д единичная нормаль поверхности Г(х). Параметр е(£) определяется следующим образом: е(£) = 1 для внутренней точки, е(£) = 0.5 для точки на гладкой границе, е(£) = а/2п для угловой точки границы (а — телесный угол, под которым видна поверхность из точки £); О (£,х) — фундаментальное решение уравнения Лапласа (1), и для плоских задач оно имеет вид

О (£,х) = -2П 1п(г) , (7)

где г — расстояние между точками £ и х.

Для численного решения интегрального уравнения (6) граница заданной области разбивается на ряд граничных элементов, в пределах которых функции ф и дф/дп изменяются линейно по длине элемента. В результате решение граничного интегрального уравнения сводится к системе линейных алгебраических уравнений (СЛАУ) для нахождения неизвестных значений потенциала и его нормальной производной в узлах элементов. Для решения СЛАУ из работы [21] взяты программы DECOMP и SOLVE, осуществляющие ее решение прямым методом Гаусса с выбором ведущего элемента. Технология решения интегрального уравнения (6), выбор шага по времени, дифференцирование функций, заданных на границе, и учет особенности при смене типа граничных условий изложены в работе [13].

4. Вычисление давления и полной энергии

Контроль за точностью численных расчетов осуществляется посредством закона сохранения полной энергии E = Ek + Ep и массы M. Значения кинетической Ek и потенциальной Ep энергий вычисляются по формулам

в

Ek = -1 / ф^з, Ep =\ f y2dx, M | (f)ydx - ф xdy\, (8)

2/ гдп' " 27 * ' 2

Г а \Г

где а и в — абсциссы точек пересечения границ Г и С\.

При исследовании взаимодействия солитонов с преградами актуальными являются задачи определения давления Р и динамических нагрузок Р3 на твердых границах. В соответствии с п. 1.13 из [12] нагрузки на вертикальную стену от воздействия прибойных волн допускается определять динамическими методами, учитывающими импульсы давления и инерционные силы. Схема расчета кинетической Ек, потенциальной Ер энергий и нагрузки Р3 приведена в [22, 17]. В настоящей работе для всех расчетов относительная ошибка полной энергии не превышает 1.5 % и массы — 0.06 %.

Для вычисления давления необходимо решить дополнительную краевую задачу:

Дфг = 0, х(х,у) € Б; (9)

1 2

Ф = -2 1УФ!2 - У, х(х,у) € Сь (10)

дф = 0, х(х,у) € Гь Г2. (11)

Решая систему линейных алгебраических уравнений для функции фt с граничными условиями (10), (11), находим ф4 на Г и Г2. После этого давление вычисляется по формуле

Р(х, у) = - (ф + 2 ^ф|2 + у) , х(х, у) € Гь Г2. (12)

5. Численные результаты

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

Таблица!

N/Ng е(фп) ) e(Vy) K (A)

100/43 6.83E-03 2.27E-03 8.55E-03 133

200/87 4.54E-03 1.11E-03 6.00E-03 273

400/175 2.51E-03 6.26E-04 3.73E-03 561

800/351 1.21E-03 6.90E-04 2.24E-03 1189

1500/659 5.90E-04 8.17E-04 1.44E-03 2395

Таблица2

N/Ng e(A), % e(M), % e(E), %

310/151 3.30 0.55 0.78

514/301 1.98 0.47 1.00

718/451 1.56 0.46 1.10

922/601 1.49 0.48 1.12

1126/751 1.35 0.50 1.17

1330/901 1.38 0.52 1.18

Для задачи о точности аппроксимации граничного интегрального уравнения (6) получены значения относительной погрешности. В табл. 1 приведены погрешности нормальной производной потенциала и компонент VX, Vy вектора скорости. В пятой колонке приводятся числа обусловленности K(A) матрицы A системы линейных алгебраических уравнений. Значения приведены в зависимости от дискретизации области (N — число узлов на всей

границе, Ng — число узлов на свободной границе области). Относительная погрешность max 1 /A - / 1

имеет вид е(/)

ции.

max | /A |

, где / — численное, /A — аналитическое значения функ-

Для задачи о распространении уединенной волны амплитуды А = 0.5 по бассейну постоянной глубины Н =1 получены погрешности основных характеристик волны. В табл. 2 приводится изменение амплитуды А, массы М и полной энергии Е на длину пробега волны в зависимости от количества точек разбиения области.

Алгоритм решения и численные результаты выполненных тестов подробно изложены в работах [22, 13, 17].

5.1. Взаимодействие уединенной волны е подводным препятствием, расположенным на дне

Ниже приведены результаты численного моделирования трансформации уединенной волны при ее движении в бассейне постоянной глубины Н = 1 над расположенным на дне препятствием в виде одиночного выступа (рис. 1). Варьируемыми параметрами задачи выбирались величины: А — амплитуда волны; к — высота выступа; а = хг — хг — ширина выступа, где хг и хг — абсциссы левой и правой вертикальных стенок выступа; Ь — расстояние между правой стенкой выступа при (х = хг) и правой границей бассейна.

Численные расчеты взаимодействия уединенной волны с подводным препятствием для диапазона варьируемых параметров А € [0.1 : 0.5], к € [0.1 : 0.9], а € [2 : 20] показали, что при движении волны над выступом она распадается на ряд отраженных и прошедших волн (рис. 2). Формы свободной поверхности при взаимодействии уединенной волны с подводным препятствием показаны для случаев:

— выступ расположен далеко от правой стенки бассейна (рис. 2, а, б) при 50 < Ь < 65;

— выступ расположен вблизи правой стенки бассейна (рис. 2, в, г), когда расстояние между правой стенкой выступа и правой границей бассейна было постоянным и равно Ь = 0.5.

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

а = 2, к = 0.2

¿=0.0

г= 53.954

X

г=о.о

109.567

-зо о 30 60

а = 15, к = 0.4 Рис. 2. Формы свободной поверхности: А = 0.3.

-60 -40 -20 0 20

X

г

а б

Рис. 3. Режим опрокидывания при А = 0.5 в случае, когда выступ расположен далеко от правой стенки бассейна: а — а = 2, Н = 0.8; б — а = 10, Н = 0.6.

?=[28.41;30.18]

г=[22.35;23.405]

Рис. 4. Режим опрокидывания для параметров А = 0.4, а = 2, Н = 0.7 (а) и А = 0.5, а = 2, Н = 0.9 (б) в случае, когда выступ расположен вблизи правой стенки бассейна.

а

волна распадается на большое количество волн достаточно большой амплитуды. Аналогичная картина течения наблюдается и для второго случая. Необходимо отметить, что на форму вторичных волн достаточно сильно влияет расположение выступа вблизи правой стенки бассейна. Из рис. 2, а, б видно, что за отраженной волной, которая распространяется в противоположном направлении к прошедшей, образуется вторая волна в виде впадины. Это так же отмечается и в работе [6].

При взаимодействии уединенной волны с подводным препятствием наблюдались режимы обрушения волн, причем увеличение ширины и высоты препятствия способствует этому. Эффект обрушения волны можно классифицировать как опрокидывание вперед (рис. 3, а, б, рис. 4, б) и опрокидывание волны назад во время отката от правой стенки бассейна, но в противоположном направлении относительно движения уединенной волны в начальный момент времени (рис. 4, а).

5.1.1. Выступ расположен далеко от правой стенки бассейна

Численные расчеты проводились для области Б £ [-40, 70]. Вершина волны для всех расчетов находилась при Ь = 0 в точке х = -15. Ширина выступа а изменялась в пределах 2 < а < 20. Левая граница выступа при (х = х\) для всех расчетов находилась в точке х = 0. Для различных значений амплитуды А распространяющейся волны изучалось влияние ширины выступа а и высоты к на величину максимального заплеска У^ах на правую стенку бассейна (рис. 5), амплитуды прошедшей волны ар (табл. 3), амплитуды отраженной волны а0 (табл. 4) и динамической нагрузки Р3 на правой стенке бассейна (рис. 6).

На рис. 5, а и 6, а линии, обозначенные цифрами 1-4, соответствуют ширине выступа а = 5, 10, 15, 20. Значение к было равно 0.4. На рис. 5, б и 6, б представлены зависимости для трех значений к при фиксированной ширине выступа а = 2. Цифры 1-3 соответствуют значениям высоты выступа к = 0.2, 0.4, 0.6. Численные расчеты показали, что при увеличении ширины выступа а и высоты к величины заплесков, нагрузка на правой стенке бассейна и амплитуда прошедшей волны убывают, а амплитуда отраженной волны возрастает.

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

ТаблицаЗ Влияние параметров А, Н и а на амплитуду прошедшей волны ар

Н А а = 5 а = 10 а = 15 а = 20

0.2 0.104 0.10925 0.11028 0.11148 0.11307

0.207 0.21884 0.21956 0.22154 0.22382

0.302 0.32210 0.32211 0.32275 0.32453

0.401 0.42888 0.42727 0.42577 0.42971

0.498 0.53611 0.53254 - -

0.4 0.104 0.10920 0.11015 0.11128 0.11276

0.207 0.21694 0.21687 0.21554 0.21487

0.302 0.31630 0.30957 0.30205 -

0.401 0.41750 - - -

0.498 - - - -

0.6 0.104 0.10772 0.11001 0.11116 0.11219

0.207 0.20967 0.20063 - -

0.302 - - - -

0.401 - - - -

0.498 - - - -

0.2 0.3 0.4 0.5

Рис. 5. Максимальный заплеск Утах на правую стенку бассейна.

0.1 0.2 0.3 0.4 0.5

а

а

Таблица4 Влияние параметров А, Н и а на амплитуду отраженной волны а0

Н А а = 5 а = 10 а = 15 а = 20

0.2 0.104 0.00532 0.00552 0.00569 0.00590

0.207 0.00920 0.00945 0.00978 0.01019

0.302 0.01191 0.01222 0.01265 0.01322

0.401 0.01402 0.01437 0.01488 0.01556

0.498 0.01552 0.01591 - -

0.4 0.104 0.01208 0.01224 0.01249 0.01276

0.207 0.02077 0.02110 0.02155 0.02215

0.302 0.02712 0.02747 0.02801 -

0.401 0.03208 - - -

0.498 - - - -

0.6 0.104 0.02110 0.02128 0.02157 0.02191

0.207 0.03648 0.03668 - -

0.302 - - - -

0.401 - - - -

0.498 - - - -

начальной амплитуды А. Однако значение нагрузки Р3 на правой стенке бассейна уменьшается (см. рис. 6), причем высота и ширина выступа способствуют этому. В табл. 4 приведены значения амплитуды отраженной волны а0. Видно, что с увеличением высоты выступа к и ширины выступа а значение амплитуды а0 увеличивается. В обеих таблицах ячейки с прочерком означают, что в момент движения волны над препятствием наблюдались режимы обрушения.

5.1.2. Выступ расположен вблизи правой вертикальной стенки бассейна

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

Численные расчеты проводились для области Б £ [-50, 20]. Расстояние Ь между правой стенкой выступа и правой стенкой бассейна было постоянным и равно Ь = 0.5. Вершина волны для всех расчетов находилась при Ь = 0 в точке х = -15. Ширина выступа а изменялась в пределах 2 < а < 20. Для различных значений амплитуды А распространяющейся волны изучалось влияние ширины а и высоты к выступа на величину максимального за-плеска У^ах на правую стенку бассейна (рис. 7), амплитуды отраженной волны а0 (табл. 5) и динамической нагрузки Р3 на правой стенке бассейна (рис. 8).

На рис. 7, а, 8, а линии, обозначенные цифрами 1-4, соответствуют ширине выступа а = 5, 10, 15, 20. Значение к было равно 0.4. На рисунках 7, б, 8, б представлены зависимости для трех значений к и фиксированной ширины выступа а = 2. Цифры 1-3 соответствуют значениям высоты выступа к = 0.2, 0.4, 0.6. Численные расчеты показали, что при увеличении ширины а и высоты к выступа величины заплесков на правой стенке бассейна, нагрузка на правой стенке бассейна и амплитуда отраженной волны возрастают.

Увеличение нагрузки на правой стенке бассейна Р3 можно объяснить тем, что в зазоре между правой стенкой выступа и правой стенкой бассейна возникает зона максимального гидростатического давления. Численное моделирование задачи о взаимодействии уединен-

Таблица5 Влияние параметров А, Н и а на амплитуду отраженной волны а0

Н А а = 5 а = 10 а = 15 а = 20

0.104 0.11028 0.11233 0.11406 0.11600

0.207 0.22000 0.22418 0.22478 0.22986

0.2 0.302 0.31995 0.32452 0.32906 0.33247

0.401 0.42400 0.42900 0.43468 0.44230

0.498 0.53058 - -

0.104 0.11015 0.11211 0.11350 0.11439

0.207 0.21254 0.21574 0.21594 0.21694

0.4 0.302 0.29295 0.29495 0.29669 0.29788

0.401 0.36454 - - -

0.498 - - -

0.104 0.11004 0.10203 0.09594 0.09180

0.207 0.20445 - - -

0.6 0.302 0.28873 - - -

0.401 - - - -

0.498 — — —

а у б

0.3 0.4 0.5

Рис. 7. Максимальный заплеск ^тах на правую стенку бассейна.

0.1 0.2 0.3 0.4

а

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

В табл. 5 показаны значения амплитуды отраженной волны ao. Видно, что при малой ширине и высоте выступа амплитуда отраженной волны увеличивается. Если ширина и высота выступа увеличиваются, то значение отраженной амплитуды ao уменьшается. Ячейки с прочерком означают, что в момент набегания волны на правую стенку бассейна наблюдались режимы обрушения.

Заключение

В настоящей работе рассмотрены результаты численного моделирования методом граничных элементов задачи о движении уединенной волны над расположенным на дне препятствием в виде одиночного выступа. Приведены результаты исследования влияния амплитуды волны A, высоты h и ширины a выступа на кинематические и динамические характеристики волны. Анализ значений максимального заплеска Ymax на правую стенку бассейна, амплитуды отраженной волны ao, амплитуды прошедшей волны ap и динамической нагрузки Ps на правой стенке бассейна показал, что выступ, расположенный далеко от правой стенки бассейна, является более эффективным, чем выступ расположенный вблизи правой стенки бассейна. В процессе взаимодействия уединенной волны с подводным препятствием наблюдались режимы обрушения волн, причем увеличение ширины и высоты выступа способствует этому.

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

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

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

[2] Synolakis C.E. The runup of solitary waves // J. Fluid Mech. 1987. Vol. 185. P. 523-545.

[3] Рузиев Р.А., Хлкимзянов Г.С. Численное исследование трансформации уединенной волны над подводным уступом // Вычисл. технологии: Сб. науч. тр. / РАН. Сиб. отд-ние. ИВТ. 1992. Т. 1, № 1. 1992. С. 5-22.

[4] Protopopov B.E. An efficient numerical method for calculation of strongly nonlinear water waves // Вычисл. технологии. 1998. Т. 3, № 3. С. 55-71.

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

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

[7] Шокин Ю.И., Рузиев Р.А., Хлкимзянов Г.С. Численное моделирование плоских потенциальных течений жидкости с поверхностными волнами. Красноярск, 1990 (Препр. ВЦ СО АН СССР; № 12).

[8] Esposito P.G., Graziani G., Orlandi P. Numerical solution of viscous flows about submerged and partly submerged bogies // 5th Intern. Conf. on Numer. Ship Hydrodinamics, Japan, 1989.

[9] Афанасьев К.Е., Стуколов С.В. Численное моделирование взаимодействия уединенной волны с препятствием // Вычисл. технологии. 1999. Т. 4, № 6. С. 3-16.

[10] Букреев В.И. Обтекание порога бурным потоком в открытом канале // ПМТФ. 2002. Т. 43, № 6. С. 54-61.

[11] Kawasaki K. Numerical simulation of breaking and post-breaking wave deformation process around a submerged breakwater // Coastal Eng. J. 1999. Vol. 41, N 3 & 4. P. 201-223.

[12] СниП 2.06.04-82. Нагрузки и воздействия на гидротехнические сооружения (волновые, ледовые и от судов). М., 1982.

[13] Афанасьев К.Е., Гудов А.М. Информационные технологии в численных расчетах: Учеб. пособие. Кемерово: КемГУ, 2001. 204 с.

[14] Clement A. The diffraction of a solitary wave by a free-surface piercing cylinder // 6th Intern. Workshop on Water Waves and Floating Bodies, Japan, 1991.

[15] Xueping G., Kunimitsu I. The characteristics of scouring and depositing in front of vertical breakwaters by broken clapotis // Coastal Eng. J. 1998. Vol. 40, N 1. P. 99-113.

[16] Kawahara M., Miwa T. Finite element analysis of wave motion // Intern. J. for Numerical Methods in Eng. 1984. Vol. 20. P. 1193-1210.

[17] Афанасьев К.Е., Березин Е.Н. Анализ динамических характеристик при взаимодействии уединенной волны с препятствием // Вычисл. технологии. 2004. Т. 9, № 3. С. 108-123.

[18] Бревия К., Теллес Ж., Вроувел Л. Методы граничных элементов. М.: Мир, 1987.

[19] Афанасьев К.Е., Стуколов С.В. О наличии трех решений при обтекании препятствий сверхкритическим установившимся потоком тяжелой жидкости // ПМТФ. 1999. Т. 35, № 1. С 27-35.

[20] Терентьев А.Г., Афанасьев К.Е. Численные методы в гидродинамике: Учеб. пособие. Чебоксары: ЧГУ, 1987. 94 c.

[21] Форсайт Дж., МАлькольм М., Моулер К. Машинные методы математических вычислений. М.: Мир, 1980.

[22] Афанасьев К.Е. Моделирование сильно нелинейных волновых течений // Вычисл. технологии. 1998. Т. 3, № 1. С. 3-12.

[23] Петров А.Г., Смолянин В.Г. Расчет нестационарных волн на поверхности тяжелой жидкости конечной глубины // ПММ. 1993. Т. 57, вып. 4. С. 137-143.

Поступила в редакцию 31 августа 2004 г., в переработанном виде —11 ноября 2004 г.