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

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

Том 8, № 4, 2003

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ВЗАИМОДЕЙСТВИЯ ПОВЕРХНОСТНЫХ ВОЛН С ПОДВОДНЫМИ ПРЕПЯТСТВИЯМИ*

М.Г. Хлжоян, Г. С. Хлкимзянов Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: pmargo@ngs.ru, khak@ict.nsc.ru

The finite-difference algorithm on movable grids is presented for calculation of the interaction between surface waves and submerged obstacles located at the bottom. The collision of a solitary wave, travelling over a horizontal bed, with a cylindrical bump of rectangular cross-section, placed parallel to the incident wave crest, is considered in the framework of the mathematical model of potential two-dimensional flows of an ideal incompressible fluid with a free boundary. The numerical results for the steady problem on fluid flow above a rectangular bump are compared with those from laboratory experiments.

Введение

Изучение воздействия поверхностных волн и течений на подводные препятствия играет важную роль в практических приложениях гидродинамики, в частности в проблеме защиты подводных сооружений. При численном моделировании процессов взаимодействия потока жидкости с подводными телами, расположенными на дне, необходимо учитывать множество факторов, например, таких, как стратификация жидкости, трение на границах, возможность отрыва слоев жидкости с поверхности тела и возникновение зон возвратного течения, возможность обрушения поверхностных волн, образование каверн за препятствиями и сложных вихревых структур. Все эти явления могут иметь место даже для препятствий простой формы, например, при обтекании порога прямоугольной формы, расположенного на горизонтальном дне [1-3]. Ясно, что для детального изучения процессов обтекания необходимо разрабатывать численные алгоритмы, основанные на сложных нелинейных математических моделях трехмерных турбулентных течений жидкости со свободной границей (см., например, [4]). Несомненно, что такие алгоритмы потребуют огромных затрат машинных ресурсов.

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант № 03-05-64108a, и Программы интеграционных фундаментальных исследований СО РАН, проекты №№ 3, 5.

© М.Г. Хажоян, Г. С. Хакимзянов, 2003.

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

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

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

1. Математическая постановка

Рассмотрим плоскопараллельное течение слоя жидкости с поверхностными волнами. Введем декартову систему координат Оху так, чтобы ось Ох лежала на невозмущенной свободной поверхности, а ось Оу была направлена вертикально вверх. Снизу слой жидкости ограничен неподвижным дном Г0, заданным функцией у = -к(х). На дне может располагаться препятствие с непроницаемыми стенками. Сверху жидкость ограничена свободной границей Г/ — подвижной кривой, описываемой в параметрическом виде формулами

х = С У = П(д>£)> (1)

где д — вещественный параметр, £ — время. В расчетах рассматривался слой жидкости, ограниченный слева и справа отрезками Г! и Г2 вертикальных прямых х = 0 и х = Ь соответственно. Область, занятую жидкостью, обозначим П(£), ее границу — Г, Г = Г0 и Г1 и Г2 и Г/ (рис. 1), ВД = ВД \ Г.

Математическая постановка задачи для потенциальных течений заключается в опре-

делении потенциала скорости, удовлетворяющего уравнению Лапласа

Д^ = 0, x 6 Q(t), (2)

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

+ = 0 д1 + = 0 xGг

dt dq dt , dt dq dt ' f

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

df |Vf|2 , ,

dt + - + 1 = const, x 6 Г/ (3)

условия. Здесь x = (x, y), u, v — декартовы компоненты вектора скорости u = Vf, dq/dt — скорость изменения параметра q для частиц жидкости, находящихся на свободной границе и движущихся вместе с ней. При t = 0 задаются начальные условия для потенциала и профиля свободной границы.

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

¥ =0, x 6 Го. (4)

дп

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

df

— = -u_TO, x 6 Г1. (5)

дп

На выходе Г2 ставится неотражающее краевое условие, которое обеспечивает свободный проход волн через эту границу. Описание используемых неотражающих условий будет дано ниже. Заметим, что в случае горизонтального дна задача будет иметь стационарное решение в виде равномерного потока u = (u_TO, 0) с невозмущенной свободной границей, если постоянную из динамического условия (3) положить равной u_TO/2.

Рассмотренная постановка задачи сформулирована в безразмерных величинах, при этом обезразмеривание проводилось по формулам

