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

Раздел I. Математическое моделирование в аэро- и гидродинамике

УДК 519.8:532

Е.В. Алексеенко, Б.В. Сидоренко, О.В. Колгунова, А.Е. Чистяков

СРАВНИТЕЛЬНЫЙ АНАЛИЗ КЛАССИЧЕСКИХ И НЕКЛАССИЧЕСКИХ

МОДЕЛЕЙ ГИДРОДИНАМИКИ ВОДОЕМОВ С ТУРБУЛЕНТНЫМ

ОБМЕНОМ

В данной работе рассматриваются и применяются на практике два подхода в моделировании гидродинамики водоемов. Один из них основан на численном решении уравнений Навье - Стокса, с помощью метода поправки к давлению. Другой подход основан на решении Lattice Boltzmann уравнения. Проводится сравнительный анализ этих двух подходов на примере водоема «Etang de Berre».

Уравнения Навье - Стокса; метод поправки к давлению; решеточное уравнение Больцмана; решеточный метод Больцмана.

E.V. Alekseenko, B.V. Sidorenko, O.V. Kolgunova, A.E. Chistyakov

THE COMPARATIVE ANALYSIS OF CLASSICAL AND NON-CLASSICAL MODELS OF FLUID DYNAMICS OF WATER BASINS WITH A TURBULENT

EXCHANGE

In this work two approaches in simulation fluid dynamics of basins are considered and applied in practice. One of them is based on numerical solution of the Navier - Stokes equations by Pressure-Correction method. Other approach based on solution of the Lattice Boltzmann equation. The comparative analysis of these two approaches on a basin example «Etang de Berre» is performed.

Navier - Stokes equations; Pressure-Correction method; Lattice Boltzmann Equation; Lattice Boltzmann Method.

Данная работа посвящена сравнению принципиально разных подходов в решении гидродинамических задач на примере мелководного водоема на юге Франции - лагуны «Etang de Berre». Первый из них является классическим и берет свое начало из уравнений Навье-Стокса, а второй более молодой и основан на кинетическом уравнении Больцмана, которое стоит выше в иерархии моделей. Кратко опишем суть каждого из методов, а затем проведем сравнительный анализ на кон.

Исходными уравнениями гидродинамики являются [1, 2]:

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

1 ' ' '

u' + uu'x + vu'y + wu' =-a'x +

б

1 ' ' '

у', + му' + уу' + ъу' =---а' +

і х у г р у \Г" х/х у)у V г)г

1 ' ' '

+ Ж; + ^'у + ^ =~ра' + 0“^ )х +

♦ уравнение неразрывности для несжимаемой жидкости:

/ / /

Мх + Уу + = 0

где V = {м, у, - компоненты вектора скорости, а - гидростатическое давление,

р - плотность, О - угловая скорость вращения земли, 8 - угол между вектором угловой скорости и вертикалью (широта местности), ц,у - горизонтальная и вертикальная составляющая коэффициента турбулентного обмена.

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

Для уравнений (1-3) выполним аппроксимацию по временной координате

u — u

1 ' ' '

+ uux + vuy + wuz = —a'x +

v — v 1 ' ' '

-— + uvx + vvy + wvz =-pa'y +

W w i / ^ /

w—w + uwx + vwy + wwz =--a +

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

Известно, что эффективным численным методом решения задач гидродинамики является МАС-метод. В работе рассматривается вариант данного метода, известный как метод поправки к давлению. Данный метод представляет собой аддитивную схему расщепления по физическим процессам и гарантирует выполне-( ), . систему уравнений (5-7) по физическим процессам. Для этого первое слагаемое каждого уравнения преобразуем следующим образом:

u — u + u — u _, 1 ,

-----------------+ uu „ + vu , , + wu , =---a „ +

x y z x

T p

/ / /

+

v v + v v + uvx + vvy + wvz =—p a'y +

w — w + w — w _ _ _ 1 ,

-----------------+ uwx + vw + ww z =---------az +

T p

+ {^w;x )x +{vwy )y +

Разобьем каждое уравнение системы (8-10) на две части следующим образом: u — u

+ uu'x + vuу + wuf =

- + u-x + w/ + wV/ =

/ f /

+ uW>X + vw/ + wW/ =

т

W — w

т

u — u 1 ,

------=-----ax

т p

-—- 1 f

--------=----------a ,

т p У

^—W = —±a'.

т p

Продифференцируем уравнения (14-16) no переменным x,y,z соответственно:

u — u' 1 „

^ = —al, <17)

