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

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

Том 1, № 3, 1996

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ УСТАНОВИВШИХСЯ ТЕЧЕНИЙ ЖИДКОСТИ В РАМКАХ МОДЕЛИ МЕЛКОЙ ВОДЫ*"»"

Г. С. Хлкимзянов, Н. Ю. ШокинА Институт вычислительных технологий СО РАН Новосибирск, Россия

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

Введение

Модели установившихся течений жидкости с поверхностными волнами имеют широкий круг приложений в практике и уже давно являются предметом интенсивных исследований с помощью методов механики сплошной среды. Для этих моделей иногда удается получить точные решения и выявить основные качественные закономерности течений с помощью аналитических методов [1, 3, 5]. Вместе с тем для большинства практически важных задач получить детальную картину установившихся волн можно лишь с помощью численных методов (см. напр., [13]).

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

*© Г. С. Хакимзянов, Н.Ю.Шокина, 1996.

^ Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, проект №97-01-00819.

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

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

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

1. Математическая модель

Пусть П — односвязная область в плоскости декартовых координат Математиче-

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

V- Ни = 0, х =(ж1,ж2) € П, (1.1)

Н 2 д

(Ни ■ V)u + gV — = дНVII - кН[ез х и] - иУ (1.2)

2 С

и краевым условиям на границе Г = дП

H u

хеГх

ai(x), u ■ n =0, Hu ■ n = v2(x), (1.3)

хеГо

хеГ2

где u = (u1,u2), ua (a = 1, 2) — компоненты вектора скорости в направлении осей Oxa, V = (d/dx1, d/dx2), n — внешняя нормаль к Г, Г = Г1 U Го U Г2, Г1 — вход в П(а1 ■ n < 0), Г0 — непроницаемая часть границы, Г2 — выход из П (v2 > 0), Г1 П Г2 = 0, H — полная глубина, H = n + h, z = — h(x) — функция, описывающая дно, g — ускорение свободного падения, e3 = (0,0,1) — базисный вектор, направленный вдоль вертикальной оси Oz, перпендикулярной плоскости x1Ox2, V = |u|, C — коэффициент Шези, k — параметр Кориолиса, плотность жидкости принята равной единице.

В безразмерных переменных система уравнений (1.1), (1.2) принимает вид

d F1 + dF = G (1 4)

dx1 dx2 ' '

где

/ Ищ \ / Ип2

F1 = Ищ? + И2/2 , F2 = Ищщ

\ Ищщ ) \Ии? + И2/2)

