CC BY
242
УДК 517.958:532.5 ББК 22.161.68
С 77
Стародумов И.О.
Аспирант кафедры математической физики института математики и компьютерных наук Уральского федерального университета, Екатеринбург, e-mail: ilyastarodumov@yahoo.ca
Кинематическое уравнение Сен-Венана. Метод решения
(Рецензирована)
Аннотация
Рассматриваются одномерная дифференциальная модель Сен-Венана и ее кинематическое приближение, описывающие свободное движение воды в открытом русле. Для кинематического приближения, являющегося квазилинейным уравнением первого порядка, в стационарном случае существует общее аналитическое решение; для нестационарного случая потребовалось разработать алгоритм поиска . .
Ключевые слова: гидравлика, уравнение Сен-Венана, кинематическое приближение, численное .
Starodumov I.O.
Post-graduate student of Mathematic Physics Department of Institute of Mathematics and Computer Sciences, Ural Federal University, Ekaterinburgh, e-mail: ilyastarodumov@yahoo.ca
The kinematic Saint-Venant equation. The solution method
Abstract
The paper examines one-dimensional differential Saint-Venant model represented with diffusion and kinematic approximation, describing free water dynamic in open bed. For the kinematic approximation, which is a quasi-linear equation of the first order, there is a general analytical solution in stationary case. For the non-stationary case we developed a numerical algorithm. In the paper, the solution results are presented.
Keywords: hydraulics, the Saint-Venant equation, kinematic approximation, numeric method.
В работах [1, 2] для моделирования гидравлического режима североамериканской Красной реки* вполне успешно применялась программа MIKE 11, разработанная Датским Гидродинамическим Институтом [http://www.dhisoftware.com]. В модуле MIKE 11 для расчета неустановившегося водного потока решается краевая задача для одномерной системы дифференциальных уравнений Сен-Венана в полном, диффузионном и кинематическом приближении.
Обратимся к модели Сен-Венана для течения воды в открытом русле [2, 3]:
dA + dQ dt dx
dQ + ±
dt dx
Q2
+gA\ —+Sf I - o.
(1)
(2)
Здесь х - протяженность по водотоку; t - время; A(x, t) - площадь живого сечения; Q (х, t) - расход воды; q (х, t) - внешний источник, отнесенный к единице длины водотока; g - гравитационное ускорение; h (х, t) - глубина; Sf (х, t) - уклон трения потока.
Глубина и уклон русла определяются следующим образом (см. рис. 1):
h = z — zъ ,
где г - высота уровня свободной поверхности (зеркала); zъ - высота русла.
(3)
Не путать с рекой Хонгха (Вьетнам, Китай), название которой также переводится как Крас-
.
5с == ї§а-дх
(4)
Рис. 1. Параметры водного потока С учетом (3), (4) уравнение движения (2) можно переписать в виде:
дЯ+А
дґ дх
г0_л
KAJ
+ gA^ + gA( - Sс ) = с.
дх
(5)
Для расчета величины уклона трения потока на практике, как правило, используют формулу [1]:
_ ем
г а2с2я
(6)
где Я - гидравлический радиус, равный отношению площади живого сечения А к смоченному периметру х этого сечения; С - коэффициент Шези, зависящий главным образом от шероховатости русла п. Часто для определения коэффициента Шези пользуются формулой Маннинга [1, 4]:
С =1 Я1/6. (7)
п
Будем считать пренебрежимо малыми и отбросим в уравнении (5) три первых члена. Тогда 5с = 5г и (6) равносильна формуле закона Шези для стационарного потока
Я
= С^ЩЁ 88И(5с).
(8)
Выразим е из (8) и подставим полученное выражение в балансовое уравнение (1). Вновь полагая А = 'мк, Я = к, получим кинематическое волновое уравнение
'#01 '
д(м>к) д
дґ дх
к5'3
Э8П (5с ) = Ч
(9)
Уравнение (9) суть квазилинейное уравнение в частных производных первого порядка:
К1 (х,ґ,к)дк(хҐ) + К2 (х,ґ,к)дк(хҐ) = К3 (х,ґ,к).
(10)
дх 2 4 ' дг
В стационарном случае (9) - это обыкновенное уравнение с разделяющимися переменными:
ёх
у/$о
w
к5'3
\
(11)
частное решение которого, удовлетворяющее условию к(х) = к, имеет вид (см.
*х=х1ф
рис. 2, 3):
h=
Г f
n
w
|qdx
+
- |qdx
W
3/5
=xieft;
(12)
Рис. 2. Характер h(x) в зависимости от знака q при постоянных w , n , S0, q Ш)
Рис. 3. Характер h(x) при постоянных w , S0 , нулевом q и n(x) = a + bsin(x)
Нестационарный случай. Будем полагать, что величины w, n, S0, q в (9) зависят только от пространственной переменной x и не зависят от времени t, кроме того, будем считать ширину w, шероховатость n и уклон S0 всюду положительными. В этом случае
к =-
З VSq"w
З n
h2^, K2 = w, КЗ = q -
1 дS0 S0 dn + S0 dw
2 dx n dx w dx
w
h3.
(1З)
Задача І.
w(х) = S0 (х) = 1, п(х) = 5/3, Ч(х) = 0, к(х,0) = 4-8ш(х), к(0,ґ) = 4 + 8Іи(ґ).
Уравнение (9): ^+*(*, Л*3 ^ = с.
дґ дх
Интегральная поверхность показана на рисунке 4. Задача 2.
х1е/і = 0, хг^ы , w, 50, п - константы, ч = 0,
h(О,t) = h1 при t > О, h(x^) = h2 при x > О.
v (q\ dh(x,t)+(З ws[S0Л
Уравнение (9): w—-—- + dt
h (x, t)
2/3 dh(x,t)
dx
= О
Аналитическое решение:
t=
3h1 2/3 —^=x + s, s > О 51
3h2~2^—^=tx + s), s < О
52 Vv
(14)
Аґ =
При ^ = 0 вслед за к функция г терпит разрыв, скачком изменяясь на величину 3 п
——^х( 2^3 -к2 2/3). Данное обстоятельство приводит к появлению в первом
квадранте плоскости {(х, г)} центрального угла (сектора), внутри которого обратная функция 5 (х, г) и, как следствие, функция к (х, г) либо не определены, либо неоднозначны (см. рис. 5).
Поскольку функция к (х. () имеет конкретный физический смысл -глубина, в указанном секторе ее естественно определить однозначно.
Функцию 5 (х, г) можно доопределить в секторе по непрерывности значением 0, в таком случае функцию к (х, г) в секторе можно считать постоянной и равной к(0) = к0. Выбрать константу к0, по-видимому, лучше всего из соответствующего условия баланса.
Представим линейный участок реки в виде цепочки из открытых резервуаров для воды со свободной поверхностью, испытывающей всюду одинаковое атмосферное давление (см. рис. 6). Каждый резервуар г имеет конечную длину Ахг, он изготовлен из непроницаемого материала, у него достаточно высокие гладкие вертикальные стенки и прямоугольное (у ), всюду одинаково шероховатое (пг) дно, имеющее нулевой крен и не-
Рис. 5. Решение системы (14)
Рис. 6. Цепочка резервуаров
большой положительный уклон (0г). В каждом резервуаре есть два входных отверстия
в левой стенке - через них в резервуар вода закачивается, а также одно выходное отверстие в правой стенке - через него из резервуара вода откачивается. Выход любого резервуара г, кроме последнего, через емкость-накопитель соединен с одним из входов (на рис. 6 он верхний) резервуара г +1, насос, установленный между этими резервуарами, может перекачивать воду из г -го резервуара в накопитель, а также из накопителя в г +1 -ый резервуар в любом заданном режиме.
Режимы работы насосов: насоса, закачивающего воду через верхний вход первого резервуара; насосов, перекачивающих воду через емкости-накопители; насосов, закачивающих воду в резервуары через их нижние входы, а также насоса, откачивающего воду из последнего резервуара, - согласуем в предположении, что для любого резервуара г справедливы условия задачи 2, то есть:
1. Движение воды в резервуаре описывается уравнением (9) и формулами (7), (8);
2. На момент времени ^ резервуар г наполнен до уровня г к . Для этого достаточно, чтобы до момента ^ в точке х{ (у левой стенки) долгое время выполнялось условие
я
у-0 + у0
Ахі =(і]к^
5/3
w,
у-0 +с У-0
Ахі - константа интенсивности притока (рас-
хода) воды в резервуар г вплоть до момента г^, а в точке х{ + Ахг (у правой стенки) было обеспечено условие свободного распространения (вода свободно стекает в накопитель);
п
3. С момента времени г. + 0 и до момента г. + АtJ■ в резервуар г закачивается вода, при этом насосы работают так, что в течение времени Аг. в точке х{ глубина воды остается постоянной, равной к. Значение к можно найти из соотношения
И 5/3 Углр0 г ■ ■ ..
— -----, где Ах. - константа интенсивности притока воды в
пг
резервуар г с момента времени г. + 0 и до момента г. + Аг.. В точке х{ + Ах. все это время обеспечивается условие свободного распространения, которому соответствует отток
воды в накопитель с интенсивностью
я+1-0у (ґ) = к( + Ах,,Ґ)5/3 —
4. Отыскав функцию к ( х, ґ) - решение задачи 2 для { хі < х < хі + Ахі ; ґу < ґ < ґу + Аґу},
ґ у +Аґ у Г7Т~ ґ у +Аґ у
1 3 3 1 W /Ол . ^ ^ І
можно найти Я'+1у =-------- І Я'+1-0у (ґ)Лґ =-——- | к(хі + Ах;.,ґ) Лґ - среднее по
Аґу Аґу п' времени значение расхода у правой стенки резервуара і . С этим постоянным расходом воды через верхний вход мы будем искать решение задачи 2 для резервуара і +1. Доба-
1 хі+2 0+1
вив к Я'+1у постоянную величину д+1уАхі+1 =-------- | | д(х,ґ)ЛґЛх - расход от независи-
у х+1 і у
мого источника, мы получим общую интенсивность притока (расхода) воды в резервуар і +1 на промежутке с ґу + 0 по ґу +Аґу. Повторим пп. 3, 4 для резервуара і +1. Подобным образом на промежутке времени с ґу + 0 по ґу + Аґу можно решить задачу 2 для всех резервуаров, то есть для всего участка реки;
1 хі+Ахі
5. Зная функцию к (х, ґ), можно найти у +1к =---- | к (х, ґу+1 )х - ее среднее в пре-
' хі
делах любого резервуара і значение на момент ґу+1, соответствующее установившемуся режиму, подобно описанному в п. 2. Мы пришли к ситуации для ґу+1, аналогичной ситуации на момент ґу ... Таким образом, повторяя пп. 3-5, можно определить гидравлическое состояние реки на любой сколь угодно большой момент времени.
Прямоугольную область Б изменения независимых переменных х и ґ покроем счетной сеткой, равномерной как по х , так и по ґ. Каждую ячейку сетки будем идентифицировать индексами ее левого нижнего узла: і - индекс по х и у - индекс по ґ; количество интервалов по х будем считать равным I, количество интервалов по ґ равным З; положение узла сетки в Б будем задавать его пространственной координатой хі , іє{0,1,...і} и его временной координатой ґу, у є {0,1,...7}.
В пределах ячейки (і, у), і є{0,1,...і -1}, у є {0,1,... З -1}, нижнее ребро будем называть ребром с индексом і , верхнее - ребром с индексом і +1 - 0, левое - ребром с индексом у, правое - ребром с индексом у +1 - 0 . Длина горизонтальных ребер ячейки равна АХ' = хі+1 - хі и не зависит от і , длина вертикальных ребер равна Аґу = ґу+1 - ґу и не зависит от у.
Будем считать исходные данные w(x, ґ), п (х, ґ), 50 (х, ґ) положительными и постоянными внутри любой ячейки (і, у) функциями, то есть кусочно-постоянными функциями в пределах области Б. Значения w, п и 50 внутри ячейки (і, у) будем обозначать wi, п и ^0'. Внешний источник д (х, ґ) будем считать сосредоточенным на верти-
[0, х Ф хг г.
кальных ребрах ячеек, так что д(х,г) = < , , где () - кусочно-постоянная
[д (г), х хг
функция, принимающая постоянные значения д на интервале г. < г < г.+1.
Рассмотрим отдельную ячейку. Будем считать глубину к (х, г) на вертикальном ребре . известной константой кг > 0, а на горизонтальном ребре г известной константой .к >0. В ситуации кг = г к внутри ячейки и на ребрах г +1 -0, . +1 -0 выполняется к = кг = .к. В ситуации же кг ф г к возможны три случая (см. рис. 7):
(*„*,) (М) <М)
I П Ш
Рис. 7. Случаи изменения глубины в пределах ячейки расчетной сетки Тестовые задачи
Ниже приведены графики численных решений тестовых задач. Все задачи решались для цепочки из 10 резервуаров длиной 665 единиц.
1. Задачи с условиями, согласованными с рисунками 2 и 3:
W г.х
Рис. 8 Рис. 9
На рисунке 8 представлен график решения задачи: wt = S0j = 1, = 5/3, qt = 0,
10h = h0j = 1. Решения во всех 10 резервуарах совпадают и равны 1 на всех ребрах во всех счетных ячейках.
На рисунке 9 представлен график решения задачи: wt = S0j = 1, = 5/3 , qt = 1/1000,
10h = h0j = 1. За несколько первых временных шагов решения h (t) во всех резервуарах устанавливаются (цветные линии, соответствующие этим решениям, становятся горизонтальными). Построив на участке установления решений h(t) график h(x) (тонкая черная линия), можно отметить совпадение характера этой кривой с характером соответствующей аналитической кривой на рисунке 2.
На рисунке 10 представлен график решения задачи, отличающейся от предыдущей лишь значением qt =-1/1000. Здесь также за несколько первых шагов решения h (t) во всех резервуарах устанавливаются, а построенная на участке установления этих решений кривая h (x) (тонкая черная линия) по характеру напоминает соответствующую кривую на рисунке 2.
На рисунке 11 представлен график решения задачи: wt = S0j = 1, ni = 5/3 + sin (Ax,-),
qt = 0, 10h = hQ] = 1. Решения h (t) в этом случае также быстро устанавливаются. Кривая h(x), построенная на участке установления всех h (t), напоминает кривую на рисунке 3.
Рис. 10
2. Нестационарная задача 1:
Рис. 11
i0
Рис. 12
На рисунке 12 представлен график решения задачи: = £0г- = 1, п1 = 5/3, = 0,
к = 4 - вш^Ах,), к0' = 4 + вт^'А. . Решения к (г) хорошо согласуются с аналитическим решением задачи 1, они относительно быстро выравнивают фазу своих колебаний и далее во всех резервуарах колебания к (г) происходят практически синфазно.
Заключение
Основным результатом настоящей работы является разработка оригинального алгоритма численного решения одномерного уравнения (9). Созданный алгоритм устойчив, поскольку для решения задачи требуется всего один раз в определенном порядке перебрать все ячейки счетной сетки. Поскольку (9) является приближением модели руслового стока (1)-(2), созданный алгоритм можно рекомендовать для поиска начального приближения при их численном решении.
Примечания:
1. Альтшуль АД. Гидравлическое сопротивление. М.: Недра, 1970. 216 с.
2. Klohn-Crippen. Red River one-dimensional unsteady flow model: final report submitted to International Joint Commission. Richmond (British Columbia), 1999. May. 88 pp.
3. Ahmad S., Simonovic S.P. Comparison of OneDimensional and Two-Dimensional Hydrodynamic Modeling Approaches For Red River Basin, final report to International Joint Commission. Winnipeg: University of Manitoba, 1999. December. 52 pp.
4. Христанович С .А. Неустановившееся движение в каналах и реках // Некоторые вопросы механики сплошной среды: сб. науч. тр. М.; Л., 1938. . 13-153.
References:
1. Altshul A.D. Hydraulic resistance. M.: Nedra, 1970. 216 pp.
2. Klohn-Crippen. Red River one-dimensional unsteady flow model: final report submitted to International Joint Commission. Richmond (British Columbia), 1999. May. 88 pp.
3. Ahmad S., Simonovic S.P. Comparison of OneDimensional and Two-Dimensional Hydrodynamic Modeling Approaches For Red River Basin, final report to International Joint Commission. Winnipeg: University of Manitoba, 1999. December. 52 pp.
4. Khristanovich S.A. Unsteady motion in channels and rivers // Some problems of mechanics of continuum: coll. of proceedings. M.; L., 1938. P. 13153.
