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

Сухинов Александр Иванович

Технологический институт федерального государственного автономного образовательного учреждения высшего профессионального образования «Южный федеральный университет» в г. Таганроге.

E-mail: [email protected].

347928, г. Таганрог, пер. Некрасовский, 44.

Тел.: 88634310599; 89281021106.

; . .- . .; .

Чистяков Александр Евгеньевич

E-mail: [email protected].

Тел.: 88634371606.

Кафедра высшей математики; ассистент.

Проценко Елена Анатольевна E-mail: rab55555@ rambler.ru.

.

Sukhinov Alexander Ivanovich

Taganrog Institute of Technology - Federal State-Owned Autonomy Educational Establishment of Higher Vocational Education “Southern Federal University”.

E-mail: [email protected]

44, Nekrasovskiy, Taganrog, 347928, Russia.

Phone: +78634310599.

The head of TIT SFedU; Dr. of Phis.-Math. Sc.; Professor.

Chistyakov Alexander Evgenjevich

E-mail: [email protected].

Phone: +78634371606.

The Department of Higher Mathematics; Assistant.

Protsenko Elena Anatol'evna

E-mail: rab55555@ rambler.ru.

Senior Lecturer.

УДК 519.6

А.В. Шишеня

ТРЕХМЕРНАЯ МОДЕЛЬ ГИДРОДИНАМИКИ И ПРОЦЕССОВ ПЕРЕНОСА

ТЕПЛА И СОЛЕЙ В АКВАТОРИИ АЗОВСКОГО МОРЯ С УЧЕТОМ СГОННО-НАГОННЫХ ЯВЛЕНИЙ

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

Гидродинамика; свободная поверхность; го/-метод; тепломассоперенос.

A.V. Shishenya

THREE-DIMENSIONAL MODEL OF HYDRODYNAMICS AND HEAT AND SALINITY TRANSFER IN THE SEA OF AZOV TAKING INTO ACCOUNT WIND-SURGE PHENOMENA

The work is devoted to the creating of the mathematical model of hydrodynamics and heat and salinity transfer in shallow water basins as well as to the development of the methods for treating with free surface of the water body for wind-surge phenomena modelling. In the article a new method for sampling the continuous model equations is suggested. It allows to take into account fullness of the sampling grid cells. The methods of the fullness of cells calculation are considered as well: based on the function of free surface and using pressure field. Numerical experiments and analysis of results were held for a model region. Simulation results coincide with real physical process.

Hydrodynamics; free surface; vof-method; heat and mass transfer.

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

[4] построена модель динамики токсичных сине-зеленых водорослей, учитывающая флуоресценцию фитопланктона, неоднородность пространственного распределения зоопланктона и биогенных веществ, а в [5] на основе модели гидродинамики проводится подробный анализ причин экологической катастрофы в Азовском море в 2001 году. В то же время ранее авторами уже рассматривалась модель распространения солей и тепла в Азовском море [6]. В настоящей работе продемонстрирована оригинальная техника аппроксимации уравнений непрерывной модели гидродинамики и тепло-массопереноса, позволяющая учитывать частичную заполненность ячеек дискретной сетки, и предложены способы расчета коэф-

.

. ,

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

В непрерывную математическую модель гидродинамики и переноса тепла и солей входят [2]:

♦ уравнения движения (Навье-Стокса):

1 f ^ f