( 0 ^

G = ИНХ1 + /сог^И - uiV/fCh У ИНх2 - fcorUlИ - u?V/f?hУ Обезразмеривание проводилось по формулам

xa = xah0, п = П^о, ИИ = И^0, h = hh0, Ua = ua\fgh0, (1.5)

где h0 — характерная глубина, символом "~"обозначены размерные величины. Через fcor и fch обозначены безразмерные значения параметра Кориолиса и коэффициента Шези:

fCor = Who/g, fch = С1/т/д.

Краевые условия для безразмерных переменных имеют вид (1.3), при этом a1 и v2 равны отношению соответствующих величин из (1.3) к h0\/gh0.

В дополнение к условиям (1.3) будем считать, что в некоторой точке M0 £ ll известна полная глубина:

И (M0) = И0. (1.6)

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

U дф U дф 4. - du1 I дщ2 (л

Ищ = —, Ищ = -ттг, ш = rotu = + ттт • (L7)

dx2 dx1 dx2 dx1

Тогда для ф и ш получим уравнения

А д /1 дф \ _

¿-dx^ яд^ -ш, (1.8)

i=1 N 7

V2 1

V ■ И иш = rot И V(h - —) - -г rot(uV)

V 2 / fC h

^ х 1 2 м ; 2

при этом граничные значения для ф однозначно определяются из условий (1.3).

(1.9)

2. Переход к криволинейным координатам

Пусть

= ж"^1,^) (2.1)

есть взаимно однозначное невырожденное отображение единичного квадрата Q, лежащего в плоскости координат д1Од2, на область П. Ради простоты будем далее предполагать, что участки входа и выхода являются связными кусками границы, причем Г1 является образом при отображении (2.1) левой стороны 71 квадрата Q, Г2 — правой стороны 72, а непроницаемая часть Г0 состоит из двух участков Г0 и Г, являющихся образами, соответственно, нижней 70 и верхней 70 сторон Q.

Уравнения (1.8) и (1.9) в новых независимых переменных да запишутся следующим образом:

д /. дф дф\ д /. дф дф\

ттт кпттт + k12T¡~2 + т;-2 к21ттт + k22 ТТ = — дд1 \ дд1 дд2/ дд2 \ дд1 дд2/

+ A J) = - §¿2 + /,

(2.2) (2.3)

где

kn =

g22

ki2 = k

HJ "21 HJ

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

1 2 2 2 1 1 2 2 1 2 2 g11 = (xq1) + (xq1) , g12 = xq2 + xq2 , g22 = (xq2) + (xq2 )

g12 , g11 k22

HJ

</ — якобиан преобразования (2.1), </ = ж^х жо2 —Хд2 , Vя контравариантные компоненты скорости, связанные с ф соотношениями

1 дф я/дд2,

v = —

1 дф HJ дд1

(2.4)

В уравнении для вихря

/ = я

дд0

h-

V2

va V

¿Ch

а

1, 2,

(2.5)

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

VI = #11^ + = #12^ + #22 V2.

При вычислении полной глубины будут использоваться уравнения движения в форме Громеки — Лэмба, записанные в криволинейных координатах:

— Я^2ш + = /1 + /осГЯ^2, Я^ш + д|2 = /2 — /соГЯМ (2.6)

где p = H2/2.

1

v

2

3. Разностные уравнения

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

/ (ки + к12 1)2) — (к21 дфт + к22 ^ =—/ / ^я1^ (з.1)

С Дс

где Дс — прямоугольник с контуром С (рис. 1, штриховая линия).

Аппроксимация соотношения (3.1) описана в работе [8]. Разностное уравнение для ф получается девятиточечным, функции ф, Я, </, предполагаются определенными в целочисленных узлах сетки, а сеточная функция ш — в центрах ячеек прямоугольной сетки с шагами Л,2, покрывающей вычислительную область

Рис. 1.

Рис. 2.

Уравнение для вихря (2.3) будем аппроксимировать схемой с центральными разностями

Лшо = (/)о, (3.2)

Л = Л1 +Л2,

где

Л^о

1 'фо — фв

2Ы1 1. Ы

1 " фо — У

2Ы,

шэ

^ - фА

-^1

) ув - у а -щ>4 +--;-щ>2

_ (/1 )м — (/1)5 (/2)д - С/2У

(/)о = ы + ы ■

Шаблон показан на рис. 2.

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

(зг + т Л)

+ Лшп = /,

(3.3)

где т — итерационный параметр, п — номер итерации, г — тождественный оператор. После приближенной факторизации [2] оператора 3Е + тЛ

зг + тЛ « ^Узг + -т=л^Узг + -т= Л,^)

уравнение (3.3) примет вид

(зг + 4зЛ-)(зг + 4зЛ,) ^ = "Л-п + /■

Если обозначить

(3.4)

(3.5)

ш

1+1/2 = (з + 4=л,)

V з

т . , ш^1 - шп

т

и ввести новую неизвестную ^ = (шп+1 — шп)/т, то разностные уравнения для функции вихря будут иметь следующий вид:

т

(л/зг + 4=Л1)шп+1/2 = —Лшп + /,

уЗ

(3.6)

(J + -J Л2)£

ÜJ

,n+l/2

(3.7)

^п+1 = шп + т£. (3.8)

Разностные уравнения систем (3.6) и (3.7) являются трехточечными, поэтому для их решения можно применять метод прогонки. Применяя его к уравнениям (3.6), находим значения ^п+1/2. Подставляя их в (3.7), находим значения £ и потом по формуле (3.8) вычисляем ^п+1. В данном случае прогонка будет устойчивой при некотором ограничении на итерационный параметр т .

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

= const.

Го

границы Г выполняются условия непротекания и ф