т p

— — V 1

-y--------У =---a'yy,

т p

Wz — Wz 1 »

—-------- =------a .

zz

т p

сложим уравнения

мх - м'% - % , й - й' 1 „ 1 „ 1 „

—---- + “---- + —------" =---а хх---а уу---агг

т т т р р р

(4) (20), :

ff.ff.ff р{~f , ~f , ~f\ /ЛП

ахх + ауу + агг = х + Уу + й ).

Уравнения (11-16), (21) и определяют метод поправки к давлению.

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

й ={хі = А,у- = А,г- = кК; -=і-^х,і=1-ыу,к=^= 1х,^уАу = 1у,= К}

В системе уравнений (11-16), (21) вместо частных производных будем использовать их конечноразностные аналоги, полученные при помощи интегро-

.

S

Первая задача представлена уравнением диффузии - конвекции, на основе которого вычисляется поле скорости на промежуточном временном шаге

+ uu0 + vu 0 + wu0 =

x y z x y z

+ uVo + vvо + wVo = )- +

T

v - v

T

w — w

+ ЫМ . + У^0 + ММ 0 = ) +

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

уравнением Пуассона

3* + Яуу + ^ = Т

с

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

Ы-Ы 1

=----------a0,

т р х

=-! а..

Т р у

М М 1

-------=------а0.

Т Р 2

Таким образом, осуществляется переход на следующий временной слой. Условием окончания перехода между временными слоями является - установление значений поля скорости. Уравнения (22-25) решаются попеременно - треугольным методом скорейшего спуска [3]. При построении дискретной модели (22-28) использовались неявные разностные схемы. Главное достоинство данных схем -большой запас устойчивости. Данную дискретную модель удобно также использовать для расчета стационарных течений.

Второй подход основан на непрерывном уравнении Больцмана - интегро-

,

функции распределения / (х, ^, /) в пространстве импульса:

^ + §.Vl f + a-V5 f = Q

д /

+ <•> ух / + а / = О

Интеграл столкновения

О

где о(О) - поперечное сечение столкновения двух частиц, которое преобразовывает скорости из (до столкновения) в {Г&] ( ). -

жение в физическом пространстве обозначено х и скорость в пространстве импульса (или скорость) обозначена ^. / (хДд) й3 хй3 £ представляет вероятность

.

нахождения частицы в объеме d x вокруг x и со скоростью между \ и \ + d,

a - , .

При численном решении уравнения Больцмана используют Lattice Boltzmann метод (LBM), использующий минимальные дискретные кинетические модели для решения задач в механике жидкости и других областях физики [5-8]. Вместо пря-

- , LBM -

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

пределе получаются течения жидкости, представляемые слабо сжимаемыми уравнениями Навье - Стокса.

Эффекты столкновения частиц при решении LBE вообще представляют собой релаксационные модели. Одна из наиболее общих - это одновременная релаксационная модель (SRT), также называемая БГК (Бхатнагар, Гросс и Крук) модель [12].

LB,

обобщенный метод LB уравнения (GLBE), основанным на использовании множе-

(MRT) .

актуально улучшенная форма квазилинейной релаксационной версии уравнения LB с матрицей столкновений [10,11,14], применяется в данной работе к водоему «Etang de Berre». Общий вид MRT моделей в двумерном и трехмерном пространствах можно найти в работах Лэлманда и Луо [13] и д'Юмиерса и др. [16].

Lattice Boltzmann

из эволюционного уравнения для функции распределения ансамбля частиц. То, как они двигаются и сталкиваются на решетке, смотрите [16, 17]:

fjx + eji,, t + S, )-fa

Здесь левая часть уравнения (31) соответствует изменению функции распределения за интервал 5t, так ансамбль частиц перемещается из положения X в его

соседнее положение X + eaSt со скоростью еа вдоль характерного направления а. Мы рассматриваем трехмерное 19-скоростное множество скоростей частиц (D3Q19), показанное на рис. 1, имеющее следующий вид

'„ =

(0,0,0) а = 0 (±Cx, 0,0),(0, ± Cy ,0), (0,0, ±Cz) a = 1,..,6 (±Cx ,± Cy ,0),(±Cx ,0, ±Cz ),(0, ± Cy ,±Cz) a = 7,..,18

Величины скоростей Сх, Су, С2 компонент скорости частицы еа определяются так с = 8 / 8, С = 5 / 8 и С = 8_ / 81, где 8 ,8 ,8 - сеточные шаги