Щ + uu'x + vu'y + wu'z =-p'x + {мК ) x +{uu'y )y + (K ) z + fx ,

p

V, + ыу х + УУ + wv

(1)

р

wt + ЫЩХ + vw у + ww ,

1 ; / \; ;

^ + ("Х Ь + "у Ь + (К Ь + Л ,

здесь ", у - микротурбулентные вязкости жидкости по горизонтали и вертикали соответственно, при этом можно считать, что " является константой;

♦ уравнение неразрывности:

/ / /

р+(р« Ь+(ру Ьу+р К =°;

♦ уравнения перен оса тепла и солей:

Т+<+уг;+м>г:=+ "т; ) + (утг;_) + /т,

$+<+< + <=") у+() + /;

♦ уравнени е состояния: р(Т, 5) = 8°°, 969°62 •1°^4 + 588,194°23 •1°^4 Т + 797, °18644 • 1°3 5

-811,465413-1°-5Т2 -325,31°441-1°5Т5 +131,71°842-1°-652 +

+476,6°°414-1°-7 Т3 + 389,187483-1°7 БТ2 +

+287,97153°-1°-8 5 2Т-611,831499-1°-1° 53.

. - : др

(2)

(3)

(4)

(5)

=V- (' )• %

Он

= 0, - = 1,2...г.

Непроницаемая граница - берега и дно: др дН

ды

= 0, РУМ^ дН

( \ ду = -Т (і^ руЦд-дН

= —у (').

- :

др

ОН

ды

ОН

( \ Оу

= -Т Л Ра^~ дН

= -т

і )•

Ро 8

На выходе - Керченский пролив, озеро Сиваш:

др ОН

= °,

дН

= о,

(6)

(7)

(8)

(9)

здесь у - участки границы, где находятся устья рек; у () - скорости течения соответствующих рек; г - количество рек; Тх (), Ту () - составляющие тангенциальных напряжений на дне и поверхности; й - интенсивность испарений; g -ускорение свободного падения; р0 - плотность пресной воды; р - плотность

воды с учетом температуры и растворенных солей.

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

(v )-

V

V

, ту = PC

(v)-

V

V

(lO)

У ' с \ / У

где рс - плотность среды, относительное движение которой вызывает тангенциальные напряжения, V - относительная скорость движения этой среды. В случае стенок и дна рс = ру, V = —V , а в случае поверхности рс = ра, V = W — V , где

ра - плотность атмосферы, W - скорость ветра или прилегающей среды. Безразмерный коэффициент определяется следующим образом:

C (И )=

о, 0088, 0,002б,

< б, б —,

М

> б, б —.

Дискретизация по времени. Схема расщепления по физическим процессам. При решении задачи гидродинамики (1)-(2) основной сложностью является нахождение поля давления. Компоненты поля скоростей входят во все уравнения (1)-(2), а поле давления содержится в уравнениях (1), поэтому не понятно, какое уравнение следует использовать для нахождения поля давления. Один из способов разрешения этой проблемы - метод поправки к давлению - используется в MAC-методе (Markers And Cells).

Будем искать поле скоростей в два этапа. Для этого введём в рассмотрение

поле скоростей ~ = (, ~, w), которое не обязательно удовлетворяет уравнению неразрывности. Требуемое поле скоростей V = (, V, w) найдем исходя из уравнения неразрывности и используя вспомогательное поля V и поле давления.

Введём следующие обозначения:

n Л П+1 n Л П+1 n Л n+1

u = u , u = u , V = V , V = V , w = w , w = w .

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

u — U ______/ _,

------+ uux + VU,, +

т

~ — V

т

w — w

■ + uv, +

■ + uw, +

т

Шу + = (/ux) x + (iUy )y + (уй') z + fx,

x + VVy + wvz = (//vx) x + /) у + (V) z + fy , x + vwy + ^ ) x + (/' ) У +((WX )z + fz .

(ll)

y z \~ x / л v y

Уравнение для расчета давления имеет вид

+ РУ2+ Р"2) = А' + (РИ Нх ^ НУ +(^wН

у А z / ‘ 4 /XV /у

Переход на следующий временной слой осуществляется по формулам

1 ,

(l2)

u — u

т

V —

т

w — w

т

=-----Px

P

l ,

=-----P У

P

l , = "P P

(l3)

с

с

Пространственная аппроксимация уравнений. Введем систе му координат: направим ось Ох на восток, ось Оу на юг, а ось 07 вертикально вниз. Покроем

расчётную область равномерной прямоугольной сеткой с шагами Нх, Ну, hz по

соответствующим трём пространственным направлениям.

Обозначим внутренние узлы сетки а :

а =|(х;, у ■, zk), = г • \, у. = ] • hy, zk = к • hz, г = 1...и1 -2, ] = 1...и2 -2, к =1...и3 -2)

где ^, П2, п3 - количество узлов по направлениям х , у , z соответственно, а границу области расчета обозначим через у.

Введем для удобства смещенную сетку с узлами £ 1 1 1 = х. ■ к . Таким

г—,1—,к — '•*'

2 2 2

образом, центрами ячеек сетки с узлами |х.. к} являются |£ . к), а |х.. к} - цен-

|£,.к }.

тры

ячеек сетки

/

5-1, у-1,к -1’5 1

V

с узлами Л /

і—, /-1,^-1 2 У V

X

5-1, /-1,к-1’5 . . 1

і-1, і—,к -1 2 У

X

Введём

5-1,/-1,к -1’5

множество

Л

1

і-1, і —1 ,к—

2 У

а заполненное жидкостью подмножество области О{ ■ к обозначим ■ к. Зададим сеточные функции и , V, w, р на сетке а, а кусочно-постоянную функцию q в ячейках к, 1 к} ^. х).

Рис. 1. Ячейки и расчетные узлы сетки Введем верхний индекс для множеств О и й по следующему правилу:

О1А0 = ок и о , Оох0 = о. . к и о , о0’0’! = о. . к и о ,

г, 1 ,к г,1,к^ 1 . , г, 1 ,к г,1 ,к ^ . 1 , г, 1 ,к г,1 ,к ^ . , 1

^ ^ г^, / ,к , ^ ^ г, /Н— ,к , ^ ^ г, / ,к^ ;

2 2 2

Д,/к = щ и ^1’°’°1 = и дод.о >

г, /+-,к 7 гЧ—, /,к

2 2

= Д,/,* и ^1’0к01 = и ,

г,/,к +- гЧ—,/,к

22

б/ = и / = б/ и Д^»,

г, /+-,к г, /,кЧ—

22

б;/і=б;/:°и б1Д" і = и б1 А\ = Д/,к и б^(.

г, і ,к +■

і , / Ч— ,к 2

і ч—, / ,к 2

Введем кусочно-постоянную функцию qi / к заполненности ячейки 5

С

іі , і ,к

ц

ц

І, і,к \

оі ,і,к '

■>1,1,1

^.1.1, 1 гЧ—, /Ч— ,кЧ—

V 2 2 2

Ґ

\

Б

. і . і , і

гЧ—, /Ч— ,кЧ—

V 2 2 2 У

8У, і ,к

кАК

где V; 1 к - объем жидкости в ячейке £ 1 к. Введем также следующие обозначения:

“і,/к= ТГь М"^, ц,',к = ТІТ №'ёа, щ,/к = Т1Г М^,

ПхПуПі д1,1,1 ПхПуПі д1,1,1 "х^у1 Д1,1,1

І ,і ,к І,/,к г,і,к

~і,/к = її"1", ~,/к = ДОИ®, *■>*= тг її™”,

пхпупі Б1,1,1 пхпупі Б1,1,1 пхпупі Б1,1,1

г ,і,к і,/,к і,/,к

',і,/,к =т^т~ Шрё”, Рі/* = т,м =^“ М™®

ККК /

І, і ,к

х У і Б1

І, і ,к

Му1; 01ДД

і, і ,к

^, і,к

1

ІхІуІі Діііі

Л! 5Й®.

Аппроксимируем теперь уравнения диффузии-конвекции. Для краткости положим ф = {и, V, ^) и проинтегрируем уравнения по заполненной части ячейки:

Шн-?+и?!+^ 1^= Щ(/х)х+(/) +У)z + Аа

,ііп 7

‘2І,і к

-.1,1,1 ‘2І,і к

Затем воспользуемся линейностью тройного интеграла:

Ш ~—~^ю+ |Ц и$а+ Щ vфydю+ |Ц ц?ф^а =

□ 1ДД ^ О1,!,1 О1,!,1 О1,!,1

°г,],к Ог,] ,к °г,] ,к Ог,] ,к

