Анализ системы уравнений Сен-Венана аналитическими и численными методами Текст научной статьи по специальности «Математика»

УДК 51-73:532.5.013

UDC 51-73:532.5.013

АНАЛИЗ СИСТЕМЫ УРАВНЕНИИ СЕН-ВЕНАНА АНАЛИТИЧЕСКИМИ И ЧИСЛЕННЫМИ МЕТОДАМИ

Семенчин Евгений Андреевич д. ф.-м. н., профессор

Кубанский государственный университет, Краснодар, Россия

Вандина Наталья Валерьевна

Ставропольский государственный университет, Ставрополь, Россия

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

Ключевые слова: СИСТЕМА УРАВНЕНИЙ СЕН-ВЕНАНА, УРАВНЕНИЕ КОНВЕКЦИИ-ДИФФУЗИИ, РАСХОД ПОТОКА, МЕТОД ПРОГОНКИ, МЕТОД РАСЩЕПЛЕНИЯ

THE ANALISIS OF SAINT-VENANT EQUATION SYSYEM ANALYTICAL AND NUMERICAL METHODS

Semenchin Evgeniy Andreevich Dr.Sci.(Phys.-Math.), professor

Kuban State University, Krasnodar, Russia

Vandina Natalia Valerievna

Stavropol State University, Stavropol, Russia

Diffusion-convection equation that has been received from Saint-Venant differential equation system describing nonstationary fluid motion in a river canal is investigated. Analytical method is considered for the solution of equation with the fixed factors and finite-difference method is considered for the solution of equation system with the float factors. The results of test calculations executed for a reaches of the river Kuban are presented

Keywords: SAINT-VENANT EQUATION SYSTEM, DIFFUSION-CONVECTION EQUATION, FLOW DISCHARGE, SWEEP METHOD, FISSION METHOD

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

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

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

дh 1 дQ

+

dt b dx

и динамического равновесия [1]:

g dt g dx dx K2

1 dv v dv dh T Q2

+---------------+ J + ^- = 0

(1.2)

при заданных начальных

Q(x,0) = a0 (x) , h(x,0) = 70 (x) ,

(1.3)

и граничных

Q(0, t) = ax(t), h(0, t) = 71 (t), Q(l, t) = a 2 (t), h(l, t) = 7 2 (t)

(1.4)

(1.5)

условиях, где t - время, 0 £ t £ T, T = const, x - пространственная координата, ориентированная по направлению движения, 0£x£l, 0, l -границы рассматриваемого бесприточного участка реки, J - уклон дна русла, g - ускорение свободного падения, b - ширина русла, K(h)-

расходная характеристика русла, h( x, t) - глубина русла, v (x, t)- средняя

скорость воды в сечении русла в точке x в момент времени t, Q(x, t)-

расход воды в указанном сечении.

Члены уравнения (1.2) при различных условиях движения жидкости

имеют различную относительную значимость [2]. Анализ паводков

показывает, что даже при значительном изменении во времени и по длине

реки скорости потока в период прохождения паводка величина двух

первых членов в уравнении (1.2) не превышает 10-5. В то же время,

Q2

величины уклонов дна J и трения для рек горных и предгорных

K2

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

dh+=0, (1.6)

dt b dx

I-7+f2=0• (1Л)

Этот факт следует учитывать при построении и изучении математических моделей потока воды в руслах горно-равнинных рек.

В работе Ж. А. Кюнжа, Ф. М. Холли, А. Вервея [2] указана возможность использования упрощенной системы уравнений Сен-Венана при исследовании и прогнозировании паводковых ситуаций, а также рассмотрено приведение системы (1.1) - (1.2) к уравнению конвекции-диффузии. Однако, в предлагаемой модели не учтена зависимость глубины потока h от времени t и не описаны методы и алгоритмы решения полученного уравнения. Поэтому цель нашей работы - предложить метод аналитического решения системы уравнений (1.6) - (1.7), указать условия применимости предлагаемого метода, разработать алгоритм численного решения рассматриваемой системы методом расщепления, выполнить расчеты характеристик потока по указанным методам на участках русла реки Кубань.

2. Аналитический метод решения системы уравнений Сен-Венана

Пусть в (1.6) b = const. Продифференцировав левую и правую части уравнения (1.6) по x, уравнения (1.7) - по t, и учитывая, что