x* = -, i = th = 1 = -1, i= -i, ua = —, f = f , (6)

ho' л/h)' h0 h) h) a \fghh0 h^VghO'

где а =1, 2, х1 = х, х2 = у, u = u, u2 = v, g = const — ускорение свободного падения, h0 - некоторая характерная глубина, символом " обозначены безразмерные величины (в последующих формулах этот символ опущен).

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

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

t

x

x(q\qV), у = y(q\qV),

(7)

которое в каждый момент времени устанавливает взаимно-однозначное соответствие между "физической" областью П(£) и неподвижной вычислительной областью < простой формы в плоскости Од!д2. В задачах с горизонтальным дном в качестве < выбирался единичный квадрат. При наличии на дне препятствия, имеющего острые кромки, например препятствия в виде порога прямоугольной формы, в качестве вычислительной области использовался единичный квадрат с вырезанным снизу прямоугольником. Свободная граница Г/ являлась образом при отображении (7) верхней стороны 7/ области < и описывалась уравнениями

х = х(д!, 1,£), у = у(д!, 1,£), (8)

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

X = £(q1,t), У = n(q1,t), 0 < q1 < 1.

(9)

Уравнение Лапласа, кинематическое и динамическое условия примут в новой системе координат следующий вид:

где

д / дф дф\ д( дф дф \

dq^ dq1 dq2/ dq^ dq1 dq2/

—+= ° дП+= °

dt dq1 ' dt dq1

дф '"'2 u2

q = (q1, q2) g Q, (10)

q G Yf,

|u|2

At +1 - Vo + ^ + n

u

2

, g22 k11 = J'

7 7 g12

k12 = k21 = - J,

q G Yf,

, gn

k22 = J,

g«e (а,в =1, 2) — компоненты метрического тензора,

(11) (12)

g11

дх \

д^у +

дУ дq1

дх дх ду ду g12 g21 дq1 дq2 дq1 дq2, g22

J — якобиан преобразования (7),

J

дх ду дх ду

дq1 дq2 дq2 дq1

дх \ / ду дq2/ ^q2

2

v — контравариантные компоненты вектора скорости

dq1 dq1 dq1 dt dx dy

dq2 dq2 dq2 dt dx dy '

dq1 1 dy dx J dq2

dq2 dx

1 дУ

J dq1

dq1 1 i dy dx dx dy dt J \dt dq2 dt dq2

dq1 dy dq2

"dt

1 dx

J dq2

1 dx dy

jydtdq1 dt dq1

dq2 1 dx dy J dq1 dy dx

дх ду

Уо= 1 + ит+ьт>

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

dp dq1 dp dq2 u =--+--.

dq1 dx dq2 dx '

dp dq1 + dp dq2

dq1 dy dq2 dy

Запишем теперь в координатах qa краевые условия на участках границы Г0 и Г1? взяв в качестве примера область, изображенную на рис. 1. Мы предполагаем, что новая система координат согласована с границей области, т.е. каждый участок границы Г лежит на координатной линии либо первого семейства (q2 = const), либо второго (q1 = const). Так, для точек, принадлежащих дну, координата q2 = const = 0, а на верхней грани препятствия — q2 = const = Условие непротекания (4) на этих участках примет в координатах qa следующий вид:

к дР + к дР дс1 дс2

На лицевой грани препятствия координата ql -и условие (4) запишется здесь как

дР дР

к11^~Т + к12^~2

дс1 дс2

q2=const

const =

0.

qt

0.

на тыльной

const

q!=const

(13)

ql

(14)

Наконец, на входе Г1 координата с^ = 0 и условие втекания (5) принимает вид

dp dp dq1 dq2

V922U-C

(15)

q1 =0

В криволинейных координатах потенциал скорости определяется уже как решение более сложного, чем в декартовых координатах, уравнения (10) эллиптического типа с переменными коэффициентами и смешанными производными. Теперь не только условия (11), (12) на свободной границе являются нелинейными, но даже и условия на участках Го и Г1 становятся нелинейными, поскольку в коэффициенты ка@ выражений (13) - (15) неявно входят искомые функции £ и п, описывающие свободную границу. Однако переход к криволинейным координатам позволяет упростить расчет около подвижных границ, в частности около подвижной криволинейной свободной поверхности, ведь в новых координатах свободная поверхность является уже неподвижной и плоской. Также упрощается логика программ расчета, так как вычисления фактически ведутся на прямоугольной равномерной сетке, покрывающей вычислительную область (( простой формы с границами, параллельными координатным осям.

1

2

v

v

v

1

q

2. Алгоритм расчета

Вычислительная область <5 покрывается равномерной прямоугольной сеткой <5^ с шагом ка и количеством узлов Ыа в направлении оси а = 1, 2. Совокупность внутренних

узлов qj сетки — мультииндекс, ^ = (^1,^2)) будем обозначать граничных — 7^. Множество узлов, лежащих на 7/, обозначим через 7/,^.

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

Остановимся на некоторых вопросах, связанных с аппроксимацией уравнений и граничных условий. Предположим, что решение задачи на п-м слое по времени известно, т.е. при £ = £п известны значения потенциала положение свободной границы ■ пп, а также координаты узлов сетки х^-. Требуется определить значения этих величин при £ = £п+1 = £п + тп, тп — шаг по времени. Динамическое условие аппроксимировалось явной схемой

- |<|2 и2

Ц_11 II гП I ' 1 ' I пп = _

тп +1 - и°1" + 2 + П = 2

+ 1 - + ■ +1 = -—г, % е 7/л \ qwl ,N2 (16)

во всех узлах, кроме углового qN1,N2, в котором, как и на всей правой границе, значение ^N+N2 определялось из неотражающего условия. В формуле (16) полагалось

1 - <1 = ■■ - , (17)

при этом скорости движения узлов брались с предыдущего временноого слоя

хп _ хп—1 11п — 11П—1

■ - ■ ■ - ■

---^__/II ---^_

= _ , = _ •

'п—1 'П—1

Отметим также, что в силу вертикальности участка границы Г1 в левой угловой точке q1)N2 выполняется равенство х = 0, поэтому в этом узле 1 - г^ = -гпу41.

Разностные уравнения для потенциала в остальных узлах сетки получены с помощью интегроинтерполяционного метода. Для этого дифференциальное уравнение (10) записывалось в виде контурного интеграла, который аппроксимировался по формуле трапеций. Контур интегрирования зависел от расположения узла на сетке. На рис. 2, а показан контур интегрирования для внутренних узлов qj е

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

П = к ^ + к ^ V = к + к П = к11ттт + к12тг^2) У = к12ттт + к22тг^2 •