X X I у у X ^ ^ ^ X ^ у ^ 2

по пространству, а 8{ - шаг по времени.

Первое слагаемое в правой части уравнения (31) представляет собой суммарный эффект столкновения частиц в эволюции функции распределения /а. вЬВБ

имеет общий вид матрицы столкновений с множеством времен релаксации, которые соответствуют основным физическим параметрам, таким как плотность, им-

пульс и тензор напряжений. Они в свою очередь представляют собой различные кинетические моменты функции распределения.

03019

Рис. 1. Набор скоростей для обычно используемой трехмерной модели D3Q19

Второй член в правой части уравнения (31) вносит изменения в развитие функции распределения обусловленное полями внешних сил Р, через источник

«а а - . -

щим образом [13,11]:

^ = ( й— м р а -

РС6-

где /ачМ (р, и) - локальное Максвелловское распределение

ам (р,й ) = ш„р\ 1 + Ьф + Й4---(4}> Ж.=

2с 2

1/3 а = 0 1/18 а = 1,..,6

1/36 а = 7,..,18

и С - С2 + С2 + С2/3 - скорость звука модели, р = {р р р } - декартовые

5 у X у 2 ] [ X ’ у’ 2 \

компоненты внешней силы.

Локальная макроскопическая плотность и скорость вычисляются следующим :

Р = ТаУа ,

] = рй = ^а/аеа+ ^ Р8С

Р =РС5 .

Р

Природа кинетического уравнения (31), и в особенности матрицы столкнове ния Лар, станет более наглядна, когда ее определим в терминах множества линей

но независимых моментов £, вместо функции распределения £ =[/" £ ... £0 ]

IV 0 ’ ^ 1 ’ ’^'18J

т.е. через £ = ^£0,/1,...,/18^ , которое получается через матрицу преобразования

Т: £ = ТГ. Элементы матрицы Т даются в работе [16]. Каждый ряд матрицы ортогонален любому другому ряду. Эта матрица такова, что матрица столкновений

Л становится диагональной матрицей Л, после преобразования Л = Т ЛТ 1.

Матрица столкновений в пространстве моментов может быть записана в сле-:

Л = diag

где 5о,51,52,...,518 релаксационные временные коэффициенты соответствую-

.

Когда есть поле внешних сил 8 = [£0,£1,S2.•••.£18], необходимо внести

соответствующий член в пространство моментов 8, то есть 8 = Т8, где

8 = [ $,, ^ ^2,..., Т .

Свойства переноса текущей жидкости, такие как объемная вязкость и кинематическая вязкость могут быть отнесены к соответствующим временам релаксации посредством анализа Чепмена-Энскога вЬББ или анализа устойчивости по Нейману линеаризованной его версии [13] следующим образом:

л

8,

1 1

V * 2

1

У=-

3

л

3 ,Р = 9,11,13,14,15.

Из уравнения (39) 59 = 5П = 513 = 5^ = 515 сохраняет изотропность тен, 51 .

релаксационных параметров не имеет гидродинамического значения, но может повысить вычислительную устойчивость при моделировании задач с большими числами Рейнольдса. Основываясь на линейном анализе устойчивости [13], берем следующие значения для других релаксационных параметров [16]: 5 = 1.19,52 = 510 = 512 = 1.4, 54 = 56 = 58 = 1.2 и 516 = 517 = ^ = 1.98. МОЖНО ,

0 < 5 а < 2.

Вычислительный эксперимент для лагуны «Etang de Вегге»

В данной работе рассматривается решение гидродинамической задачи, а именно расчет параметров течения водной среды в определенные моменты времени в лагуне «Б1а^ de Вегге» с заданными внешними силами, начальными и гра. ,

G = {(х, у, 2): Н (х, у) < 2 < 0} - мелководном водоеме «Ей^ de Вегге» (рис. 2,3), требуется найти компоненты скорости течения й = й (X, у, 2, X), V = у(X, у, 2, X), w = w(X, у, 2, X), давление водной среды Р = Р(X, у, 2, X), где

Н (X, у) известная функция рельефа дна. Внешними силами в данной модели являются: сила ветра, сила трения о дно и сила Кориолиса, вычисляемые по формулам (42), (46) и (49).

Рис. 2. Снимок «Etang de Berre» Рис. 3. Область решения задачи с г.у

В начальный момент времени внутри области G жидкость покоится V = 0, ЧТО соответствует ДЛЯ функции распределения fa = (Оар . Изменение уровня свободной поверхности в модели не учитывается, считается, что она всегда горизон-.

Граничные условия на свободной поверхности определяем таким образом: ди ^ , dv ^ . . др л

w = dz = F-

где Fsx, Fy - компоненты силы трения ветра о морскую поверхность

F = FJ + Fsyj = 0Pa |W|'W ,

где 0 = 8.8-10-3 - коэффициент трения о морскую поверхность,

ра = 1,25 (кг / ж3) - плотность атмосферы, w - скорость ветра (закон Ван-

Дорна). Этим условиям для скорости (41) соответствует граничное условие для функции распределения частиц:

fa

vc.

На дне при z = H(x, y) и на твердых частях границы для скорости применяются условия прилипания U = V = W = 0 и для давление условие др = о,

дп

которым соответствуют граничные условие отскока (Bounce-Back BC) для функ:

fa = fа.

Также для дна могут использоваться такие условия:

Эи „ . ч Эу Эп

^^=-ръх

где ЕЬх, ¥ьу>, Рь2 - компоненты силы трения воды о поверхность дна

ГЬ = ^Ьх1 + ^Ьу} + РЫк = ЗРо М ™ ’

где ^1 = 1.3 • 10-6 - коэффициент трения о дно, Ро - плотность воды, П - скорость течения. Соответствующее граничное условие для функции распределения частиц:

fa <Л /, Н(X. У).t) = fa (X. /, Н(X. У) +dz, t) + дгаР

На жидких частях границ и в источниках (в устьях рек) U

VC

U

CGnst

задаются

известные потоки жидкости. Для данного типа граничных условий для компонент скорости на границе, часто используют ВВ ВС [8,9] для частичной функции распределения /а, которые на входе имеют вид

,

fa= fa + 2aP-

где ea и ea обозначают направления друг напротив друга: ea — —ea .

В оипскгу по с1е 1и1ег Ь оипЛагу

Рис. 4. Расположение границы на входе

. 3. :

1-2) йп = й/(*), (|йп| ~ 1 т /5) -в общем случае функция от времени (направление скорости может еще зависит от того прилив или отлив в данный период времени для 1;

3-4) йп = й0, устья рек 1’Агс и 1а Тои1оиЬге, где задаются известные потоки ;

5-8) ^ = о, где п - нормаль к боковой стенке «Е1а^ de Вегге».