dK dK dh dK ( 1 dQ Л

b dx

дt dh дt dh

\ /

систему уравнений Сен-Венана (1.6), (1.7) сведем к дифференциальному уравнению в частных производных второго порядка [2]:

dQ

dt

+

Q dK bK dh

dQ K2 d 2Q

dx 2bQ dx2

Уравнение (2.1) представляет собой дифференциальное уравнение

K2

конвекции-диффузии с коэффициентом диффузии --------------= а и скоростью

2bQ

Q dK т

конвекции для расхода —--------= b .

bK dh

Если выполняется одно из условий:

K2 d 1nK = C, C = const; K2 = C1h + C2,C12 = const, (2.2)

dh 1 2 I’2 V 7

то при малых h, h £ 1 м, уравнение (2.1) преобразуется к уравнению с постоянными коэффициентами (а = const, b = const):

= а ^ - b . (2.3)

dt dx2 dx

В этом случае может быть построено аналитическое решение уравнения (2.3), удовлетворяющее начальным и граничным условиям для расхода воды в русле реки (1.3) - (1.5).

Введем в рассмотрение функцию u( x, t), удовлетворяющую условию

Q( x, t) = q( x, t) + u (x, t),

x

где q(x,t) = a1(t) + у (a2(t) - a1(t)). Тогда задача построения Q(x,t)

сводится к задаче построения решения уравнения второго порядка с

постоянными коэффициентами

du _d2u 7-du N .N

^ = а ^7 - b — + f (x, t), (2.4)

dt dx dx

f (x, t)=a - b -®q

7 v ^ dx2 dx dt

и однородными граничными условиями

u (x,0) = j (x) = a0 (x), (2.5)

u (0, t) = 0, (2.6)

u (l, t) = 0, (2.7)

где о0 (х) - заданная функция.

Решение задачи (2.4) - (2.7) будем искать методом Фурье. Представим и ( х, t), / ( х, t), р ( х) в виде тригонометрических рядов:

и(х, t) = Х ип (t )• е2а Бт — х,

п=1

„ . ^ \ ■ рп „/л 2 '[/(%, t) . лп,.,

/ (x,t) =Е /п tt )• е2а ЯП — х, где /, (?) = Т^^Н1Ш;

п=1 I I 0 ^=х I

е

/ ч V- ^ . лп 21р(Х) • ш-„г

р( х) = Х Рп ■ е а 81ПТ" х, где рп =-т]^т^81П — .

п=1 I I 0 ^=х I

0 е2а

Вычислив значения частных производных их , и;, ихх и подставив в (2.4), после преобразований найдем:

и (

(х, () = X / е-4 (,-т (т )о?т + р(х) • е

\ Ь

— х . Лп

2 а

У

Бт —х, (2.8)

Ь

2

где к = ^^ + а 2а

3. Численный метод решения системы уравнений Сен-Венана

Рассмотрим начально-граничную задачу, представляющую собой объединений уравнений конвекции-диффузии (2.1) и уравнения неразрывности (1.1) при условиях (1.3) - (1.5).

Численное решение рассматриваемой задачи целесообразно находить с помощью неявной разностной схемы. Использовать явную схему нецелесообразно. При использовании явной схемы необходимо потребовать выполнения условия Куранта-Фридрихса-Леви [1], определяемого следующим соотношением между расчетным шагом по времени Дt и шагом по длине Д :

Ь

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

Построим решение рассматриваемой задачи методом расщепления ее по физическим процессам [3].

Разобьем интервал [0,Т ] на частичные интервалы (Г, Г+1], у = 0,1,..., п

точками t0 < t. <...< t , t0 = 0, t = Т. Учитывая, что

0 1 п ? 0 ? п ?

п

где п - коэффициент шероховатости, построим решение на интервале (Г, t]+1] задачи

Ь

К (И) = -к5/3,

дк 1 дQ Л — + —— = 0,

дt Ь дх

(3.1)

дt 3 Ьк дх ’

00 = о0, й0 = у0 при t = Г, у = 0;

31 = <31, И1 = к1 при t = tJ, у = 1, п ;

31 = о1, к} = у1 при х = 0;

3] = о 2, к] = у2 при х = I;

(3.2)

(3.3)

(3.4)

(3.5)