дд1 дд2 дд1 дд2

Тогда уравнение для потенциала можно переписать так:

дП дУ

^ + а? = ° qе < (18)

а краевые условия (13) - (15) примут такой вид

V = 0,

д2=сопв1

и

= 0, (20)

д1=сопв1

и = л/д^и-ж- (21)

д1=0

Для контура, изображенного на рис. 2, а, интегральное соотношение, соответствующее уравнению (18), имеет вид

J ивсс2 - J ивсс2 + J УвсС1 - J УвсС1 = 0. (22)

(ВС) (ЛБ) (БС) (ЛВ)

Будем предполагать, что сеточные функции и и У определены в центрах ячеек, т.е. в точках = ^+1/2^2+1/2. Для интегралов в (22) по сторонам прямоугольника ЛБСР

применяется формула трапеций. В результате вместо (22) получим разностное уравнение

и (Б) + и (С) _ и (Л) + и (Р) + У (Р) + У (С) _ У (Л) + У (Б)

2 2 2 2 + 2 1 2 1 '

или

и,1 + Уд2 = 0, % € (23)

Здесь использованы разностные производные в целых узлах %, которые для произвольной сеточной функции ф, определенной в центрах ячеек, задаются по формулам

ф(Б)+ ф(С) - ф(Л) - ф(Р) ф(Р)+ ф(С) - ф(Л) - ф(Б) ф^ =-Шу-' ф2'1 =-2Н2-•

Чтобы записать разностное уравнение (23) в виде уравнения относительно потенциала, необходимо указать способ вычисления величин и и У в центрах ячеек. Предполагается, что сеточные функции р и ха определены в целых узлах q3, а коэффициенты ка@ — в узлах ц_2+1/2 с полуцелыми индексами. Тогда, например, для точки С имеем

и (С) = (к11 Рд1 + к12Рд2 ) (С)' У (С) = (к12Рд1 + к22Р д2 ) (С)'

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

