CC BY
45
6. Самарский А.А. Теория разностных схем. - М.: Наука, 1989.
7. Сухинов А.И., Тимофеева Е.Ф. Чистяков А.Е. Построение и исследование дискретной математической модели расчета прибрежных волновых процессов // Известия ЮФУ. Технические науки. - 2011. - № 8 (121). - С. 22-32.
8. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. - М.: Наука, 1978.
9. Коновалов А.Н. К теории попеременно-треугольного итерационного метода // Сибирский математический журнал. - 2002. - № 43:3. - С. 552-572.
10. Чистяков А.Е. Теоретические оценки ускорения и эффективности параллельной реализации ПТМ скорейшего спуска // Известия ЮФУ. Технические науки. - 2010. - № 6 (107). - С. 237-249.
Статью рекомендовал к опубликованию д.т.н., профессор Я.Е. Ромм.
Чистяков Александр Евгеньевич - Федеральное государственное автономное образовательное учреждение высшего профессионального образования «Южный федеральный университет»; e-mail: [email protected]; 347928, г. Таганрог, пер. Некрасовский, 44; тел.: 88634371606; кафедра высшей математики; доцент.
Семенякина Алена Александровна - e-mail: [email protected]; кафедра высшей математики; студентка.
Chistyakov Alexander Evgenjevich - Federal State-Owned Autonomy Educational Establishment of Higher Vocational Education "Southern Federal University"; e-mail: [email protected]; 44, Nekrasovskiy, Taganrog, 347928, Russia; phone: +78634371606; the department of higher mathematics; associate professor.
Semenyakina Alena Alexandrovna - e-mail: [email protected]; the department of higher mathematics; student.
УДК 532.5.031
А.Е. Чистяков, И.Ю. Кузнецова
ЗАДАЧА РАСЧЕТА И ПРОГНОЗИРОВАНИЯ ПРОЦЕССА ПОДЪЕМА УРОВНЯ ВОДЫ В МЕЛКОВОДНЫХ ВОДОЕМАХ
Для реконструкции сценария экологической катастрофы, а также для моделирования возможных вариантов биологической реабилитации водоема был создан ряд высокоточных математических моделей гидрофизических процессов в мелководных водоемах. Данные модели описывают движение водной среды с учетом различных факторов. Другое актуальное назначение высокоточных математических моделей гидрофизических процессов связано со своевременным предсказанием различных природных катаклизмов, связанных с изменением уровня воды: затоплением прибрежных районов, обмелением судоходных каналов и др. Прогнозирование процесса подъема уровня воды в Азовском море является актуальной задачей на данный момент.
Математическое моделирование; рельеф дна; прогнозирование.
A.E. Chistyakov, I.Y. Kuznetsova
PROBLEM OF CALCULATION AND PREDICTIN PROCESS OF LIFTING OF WATER LEVEL IN SHALLOW WATERS
For reconstruction of the scenario of ecological disaster, and also for modeling of possible options of biological rehabilitation of a reservoir was created high-precision mathematical models of hydrophysical processes in shallow waters. These models describe motion of the water environment taking into account various factors. Other actual purpose of high-precision mathematical
models of hydrophysical processes is connected with a timely prediction of the various natural cataclysms what change water level: coastal inundation, shoaling of navigation channels, etc. Prediction of process-level rise in the Sea of Azov is important task at the moment.
Mathematical modeling; underwater topography; prediction.
Введение. Для реконструкции сценария экологической катастрофы, а также для моделирования возможных вариантов биологической реабилитации водоема был создан ряд высокоточных математических моделей гидрофизических процессов в мелководных водоемах [1-2]. Данные модели описывают движение водной среды с учетом следующих факторов: ветровые течения и трение о дно, стоки рек, испарение, сила Кориолиса, турбулентный обмен, сгонно-нагонные явления, сложная геометрия дна и береговой линии. Другое актуальное назначение высокоточных математических моделей гидрофизических процессов связано со своевременным предсказанием различных природных катаклизмов, связанных с изменением уровня воды: затоплением прибрежных районов, обмелением судоходных каналов и др. Прогнозирование процесса подъема уровня воды в Азовском море является актуальной задачей на данный момент.
Постановка задачи. Исходными уравнениями гидродинамики являются [3-4]:
- уравнение движения (Навье-Стокса):
1 ' ' '
u't + uux + vuy + wu[ =--p'x + (juu'x )x +(ми'у) +{vu[ )z + 2Q(v sine-w cose),
1
v't + uv'x + vv'y + wv'z =-— p'y +(juv[ )x +{juv'y ) + (vv'z ) - 2Qu sine, (1)
1
w't + uw'x + vw'y + ww' =--p' + (juw'x )x + (juw'y) +(vw'z ) + 2Qu cos в ;
- уравнение неразрывности в случае переменной плотности:
P, +(pu )x +(pv )y +(pw ) = 0, (2)
где V = \u,v, w} - компоненты вектора скорости; p - превышение давления над
гидростатическим давлением невозмущенной жидкости; р - плотность; Q - угловая скорость вращения земли; в - угол между вектором угловой скорости и вертикалью; - горизонтальная и вертикальная составляющие коэффициента турбулентного обмена.
Система уравнений (1)-(2) рассматривается при следующих граничных условиях:
- на входе (устье рек Дон и Кубань):
u(x, y, z, t) = u(t), v(x, y, z, t) = v(t), p'n (x, y, z, t) = 0, V'n (x, y, z, t) = 0,
- боковая граница (берег и дно):
pj(u )n (x, y, z, t) = -тх (t), pj(v')n (x, y, z, t) = -ту (t),
Vn (x, y, z, t) = 0, p'n (x, y, z, t) = 0,
- верхняя граница:
pj(u\ (x, y, z, t) = -тх (t), pu(v )n (x, y, z, t) = -ту (t), (3)
w(x, y, t) = -ю- p / pg, p'n (x, y, t) = 0,
- на выходе (Керченский пролив):
p'n (x, y, z, t) = 0, Vn (x, y, z, t) = 0, где ю - интенсивность испарения жидкости; Tx,Ty - составляющие тангенциального напряжения (закон Ван-Дорна).
Составляющие тангенциального напряжения для свободной поверхности равны
тх =РаСр (И )мх\Щ , ту =раСр ,
где И - вектор скорости ветра относительно воды; ра - плотность атмосферы, 10,0088, х < 6,6 м / с
с, ). {0;
- безразмерный коэффициент. 0026, х > 6,6 м/с
Составляющие тангенциального напряжения для дна с учетом введенных обозначений могут быть записаны следующим образом:
^ =рср (И ж, Гу =РСР (V ж,
где ру - плотность донных отложений.
Рассмотренная ниже аппроксимация позволяет на основании измеренных пульсаций скоростей строить коэффициент вертикального турбулентного обмена, неоднородный по глубине [5]:
v = С2 А21, s 2 ^
( дй V
V dz у
(dV v
V у
(4)
где и ,И - осредненные по времени пульсации горизонтальных компонент скорости; А - характерный масштаб сетки; С5 - безразмерная эмпирическая константа,
значение которой обычно определяется на основе расчета процесса затухания однородной изотропной турбулентности.
Дискретная модель гидродинамики мелководных водоемов. Расчетная область вписана в параллелепипед. Для численной реализации дискретной математической модели поставленной задачи гидродинамики вводится равномерная сетка:
Ик ={гп = пт, х. = 1кх, у. = ]ку, гк = ккг; и = О^, I = 0.. Ых, ] = 0..Ыу, к = 0..Ыг; Игт = Т, N к = I , N к = I , N к = ¡} ,
t 'хх х^ у у у' г г г '
где т - шаг по времени; кх, ку кг - шаги по пространству; N - количество временных слоев; Т - верхняя граница по временной координате; Nx, Ny Nг - количество узлов по пространственным координатам; 1х, 1у 1г - границы по параллелепипеду в направлении осей Ох, Оу и Ог соответственно.
Для решения задачи гидродинамики использовался метод поправки к давлению [5]. Вариант данного метода в случае переменной плотности примет вид и - и
+ ии\ + Уи,, + Ии\ = ( /1и, ) + ( ци1
т
w - w
+ uux+ vuy + wuz =(juux )x + (juy) +(vu'z)i + 2Q.(v sinö- w cosö),
x x x
+ uvx + vvy + wvz=(u~Vx)x +(u~Vy) +(vv'x)i - 2Q.u sinö,
v - v
т
uw'x + vw'y + ww'z=(uw'x )x +(jw'y) +(vw[) + 2Q.u cosö, (5)
++=t-p+(M,+(pL+M,
Pxx +Pyy Pzz т т т т '
u - u 1 v - v 1 w - w 1
-=--IPx ; -=--Py ; ^— = -~Px ;
т p т p т p
где V = {и,V, и>} - компоненты вектора скорости; {й, V, н>},{й, V, м>} - компоненты полей вектора скорости на «новом» и промежуточном временных слоях соответственно; й = (й + и) / 2 ; р и р - распределение плотности водной среды на новом и
предыдущем временных слоях соответственно.
При построении дискретных математических моделей гидродинамики учитывалась «заполненность» контрольных ячеек, что позволяет повысить реальную точность решения в случае сложной геометрии исследуемой области за счет улучшения аппроксимации границы.
Через о; ] к обозначена «заполненность» ячейки (г, ], к). Степень «заполненности» ячейки определяется давлением столба жидкости внутри данной ячейки. Если среднее давление в узлах, которые относятся к вершинам рассматриваемой ячейки, больше давления столба жидкости внутри ячейки, то ячейка считается заполненной полностью (ог -к = 1). В общем случае «заполненность» ячеек можно вычислить по следующей формуле:
р + р + р + р
11, 1, к т 11 -1,1 к 1,1 -1, к^ -1,1 -1, к
\рфг
(6)
где Р = р + рёг - давление.
Вводятся коэффициенты 90, д1, 92, 93, 94, 95, 96, описывающие «заполненность» областей, находящихся в окрестности ячейки (контрольных областей). Значение 90 характеризует «заполненность» области Д0: { хе (х, -1, х1+1), у е( У;-1, +1),
г гк+1)}, 9 - Д: {х е(х, хт), У е(yi-l, Ум), г гк+1)}, 42 - Д2:
{хе (х,^х,), уе^-рУ1+1), ге ()}, 9з - Дз: { х^х^хт) уе(,у^) г гк+1 ) 94 - Д4: {х ^-п х+1), У Цу^ У.), г е^ гк+1)}, 95 - Д5: {х е (х-l,х{+1), у е(у.^ у.+1), г е(, гк+1) }, 9б - А: { хе(хх+1), у е(у.-1, у.+1), г е (гк-1, гк) }. Заполненные части областей Дт будем называть 0.т
, где т = 0...6. В соответствии с этим коэффициенты 9т можно вычислить по формулам
5Г
9т 1 к =, (91 \
ог+1,1, к + ог+1,1+1, к + ог+1,1, к+1 + ог+1, 1+1,к+1
92),,,,, = ^
к + ог, 1 + 1, к + ог, 1 ,к + 1 + ог, 1+1, к + 1
9з, 1, к =
ог+1,1+1,к + ог, 1+1, к + ог+1,1+1,к+1 + 1+1,к+1
\ и1,1 ,к 94, 1, к =-
4
', 1 ,к+1
95 )и.к =
ог, 1 ,к+1 + ог+1,1,к+1 + ог+1,1+1,к+1 + ог, 1+1,к+1
(96),
оц,к + ог+1,1,.к + ог+1,1+1,к + оц+1,к 4
, (9о , 1, к =1 ((91 ),,1,к +(92 ).
В случае граничных условий третьего рода
4
4
4
4
с'я (х, у, х) = апс + ря
дискретные аналоги операторов конвективного ис'х и диффузионного (/ис'х) переноса, полученные при помощи интегро-интерполяционного метода [7], учитывающие частичную «заполненность» ячеек, могут быть записаны в следующем виде [9-10]:
с , - с
(Чо ) ¡UC', = (l ) jui
+1/2,j
2h
-(Ч2 )ijui-
,j
-1, j
1/2, j
(Чо L [wX)'x = (1 L M
2hx
c, , - c.
+1/2, j
( Ч2 hj M-1/2,j"
,j
1, j
h2
(Ч1 L -(Ч2 )::
Mi,
a c. + ß
x i, j f x
Аналогичным образом запишутся аппроксимации по оставшимся координатным направлениям.
Погрешность аппроксимации математической модели равна О (г +1|й||21, где
\\Н\\ = ^1 Нх + Н2у + Н2г . Доказано сохранение потока на дискретном уровне разработанной гидродинамической модели, а также отсутствие неконсервативных дисси-пативных слагаемых, полученных в результате дискретизации системы уравнений. Достаточное условие устойчивости и монотонности разработанной модели определяется на основе принципа максимума [7] при ограничениях на шаг по пространственным координатам:
кх < \2р/ и| , ку < , кг < / или Яе < 2N ,
где Яе = |у| • I / /и - числа Рейнольдса; I - характерный размер области, N = тах ^х, Nу, Nг} .
Дискретные аналоги системы уравнений (5) решаются адаптивным модифицированным попеременно-треугольным итерационным методом вариационного типа [11-13].
На рис. 1 темным цветом показана функция возвышения уровня при западном ветре 5 м/с, рассчитанная на основе полной трехмерной модели гидродинамики.
Рис. 1. Функция возвышения уровня
x
Следует отметить, что расчет функции возвышения уровня на основе полной трехмерной модели гидродинамики весьма вычислительно трудоемок и возникает задача расчета функции возвышения уровня на основе сравнительно более простых и менее трудоемких алгоритмов.
Моделирование и прогнозирование динамики возвышения уровня на двумерной модельной задаче. Решение задачи гидродинамики в двумерном случае сводится к решению системы уравнений мелкой воды [5]. Согласно методу поправки к давлению [6], система уравнений мелкой воды решается согласно следующему алгоритму:
- вычисляется поле скорости без учета давления (функции возвышения уровня)
п+а п
(я + )))L_u_+(я + ф<+(я + фи; =
= ((я +))циХ) + (( + ))К ) +р-р,
1;п+ст - Vп
(я + + (я + +(я + = (7)
= ((я + + ((я + ж)' у -р,
- рассчитывается функция возвышения уровня
) - к, ((я+))))х - \ ((я+))))) =
= -* ((я + )ип+а) - * ((я + ))-)) , (8)
- уточняется поле скорости
п+1 _ п+а п+1 _ п+а
(Я + )) = -8 (Я + ))), (Я + )) = -* (Я + )))), (9)
к, к,
где ) - функция подъема уровня; V = {и, V} - вектор скорости движения водной среды; ц , 1] - коэффициенты турбулентного обмена по горизонтальному и вертикальному направлениям соответственно; * - ускорение свободного падения; р - плотность жидкости; тх ,т - тангенциальное напряжение на дне жидкости.
Система уравнений (7)-(9) рассматривается при следующих граничных условиях: и'п = 0, уп= 0, )) = 0.
Будем предполагать, что основной вклад в изменение функции возвышения уровня вносят ветровые тангенциальные напряжения и слагаемыми, описывающими турбулентный обмен и нелинейность движения водной среды, можно пренебречь. Таким образом, система уравнений (7) сводится так:
(Я + ) (ип+а - ип) = к, ТР , (Я + ) (а - V") = к, ТР .
Подставляя полученное выражение в (8) и принимая во внимание, что поле вектора скорости на п-м слое подчинено уравнению неразрывности, получим:
) - к, ((Я + ) - к, ((Я + ))))'у = -р(^ )) - *рр(у,р)у. (10)
Дискретный аналог уравнения (10) запишется следующим образом:
£n+1 _ £n Zn+a _ Zn+a
(Чо )i, j - h (Ч1 )i, j (^i+1/2,j + 5n+1/2, j )+
(11)
Zn+a _ pn+a pn+a _ pn+a
+h (Ч2 )i,j ((-1/2,j + £-1/2,j) ¡2 -1,j -h (Ч3 )i,j (Hi,j+1/2 + £+1/2 +
+h (44 )i,j (( -1/2+£j -1/2)
pn+a pn+a bi,j -£,j -1
hi
p
i \ x,p и+1/2,
-(Ч2
\j \ x,p !i-1/2,
(Ч3
',П y,p h, j+1/2
,-(Ч4
j \ у, p h, j-1/2
h
h
где а - вес схемы [8].
В случае постоянного тангенциального напряжения на поверхности водоема решение задачи (10) сводится к решению однородного уравнения с граничными условиями второго рода.
Программная реализация. На основе разработанных алгоритмов построен программный комплекс, предназначенный для моделирования и прогнозирования процесса подъема уровня. На рис. 2 цветом показана функция возвышения уровня, рассчитанная на основе уравнения (10).
Рис. 2. Функция возвышения уровня
x
На рис. 3 приведена зависимость установившейся функции возвышения уровня на боковой границе от тангенциального напряжения, полученная экспериментальным путем. Из приведенного графика видно, что значение подъема уровня на боковой стенке зависит от тангенциального напряжения линейно.
Будем предполагать, что поведение функции возвышения уровня на боковой поверхности описывается линейным дифференциальным уравнением первого порядка вида = + Ь. Решением данного уравнения является функция
= C0ekk + C1 . На рис. 4 приведена зависимость функции возвышения уровня от времени, полученная на основе уравнения (10) и основе линейного дифференциального уравнения вида = k£ + Ь .
Рис. 3. Зависимость установившейся функции возвышения уровня на боковой границе от тангенциального напряжения
1-/7 ¡1 2
Рис. 4. Зависимость функции возвышения уровня от времени: 1 - полученная на основе уравнения (10); 2 - полученная на основе линейного дифференциального уравнения
Из приведенного рис. 4 видно, что зависимости функции возвышения уровня достаточно близко проходят друг от друга. Показана возможность прогнозирования поведения функции возвышения уровня на основе модели, базирующейся на дифференциальном уравнении первого порядка. К достоинствам данной математической модели можно отнести ее малую вычислительную трудозатратность. К недостаткам будем относить низкую точность и ограниченную область применимости. Для повышения точности подхода следует использовать модели, основанные на дифференциальных уравнениях более высоких порядков.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Якушев Е.В., Сухинов А.И. Комплексные океанологические исследования Азовского моря в 28-м рейсе научно-исследовательского судна «Акванавт» // Океанология. - 2003. - Т. 43, № 1. - С. 44-53.
2. Сухинов А.И. Прецизионные модели гидродинамики и опыт применения в предсказании и реконструкции чрезвычайных ситуаций в Азовском море // Известия ТРТУ. - 2006.
- № 3 (58). - С. 228-235.
3. Сухинов А.И., Чистяков А.Е., АлексеенкоЕ.В. Численная реализация трехмерной модели гидродинамики для мелководных водоемов на супервычислительной системе // Математическое моделирование. - 2011. - Т. 23, № 3. - С. 3-21.
4. Чистяков А.Е. Трехмерная модель движения водной среды в Азовском море с учетом транспорта солей и тепла // Известия ЮФУ. Технические науки. - 2009. - № 8 (97).
- С. 75-82.
5. Белоцерковский О.М. Турбулентность: новые подходы. - М.: Наука, 2003.
6. Белоцерковский О. М., Гущин В. А., Щенников В. В. Метод расщепления в применении к решению задач динамики вязкой несжимаемой жидкости // Вычислительная математика и математическая физика. - 1975. - 15:1. - С. 197-207.
7. Самарский А.А. Теория разностных схем. - М.: Наука, 1989.
8. Сухинов А.И., Чистяков А.Е., Бондаренко Ю.С. Оценка погрешности решения уравнения диффузии на основе схем с весами // Известия ЮФУ. Технические науки.- 2011. - № 8 (121). - С. 6-13.
9. Сухинов А.И., Тимофеева Е.Ф. Чистяков А.Е. Построение и исследование дискретной математической модели расчета прибрежных волновых процессов // Известия ЮФУ. Технические науки. - 2011. - № 8 (121). - С. 22-32.
10. Сухинов А.И., Чистяков А.Е., Тимофеева Е.Ф., Шишеня А.В. Математическая модель расчета прибрежных волновых процессов // Математическое моделирование. - 2012.
- Т. 24, № 8. - С. 32-44.
11. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. - М.: Наука, 1978.
12. Коновалов А.Н. К теории попеременно-треугольного итерационного метода // Сибирский математический журнал. - 2002. - 43:3. - С. 552-572.
13. Чистяков А.Е. Теоретические оценки ускорения и эффективности параллельной реализации ПТМ скорейшего спуска // Известия ЮФУ. Технические науки. - 2010. - № 6 (107). - С. 237-249.
14. Сухинов А.И., Никитина А.В. Математическое моделирование и экспедиционные исследования качества вод в Азовском море // Известия ЮФУ. Технические науки. - 2011.
- № 8 (121). - С. 62-73.
Статью рекомендовал к опубликованию д.т.н., профессор Я.Е. Ромм.
Чистяков Александр Евгеньевич - Федеральное государственное автономное образовательное учреждение высшего профессионального образования «Южный федеральный университет»; e-mail: [email protected]; 347928, г. Таганрог, пер. Некрасовский, 44; тел.: 88634371606; кафедра высшей математики; доцент.
Кузнецова Инна Юрьевна - e-mail: [email protected]; тел.: 89185954091; кафедра высшей математики; аспирантка.
Chistyakov Alexander Evgenjevich - Federal State-Owned Autonomy Educational Establishment of Higher Vocational Education "Southern Federal University"; e-mail: [email protected]; 44, Nekrasovskiy, Taganrog, 347928, Russia; phone: +78634371606; the department of higher mathematics; associate professor.
Kuznetsova Inna Yurievna - e-mail: [email protected]; phone: +79185954091; the department of higher mathematics; postgraduate student.