(3.6)

описывающей процесс конвекции для расхода [2], и задачи

дt Ь дх

(3.7)

10

дt 2 п~з дх

3 = 0\ ~0 = й1 при t = tJ, у = 0;

(3.8)

(~1 = QJ+1, И1 = И111 при ї = Г, і = 1,п; (3.10)

Q1 = а1, И1 = 71 при х = 0; (3.11)

Q1 = с2, И1 = 72 при х = I, (3.12)

описывающей процесс диффузии для расхода.

Учитывая аддитивность процессов конвекции и диффузии на достаточно малом интервале (ї , ї+1], найдем приближенное решение

общей задачи (1.1), (1.3) - (1.5), (2.1): Q = Q1+1, И = И1+1.

Также, решения задач (3.1) - (3.6) и (3.7) - (3.12) находим на каждом из интервалов (ї , ї+1 ].

Построим решения каждой из рассматриваемых задач (3.1) - (3.6) и

(3.7) - (3.12) конечно-разностным методом. Введем в области

В(0 < х < 1,0 < ї < Т) сетку Шт = Щ хЩ ={(х = ІЦ. = 1т),і=0,1,...Д, 1=0,1,...К),

= 1_ = Т *=N ’т = к ‘

Обратимся к задаче (3.1) - (3.6). Аппроксимируем производные в (3.1) - (3.2) разностными отношениями:

зе»а111 - д/; ЭН » иг - и ; эе»Щ+1 - аг.

Э? т Э т Эх 2^

Тогда, неявная разностная схема построения решения (3.1) - (3.6), определенная на четырехточечном шаблоне, имеет вид:

И1+1 - И/ 1 Q.+]+I - Q.;]+I А г г | -£^г+1 -£^г-1 0 т Ь 28 (3.13)

oiі+I - 0/ + 5 ё.1 шг - и»1 = 0 т 3 ЬН1 2^ ’ г (3.14)

г = 1,2,...,N -1, 1= 0,1,...,К -1; (3.15)

Ц0 = а 0 (х.), ъ10 = 70 (х.), І = 0,1,..., N; (3.16)

°+1 = а (1), И01+1 = 71 (ї1+1), (3.17)

а

]+1

N

О

І+1 ), *»+' = Уг1 ), 1 = 0,1,...,К - 1.

(3.18)

Перепишем (3.13) - (3.18) в матричном виде. Введем векторы

21 =

к1 V 1

3Ьк1О1

11

Тогда систему (3.13) - (3.18) можно записать в виде:

— ■ т ' 1 0" ' 0 Ь"

Л1 = — , В1 =

1 2 я , 5О1 0 у 5 г 3Ьк 1 0 V г у

2(

Л] 2 .]11 - ВІ2ІЛЛ - А1 2 ++1 =-FJ 1 1 -1 1 1 1 1 +1 1

\

О 0 (Х ) ЧУ0 (х1 ),

Г1 7 1+1

.

/

2,

+ 1

о 1 (|1>1) Г. (1)

121+1 \

2

+1

N

О2 (і1>1 )' УУ2 (|1>1 ),

(3.19)

(3.20)

Система (3.19) линейна относительно искомого вектора 2/+1. Поэтому ее решение на каждом из интервалов (і і ] можно искать методом линейной факторизации [4]:

2/+1 = а2,++1 + Д, 1 = 1,..., N -1, (3.21)

где а. - матрица размерности 2 х 2, Д - двухкомпонентный вектор. Прогоночные коэффициенты вычисляются по формулам:

а+1 =(А+д-в+1 )-1 • А+1, Д+1 =(А+д-В+1 )-1 • (А+1Ь + , (3.22)

1 = 1,2,...,N -1.

На первом этапе выполним прямую прогонку по рекуррентным формулам (3.22) и определим значения прогоночных коэффициентов

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

условиями а1 =

00

00

V У

, А, =

О

,Ы

уУ11+1 ),

. Затем, с помощью краевого условия

2

1+1

, выполним обратную прогонку и вычислим 2, 1+1,

(Г„ )Л У 2(^. )

V /

I = Ж -1,Ж - 2,...,1, по формуле (3.21).

Метод матричной прогонки устойчив по отношению к случайной ошибке, если для рассматриваемой задачи выполнены условия