/ ф31 + 1,32 + ф31 + 1,32 + 1 - ф31,32 - ф31,32 + 1

ф1,з+1/2 =-^-'

фд2,3+1/2 2^2 •

Эти формулы используются и при вычислении коэффициентов кар. Например,

к11,3+1/2 = (3+1/2 ' д22,3+1/2 = (хд2 + У^) 3+1/2 ' 33+1/2 = (хд1 Уд2 - хд2Уд1)3+1/2 •

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

(к11Рд1 + к12Рд2)д1 + (к12Рд1 + к22Рд2)д2 = 0' % € (н- (24)

ф31,32+1 + ф31+1,32+1 - ф31,32 - ф31+1

32

32-

14

.71

|7

•- £>.- -1 -1 |- - ..С 1

1 1 Ю 1 |

1

1 1 "В

5 •— \2 -•- -

б

Рис. 2. Шаблон и контур интегрирования для внутреннего узла сетки (а), граничного узла на входе (б), граничного узла, совпадающего с точкой пересечения лицевой и верхней грани порога (в).

а

в

Шаблон этого уравнения, изображенный на рис. 2, а, является девятиточечным.

Интересно заметить, что для уравнения Лапласа (2), аппроксимируемого на квадратной сетке, т. е. при

кп = к22 = 1, к 12 = о, Ьл = Ь,2 = к, (25)

построенная схема (24) переходит в известную схему "косой крест" [14]:

<7 - 2< + <5 + < - 2<0 + < = о

Перейдем теперь к описанию способа аппроксимации краевых условий, рассмотрев в качестве примера условие на входе (15) (или (21)). В настоящей работе мы вместо аппроксимации этого условия аппроксимируем в граничных узлах уравнение для потенциала (18), предполагая, что оно выполняется вплоть до границы 71, являющейся прообразом входа Г1. При этом, так же как и во внутренних узлах, используется интегроинтерполяци-онный метод, но теперь с контуром интегрирования, зависящим от типа граничного узла. Для узлов, лежащих на 71 и не совпадающих с угловой точкой q11, контур интегрирования изображен на рис. 2, б. На его стороне АБ выполняется условие (21), поэтому разностное уравнение выглядит так:

ЧЮ+тк2 - /Ши-^2 + V(С)£ - V(В)| = о,

или

2_____ _ _ 2

— и (Е) + (Е) = — л/022 (0 )и-ж. (26)

к1 к1

Нетрудно проверить, что это шеститочечное разностное уравнение аппроксимирует краевое условие (15) на решениях уравнения (10) со вторым порядком по ка.

Интересно отметить, что при условиях (25) разностное уравнение (26) превращается в уравнение с односторонней разностью

1 ( <6 + <7 \

<о = и-ж,

Н\ 2

которое тем не менее аппроксимирует условие Неймана (5) дР/дх = и-ж при х = 0 со вторым порядком на решениях уравнения Лапласа (2).

Таким же образом получаются разностные аппроксимации уравнения (18) во всех других граничных узлах, за исключением узлов, лежащих на и на выходе из области,

п + 1

поскольку для этих узлов значения потенциала Р3 определяются из динамического и неотражающего условий соответственно. К примеру, для граничного узла, совпадающего с точкой пересечения лицевой и верхней граней порога, и контура интегрирования, изображенного на рис. 2, в, с учетом условий (19), (20) на отрезках ОЕ и БО получим восьмиточечное разностное уравнение

2 2 2 2

- зь и (Л) + 3 ид1 (К) + ^ У (С) + 3 Уд2 №) = 0.

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

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

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

Разработанный алгоритм вначале тестировался на задаче о взаимодействии уединенной волны с подводной ступенькой. В этой задаче определяющими параметрами являются высота ступеньки Ь и амплитуда а набегающей волны. Теоретические и экспериментальные исследования [9] показывают, что уединенная волна, проходя над ступенькой, распадается на отраженную волну и ряд прошедших волн, количество которых зависит от а и Ь. Численное моделирование выполнено в работах [9, 16-18]. В статье [9] использовалась неявная конечно-разностная схема для нелинейно-дисперсионной модели мелкой воды, в [16] — дискретная модель несжимаемой жидкости, в [17] — модель потенциальных течений, но с другой, чем в настоящей работе, аппроксимацией уравнения (10) и использованием более грубых сеток, в работе [18] проведен сравнительный анализ результатов решения рассматриваемой задачи в рамках различных моделей мелкой воды. Таким образом, эта задача достаточно хорошо изучена и, действительно, может рассматриваться как тестовая.

Начальные условия для уединенной волны, распространяющейся по каналу единичной глубины, задавались по формулам Л.В. Овсянникова [19] (см. также [20], формула (4.3.36)). Неотражающие условия на правой границе имели следующий вид [20]:

¿п I „- („, „,п , /й~„,п\ ип

Рп+1 = ( Рп + Тп ^,3УП -л/¥*и3) ип > 0' X € Г2' (27)

3 у Рп ип < 0' 3

где к3 = 1 - Ь — глубина над ступенькой. В качестве примера на рис. 3 показана свободная граница жидкости в различные моменты времени. Видно, что при а = 0.148 и Ь = 0.45

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

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

Далее представлены новые результаты численного моделирования трансформации уединенной волны при ее движении над подводным выступом — порогом прямоугольной формы высоты Ь и ширины I. В расчетах использовались сетки, адаптирующиеся лишь к границам области. Сетка строилась на основе простого "спрямляющего" преобразования координат [20]. При этом преобразовании первое семейство координатных линий состоит из кривых, расположенных между дном и свободной поверхностью, а второе — из вертикальных прямолинейных отрезков (рис. 4). На правой границе задавались условия (27) с К8 = 1. В приведенных ниже примерах расчетов мы ограничились рассмотрением необрушающихся волн. В таких случаях вместо (1) допустимо использование представле-

Рис. 3. Динамика свободной границы при накате уединенной волны амплитуды а = 0.148 на подводную ступеньку высотой Ь = 0.45.

Рис. 4. Сетка в задаче о набегании уединенной волны на порог.

ния свободной границы в виде однозначной функции

У = (28)

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

Рис. 5. Картина взаимодействия уединенной волны амплитуды а = 0.148 с порогом высоты Ь = 0.5 и ширины I = 2 а, I = 20 б.

Таблица 1

Влияние параметров а, Ь и I на амплитуду аг,1 первой прошедшей волны в задаче о взаимодействии уединенной волны с подводным порогом

Ь а аг, 1

1 = 2 1 = 7 1 = 15 1 = 20

0.334 0.199 0.195 0.194 0.193 0.192

0.334 0.287 0.277 0.272 0.258 0.25

0.334 0.41 0.385 0.361 0.303 0.277

0.334 0.487 0.444 0.391 0.296 0.266

0.5 0.21 0.204 0.198 0.177 0.163

0.5 0.293 0.279 0.257 0.196 0.169

0.5 0.4 0.368 0.298 0.2 0.175

0.5 0.495 0.432 0.31 0.211 0.186

0.667 0.197 0.186 0.167 0.117 0.1

0.667 0.29 0.266 0.201 0.134 0.115

0.667 0.387 — 0.255 0.148 0.123

0.667 0.437 — — 0.152 0.126

В табл. 1 приведены результаты расчетов амплитуды ац первой прошедшей волны для различных значений аг, Ь и I. Из таблицы видно, что при невысоком пороге и небольших значениях амплитуды набегающей волны величина ац мало отличается от аг, но при увеличении амплитуды набегающей волны разница увеличивается, причем и ширина порога способствует этому. В табл. 2 приведены значения амплитуды агд первой отраженной волны. Видно, что при I > 7 амплитуда отраженной волны при прочих равных условиях не зависит от ширины порога.

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

Отметим, что при набегании уединенной волны на порог квадратного поперечного сечения (Ь = I) волновые картины взаимодействия весьма схожи с теми, которые возникают, когда подводным препятствием служит полукруговой цилиндр [12, 13]. Так, например,

Таблица 2

Влияние параметров аг, Ь и I на амплитуду агд первой отраженной волны в задаче о взаимодействии уединенной волны с подводным порогом

Ь аг аг,1

1 = 2 1 = 7 1 = 15 1 = 20

0.5 0.21 0.021 0.026 0.026 0.026

0.5 0.293 0.027 0.031 0.031 0.031

0.5 0.4 0.032 0.037 0.037 0.037

0.5 0.495 0.034 0.04 0.04 0.04

0.667 0.197 0.036 0.038 0.038 0.038

0.667 0.29 0.047 0.05 0.05 0.05

0.667 0.387 — 0.058 0.058 0.058

0.667 0.437 — — 0.062 0.062

Рис. 6. Взаимодействие уединенной волны с порогом квадратного сечения: а — аг = 0.7, Ь = I = 0.4; б - аг = 0.462, Ь = I = 0.7.

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

4. Взаимодействие потока жидкости с подводным препятствием

Задача о взаимодействии равномерного потока жидкости с подводным препятствием решалась методом установления. Свободная граница при Ь = 0 задавалась невозмущенной. Начальные значения компонент скорости определялись по формулам

дф

Щ (Х,У) = ду> ^0(х,У) = - дх>

при этом функция тока ф являлась решением задачи о потенциальном течении несжимаемой жидкости в области, показанной на рис. 1, с непроницаемой верхней стенкой у = 0:

Дф = 0, х € П(0),

I I Щ—

ф |х=0= Щ_те(у + 1), -1 < у < 0, ф |л=ь= (у + Н3), -Н3 < у < 0, ф |Го =0 ф |у=0 = 0 < х <

Таким образом, при Ь = 0 задавалось стационарное течение в канале с твердой верхней крышкой у = 0. При Ь > 0 крышка убиралась, возникало существенно нестационарное течение со свободной границей, которое впоследствии устанавливалось.

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

= 2^1-1,М2 - > ^2 = 0, ¿2 <^2. (29)

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

В табл. 3 для задачи о взаимодействии набегающего потока с подводной ступенькой приведены значения разностей уровней установившейся поверхности до ступеньки (Л._) и над ней (Л-+). Видно, что увеличение скорости набегающего потока приводит к росту разности уровней = Н_ — Л.+ , причем при достижении определенной скорости щ_ установившегося режима уже не существует.

Рис. 7. Выход на стационарный режим процесса взаимодействия равномерного потока со ступенькой (а), с широким порогом (б): п-1Х) = 0.3, Ь = 0.334, I = 15.

Таблица 3

Влияние параметров п-1Х) и Ь на разность уровней жидкости до ступеньки и над ней