Функция р определяется на основе следующего выражения, вытекающего из уравнений (2.6):

p(qj) = Ро + J (/i + HJv2(w + /cor)) dqi + (/2 - HJvV + /cor)) dq2

7(90 ,qj)

(3.9)

где p0 = H2/2, q0 — прообраз точки M0 из (1.6). В качестве кривой y, соединяющей точку qj, в которой требуется определить p, с точкой q0, в которой задано p0, будем брать ломаную, звенья которой параллельны осям Oqa и проходят по сторонам ячеек.

Покажем независимость вычисления p(qj) от пути интегрирования. Рассмотрим произвольную ячейку вычислительной области Q (см. рис. 2). Покажем, что при указанной ниже аппроксимации интегралов из (3.9) вычисление p в точке C по разным контурам ABC и ADC даст одно и то же значение. Имеем

pabc(c) = pa + hi(/i)s - (фб - Фа)/сог - (фб - Фа)+ ш +

2

+М/2)е - (Фс - Фб)/cor - (Фс - Фб)

Ш0 + Ш3 2 ;

padc(C) = pa + M/2V - (фя - Фа)/сог - (фя - Фа)Ш + Ш +

2

+hi(/i)w - (Фс - Фя)/cor - (Фс - Фя)

Ш +

2 '

(3.10)

(3.11)

Если из (3.10) вычесть (3.11), то получим pabc (C) - padc(C) = -hih2

1 f Фс - фя Фб - Фа

2h

hi

-ш4 -

1 / фс - Фб _ фя - Фа 2hi V h2 h2

(/i)N - (/i)s (/2)e - (/2)

шН +

W

hi J h2

= -hih2 [Л1Ш0 + Л2Ш0 - (/)0] .

hi

(3.12)

Предполагая, что итерационный процесс (3.6) - (3.8) сошелся и, следовательно, выполняется стационарное разностное уравнение (3.2) для вихря, получим, что правая часть

равенства (3.12) равна нулю, то есть значение р действительно не зависит от контура интегрирования. Таким образом, для вычисления функции р используется метод согласованной аппроксимации. Для случая прямоугольных сеток он рассматривался в работе [9]. Для течений вязкой несжимаемой жидкости он известен как маршевый метод расчета давления [12].

После вычисления р полная глубина определяется по формуле

Н = фРр. (3.13)

4. Итерационный процесс

Задача об установившемся течении жидкости с поверхностными волнами решается в три этапа.

1. Вначале решается задача о потенциальном течении "под крышкой". Для этого полагается п(х) = 0, Н(ж) = Л,(ж), й(ж) = 0 и методом последовательной верхней релаксации решается система разностных уравнений для —.

2. Решается задача о вихревом течении "под крышкой". Для этого в качестве начального приближения берется решение задачи о потенциальном течении "под крышкой", полагается Н(ж) = Л,(ж) и организуется итерационный процесс, на п + 1-м шаге которого производятся следующие вычисления.

1. Решается система разностных уравнений для —™+1.

2. Вычисляются граничные значения вихря на входе Г1. Для этого используется выражение для вихря в криволинейной системе координат

1 ( дг1

й = 7 I-д? + а?) . (41)

Ковариантные компоненты скорости вычисляются по заданному на входе расходу (1.3). При вычислении касательной производной дг^/д^2 используются компоненты вектора а1 и значения полной глубины Н на Г1. При вычислении производной дг2/д^1 используется г2, определенное в центре приграничной ячейки по значениям — П+1. Разностный аналог формулы (4.1) применяется для определения предварительного значения й, а окончательные значения вихря на входе получаются релаксацией й и граничных значений вихря с предыдущей итерации:

' 5Шй + (1 - 5Ш)йп , (4.2)

Гх

йп+1

Гх

где положительный параметр релаксации [6] 8Ш = 0(^1).

3. Методом стабилизирующей поправки решается система уравнений (3.2) для

п + 1

определения вихря й^.+1/2 во внутренних узлах.

Итерации — — й проводятся до тех пор, пока йп+1 и йп различаются на величину, большую некоторого заданного числа £ш > 0.