в;1 а + в- (- А )||< 1, (3.23)

где ЫЦ = вир-

Лх

хф0 X

матричная норма, индуцированная векторной нормой

, введенной в линейном «-мерном пространстве. Для случая евклидовой

имеем ||А|| = ^1, где /лг - сингулярные числа

( ( и . 2 і 1/2 \

нормы X = X х1

11= 0 0

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

(л, > Л2 > .. > Л„ > 0).

Для рассматриваемой задачи (3.13) - (3.18) условие (3.23) устойчивости и корректности метода прогонки имеет вид:

3к/яЬ > 5tQlJ.

Аналогичным образом построим решение задачи (3.7) - (3.12). Производные в (3.7) - (3.8) аппроксимируем конечными разностями:

Эй » Q/+1 - й!. Э~ » ~1+1 - ~1 .

Эt Т ’ Эt Т ’

Эй»й,++1 - й/_г. »й;:: - 2~г+а,-г.

Эх 2^ ’ Эх2 я2

Тогда (3.7) - (3.12) примет вид:

к;+' - к1 1 а ++1 -

г______+_■£-'1 +1_-г^г-1

Т Ь 2я

0

аг1 - а; -1 ь(к/ )3 _ - 2.У1 + о„-+1 = 0

т 2 и О, 1

1 = 1,2,...,N -1, 1 = 0,1,...,К -1;

О: = О (х,), ~ = к(х,), 1 = 0,1,...,N;

Й+1 = 01 (|1+1 )> ~1+1 = У1 (|1+1) >

§Г = О,(І1+1), ~'+‘ = у21), 1 = 0,1,...,К -1.

(3.24)

(3.25)

(3.26)

(3.27)

(3.28)

(3.29)

Перепишем систему (3.24) - (3.29) в матричном виде:

А]2 '+1 - В]21+1 + С]2.І+1 =-К1.

г г-1 г г г г+1 г

20 =

у1+1 _

0

а

^+,) Г,((,+1)

у І+1 _

а

2 ((, +1 ) 72 ((, +1 )

(3.30)

(3.31)

где 21. =

' 1 0" /

аі И.1 , А 1 = — г 2л ьИ і )- 0 , І?1' = “ г

V г Vі 1 Л 0 V

0

Ь

п 2аі +—

ьИ і)3 0

С] = — г 2л

-1 0 ь(0 0

К1

Г Ь.1 Л

п

а,1 )2

Решение (3.30) как и (3.19) будем искать методом линейной факторизации:

2/+1 = а, 2% + Д, ,= 1,..., N -1, (3.32)

сД - матрица размерности 2 х 2, Д - двухкомпонентный вектор.

Рекуррентные формулы для прогоночных коэффициентов имеют

вид:

=(2: - А:д )Т1 • с,:,, 2+ = (2+ - А:д Г • (/?^, + Д+,Д), (3.33)

,= 1,2,...,N -1.

а

+1

Найдем (Х,п, Дп, с начальными условиями с~ =

00

00

V У

71

1 ((,+1) 1 ((1.1 )

Используя краевое условие 2+1 =

а

2((і+1 ) У72 ((і+1 )/

выполним обратную прогонку

и вычислим 2/+1, , = N -1,N - 2,...,1, во внутренних точках

рассматриваемой области по формуле (3.32).

Условие устойчивости и корректности метода матричной прогонки по отношению к случайной ошибке

В-1 А + В-1 (С . )||< 1

г г г \ г /||

в данном случае можно свести к выполнению неравенства Ъ8 > т.

Л'

4. Результаты расчетов, выполненных для некоторых участков русла реки Кубань

Применим предложенные выше методы решения системы уравнений Сен-Венана для расчета основных параметров паводковой волны в русле реки Кубань.

Вычислим эти характеристики с помощью аналитических методов, указанных в пункте 2, на участке, расположенном между гидропостами на хуторе Дегтяревском и в селе Успенском. Данный участок находится в предгорной зоне реки Кубань. Средняя величина уклона дна на

3 3

рассматриваемом участке изменяется от 2,7-10" до 1,8-10" , что позволяет использовать для описания движения паводковой волны уравнение (2.3).

В качестве исходных данных для определения Q(х,0) = о0(х), Q(0, г) = о1(г), Q(l, г) = о2(г) используем результаты измерения расхода потока, проведенные в июне 1985 года [5]. Для их обработки была использована прикладная статистическая программа Ва1аБ11:, позволяющая построить аппроксимирующую функцию с наибольшим значением коэффициента детерминации. Согласно проведенным расчетам о0 (х) = 6,0606 • 10-10 х2 - 2,739 • 10-4х + 46; о1(г) = 2,0552• 10-15г3 -7,663• 10-10г2 +1,4• 10-4г + 45,779; о2(г) = 1,956 • 10-15 г3 - 9,285 • 10-10 г2 + 2,294 • 10-4 г + 24,746.

Длина рассматриваемого участка по руслу реки Ь = 90000 м, ширина русла Ъ = 70 м. Коэффициент шероховатости находим по справочным таблицам [6]: п = 0,04.

Полученное по формуле (2.8) аналитическое решение позволяет найти значение расхода потока воды в произвольный момент времени г в любом интересующем створе рассматриваемого участка русла реки. На рисунке 1 представлен график решения задачи (2.3), (1.3) - (1.5) (период

Научный журнал КубГАУ, №64(10), 2010 года времени равен пяти дням).

Рис. 1. Зависимость расхода потока Q от времени t и координаты х

Алгоритм численного решения задачи (1.1), (1.3) - (1.5), (2.1), описанного в п.3, реализован в программном продукте, написанном в среде Мар1е. Расчетным выбран участок реки Кубань, длиной Ь = 100000м, расположенный между городом Армавиром и станицей Темижбекской. Ширина русла на рассматриваемом участке Ь = 140м, коэффициент шероховатости п = 0,04.

Для задания начально-граничных условий (1.3) - (1.5) были использованы результаты измерений расхода Q. Значения к были определены с помощью соотношения:

к = У0 + У - УЬ,

где у0 (х) - ноль поста относительно уровня моря, У(х, t) - функция, значение которой в момент t в точке х совпадает с уровнем воды в русле, уЬ (х) - отметка дна относительно уровня моря [5].

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

времени ? в различных точках х можно определить с помощью таблиц [5]. Отметку дна можно определить по результатам нивелирования, осуществляемого с помощью картографической системы Ооо§1еЕагШ [6].

Для рассматриваемого участка русла при т = 21600 с, ^ = 10000 м, краевые условия имеют вид:

^110 У108,7 У107,4 У106ЛУ104,8 У103,6 У102,3 У101 У 99,7 У98,5 лЛ

1,25 УЛ А

1,23

1,217

А ’

1,2

1,18

1,17

1,15

1,13

А ’ А

А А А А /V і іа VI іоЛАюоЛАї/тЛАї VI т V осп \\

1,117

А

1,1

У У

((

110

1,25

118

1,37

128

1,5

140

1,68

160

1,95

А А

А ’ А ’ А 98,5 У107 У116 У126 У140

193

2,4

250

3,19

1,1

1,22

1,34

А

1,49

А /У V 1 V

А

1,7

А

160

2,01

А

200

2,51

\\

/У

А ’ А

Результаты расчетов Q, ^, выполненные по схеме (3.1) - (3.12), представлены на рисунках 2, 3.

Рис. 2. Гидрографы потока Q(t) (на участке реки Кубань)

3,0 -2,8 -2,6 -2,4S 2’2'

-=! 2,0 -

1,8 -1.6 -1,4 -1,2 -О

Рис. 3. Зависимость глубины потока h от времени t (на участке реки Кубань)

Заключение

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Грушевский М. С. Неустановившееся движение воды в реках и каналах. - Л.: Гидрометеоиздат, 1982. - 288 с.

2. Кюнж Ж. А., Холли Ф. М., Вервей А. Численные методы в задачах речной гидравлики. - М.: Энергоатомиздат, 1985. - 256 с.

3. Марчук Г. И. Математическое моделирование в проблеме окружающей среды. - М.: Наука, 1982. - 320 с.

4. Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. - М.: Наука, 1978. - 592 с.

5. Государственный водный кадастр. Ежегодные данные о режиме и ресурсах поверхностных вод суши 1985 г. Выпуск 1. Бассейны рек северо-восточного побережья Черного моря, бассейн Кубани. - Обнинск, 1987. - Т. 1. - 204 с.

6. www.googleearth.ru