Ь 5Н Ь 5Н Ь 5Н

0.334 0.1 0.009 0.5 0.1 0.021 0.667 0.1 0.04

0.334 0.22 0.044 0.5 0.22 0.08 0.667 0.22 0.2

0.334 0.3 0.085 0.5 0.3 0.189 0.667 0.28 0.265

0.334 0.34 0.107 0.5 0.34 0.286 0.667 0.3 —

0.334 0.38 0.14 0.5 0.38 0.331

0.334 0.42 — 0.5 0.42 —

Результаты решения задачи о взаимодействии равномерного потока жидкости с порогом сведены в табл. 4. В настоящей работе рассчитывался только подтопленный режим сопряжения бьефов, т.е. когда степень подтопления 5 = (Н+ — Ъ)/Н* превышала критическое значение 5* (в справочнике [5] рекомендуется использовать значение 5* = 1.25, в экспериментах [1] получено уточненное значение 5* = 1.33). Здесь Н* — критическая глубина, которая в безразмерных переменных (6) выражается через расход по формуле Н* = (п-жН-)2/3. Напомним, что в наших численных экспериментах задавалась скорость потока п-Х), а расход определялся в результате решения задачи. Отметим также, что здесь рассмотрены только докритические течения над порогом, поскольку во всех расчетах НтП/Н* > 1. Из табл. 4 видно, что для рассмотренных здесь режимов течений соответствие результатов расчетов и экспериментальных данных [1] вполне удовлетворительное. Это дает основание сделать вывод о возможности использования разработанного алгоритма для решения прикладных задач.

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

