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

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

Том 9, № 2, 2004

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ СТРАТИФИЦИРОВАННЫХ ТЕЧЕНИЙ В СИСТЕМАХ ОТКРЫТЫХ РУСЕЛ И ВОДОЕМАХ РАЗВЕТВЛЕННОЙ ФОРМЫ*

О. Ф. Васильев Институт водных и экологических проблем СО РАН,

Новосибирск, Россия e-mail: vasiliev@ad-sbras.nsc.ru

А. Ф. Воеводин, В. С. Никифоровская Институт гидродинамики СО РАН, Новосибирск, Россия e-mail: voevodin@hydro.nsc.ru

The paper is concerned with 2D hydrodynamic models of flows and related hydrophysical processes in river channels, lakes and reservoirs in particular for the case of stratified flows. Such models are used in the physical limnology and river hydrology for a study of flow patterns, e.g. in the branched areas, when the problem under certain assumptions may be reduced to the 2-dimensional case. Special attention is given to the modelling of stratified bodies of water with a complex configuration, which can be represented (or composed) as a system of longitudinal and meridional arms. On the basis of the fractional step method and using implicit schemes we have developed a numerical model of the hydrophysical processes in bodies of water with a complex topological structure, which is described by a tree graph. Test calculations have shown the effectiveness of our mathematical model and numerical method.

1. Поперечно-осредненные двумерные модели течений (основные уравнения)

Если поперечные размеры русла (или водоема) достаточно малы по сравнению с его протяженностью, трехмерная (3D) гидродинамическая модель может быть аппроксимирована двумерной (2D продольно-вертикальной моделью). В частности, такой тип модели может применяться к стратифицированным проточным водоемам вытянутой формы, таким как вытянутые водохранилища или горные озера. Другая область их применения — это гидродинамика эстуариев. Основные уравнения двумерной продольно-вертикальной модели

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 03-05-65399) и Совета Программы поддержки ведущих научных школ (грант № 22.2003.5).

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

(ЬУМ) вертикального и продольного распределения плотности и скорости могут быть получены путем поперечного осреднения трехмерных динамических уравнений и уравнения неразрывности [1-3]. После осреднения основные уравнения в гидростатическом приближении выглядят следующим образом:

du du du д I . 1

t¡7 + + w— = -g— К + —

ot ox oz dx \ po

z