Эп

На движущуюся морскую среду все время действует сила Кориолиса, вычисляемая следующим образом:

F = 2р

CGr Г

Qx U

где Q - вектор угловой скорости вращения земли.

Численные результаты моделирования с помощью Lattice Boltzmann метода представлены на рис. 5, для следующего набора параметров: Ма = 0,002, Re = 10, западный ветер (5 м/с). Аналогичные результаты для данного водоема при таком , ,

рис. 5.

C

S

Для сравнения программой Mars 3D [18] получены результаты, одна из реализаций которых изображена на рис. 5. Для северо-западного ветра (около 5 м/c) с учетом эффекта приливов-отливов в канале Caronte при прочих равных граничных условиях, изображенных на рис. 3.

5гОО°Е 5,0+пЕ 5,0SDE 5.12°Е 5,1 6"Е 5,20°Е 5,24“Е

и , v-> 0.1Q+ LONGITUDE

U-2+V-2

. З.

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

Для параметризации коэффициента вертикального турбулентного обмена была выбрана прямая численная модель Белоцерковского [4]. Она показала себя с наилучшей стороны в вычислениях гидродинамики в мелких водоемах.

K = (cAz )2

dz

+

dz

ч / V /

Ниже на рис. 6 показаны распределения коэффициента вертикального турбулентного обмена неоднородного по глубине, в численном моделировании, в сравнении с данными экспедиционных исследований в лагуне Ейщ de Вегге (Франция).

Рис. 6. Коэффициент вертикального турбулентного обмена на различных уровнях глубины

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

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

Lattice Boltzmann , разностные схемы и обладает меньшими вычислительными затратами. Однако более точным является классический метод поправки к давлению. Совместное использование этих двух отличных по постановке, методов может служить хорошей проверкой на правильность полученных с их помощью результатов.

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

1. Марчук Г.И., Саркисян А.С. Математическое моделирование циркуляции океана. - М.: Наука, 1988. - 304 с.