/ ' /

= (/х )хС Ш(<Ру ШУ )zdа+ 111

І, / ,к

І, / ,к

І, / ,к

І, / ,к

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

ф-ф^ Фі, і ,к - Фі, і ,к , , ,

-----— q. і. і, і КМ.

Т г--2,/--2,к--2

Рассмотрим аппроксимацию конвективного слагаемого:

щ>'$ю= и0хdm + и0хdm = q 1 1 и0хdm + q 1 1 и0хdm.

О-дд О»дд О»дд 1,1 2’ 2 О»д,д 1 2’ 2 О»дд

Для каждого из полученных интегралов воспользуемся теоремой о среднем:

q 1 1 [[[ и¥Аа+q 1 1 [[[ и¥Аа=

—,к— 1 -1,1—,к—

г, і—,<-2 2

Б-

;,і-2- 2

Б"

— q.. і, і “ і., Ш ф^^®+q... і, і “ .., Ш

г, / —,к — гЧ—, /,к г-1, / —,к — г—, /,к

2 2 2

Б"'

-1, / — ,к — г—, / ,к 2 2 2

Б-

Теперь вычислим приближенно эти интегралы, используя формулу Ньютона-Лейбница по переменной х и пользуясь теоремой о среднем для интегралов по переменным у и z :

х 1 У. 1 г 1

і+— /+— к+— 2 2 2

і. і. “ і. [[[ Фхії® = q . . “ . [ [ [ фхйхйуйг —

—,к — ІЧ—, і к Ш і, і — ,к — ІЧ—, і к ** *

і, і—,к — ІЧ—, і,к

2 2 2

Б"

і,і —,к — ІЧ—, і,к

2 2 2 х{ у і г і

і— *— 2 2

— q і і “ і

і,і—,*— і^-,і ,*

2’ 2 2 V 2

Ф. 1 - Ф, і ,к

г+- ] ,к

КА,

у. і 1 і

x /+— *+—

xi ^ 2 2

Б"

1, / — ,к — г—, / ,к

2 2 2 х і у і г і

г— і— к—

2 2 2

Ґ

— <І. . . 1 1 “. 1 .

г -1,/—,к— г—, / ,к 2 2 2

Фі, і ,к -Ф 1

г—, / ,к 2

КК.

Имеем следующую аппроксимацию рассматриваемого слагаемого:

(ГГ _________ Ф-1+1,1 ,к 1 ,к 1 , , Ф-1,1 к ],к , ,

III ифхаа= q 1 1 и 1 ---- --— hyhz + q 1 1 и 1 —-----------— Нп,.

О,»

2. . 1 , 1^. 1 . ,

і, /—,к— і+—, /,к

2 2 2

-і,/ —к— І—* 2 2 2

Аналогично производится аппроксимация остальных конвективных слагае-

мых:

Ж -/7____ ФІ,-' +1,к ФІ,Ф ,к 7 7 ФІ,/,к ФІ,/-1,к 7 7

Ф.б® — q і і V і —1—------- КхКг + q і іу і —1--------------—КК

і—, і,к— і, і +—,к ..і;-..

2 2 2

і—, і-1,к— і, і —,к 2 2 2

ФІ,-,к+1 ФІ,-,к , , ФІ,-,к ФІ,-,к-і , ,

—------------— К К + q . . ^ . ——-----------------—— К К

-~к і,-,*+± 2 у х і-і,--і,*-і і,ік-і 2 у у

Ш™¥^0)—q. і. і іФ

*>*>*> і—, і—,к і, і,к-—

Рассмотрим теперь аппроксимацию диффузионного слагаемого:

/ / /

Ш (/&)^а= ш (мФх)^а+ /л (м?х)^а=

/ /

—qii-ік-і Ш (ц^х) х d® + q . .

2 2

і-1, - — ,к 2 2

Б",

Теперь вычислим приближенно эти интегралы, используя формулу Ньютона-Лейбница по переменной х и пользуясь теоремой о среднем для интегралов по переменным у и z :

qi,j -! к-1 III (/х) х а+ qi -1,,-1к-1 Ш (м) х =

Б",

х 1 у 1 г 1

-+—

2 2 2

% ■ і* І I I I Ц) х dxdydz + 9 і і І І I (<Ф) х dxdydz

2 2 4 2

ц 1 .Л 1 -Ц-Ф--

г +-, - ,к г +-, - ,к ^ ^

2 2 2 /

КуК +

-1,--і'k-і 2 2 4

Ц, - Ф - --Ц 1 .(Ф,1 .

- - г--, - -

2 2 У

КуК .

Для аппроксимации оставшихся производных используем центральные ко:

^ . 1 , іЦі.Ді., КуКг + q. . . 1, іЦ 1 .Ух 1 КК + . . 1, 1 -^ . 1 , 1 ІЦі.І.кФі.І.кКК =

1 Г--&•-- 1-І-Г &• 1-І-Г &• У 1—1 Г-&•---1----Г &• 1--- ! Ъ- ■' І 1—1 Г--&•-- 1 г--&•-- І •> •> у

9 . і іц і .

г, і—,- — г-—, і,к 2 2 2

ФІ+1,і,к ФІ,і,к , , , ,, ФІ,/,* Фг-1,/,*

—1-----------— КК + ^ . к іЦ і —1----------------------------------------—

К ‘Л

Ч ^ :-і*-і- 9і:-1*-і \ЦІ,і,кФх,і,кКК . ,- 2, 2 ,- 2, 2 ,

Имеем следующую аппроксимацию рассматриваемого слагаемого:

/

Ш(цф)xd®—q . і, іЦ і.

х 1, 1—,к— і +—, І,

Фі +1,-,*_____________________ФІ2І2^и и + ^ Ц Фі, і,к Фі-1, і,к К К

-КК + q

Пі,

1 ■ 1 !• іЦ 1 ■

-1, і—,-— г—, і,к

2 2 2

q і.-q і. і-1, і-1 ,*-- і, і--,*--

V 2 2 2 2 У

Ц . Ф. ., К К .

” г,-,*' х,-,* у г

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

Ш (ц^у )у6ю — 4 - ■ Фг'-'- кхК + 9. - . . .-м.. - * Фг'-' * Ф—1' * КА +

4 ' у г—,-—— г,--—,к К г—, і-і,к — г,:—, к К

--1,*-— і,-—к К

2 2 2 Пу

Ці.-Ж-АК.

Ш (тМ®—9і-і і-1

Ф,:, *+і-Ф

■'--і' * І, -, к +1 К

2 2 2 Пг

і ,- +1 П-, к _ „ Ф-,* Ф1,— -1 к к +

"" КуКх + 1. 1 Г

-1,--1 ,*-1 І, - „к-- К

2 у 2 ^ 2 "г

п-Ж-АКх.

Найдем аппроксимацию уравнения Пуассона для давления.Для этого проинтегрируем его по заполненной части ячейки:

Ш7 ■ ( + р}+р"72 ) а® = (+() х+() у+р) г) 6?®

Пі,!,! Й1"

Воспользуемся свойством линейности интеграла:

^ Р/'Х2 dа + 7• Ш Р/ Ш Р*' ^ =

Он О1' ЙЧ'1

/ / /

= Ш ^С + Ш (^“ МС + Ш (^ )ydа+ 111 (^

а1" йч1 а1" а1"

Рассмотрим аппроксимацию каждого слагаемого. Ранее были рассмотрены аппроксимации интегралов, аналогичных первым четырём:

Ш Рх-2

рі+1,і,- рі,-■,- 7 7 .

---------:----------КА + q... ііЦ і

г'-'k i-1'і'k КД +

2 2 2 '

. . і, і-ч . і, і і -1, і —,-— і, і ———

V 2 2 2 2 У

111 Р/^ — .. і . іЦ. і -і-1*--------------------— Нх\ + ^ - іЦ -

у і+1, /+—, к-— і+—і к и

2 2 2

-і, -- Рі/-і-

КхКг +

Р.^'^'kКxКг ,

г, ;—,к— гН—, ;,к

2 2 2

'кЛ+4 . к 1Л1

Рш~Ри,к -1

-1, ;—,к— г—, /,к 2 2 2

кД +

4. , . 1. 1 -4. . 1 . 1

г-1,.—, к— 1,1—, к —

V 2 2 2 2 У

р' . ,к к ,

Г х.,к у х’

г 2 ’3 2 ’к 2

Аппроксимация правой части уравнения для расчета давления имеет вид

ГГГ/ -) , — )г+1,3 ,к -( — );,з,к,, , (Рй ),3, к -( — ),■- 1,з,к,,

.Ш (-и) х d® = 4 . 1--------—--------------------— кукг + 4 _ 1 ^ 1-3----------к.к

^1-1-1 “,' ,

г, 3—,к— 2 2

- 1,3-—,к--2 2

/

Л!(—4 1

( —Ч,3+1,к -( —4

О»

—, 3,к— 2 2

^-1-кхкг + 4, 1 3 1

'- ^, 3 -1,к-т

(—Ч,а -(—Ч,М,

— кк,

ГГГ ( — )<,3-,к+1-(-ЩЦ,, , (-Щ )цл-1, ,

_)_)] (-щ) г «;« = 4 --------------------кк + 4 --------------------кк.

□1.1,1 2, 2, 2 2,

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

Расчет заполненности ячеек. При расчете заполненности будем считать, это колебания плотности достаточно малы и слабо влияют на форму поверхности водоема. Для учета сгонно-нагонных явлений можно считать, что свободная поверхность водоема задана уравнением г = £(х, у). Зная функцию £(х, у), можно найти объем жидкости в ячейках по формуле

Г У3+1 Х;+1 (к(х,у)гм) ^

4г,3,к =тах II | dzdxdy, 0 . (14)

V У3 х1 тах(#(х,у)г) у

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

dz & . I & . I &

ще=—. = & + и\& + Н1 &.

dt

а

В то же время щ р = щ£Н-------------, тогда

# ? -о

£' + и 1^ + Н Л’у -Щ; =

а

-0

Аналогично на дне:

dh

dt

к I-к х I-к у

(15)

(16)

Проинтегрируем уравнение неразрывности для постоянной плотности по переменной z на интервале ((x, у); h(x, у)):

h , k ,

Ju xdz + Jv ydz + w\h - w\#=0. (17)

, :

J f' (x, у )dx = dy J f (x, у )dx - f ((у), у )-p'(y) + f ((у), у ))'(y). (18)

v(y) у р(у)

Обозначим

h h U = Judz, V = Jvdz, (19)

(17)

(h + £)t + U x + V у =< - u |# • - V| # W # W + U| К + V| - h к.

(20)

С учетом соотношений (15), (16) имеем:

J ' ' Ю

• +U x +V у =—. (21)

у Ро

Проинтегрируем теперь выражения (13) с учетом того, что plh = pgн+• р\4 =0:

h '

U - U = -— Jp'xdz = -*Hg (h + 4)x , p •

—h

V- V =-------Jp\dz = -THg ( + •) у.

Р •

(21), :

• -—g(h(h+i)'x) -—g(h(h+ •) у) =РЮ+U' + Vy. (22)

, -

, . приповерхностной ячейке возвышение уровня кусочно-линейной функцией и найдем заполненность по формуле (14), тогда после интегрирования получим:

9,j ,k = КК •j + •" ■ +4•j+■+•+■. (23)

,

:

1. Численно решаются уравнения диффузии-конвекции (11) и находятся промежуточные значения скоростей.

2. Решается уравнение (12) для нахождения поля давления.

3. Производится коррекция поля скорости согласно соотношениям (13).

4. -ности на новом временном слое.

5. Вычисляются интегралы скоростей на промежуточном временном слое (19) и из уравнения (22) находится функция возвышения уровня на новом слое.

6. Исходя из функции возвышения уровня по формуле (23) рассчитываются

.

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

, , ячейка считается заполненной полностью. В общем случае используется следующая формула для расчета заполненности ячейки:

0.1, ],к

р + . + р + + р + + р + . + р + . + р + + р + . + р +

і— 1, ] ,к і, ] —1,к і, ] ,к — 1 і— 1, ]—\,к і— 1, ] ,к — 1 і, ]—\,к — 1 і—1,І ,к — 1 і—1, І _1,к —1

8gph

у

ри,к = р,і,к -л{Рі,і,к ^ Л - Функция Хэвисайда.

где Рик = Р,;,к -ЛЩ,;,к ), V

Численные эксперименты. Проверка корректности построенной модели производилась на распространенной модельной задаче. Область расчета представляет собой параллелепипед, вдоль одной из границ которого течет жидкость. Для расчета свободной поверхности был применен описанный выше эмпирический . . 2,3 -

ным физическим процессом.

Рис. 2. Трехмерная визуализация функции возвышения уровня

Рис. 3. Решение модельной задачи. Поле скорости в приповерхностном слое и функция возвышение уровня (вид сверху)

Заключение. В работе построена непрерывная модель гидродинамики и переноса тепла и солей в мелководных водоемах. Для учета формы свободной поверхности вводится функция заполненности ячеек сетки. Для расчета заполненности предложены две модели: на основе функции свободной поверхности и на основе поля давления. Метод на основе функции свободной поверхности требует, чтобы область расчета была простой по оси аппликат, в то время как эвристический метод не налагает никаких ограничений. Рассматривается оригинальный метод аппроксимации непрерывных уравнений, позволяющий учесть частичную заполненность ячеек расчетной сетки. Такой подход имеет ряд преимуществ по сравнению с 7 -сетками [7]. Во-первых, предлагаемый в статье метод позволяет проводить аппроксимацию на структурированных сетках, для которых расчетные узлы расположены в центрах контрольного объема, что позволяет получить более точную аппроксимацию. Во-вторых, данный подход более прост в реализации по сравнению с 7 -сетками при тех же условиях применимости, что позволяет строить более сложные модели. Проверка адекватности построенной модели проводилась на модельной задаче. Результаты численного моделирования являются фи-зичными и соответствуют ожидаемым.

Автор выражает благодарность АЛ. Сухинову и А.Е. Чистякову за полезное обсуждение материалов статьи.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Сух иное А.И., Чистяков А.Е. Двумерная модель турбулентного движения водной среды в Миусском лимане. Математическое моделирование и информационные технологии / Юж.-Рос. гос. техн. ун-т (НПИ). г. Новочеркасск: Ред. журн. «Изв. вузов. Электромеханика», 2007. (Приложение к журналу). - С. 43-47.

2. Чистяков А.Е. Трехмерная модель движения водной среды в Азовском море с учетом транспорта солей и тепла // Известия ЮФУ. Технические науки. - 2009. - № 8 (97). - С. 75-82.

3. Сух иное А.И., Чистяков А.Е., Алексеенко ЕМ. Численная реализация трехмерной модели

// .

моделирование. - 2011. - № 23. - C. 3-21.

4. . . -

ском заливе // Известия ЮФУ. Технические науки. - 2010. - № 6 (107). - С. 113-117.

5. Sukhinov A.I. Reconstruction of 2001 Ecological Disaster in the Azov Sea on the Basis of Precise Hydrophysics Models.// Parallel Computational Fluid Dynamics, Mutidisciplinary Applications, Prcoceedings of Parallel CFD 2004 Conference, Las Palmas de Gran Canaria, Spain, ELSEVIER, Amsterdam-Berlin-London-New York-Tokyo, 2005. - P. 231-238.

6. .. ..

мелководных водоёмах на примере Азовского моря. КРЭС, 2008. - C. 276-277.

7. Ezer T., Mellor G.L. Sensitivity studies with the North Atlantic sigma coordinate Princeton Ocean Model. Dynamics of Atmospheres and Oceans. - 2000. - Vol. 32. - C. 155-208.

Статью рекомендовал к опубликованию к.т.н., доцент В.Е. Мольдерф.

Шишеня Александр Владимирович

Технологический институт федерального государственного автономного образовательного учреждения высшего профессионального образования «Южный федеральный университет» в г. Таганроге.

E-mail: [email protected].

347928, г. Таганрог, пер. Некрасовский, 44.

Тел.: +7928322282; +79081761837.

Кафедра высшей математики; аспирант.

Shishenya Alexander Vladimirovich

Taganrog Institute of Technology - Federal State-Owned Autonomy Educational Establishment of Higher Vocational Education “Southern Federal University”.

E-mail: [email protected].

44, Nekrasovskiy, Taganrog, 347928, Russia.

Phone: +7928322282; +79081761837.

The Department of Higher Mathematics; Postgraduate Student.

УДК 519.6

. . , . .

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ДЕФОРМАЦИИ СТЕРЖНЯ

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

Уравнение балки Эйлера; ряд Фурье; метод баланса; численный эксперимент.

A.E. Chistyakov, E.A. Kostyrko

MATHEMATICAL MODELLING OF DEFORMATION OF THE BAR

This work deals with the development and study of algorithms for solving problems of elasticity theory that is an actual problem of mathematical physics. Euler's theory about beam is used for the description of elastic properties of a beam. The problem dared in two ways: the decomposition method on Fourier series had been received the analytical decision of a problem; the numerical decision is received with the help finite-difference approximations. Analytical research of the proposed discrete model was completed in this work. Quantitative coincidence of results of numerical and analytical calculations of deformation of a beam is received in case of stationary pressure on time (the distributed loading).

The equation of a beam of Euler; Fourier series; the balance method; a numerical experiment.

. .

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

.

( ).

.

нагрузки стержня решается в линейной постановке. Колебание стержня описывается уравнением балки Эйлера [1]:

, Э 2W Э 4W

Ph^~r + EI^t = P(x), (1)

Эг Эх