, . 1 д (1 du\ к

pdz\ + bdZ [bUt0~z — bT;

д (bu) + д (bw)

где

к =

E

i= 1

1 +

dx

0Ы dx

+

dz

дЫ

dz

(1.1)

(1.2)

1/2

A l l T = 8lulu-

Здесь х— продольная (горизонтальная) и вертикальная декартовы координаты; Ь — время; и(х, г, Ь), /ш(х, г, Ь) — горизонтальная и вертикальная составляющие скорости соответственно; Т(х, г, Ь) — температура; Б(х, г, Ь) — соленость; р(Т, Б) — плотность жидкости; р0 — характерное значение плотности; g — гравитационное ускорение; Ь = Ь1 + Ь2 — ширина русла, Ь1(х, г), Ь2(х, г) характеризуют положение и контур средней поверхности канала; ^ (х, Ь) — колебания свободной поверхности; V — коэффициент турбулентной вязкости; т — сопротивление трения (касательное напряжение) на боковых поверхностях; Л — коэффициент трения.

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

дТ дТ дТ дt дx дz

дБ дБ дБ — + u— + w— : дt дx дz

1L (ЬКт дТ

bдz \ дz

1A (bKs дБ

bJZ \ дz

+

1

poco

(Qt -

к

-qr

(1.3)

(1.4)

где Кт ,Кв — коэффициенты турбулентной температуропроводности и диффузии; р0с0 — произведение средней плотности и коэффициента теплопроводности; — поток тепла через боковую поверхность; Qт — поступление тепла от источника, или объемный приток тепла.

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

4 + ¥ = о- Q

CJt CJx

z

budz, ( — z0 = h,

(1.5)

zo

где В(х) = Ь(х, () — ширина русла на свободной поверхности; Q — общая скорость течения или расход; (х) — отметка дна; Н — глубина потока.

Используя уравнение (1.5), мы удовлетворяем кинематическому условию на свободной поверхности и на дне. Поэтому уравнение (1.5) часто называют кинематическим условием при решении системы (1.1), (1.2).

Для решения системы (1.1)—(1.5) необходимо определить коэффициент турбулентной вязкости V и коэффициенты диффузии Кт, К в. Самый простой способ — применение

0

2

2

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

1.1. Граничные условия

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

1) при г = (

Ои ОС ОС

Ч^т- = 0, ш = ^г + и—, Ог ОЬ Ох

г, ОТ ^ ОБ

Кт— = д, = 0;

Ог Ог

2) при г = го

Ои Ог о

V Ог = Кь|и|и, ш = и—,

ОТ ОБ

Кт— = 0, К^ = 0.

Ог Ог

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

В обоих случаях обычно рекомендуются следующие условия. Например, в первом случае на входе в водохранилище ставятся скорость притока (расход) Qin(í), температура воды Т(Ь) или плотность, и в предположении, что глубина к здесь известна, может быть задано также вертикальное распределение продольной скорости и;п(г,Ь). В случае эстуария такие условия могут быть заданы в сечении реки, расположенном далеко вверх по течению от устья. Плотность (или соленость) речной воды также должна быть задана.

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

и используемым типом граничных условий вниз по потоку. Даже для неприливных эстуариев при решении уравнения (1.1) и уравнения солевого баланса (1.4) на границе с морем (при x = L), где морские воды проникают в рассматриваемую область, необходимо задать значения солености S = Sin и, вообще говоря, скорости u = uin, заранее неизвестные [3].

Здесь представляется возможным ввести простое и естественное предположение, касающееся кинематического характера течения на выходе из эстуария в море (x = L), а именно допустить, что в этом случае имеет место преимущественно горизонтальное направление, и задать приблизительное значение вертикальной компоненты скорости вместо горизонтальной на этой открытой границе. Вертикальная компонента значительно меньше горизонтальной (w ^ u) . Таким образом, в простейшем случае можно предположить, что w = 0.

С точки зрения кинематики более точным было бы предположение о том, что на грат „ dzo нице x = L вертикальная компонента скорости w линейно изменяется от w = u—— на дне

dx

до w = dC/dt на свободной поверхности. Это предположение в сочетании с уравнением (1.2) позволяет проинтегрировать уравнение (1.1) без дополнительных предположений. Таким образом, условия на границе с морем при x = L будут следующими:

пт л г (л dZ¡n (z - ^ , dzo (С - z) п йч

С(L,t) = C¡n(t), w = + u; (16)

S = í Sin(t) j если u< 0j (1 7)

[не задается, если u > 0,

где Zin(t) — функция, определяющая колебания уровня моря в зависимости от времени.

Вернемся к постановке граничных условий вверх по потоку. При постановке условия на горизонтальную компоненту скорости uin(z,t), описанной ранее, в некоторых типичных для стратифицированных течений случаях (например, вторжение соленого клина) нельзя быть уверенным, что возмущения в нижнем течении не достигнут сечения, в котором вводятся условия. Имея дело со стратифицированными течениями, предпочтительно не задавать заранее профиль скорости. В некоторых случаях возвратные течения могут возникать даже в области сечения. Поэтому, используя в этом случае тот же подход, что и при постановке условий на границе вниз по потоку в задаче для эстуария, мы предлагаем следующие условия в сечении вверх по потоку:

Q = Qin(t); (1.8)

S =íPin(t)j если u> 0j (1 9)

не задается, если u < 0,

и "мягкое" кинематическое условие с уравнением неразрывности

w = dCin (z - zo) + ud£o (С - z) (1 10)

dt h dx h ' '

Возможно, в случае очень малого уклона русла последним членом можно пренебречь.

1.2. Условия сопряжения

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

Рис. 1.

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

Применяя основные законы сохранения и делая некоторые предположения о динамике и кинематике, мы получим следующие условия, которые должны выполняться в каждой точке узловой вертикали (или вертикали ветвления) [4]:

1. Баланс водных потоков:

^ = 0, г =1, 2, 3. (1.11)

2. Баланс потоков тепла и солености:

^biTiUi = 0, i =1, 2, 3; (1.12)

J^biSiUi = 0, i =1, 2, 3. (1.13)

3. Равенство давлений в рукавах в каждой точке узловой вертикали:

pi(z,t) = idem. (1.14)

4. Равенство уровней свободной поверхности на узловой вертикали для всех рукавов:

Сi(t) = idem. (1.15)

5. Мы будем предполагать, что вертикальное распределение вертикальной же компоненты скорости w линейно во всех ветвях на узловой вертикали (так называемое "мягкое" кинематическое условие). Пренебрегая уклоном дна в узле, мы можем записать

z - Zo

w = Wa---, (1.16)

h

где Wa = d(/dt.

В дополнение к этим условиям, из уравнения (1.5) будем иметь

1 ддг

Bi dxi

idem. (1-17)

Аналогичные условия можно получить из уравнения (1.12), используя предположение (1.16).

2. Численный метод

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

2.1. Построение сетки

Дискретная математическая модель одним из необходимых элементов включает построение разностной сетки. Исходная область для каждой секции при £ > 0

0 < х < Ь(т), 0 < г < ((т)(х,£),

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

0 < х < Ь(т), 0 < ц < 1

с помощью преобразования

z -

Z - z0m)'

(2.1)

Такое преобразование позволяет построить сетку с одинаковым шагом по оси п и шагом Дх(т) = Ь(т)/Ы(т), различными для каждой области. В результате для каждой секции получим сеточную область

п(т) Г х, = гДх(т), г = 0,..., М(т) н \ Пз = 3Дп, 3 = 0,... М

Шаг по оси £ является общим для всех секций и равняется т. Считаем, что при £ = 0 решение задано, т. е. для каждой секции известны

и(х, п, 0) = и0(х, п), ~ш(х, п, 0) = т°(х, п), ((х, 0) = (0(х).

Решение ищется на момент времени ¿к+1 = ¿к + т, к = 0,...

п

2.2. Расчет свободной поверхности

Для расчета ((т)(х.£) кроме уравнения (1.5) привлекается уравнение (1.1), проинтегрированное по глубине, которое перепишем в виде [4]

= * (2.2)

где

F = — ±

dx

Z

buudz + ( bv— dz

z о

z=Z

z=zo

A,

z z

— к ^ |—|—dz — g b 8pdz,

z о

zo

U

bd(, 5p

P — Po Po

Здесь и далее индекс т иногда будем опускать.

Уравнения (1.5) и (2.2) образуют замкнутую систему относительно Q(xi,tk) = Q/k и С) = ^, к = 1,... Разностная аппроксимация этих уравнений строится следующим образом. Вначале уравнение (1.5) аппроксимируем в виде

Z k+1(x) = Zk — т

0Q\

k+i

dx J

(2.3)

Используя соотношение (2.3), в уравнении (2.2) исключаем д(/дх и аппроксимируем производные разностными соотношениями. Из (2.2) получим следующие уравнения для определения неизвестных Q/¡+1 во внутренних точках области:

Qk+1 — Qk k —-— + guik

Zk — Zk Si+1 Si—1

Qk+1 — 2Q1+i + Qi-1 \

Ax2 J

k+1

k+1

2Ax

— т

Fk

(2.4)

i = 1,

N — 1.

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

После определения Q/¡+1 всех секций неизвестные вычисляются из соответствующей аппроксимации соотношений (2.3).

2.3. Расчет поля скоростей

Уравнения (1.1), (1.2) после преобразования (2.1) будут иметь вид

д— д— д— д bvt д—

bh— + bh—- + bw— = -—— ---hf;

ot ox on on h on

dbh— + dbw + bdZ о dx dn dt '

где

dh dz0

w = w — n—t;--;

dx dx

(2.5)

(2.6)

(2.7)

1

z

w=w—ndi;f = gbdx[Z+po 16Pdz

z

z

Для расчета поля горизонтальных скоростей u(xi, nj, tfc+i) = uk+1 (k = 0,... ,) уравнения (2.5) аппроксимируются с помощью следующего алгоритма расщепления. Конвективные слагаемые

du ди Ku = bhu——+ bw— = K1u + K2u дх дп

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

1 / du dbhu2\

Kiu =2 [bhu8X + —) ;

1 Л _du dbwu\ д(

Ku = 2 дП + + buW

(2.8) (2.9)

Представление конвективных слагаемых в виде (2.8), (2.9) позволяет построить так называемые энергетически нейтральные схемы без аппроксимационной вязкости [5, 6].

Далее расчет одного шага по времени разбивается на два этапа. На первом этапе на момент времени Ьк+1/2 = tk + т/2 уравнения (2.5) аппроксимируются во всех внутренних точках сетки П^ следующим образом:

k+1/2 k

bhkl Uij - uij | + Kfuk+1/2 = -Khukj + Dhukj - (hf)kj,

т/2

здесь обозначено

Khuu

i j

k+1/2 _ (bhu)k+l/2,juk+1 ,j (bhu)k_i/2 juk_i j

2Ax

K2huik j

(bw)k,j+1/2uk,j+1 - (bw)k,j_1/2ulj_1

2Ax

k+1 k т

D>wj = ( 'T

ui,j+1 ui,j \ I uut

i,j+1/2

(An)2

(Г

i,j-1/2

kk ui,j — ui,j_1

(An)2

fi+1 = (gv)k,j

k I Si+1

Z k+1 — z k+A

Si+1 Si_1

2Ax

+(gbh)k ±Apj,

p0 Ax

Pi,j Ax

zk+1 — zk+A м /, (x ) zk+1 — zk+1\

Si+1 Si_1 1 r k \д I x k / dz0(xi) Si+1 Si_1

^Ax— )An + ^ '

/ i=j

k i,j

dx

+ nj

2Ax

J +

hk+

1м

м

Z^+V <^_1,г An

2Ax\^ ~ri+1'1

i=j i=j

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

к+1/2 £ неизвестных , которая вместе с условиями сопряжения и граничными условиями об-

разует замкнутую систему. Метод решения таких систем уравнений (для фиксированных ] = 1,... ,М — 1) аналогичен алгоритму решения осредненных уравнений, рассмотренных выше.

k

Граничные условия для входных и выходных створов рассчитываются путем решения специального уравнения, которое является следствием уравнений (2.5), (2.6) и мягкого кинематического условия:

ЬЛ

ди

д_

дп

Ь^ ди Л дп

дЛ д - Л/ - ( Ь — + и^х^Л 1 и.

(2.10)

Особенность уравнения (2.10) заключается в том, что оно не содержит производных по х от искомой функции и при известных Лк, Лк+1 может быть проинтегрировано по вертикальной координате [3].

На втором этапе, т.е. для момента времени ¿к+1 = ¿к + т, уравнение (2.5) аппроксимируется во всех внутренних точках сеточной области с помощью неявной схемы по направлению п (индексу 3 ) следующего вида:

ьькк

ик+1 - ик+1/2*

т/2

+ К2\к+1 - Оик+1 = -К?ик+1/2 - (Л/)

(2.11)

Уравнения (2.11) решаются совместно с граничными условиями на дне и свободной поверхности

и

к+1

г,1

-и

к+1 г,0

Дп

ЙКЛ К'1 +1

и

к+1 г 0

+ ик

к+1

г 1

к+1

к+1

иг,м иг,м-1

Дп

методом прогонки.

На пограничных вертикалях г = 0 и г = N для вычисления ик+1, и^/,1 используются

N,3

аппроксимация уравнения (2.10) или балансовые соотношения (1.1) в зависимости от на-

(1) / , \ (2)

правления течения. Так, например, если на первой или второй секции (и^

к+1

и

N,3

неотрицательны, то для вычисления граничных значений и^3- используется уравнение типа (2.10), а граничное значение ^ик+^ для третьей секции вычисляется из балансовых соотношений. Аналогичная процедура используется и в других ситуациях.

После расчета горизонтальных скоростей во всей расчетной области рассчитывается поле вертикальных скоростей 7к+1. На пограничных вертикалях каждой секции используется кинематическое условие

7о,

п3

л к+1 40

- С,

т

WN^

/■ к+1 ZN

С к

Во внутренних точках аппроксимируется уравнение неразрывности (2.6) и значения 7гк+1 вычисляются по формулам

7?,к+1 = 7?,к+1 - Дп

= 7г'3-1 4Дх

(ЬЛи)

Л к+1 г+1,3+1

+ (ьли к+11>з+1- (ьли г-1,3+1- (ьли г-1

\ к+1

\ к+1

Значения поля скоростей 1 восстанавливаются из соотношений (2.7). На этом заканчивается расчет одного шага времени, следующие шаги рассчитываются аналогично.

к

0

2

п

3

2.4. Расчет поля плотности

Уравнения переноса тепла и солености после преобразования координат (2.1) переписываются в виде

д д д д ЬК д

Ж+ дХ ^ + апЬ** = вП-й^ + Г (2'12)

здесь под у(х,'ц,Ъ) понимается Т(х,гц,Ъ) или Б(х,'ц,Ъ), а правая часть Г(у,х,гц,^) вычисляется в соответствии с правой частью уравнения (1.3) или (1.4).

Граничные условия во входных створах, т. е. при х = 0, для первой или второй секций задаются обычно в виде

¡(¿'о('П^), если и(0,ц,Ь) > 0,

I не задается, если и(0,п,ь) < 0.

В выходном створе третьей секции, т.е. на границе х = Ь (река-море), соответственно

\ (¿Чп(п^), если и(Ь,ц,г) < 0, Ф(Ь,П,'Ч=\

I не задается, если и(Ь,п,ъ) > 0.

В узловых вертикалях для вычисления субстанции (температуры или солености) привлекаются балансовые соотношения (1.12) или (1.13)

3

^2(Ьи)трт = 0.

т=1

Для аппроксимации уравнений (2.12) используются монотонные разностные схемы с

направленными разностями. При этом расчет одного шага по времени разбивается на два

этапа. На первом этапе для момента времени ¿^+1/2 = ¿к + т/2 вычисляются величины к+1/2

, с помощью разностных уравнений

к+1/2 к — ,,

Т/2 I ' " ,3 ,3 1 " Тг ,3 1 "г ,,

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

'к+1 |ик+1Л /,/+1/2 тк+1/2

(Ь^)к+М /0 3 + К1ук+1/2 = -К2& + + г,, (2.13)

К ук,+1/2 = (Ь^Г1

"%,] I -г,з I Гг+1,3 +

2 )\ Ах 1

ик+1 + 1ик+1Л /тк+1/2 Ук+1/2)'

+ .

Ах

2

К* = Ьг

1 — 1|

г,3

Ах

+

+

IV

к+1

+ 1|

к+1/2

к+1/2

Ах

к

2

3

2

В <Рк,3 =

ьк к

к+1

3,3+1/2

1,3

(Ап)2

ьк, к

к+1

р.

1,3

3

(Ап)2

3,3-1/2

Уравнения (2.13) рассматриваются для всех ] = 1,... , М — 1; г = 1,... , N — 1.

На втором этапе рассчитываются значения 1 для момента времени ¿к+1 помощью разностных уравнений

(Ьк)

к+1 г

,.к+1 — ^к+1/2\

Ухту \ т^Ьгк+1

-Т/2- )+ К ^

— ^ ^

г,3

к+1/2 ^к

¿к + Т с

(2.14)

М — 1; г = 1,

N — 1.

Уравнения (2.14) рассматриваются для всех ] = 1,

Для вычисления граничных значений при ] = 0 и ] = М используется аппроксимация граничных условий на дне и свободной поверхности (1.11), (1.12)

— к,

^к+1 гг,0

Ап

як+1 Уг,0 ,

—к,

^к+1 гг,Ы

^к+1 гг,Ы -1

Ап

„к+1 %Ы ,

здесь потоки у^1 и ^^„М1 вычисляются в зависимости от того, какая из субстанций рассматривается — Т или

Для замыкания системы разностных уравнений (2.13), (2.14) и вычисления граничных значений при г = 0 и г = N для каждой секции, т. е. во входных и выходных створах и узлах сопряжения, используются граничные условия или условия сопряжения. При этом учитывается направление горизонтальной скорости потока.

Так, например, для первой секции:

1) при х = 0 (г = 0) имеем

{задается значение £0+1 или Т0к+1, если «к+1 > 0,

используется уравнение (2.13) или (2.14), если «к+1 < 0;

2) при х = Ь (г = N) имеем

I вычисляется б^/ или из условий сопряжения, если < 0,

(используется уравнение (2.13) или (2.14),

если > 0.

Аналогично вычисляются граничные значения для второй секции. При рассмотрении третьей секции аналогичная ситуация возникает при г = 0 и при г = N.

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

к

к

3. Алгоритм решения систем разностных уравнений

При расчете свободной поверхности, т. е. при решении проинтегрированных уравнений (1.1), (1.2), при расчете поля горизонтальных скоростей или субстанций на момент времени ¿к+1/2 возникает необходимость решения системы неявных разностных уравнений,

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

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

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

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

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

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

4. Примеры расчетов

Все эксперименты проводились с использованием трех секций (в дальнейшем будем называть это сопряжение "тройником"), причем считалось, что третья из них выходит в водоем с, вообще говоря, солеными водами (р < рв). Дно считалось горизонтальным, все секции имеют одинаковую длину 30 м, а глубина в различных экспериментах была порядка 0.15... 0.25 м. Для аппроксимации уравнений использовалась равномерная сетка по оси х с шагом Ах = 3 м, т. е. всего 21 узел, а по вертикальной оси п — равномерная сетка с шагом Ап = 0.1 (всего 11 узлов). Уравнение состояния включало зависимость плотности от солености р = р(Б). Кроме того, в большинстве расчетов (за исключением случая одиночного канала) полагалось Ь1 = 0.4 м, Ь2 = 0.6 м, Ь3 = 1 м, а для одиночного канала — Ь1 = Ь3 = 1 м.

z/h

1.0

0.8

0.6

III

0.4

0.2

1 U, м/с

0

0.02 0.04 0.06 0.08

Рис. 2. Профили горизонтальной скорости вблизи вертикали сопряжения для трех стационарных режимов: секции I, II — x = 27 м, секция III — x = 33 м.

Пример 1. Тестирование условий сопряжения в стационарных условиях.

Этот эксперимент был проведен для однородной жидкости в двух вариантах. Первый вариант — это расчет для одиночного русла, т. е. когда вторая секция исключалась, второй — для сопряжения трех секций. В первом случае в качестве граничных условий при x = 0 задан расход (м/с) Q1 = 1 • 10-3, во втором — заданы оба расхода Q1 = 0.4 • 10-4, Q2 = 0.6 • 10-3, т. е. соблюдалось равенство удельных расходов в первой и второй секциях. Как на входе, так и на выходе ставились мягкие кинематические условия. Начальное состояние — покой, во всех секциях Z(x) = 0.2 м. Расчеты выполнялись по двум сценариям, в каждом из которых на этапе I достигался стационарный режим при заданном расходе на входе в первую и вторую секции (либо только в первую). Далее задавалось изменение уровня на правой границе третьей секции с водоемом, причем в первом сценарии уровень поднимался за время At = 5 мин от отметки h\ = 0.2 м до h2 = 0.25 м, а во втором — опускался за такое же время от отметки h\ до h3 = 0.15 м.

Таким образом, в данном эксперименте мы имеем дело с тремя различными стационарными режимами. Эксперименты показали, что при равном суммарном расходе для "тройника" и одиночного русла картины течений оказались практически идентичными (различие между горизонтальными компонентами скоростей не превышало 10-4, а отличие уровней свободной поверхности имело порядок 10-5). На рис. 2 приведены профили горизонтальной скорости после выхода на стационарный режим для всех трех случаев. Уровень свободной поверхности при этом почти линейно убывал от правой границы к левой.

П р и м е р 2. Тестирование сопряженной модели при различных нестационарных режимах однородной жидкости. В этой серии экспериментов мы хотели проверить случай слияния секций с различными расходами и скоростями. Как и в предыдущем эксперименте, на левой границе первой и второй секций задавался расход, на правой границе с момента времени t\ до t2 (Ati = t2 —1\ = 5 мин) уровень линейно возрастал от отметки h1 = 0.2 м до h2 = 0.25 м, затем держался постоянным. В первом расчете в качестве граничных условий задавались Q1 = 0.4 • 10-3, Q2 = 0.6 • 10-3, т. е. соблюдалось равенство удельных расходов в первой и второй секциях (q1 = q2). Во втором расчете расходы "поменялись местами", т.е. теперь Q1 = 0.6 • 10-3, Q2 = 0.4 • 10-3 и равенство удельных расходов в первой и второй секциях уже не наблюдалось. Здесь мы не стреми-

Рис. 3. Профили горизонтальной скорости в случае Ь = 1.75, 5, и 10 мин (а — в соответственно).

Рис. 4. Профили горизонтальной скорости в случае Ь = 1.75, 5, и 10 мин (а — в соответственно).

лись выйти на стационарный режим. До начала подъема уровня расчет был выполнен на период 10 мин. На рис. 3, 4 приведены профили горизонтальной скорости для первого и второго расчетов соответственно на расстоянии одного шага сетки от вертикали сопряжения в три момента времени: на подъеме уровня, в момент окончания подъема и через 5 мин после окончания подъема.

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

(Б = 1).

Рис. 5. Картина изохалин.

Рис. 6. Профили горизонтальной скорости.

На рис. 5 приведены изохалины на момент времени £ =150 мин (установившийся режим), на рис. 6 — профили скоростей на тот же момент времени.

Анализ расчетов показывает, что при £ =150 мин устанавливается стационарный режим течения, при этом клин соленой воды проходит весь второй и третий участки, а на первом участке устанавливается в районе концевого створа.

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

[1] Blumberg A.F. Numerical model of estuarial circulation //J. Hydraul. Div., ASCE. 1977. Vol. 103, N 3. P. 295-310.

[2] Васильев О.Ф., Квон В.И., Чернышева Р.Т. Температурно-стратифицированное течение в водоеме вытянутой формы // Гидротехническое строительство. 1974. № 4. С. 35-43.

[3] VASILIEV O.F., Dumnov S.V. A two-dimensional mathematical model for salt water intrusion in an estuary // Proc. of the 20th IAMR Congress, Moscow, USSR. 1983. Vol. 2. P. 10-19.

[4] VASILIEV O.F. Vertical two-dimensional hydrodynamic models of water bodies: the state of the art and current issues // Proc. of the 5th Intern. Conf. on Hydro-Science and Eng., Warsaw, Poland, 2002. Vol. 5. P. 3.

[5] Астраханцев Г.П., Руховец Л.А. Дискретная гидродинамическая модель климатической циркуляции глубокого озера // Вычисл. процессы и системы. 1986. Т. 4. С. 135-178.

[6] Марчук Г.И., Дымников В.П., Залесный В.Б. и др. Математическое моделирование общей циркуляции атмосферы и океана. Л.: Гидрометеоиздат, 1984.

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

[8] Васильев О.Ф., Атавин А.А., Воеводин А.Ф. Методы расчета неустановившихся течений в системах открытых русел и каналов // Числ. методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ, ИТПМ. 1975. Т. 6, № 4. С. 21-30.

[9] Воеводин А.Ф., Шугрин С.М. Численные методы решения одномерных систем. Новосибирск: Наука, 1981.

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