3. Решается задача о вихревом течении жидкости со свободной границей. В качестве начального приближения берется решение задачи о ви'хревом течении "под крышкой". В итерационном процессе кроме определения —п+1, йп+1 и йп+11/2 производится вычисление полной глубины НП+1. При этом на основе формулы (3.13) находятся предварительные значения Н, а в качестве окончательных значений полной глубины берутся интерполированные по Н и Нп:

НП+1 = 6нН + (1 - 6и)НП, (4.3)

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

Итерации в этой задаче продолжаются до тех пор, пока Нп+1 и Нп различаются на величину, большую некоторого заданного > 0.

5. Описание геометрии водоема

В рассматриваемом ниже примере моделирования течения жидкости с поверхностными волнами на заданном участке речного русла предполагается, что боковые границы являются вертикальными стенками. Геометрия границ и поверхности дна задается набором опорных сечений и опорных точек. Пусть П является проекцией на плоскость г = 0. Тогда опорное сечение — это отрезок А, В, прямой, пересекающей ось водоема и две противоположные части Г0 и Г0 границы Г области П. Здесь ] — номер сечения, ] = 1, ..., N, N — количество сечений. В качестве оси водоема можно взять произвольную прямую или кривую, располагающуюся в направлении вытянутости водоема и проходящую примерно посередине области П.

Будем предполагать, что П представляет собой плоский криволинейный четырехугольник со сторонами Г1, Г2, Г0, Г0, причем Г1, Г2 являются прямолинейными отрезками, совпадающими с первым А1В1 и последним А^В^ сечениями соответственно. Геометрия двух других сторон Г0, Г0 определяется местоположением точек А, и В, соответственно.

С каждым сечением А, В, связаны следующие вводимые параметры:

ж/, у0 — координаты точки г0 пересечения отрезка А, В, с осью водоема;

а, — угол наклона вектора А, В, к оси Ож;

, — знакоопределенное расстояние к-й опорной точки ]-го сечения от точки г/, при этом для опорных точек, располагающихся на отрезке А, В, между г/ и В,, полагаем > 0, а между А, и г/ полагаем , < 0, к = 1, ..., NAjв,, NAjв, — количество опорных точек в сечении А, В;

— вертикальная отметка уровня дна в к-й опорной точке сечения А, В,, отсчитываемая от уровня моря.

Вкратце остановимся на ограничениях для указанных выше вводимых данных. Глубины в опорных точках = — должны отличаться от нуля, > 0. Здесь — вертикальная отметка уровня невозмущенной поверхности жидкости (плоскости г = 0). Угол наклона а, отсчитывается от положительного направления оси Ож, при этом —180° < а, < 180°. Точки А, и В, обязательно входят в число опорных точек, они имеют номера к =1 и к = NAj в,. Другие опорные точки пронумерованы в порядке их расположения от А, к В,. Если в некотором сечении А,В, глубина постоянна, то при вводе данных для этого сечения достаточно взять Nа^ в, = 2 и задать расстояния /,;1 < 0 и /,,2 > 0 точек А, и В, от г/.

Для поверхности используется кусочная параметризация, причем для куска, проектирующегося на А,В,В,+1А,+1, используется своя, независимая от других четырехугольников, параметризация параметрами р, д, 0 < р, д < 1. Для координат ж1 = ж, ж2 = у произвольной точки из А, В, В,+1А,+1 будем использовать билинейную интерполяцию

