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

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

Том 12, № 1, 2007

РАСЧЕТ ВОЛНОВЫХ ДВИЖЕНИЙ ЖИДКОСТИ НА ОСНОВЕ УРАВНЕНИЙ ЭЙЛЕРА*

Б. Е. Протопопов Институт гидродинамики им. М.А. Лаврентьева, Новосибирск, Россия

e-mail: boris@hydro.nsc.ru

A numerical method for computation of strongly non-linear water waves is constructed on the basis of the Euler's equations. In order to test the method and to demonstrate its strength, the problem of stretching of a gravity-free liquid ellipse and the problem on overturning of a smooth bore have been solved numerically. A comparison of numerical results with an exact solution (in the case of the former problem), as well as with the similar results obtained in the frames of potential flow model, makes it possible to conclude that the developed method is highly accurate and capable.

Введение

Данная работа посвящена построению численного метода расчета волн на воде на базе уравнений Эйлера. Идеология такого построения в основном аналогична использованной в численном методе, разработанном ранее [1] в рамках модели потенциального движения идеальной несжимаемой жидкости.

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

Среди других методов, напоминающих данный, следует отметить так называемые "бессеточные" методы, в частности метод MPS (moving particle semi-implicit method) [2-4]. Предлагаемый в настоящей работе метод с таким же успехом можно отнести к бессеточным, поскольку в качестве лагранжевых координат жидких частиц здесь просто используются их номера (при этом номера могут быть и дробными), а эйлеровы координаты

*Работа выполнена при финансовой поддержке Интеграционного проекта СО РАН (№ 2.12) и частичном финансировании Совета поддержки ведущих научных школ (грант № НШ-5873.2006.1).

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

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

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

Уравнения Эйлера имеют вид

Ux + wz = 0;

+ i *Рx I I IP z

uux + wuz = —, wt + uwx + wwz =--g,

p p

(1) (2)

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

Из системы (1), (2) вытекает уравнение Пуассона для давления:

Рхх + 'Рхх = —2р (их'Шх — ихШх) ,

которое и будет использовано ниже вместо условия (1).

Переходя от эйлеровых координат (£,х,г) к лагранжевым (£, ), получаем

xt = u, zt = w;

xt a + xz c zt a + zz c ut =-f-, wt =----g;

5

5

at + cz = — 2e, 5pt = p(aa — fic), 5pz = p(yc — в a);

(3)

(4)

(5)

а = + в = — (х^х^ + ), 7 = х2 + ,

(6)

8 = — х^ , £ = и — — и£ —.

Здесь уравнения (4) получены из (2), к ним добавлены уравнения траекторий (3). Уравнение Пуассона для давления принимает общий вид эллиптического уравнения. Последнее удобнее записать в виде системы (5) из трех уравнений первого порядка, в которой искомыми функциями, наряду с давлением, являются а и с — проекции ускорения жидких частиц на криволинейные координатные линии. Коэффициенты, используемые в (4), (5), определены формулами (6).

Эволюционную часть задачи — уравнения (3), (4) — следует дополнить начальными данными, а эллиптическую (уравнения (5)) — граничными условиями.

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

х = хо(£,С), 2 = 2о(£,С ), Ь = 0; (7)

и = ио(£,С), - = -о (£,С), Ь = 0. (8)

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

лагранжевых координатах.

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

С(Ь,х,г ) = 0 (9)

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

с(ь,х(и,с ),г(а,с)) = о,

где х(Ь, £), £) — координаты жидких частиц. Дифференцирование этого равенства

два раза по Ь дает условие на ускорения частиц:

Сли4 + С*-г = -Сй, (10)

Сй = Си + 2С4жи + 2С4* - + Схх и2 + 2СЛ* и- + С**(11) Подстановка в (10) выражений (4) для компонент ускорений приводит к выражению

С^ а + С^ с = Сгг,

где использованы обозначения:

С? = Сжх? + С*, Сс = Сжхс + С*2С, С« = -8 (Сй - С*(12)

Пусть твердая граница состоит из трех участков — левой и правой стенок и дна. Уравнения этих участков в эйлеровых координатах имеют вид (9), но с функциями Ь(£,х,г), Д(£, х, г), В(£, х, г) вместо х, г), а в лагранжевых координатах — £ = £ = £д, £ = Св соответственно. Тогда, очевидно, Ь = Д^ = В« = 0 и граничные условия принимают следующий окончательный вид:

Ь*

Ь д**

Д«

Ви Вс

-Г-, £ =

а = ъ", £ = £я;

С = Св.

(13)

(14)

(15)

Здесь производные в правых частях вычисляются по формулам (11), (12), в которых функцию О следует заменить на Ь, Я или В.

Граничное условие на свободной поверхности — нулевое давление:

р = 0, С = (з.

(16)

а

с

2. Метод расчета

Метод расчета задачи (3)-(8), (13)-(16) построен по аналогии с [1], поэтому излагается здесь достаточно кратко.

Эволюционные уравнения (3), (4) рассчитываются по неявной двуслойной схеме, например:

Д£

хп+1,т+1 = хп + (ига+1'т + ип) . (17)

Здесь Д£ — шаг по времени (постоянный), первый верхний индекс означает номер временного слоя, второй — номер итерации. Аналогично рассчитываются остальные эволюционные уравнения. Итерационный процесс необходим в силу нелинейной и неявной зависимости правых частей этих уравнений от искомых функций.

Решение эллиптической задачи (5), (13)—(16) находится установлением (на каждом шаге по времени) решения соответствующей параболической задачи. Последняя рассчитывается методом дробных шагов, т. е. решение двумерной задачи находится последовательным решением двух одномерных задач:

гт дт

рк+1/2 — гак+1/2 = рк + т (ск + 2ет) , ак+1/2 — — рк+1/2 = ^ск; (18)

« ^ V с у рат « ат

т т Я т

ак+1/2 = Ьт, £ = £ь, ак+1/2 = Д^, £ = £д (19)

Ьт' ' ' ' Дт

(горизонтальная прогонка);

/ \ Лт Дт

0к+1 _ Тск+1 = рк + т (ак+1/2 + 2ст\ ск+1 _ _рк+1 = в ак+1/2.

— тсс = / + т (аГ'' + 2^) , ск- — —'Г = ^^(20)

Вт т

ск+1 = Вт, с = Св, рк+1 = 0, с = С* (21)

Вс

(вертикальная прогонка). Здесь т есть итерационный параметр (фиктивный шаг по времени); верхние индексы означают номера итераций. Итерационный процесс, необходимый для установления решения параболической задачи (помеченный индексом к), может быть организован независимо от итерационного процесса, используемого в эволюционных уравнениях (индекс т), а может с ним совпадать (к = т).

Численное дифференцирование, необходимое как для расчета коэффициентов (6), так и для решения написанных выше одномерных краевых задач, осуществляется с применением сплайнов (квадратичных или кубических). При этом использование квадратичных сплайнов дает значения производных гладких функций с погрешностью О (Ж-2), а решения краевых задач вида (18), (19) или (20), (21) — с погрешностью О (Ж-3), N — число ячеек сетки. В случае применения кубических сплайнов порядки соответствующих погрешностей получаются О (Ж-3) и О (Ж-4). Это установлено экспериментальным путем с помощью несложных тестов. Продвижение по времени согласно схеме (17) сопровождается погрешностью О (АЬ2).

3. Результаты расчета

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

3.1. Вытягивание жидкого невесомого эллипса

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

^ +(г, =1

2 / 2 N 2

X) + Ь

с горизонтальной X(Ь) и вертикальной Z(Ь) полуосями, связанными зависимостью Х^ = 1. Из-за отсутствия массовых и поверхностных сил здесь волн как таковых не образуется (поверхность жидкости остается монотонной). Однако достоинством этой задачи является наличие точного решения (в нелинейной постановке) [5]. Более строго, полная задача сведена к задаче Коши для одного обыкновенного дифференциального уравнения, описывающего эволюцию полуоси эллипса (вертикальной, например):

* = Т^Ь • * <°) = *. (22)

Здесь точка сверху означает производную по времени; о — произвольный параметр, с помощью которого можно регулировать скорость вытягивания эллипса вдоль своей вертикальной (о > 0) или горизонтальной (о < 0) оси.

Расчеты выполнены не во всем эллипсе, а лишь в его четырехугольном фрагменте, а именно задаются следующие "твердые" границы:

Ь(Ь,х,2) = —х, Я(Ь,х,г) = х — ^Х (Ь), В (Ь,х,2) = —2.

Таким образом, левая стенка "установлена" на вертикальной оси эллипса, дно — на горизонтальной, правая стенка также вертикальна, но подвижна, и закон ее движения согласован с точным решением задачи (параметр ^ следует брать в пределах от 0 до 1).

В начальный момент времени должны быть заданы форма свободной поверхности и значение потенциала скорости ф(г, х,г) на ней [5, 1]:

1 — (I)'

*=0

х

1 — 1 Х 1 = 25 (х)

Ф (22 — х2

4=0,^=^5 (ж)

2До

(1 — Д^х2) , Хо = X(0), До = л/Хо2 + ^о2

0

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

Р

1,

9

0,

а

А

1, ^ 2'

1 (^ Хо = 1)

на сетках 10 х 10 и 20 х 20 с шагом по времени Дг = 0.1Дх в случае квадратичных сплайнов и Дг = (Дх)2 в случае кубических (Дх = ^Хо/М, М — количество ячеек сетки по оси £).

Результаты расчетов представлены в табл. 1. Здесь значения X(г) вычислены как Z-1(г), а значения Z(г) получены численным решением задачи (17) по схеме, описанной выше. ¿х и ¿2 есть относительные погрешности координат узлов сетки (жидких частиц) на свободной поверхности:

¿х =

х — х Из

X ;

¿2 =

2 — г ||5

Z '

Под нормой функции понимается ее максимальное абсолютное значение, индекс Б означает, что максимум ищется по узлам сетки не во всей расчетной области, а лишь на поверхности (£ = Сз), "тильдой" помечено точное решение, которое можно посмотреть в [1].

2

2

Таблица 1. Результаты расчетов задачи

Квадратичные сплайны Кубические сплайны

г X 2 ¿х 103 ■ 103 6х 104 ■ 104

10 х 10 20 х 20 10 х 10 20 х 20 10 х 10 20 х 20 10 х 10 20 х 20

0 1.000 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

1 0.439 2.28 0.11 0.04 0.21 0.05 0.04 0.00 0.03 0.00

2 0.272 3.68 0.60 0.24 0.37 0.14 0.14 0.01 0.11 0.01

3 0.196 5.09 1.43 0.58 1.05 0.38 0.32 0.02 0.30 0.02

4 0.154 6.50 2.54 1.03 1.83 0.67 0.60 0.04 0.66 0.04

5 0.126 7.92 3.92 1.55 2.64 0.97 1.03 0.06 1.26 0.07

6 0.107 9.33 5.58 2.14 3.41 1.27 1.66 0.09 2.15 0.13

7 0.093 10.75 7.54 2.79 4.14 1.58 2.56 0.13 3.42 0.20

8 0.082 12.16 9.81 3.49 4.79 1.88 3.80 0.17 5.11 0.31

9 0.074 13.58 12.43 4.23 5.89 2.18 5.49 0.23 7.31 0.44

10 0.067 14.99 15.43 5.01 7.51 2.48 7.71 0.32 10.06 0.61

Таблица 2. Результаты расчетов задачи

Квадратичные сплайны Кубические сплайны

г X 2 ¿х 104 ■ 104 6х 105 ■ 105

10 х 10 20 х 20 10 х 10 20 х 20 10 х 10 20 х 20 10 х 10 20 х 20

0 1.000 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

1 0.439 2.28 0.30 0.08 3.18 0.80 0.33 0.03 0.19 0.02

2 0.272 3.68 0.68 0.21 5.06 1.27 0.72 0.08 0.34 0.03

3 0.196 5.09 1.04 0.34 6.34 1.59 1.10 0.12 0.48 0.06

4 0.154 6.50 1.39 0.45 7.31 1.84 1.46 0.16 0.59 0.08

5 0.126 7.92 1.74 0.57 8.10 2.04 1.80 0.20 0.70 0.11

6 0.107 9.33 2.10 0.71 8.76 2.21 2.13 0.27 0.80 0.13

7 0.093 10.75 2.48 0.88 9.34 2.36 2.46 0.35 0.90 0.15

8 0.082 12.16 2.89 1.08 9.85 2.49 2.79 0.48 0.99 0.17

9 0.074 13.58 3.34 1.34 10.31 2.61 3.13 0.70 1.08 0.18

10 0.067 14.99 3.85 1.65 10.73 2.72 3.50 1.02 1.16 0.19

К моменту времени Ь =10 жидкий эллипс вытягивается по вертикали практически в 15 раз и во столько же раз уменьшается его горизонтальный размер, так что отношение полуосей эллипса изменяется примерно в 225 раз. Столь существенные деформации жидкой области приводят к росту погрешности, однако даже в худшем случае (квадратичные сплайны, сетка 10 х 10) она не превышает 1.55%. Применение кубических сплайнов значительно (на порядок) уменьшает вычислительную погрешность в сравнении с квадратичными сплайнами. Двукратное уменьшение шагов сетки влечет существенное уменьшение относительных погрешностей — примерно трехкратное для квадратичных сплайнов и 15-24-кратное для кубических, что в принципе не противоречит порядкам аппроксимации, указанным в конце предыдущего параграфа, если принять во внимание конечные размеры использованных сеток и нелинейность полной задачи.

Для сравнения выполнены расчеты этой же задачи, но по модели потенциального движения идеальной жидкости. Результаты расчетов представлены в табл. 2. Можно видеть, что использование потенциальной модели приводит к существенно меньшим (на порядок) погрешностям. Однако с измельчением расчетной сетки различие двух моделей ослабевает.

3.2. Опрокидывание плавного бора

В этой задаче твердые границы задаются с помощью функций [1]:

£(Ь,х,г) = —х, Я(Ь,х,г) = х — (/ + !Л), В(Ь,х,2) = —(2 + й). Начальная форма свободной поверхности выглядит следующим образом [6]:

2 = — (1 + 1^^ Ь = 0,

^(х) = 1 и (1+1апЬ х — х*

2 V Л

Необходимо также задать начальное значение потенциала скорости на свободной поверхности [6]:

1 -г т I ч т х х*

ф = ^(х)йх = 2и ( х + Л 1п еовЬ —) , Ь = 0, < = Сз.

Расчеты выполнены при значениях параметров задачи р = 9 = Н = 1, I = 36, и = —1.7, х* = 20, Л = 2 на сетках 180 х 10 и 360 х 20 с шагом по времени Дг = 0.05Дх в случае применения квадратичных сплайнов и Дг = 0.25(Дх)2 в случае кубических сплайнов (Дх = 1/М).

Волна движется справа налево, ее профили в моменты времени г = 0.00, 0.50,1.00,..., 5.00, 5.50, 6.00 показаны на рис. 1, а, в моменты г = 0.00, 0.50,1.00,... ,6.00, 6.50, 6.915 — на рис. 1, б, 360 х 20); в моменты г = 0.00, 0.50,1.00,... ,6.00, 6.50, 6.87 — на рис. 1, в, в моменты г = 0.00, 0.50,1.00,... ,6.00, 6.50, 7.00 — на рис. 1, г. Точки соответствуют узлам сетки (жидким частицам). Для большей наглядности соседние точки соединены отрезками прямых. Показана не вся жидкая область, а лишь левая ее часть, где и происходит опрокидывание волны.

Рис. 1. Опрокидывание плавного бора: квадратичные сплайны, сетка 180 х 10 (а) и 360 х 20 (б); кубические сплайны, сетка 180 х 10 (в) и 360 х 20 (г).

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

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

По результатам расчета в рамках потенциальной модели построен рис. 2, который аналогичен рис. 1. Волновые профили соответствуют моментам времени £ = 0.00, 0.50, 1.00, ..., 5.00, 5.50, 5.93 (а); £ = 0.00, 0.50, 1.00, ..., 5.00, 5.50, 5.93 (б); £ = 0.00, 0.50, 1.00,..., 5.00, 5.50, 5.74 (в); £ = 0.00, 0.50, 1.00, ..., 5.00, 5.50, 5.67 (г).

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

Рис. 2. Опрокидывание плавного бора (расчет по потенциальной модели): квадратичные сплайны, сетка 180 х 10 (а) и 360 х 20 (б); кубические сплайны, сетка 180 х 10 (в) и 360 х 20 (г).

Рис. 3. Положение жидких частиц в последний момент времени t = 8.145 (уравнения Эйлера, кубические сплайны, сетка 360 х 20).

360 х 20), оказался уникальным по степени деформации жидкой области — вытягивании волны в тонкую струю. Вид волны (струи) в последний момент времени t = 8.145, до которого удается досчитать (на следующем временном шаге итерационный процесс расходится), показан на рис. 3. В отличие от предыдущих рисунков здесь изображены не только поверхностные, но и внутренние жидкие частицы. Для большей наглядности соседние (по направлению £) жидкие частицы соединены отрезками прямых. Получающиеся кривые (ломаные) представляют собой координатные линии £(t, x, z) = const. Видно, что некоторые внутренние частицы оказались за пределами волны. Другими словами, произошел перехлест соседних координатных линий. Такой перехлест приводит к тому, что в каких-то точках якобиан S становится отрицательным, что и делает невозможным дальнейший счет. Заметим, что значения якобиана есть не что иное, как площади ячеек сетки в физических координатах. В силу несжимаемости жидкости эти площади должны сохраняться, однако при сильных деформациях жидкой области погрешность сохранения локальной массы (якобиана) может накапливаться вплоть до появления ячеек с отрицательными значениями площади. При этом интегральная масса жидкости сохраняется хорошо. В частности, для момента времени, указанного на рис. 3, отклонение полной массы жидкости от своего начального значения не превысило 0.1 %.

Заключение

На основе уравнений Эйлера построен метод расчета сильно нелинейных волн на воде. Выполнены расчеты двух модельных задач: о вытягивании жидкого невесомого эллипса и об опрокидывании плавного бора. Результаты расчетов демонстрируют работоспособность построенного численного метода при весьма значительных деформациях физической (жидкой) области.

Как показывают расчеты первой задачи, предлагаемый метод может заметно уступать

по точности аналогичному методу, базирующемуся на потенциальной модели. Однако с измельчением расчетной сетки разница в точности двух методов ослабевает. В некоторых же случаях (как показывают расчеты второй задачи) построенный метод может давать более "качественное" (т. е. более "долгоживущее") численное решение в сравнении с той же потенциальной моделью. К числу достоинств предлагаемого метода, несомненно, относятся его универсальность (он может быть применен к расчету не только потенциальных, но и вихревых движений) и возможность дальнейшего обобщения на случай стратифицированной и/или вязкой жидкости.

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

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

[2] Koshizuka S., Tamako H., Oka Y. A partical method for incompressible viscous flow with fluid fragmentation // Comput. Fluid Dyn. J. 1995. Vol. 4. P. 29-46.

[3] Koshizuka S., Nobe A., Oka Y. Numerical analysis of breaking waves using the moving particle semi-implicit method // Intern. J. Numer. Meth. Fluids. 1998. Vol. 26. P. 751-769.

[4] Ataie-Ashtiani B., Farhadi L. A stable moving-particle semi-implicit method for free surface flows // Fluid Dyn. Res. 2006. Vol. 38, N 4. P. 241-256.

[5] Овсянников Л.В. Задача о нестационарном движении жидкости со свободной границей. Новосибирск: Наука, 1967. 108 с.

[6] Cooker M.J. The Interaction of Steep Water Waves and Coastal Structures: Ph. D. Thesis. University of Bristol, 1990.

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