2. Роуч П. Вычислительная гидродинамика: Пер с англ. - М.: Мир,1980. - 616 с.

3. . . .

журнал, 2002. 43:3. - 552-572 с.

4. . ., . ., . . . .

- М.: Наука, 2002.

5. S. Chen and G. Doolen, Ann. Rev. Fluid Mech. 8, 2527 (1998).

6. S. Succi. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond (Clarendon Press, Oxford, 2001).

7. S. Succi, I. Karlin, and H. Chen, Rev. Mod. Phys. 74, 1203 (2002).

8. D. Yu, R. Mei, L.-S. Luo, and W. Shyу Prog. Aero. Sci. 39, 329 (2003).

9. G. McNamara and G. Zanetti, Phys. Rev. Lett. 61, 2332 (1988).

10. F. Higuera and J. Jimenez, Europhys. Lett. 9, 663 (1989).

11. F. Higuera, S. Succi, and R. Benzi, Europhys. Lett. 9, 345 (1989).

12. P. Bhatnagar, E. Gross, andM. Krook, Phys. Rev. 94, 511 (1954).

13. P. Lallemand and L.-S. Luo, Phys. Rev. E 61, 6546 (2000).

14. R. Benzi, S. Succi, andM. Vergassola, Phys. Rept. 222, 145 (1992).

15. P. Resibois and M. D. Leener, Classical Kinetic Theory of Fluids (John Wiley and Sons, New York, 1977).

16. D. d'Humieres, I. Ginzburg, M. Krafcz/k, P. Lallemand, and L.-S. Luo, Phil. Trans. R. Soc. Lond. A 360, 437 (2002).

17. K.N. Premnath and J. Abraham, J. Comput. Phys. 224, 539 (2007).

18. User manual Mars 3D v.7.34, IFREMER , 2007.

Алексеенко Елена Викторовна

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

E-mail: lena.alekseenko @gmail.com.

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

Тел.: 8(8634)371-606.

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

Alekseenko Elena Viktorovna

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

E-mail: lena.alekseenko @gmail.com.

44, Nekrasovskiy, Taganrog, 347928, Russia.

Phone: 8(8634)371-606.

The Department of Higher Mathematics; post-graduate student.

Сидоренко Борис Владимирович

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

E-mail: sidorenkobv@mail.ru.

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

.: 8(8634)371-606.

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

Sidorenko Boris Vladimirovich

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

E-mail: sidorenkobv@mail.ru.

44, Nekrasovskiy, Taganrog, 347928, Russia.

Phone: 8(8634)371-606.

The Department of Higher Mathematics; post-graduate student and assistant.

Колгунова Олеся Владимировна

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

E-mail: lena.alekseenko @gmail.com.

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

.: 8(8634)371-606.

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

Kolgunova Olesya Vladimirovna

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

E-mail: lena.alekseenko @gmail.com.

44, Nekrasovskiy, Taganrog, 347928, Russia.

Phone: 8(8634)371-606.

The Department of Higher Mathematics; post-graduate student.

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

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

E-mail: cheese_05@mail.ru.

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

.: 8(8634)371-606.

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

Chistyakov Alexander Evgenievich

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

E-mail: cheese_05@mail.ru.

44, Nekrasovskiy, Taganrog, 347928, Russia.

Phone: 8(8634)371-606.

The Department of Higher Mathematics; post-graduate student and assistant.

УДК 519.8:532

Б.В. Сидоренко

МОДЕЛИРОВАНИЕ ГИДРОДИНАМИКИ МЕЛКОВОДНЫХ ВОДОЕМОВ НА ОСНОВЕ SRT LATTICE BOLTZMANN МЕТОДА

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

Решеточное уравнение Больцмана; решеточный метод Больцмана; одиночная релаксационная модель; многовременная релаксационная модель; граничные условия отскока.

B.V. Sidorenko

MODELING OF FLUID DYNAMICS FOR SHALLOW WATER BASINS BASED ON SRT LATTICE BOLTZMANN METHOD

In the work SRT Lattice Boltzmann model in CFD are considered. With their help numerical experiments for real basins are made. At numerical modeling used D3Q19 model which has been modified for not cubic grids, that is with prevalence of one of step in space. SRT models have shown high speed of calculations and simplicity of using in numerical computations.

LBE; LBM; LB-BGK; Lattice Boltzmann; MRT; Bounce Back BC.

В настоящее время все больше внимания в вычислительной гидродинамике , Lattice Boltzmann(LB) ,