Таблица 4

Влияние параметров Ь, I и п-(Х на уровень жидкости перед порогом и над ним

b 1 Расчет Эксперимент [1]

h- h • hmin hmin/h* h+ S h- h hmin

0.334 2.08 0.1 0.99 0.651 3.04 0.995 3.11 0.99 0.62

0.334 2.08 0.18 1.012 0.635 1.98 0.99 2.04 1.01 0.65

0.334 2.08 0.22 1.029 0.626 1.69 0.991 1.77 1.08 0.60

0.334 2.08 0.25 1.038 0.598 1.47 0.995 1.63 1.05 0.60

0.334 2.08 0.3 1.071 0.577 1.23 1.012 1.44 1.09 0.56

0.334 2.08 0.33 1.104 0.598 1.17 1.028 1.36 1.13 0.55

0.334 2.08 0.35 —

0.5 3.12 0.1 0.988 0.461 2.15 0.989 2.29 0.99 0.48

0.5 3.12 0.15 1 0.415 1.47 0.979 1.70 1.01 0.44

0.5 3.12 0.18 1.002 0.354 1.109 0.981 1.51 1.03 0.40

0.5 3.12 0.2 1.04 0.367 1.046 1.005 1.44 1.07 0.41

0.5 3.12 0.21 —

0.667 4.16 0.05 0.989 0.298 2.21 0.985 2.36 0.99 0.31

0.667 4.16 0.1 0.99 0.23 1.079 0.988 1.50 1.02 0.27

0.667 4.16 0.13 —

го. В дальнейшем предполагается, оставаясь в рамках конечно-разностного алгоритма на адаптивных сетках, применить для решения рассмотренной стационарной задачи итерационный метод (см., например, [21]).

В заключение авторы выражают благодарность Г. А. Тарнавскому за полезные обсуждения и постоянное внимание к работе.

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

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

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

[3] Букреев В.И., Гусев А.В. Каверны за водосливом с широким порогом // ПМТФ. 2002. Т. 43, № 2. С. 129-135.

[4] HüDGes B.R., Street R.L. On simulation of turbulent nonlinear free-surface flows // J. Comput. Phys. 1999. Vol. 151, No. 2. P. 425-457.

[5] Справочник по гидравлическим расчетам / Под ред. П.Г. Киселева. М.: Энергия, 1972.

[6] Букреев В.И., Гусев А.В., Ляпидевский В.Ю. Транскритическое течение над порогом в открытом канале // Изв. РАН. Механика жидкости и газа. 2002. № 6. С. 5562.

[7] Ляпидевский В.Ю., Тешуков В.М. Математические модели распространения длинных волн в неоднородной жидкости. Новосибирск: Изд-во СО РАН, 2000.

[8] Железняк М.И. Воздействие длинных волн на сплошные вертикальные преграды // Накат цунами на берег. Горький: ИПФ АН СССР, 1985. С. 122-140.

[9] SEABRA-SANTos F.J., Renouard D.P., Temperville A.M. Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle // J. Fluid Mech. 1987. Vol. 176. P. 117-134.

[10] Железняк М.И., ПЕлиновский Е.Н. Физико-математические модели наката цунами на берег // Накат цунами на берег. Горький: ИПФ АН СССР, 1985. С. 8-33.

[11] Марчук А.Г., Чубаров Л.Б., Шокин Ю.И. Численное моделирование волн цунами. Новосибирск: Наука, Сиб. отд-ние, 1983.

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

[13] Cooker M.J., PEREGRINE D.H., Vidal C., Dold J.W. The interaction between a solitary wave and a submerged semicircular cylinder //J. Fluid Mech. 1990. Vol. 215. P. 1-22.

[14] Самарский А.А., Андреев В.Б. Разностные методы для эллиптических уравнений. М.: Наука, 1976.

[15] ХАжоян М.Г. Численное моделирование взаимодействия волн с подводными препятствиями: Дипл. работа. НГУ, Новосибирск, 2002.

[16] Паутов В.Н., Франк А.М. Дискретная модель несжимаемой жидкости с частицами переменной массы // Тр. Всесоюз. совещ. по численным методам в задачах волновой гидродинамики. ВЦ СО АН СССР, Красноярск, 1991. С. 80-86.

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

[18] ChubARov L.B., FedotovA Z.I., Shokin Yu.I., ElNARssoN B.G. Comparative analysis of nonliner dispersive shallow water models // Intern. J. of Comp. Fluid Dynamics. 2000. Vol. 14. P. 55-73.

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

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

[21] Brummelen E.H., Raven H.C., Koren B. Efficient numerical solution of steady free-surface Navier-Stokes flow // J. Comput. Phys. 2001. Vol. 174, No. 1. P. 120-137.

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