(р,д) = (1 — д — р + + д(1 — р)жВ, + р(1 — + Р?хВ,.+1, (5.1)

при этом координаты х^. , х'В. точек А^- и Б^- вычисляются с помощью формул для координат опорных точек:

Xj,k = lj,k cos aj + x0, yj¡k = j sin aj + yj0, k = 1, ..., N^b Для третьей координаты используется следующая интерполяция:

(5.2)

z (p,q) = (1 - p)zj (q) + pzj+i(q),

(5.3)

где Zj (q) — кусочно-линейная функция, задаваемая на интервалах qj,k < q < ?j,k+i фор-

j

мулой

zj (q) = —j +

q- qj,k qj,k+1 — qj,k

( —hj,k+1 + hj,k)

k = 1,

n4..r.. - 1,

а qj,k — значения параметра q, соответствующие опорным точкам:

lj,k lj,1

qj,k = l l j.N* o ~ l

j,1

(5.4)

(5.5)

Таким образом, участок поверхности который проектируется на четырехугольник А^- Б^- Б^+1А^+1 плоскости г = 0, взаимно однозначно отображается на квадрат (( = [0; 1] х [0; 1], и для произвольных параметров (р, д) £ (( координаты соответствующей точки г(р, д) £ могут быть вычислены по формулам (5.1), (5.3).

На рис. 3 показана геометрия участка реки и используемые опорные сечения и опорные точки (последние помечены маркерами). Пример взят из работы [10]. При обезразмерива-нии характерная глубина Л,0 равнялась трем метрам.

Рис. 3. Геометрия водоема, опорные сечения и опорные точки.

6. Построение криволинейной сетки

Вначале с помощью одномерного метода эквираспределения [7] расставляются по N1 узлов на Г0 и Г0 и по N2 узлов на Г1 и Г2. Далее генерируется сетка внутри области О. Координаты ее узлов рассчитывались на основе двумерного аналога метода эквираспределения

[7, 11]:

д ( дха \ д ( дха \

ддг ^22+ ^ =°, а = 1,2, (6-1)

где т — управляющая функция. При т = 1 получается криволинейная "квазиравномер-ная"сетка. При переменном т сетка измельчается в тех подобластях, где т принимает большие значения. При расстановке узлов на границе Г использовалась та же управляющая функция, что и в уравнении (6.1).

Уравнение (6.1) аппроксимировалось на 5-точечном шаблоне. Полученная система разностных уравнений решалась итерационным методом переменных направлений.

На рис. 4 показана сетка, покрывающая О, с числом узлов N1 = 61, N = 23, а = 1000, где а — параметр в управляющей функции, которая зависела от ширины русла В:

т(х1,х2) = 1 + а/В (х1). (6.2)

Рис. 4. Расчетная сетка.

7. Определение глубины в узлах сетки

После того как сетка построена, необходимо вычислить глубину во всех ее узлах. Пусть Р(х*, у*) — один из узлов сетки. Так как О представляет собой объединение выпуклых четырехугольников, то вначале надо определить, какому из этих четырехугольников принадлежит Р. Для этого производится разбиение четырехугольника Л, В, В,+1Л,+1 на два треугольника и проверяются условия принадлежности точки Р каждому из них. Точка Р будет, например, лежать в треугольнике Л,В,Л,+1, если одновременно будут выполняться три неравенства:

Ьл,в, (Р) ■ ЬА,в, (Л,-+1) > 0, Ьв,, (Р) ■ Ьв,, (Л,-) > 0,

Ь, л, (Р) ■ Ьл,+1 л, (В,) > 0, (7.1)

где через Ь обозначены выражения вида

Ьл, в,(Р) = -(у* - у л, )(хв, - хл,) + (х* - хл, )(Ув, - у л,).

Пусть, например, Р £ Л,В,В,+1Л,+1. Тогда точке Р соответствует единственная пара параметров (р*, д*). Знание значений этих параметров необходимо для вычисления аппликаты (и тем самым глубины) точки Р по формуле (5.3). Определить эти параметры

можно из соотношений (5.1), в левой части которых следует использовать известные значения жа Исключая из (5.1) параметр р, получим квадратное уравнение для определения

д*:

ад2 + Ьд + с = 0, (7.2)

где _^ _^ _^ _^ _^ _

а = [——В, х Л—1 —я^], Ь = ,+1рР х —В,] — [——-Р х —

с = [—Р х ——Р].

Легко показать, что если Р £ А,В,В,+1А,+1, то уравнение (7.2) имеет единственный корень д*, удовлетворяющий условию 0 < д* < 1.

Пусть М, £ А,В,, М,+1 £ А,+1В,+1 — точки, соответствующие при отображении (5.1) параметрам д = д*, р = 0 и д = д*, р = 1. Тогда точке Р будет соответствовать параметр р*, определяемый по формуле

р* = |—МI / |М—М,|. (7.3)

Далее на основе полученных значений р*, д* и используя формулу (5.3) вычисляем глубину в узле Р.

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

Рис. 5. Изолинии глубин.

8. Результаты расчетов

На рисунках 6 и 7 представлены некоторые результаты расчетов для водоема с неровным дном при тех же значениях параметров, что и в статье [10]:

к = 1.19 ■ 10-4, С = 45 м1/2с-1, ща = 0.5 м ■ с-1, (8.1)

где ип — скорость втекания жидкости через участок Г1. Вектор потока массы а1 и расход из (1.3) задавались следующим образом:

а1(ж) = (Л(ж)ип; 0), ^(ж) = Л(ж)и,оиъ

хеГ1 хеГ2

Рис. 6. Линии тока установившегося течения.

Рис. 7. Поле вектора скорости установившегося течения.

где постоянная иои определялась из условия равенства прихода жидкости через Г1 и ее расхода через Г2:

I щ^Г + I(ах ■ п)^Г = 0.

Г2 Г1

В качестве точки М0, в которой по формуле (1.6) задается полная глубина, была взята точка пересечения входа Г1 и стенки Г^. В этой точке полагали Н0 = Н(М0). Глубина дна около боковых стенок Г^ и Г^' была равна 0.5 м, а в средней части водоема — около 3 м.

Как видно из рис. 7, скорость жидкости велика в местах сужения потока, в особенности около входа и выхода. В других частях водоема поток замедляется, около левого берега образуются застойные зоны. Отметим, что в отличие от работы [10] в наших расчетах по модели идеальной жидкости замкнутые линии тока не возникали. Заметим также, что варьирование параметров (8.1) в довольно широких пределах не ведет к качественным изменениям в поведении линий тока. В частности, при отсутствии трения о дно характер линий тока остается таким же, как на рис. 6. С другой стороны, изменение функции г = -Н(х), описывающей геометрию дна, приводит к существенной перестройке течения.

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

[1] Андреев В. К., КАпцов О. Б., Пухначев В. В., Родионов А. А. Применение теории групповых методов в гидродинамике. Наука, Новосибирск, 1994.

[2] МАРЧУК Г. И. Методы вычислительной математики. Наука, М., 1980.

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

[4] РоУЧ П. Вычислительная гидродинамика. Мир, М., 1980.

[5] Сретенский Л. Н. Теория волновых движений жидкости. Наука, М., 1977.

[6] ТАРУнин Е. Л. О выборе аппроксимационной формулы для вихря скорости на твердой границе при решении задач динамики вязкой жидкости. Числен, методы мех. сплошной среды, 9, №7, 1978, 97-111.

[7] ХАКимЗЯнов Г. С., ШокинА Н. Ю. О методе эквираспределения для построения двумерных адаптивных сеток. В "Вычислительные технологии", ИВТ СО РАН, Новосибирск, 4, №13, 1995, 271-282.

[8] ХАКимЗЯнов Г. С., ШокинА Н. Ю. Численное моделирование на адаптивных сетках двумерных установившихся внутриканаловых течений идеальной жидкости. Там же, 283-294.

[9] ХАкимзянов Г. С., Яушев И. К. О расчете давления в двумерных стационарных задачах динамики идеальной жидкости. "Журн. вычис. мат. и мат. физики", 24, №10, 1984, 1557-1564.

[10] BORTHWICK A. G.L., Barber R. W. River and reservoir flow modelling using the transformed shallow water equations. In "Int. J. Numer. Meth. in Fluids", 14, 1992, 11931217.

[11] CHRISTOV C. I. Orthogonal coordinate meshes with managable jacobian. In "Numer. Grid Generation, Appl. Math, and Computation", 10/11, 1982, 885-894.

[12] Richards C.W., Crane C.M. Pressure marching schemes that work. In "Int. J. Numeric. Meth. Eng.", 15, 1980, 599-610.

[13] Yeung R. W. Numerical methods in free-surface flows. In "Ann. Rev. Fluid Mech.", 14, 1982, 395-442.

Поступила в редакцию 15 сентября 1996 г.