CC BY
40
..
МОДЕЛИРОВАНИЕ ТРЕХМЕРНОЙ ГИДРОДИНАМИКИ МЕЛКОГО МОРЯ В ГИДРОСТАТИЧЕСКОМ ПРИБЛИЖЕНИИ
Безразмерные переменные. Введем характерные масштабы переменных Кр , К , К , К , К , К , Ку , , Кр , Кф , Кц , Кц' , . Физические ПОСТОЯННЬШ
(g и, возможно, ц и ц') примут в выбранных масштабах соответствующие числовые значения. Потребуем от масштабов равенства чисел Струхала 8ИХ 8И у = ЬуЬ-1Ц1 и 8И2 = Ьг17^Ц1 между собой 8И х = 8И у = 8И г = 8И и равенство их единице 8И = 1. Потребуем также выполнения для чисел Эйлера
Fry = LvLl, Frz = LvLl соотношений Eux = Frx Ч Euy = Fry\ Euz = Frz !. Тогда
[1]
вращающейся с постоянной угловой скоростью системе координат примет вид
Кориолиса (при выборе направлений х, у , г на восток, на север и вертикально вверх соответственно), Ф - северная широта, X и X' - коэффициенты турбулентного обмена импульсом.
Если потребовать равенства единице чисел Эйлера и Фруда в уравнениях для горизонтальных компонент вектора скорости Ей х = Ей у = Ей = 1,
Frx = Fry = Fr = І, то Euz = R2Eu = R2, Frz = R_2Fr = R~2 и R = diag{UR}, где
R = Lh/Lz , Lh = Lx = Ly .
Вертикальное поле. Пусть, аналогично [2], обобщенное время т синхронно во всех точках пространства (т'х = т'у = т'г = 0). Пусть координатные линии направлений z и Z совпадают с силовыми линиями (вертикалями) гравитационного поля grad ф= -g = gk (следовательно, £'z = nZ = 0) и остаются теми же самыми (растягиваясь или сжимаясь по вертикали, но не смещаясь по горизонтали) в любой момент времени (следовательно, %t = П = 0). Тогда портрет матрицы Якоби преобразования координат (,n, Z)^(x, y, z) будет следующим:
(І)
(2)
' 0 sin Ф - cos Ф'
- sin Ф 0 0 - матрица ускорения
v cos Ф 0 0 у
Физические компоненты вектора скорости содержит вектор R 1v .
п; СП 0 0 СП
Т'х ъх П'х сх 0 ^>х п'х сх
ТУ ъу пУ С' У 0 У п У С У
ъ; п; С' Ъ; ) ,0 0 0 С' ;
(3)
Аналогичным будет портрет матрицы Якоби преобразования координат (у y, г п С):
(4)
Это дает следующее выражение элементов матрицы (3) через элементы матрицы (4):
К г хт уТ (< 2 ; п; СР -1 (Ґ 0 0 ;;1
4 / х^ Уъ Т'х ъ х п х сх 0 / х^ Уъ ;ъ
4 хП уП ;П ТУ ъ У п'у су 0 хП уП ;П
Л х^ у? Vх'; ^ ; п; СП I0 0 0 ;^у
К а п СП ГУ 4 0 0 - ;Т/( /?) ^
Т'х ъ / / х п х сх 0 уЩ у(М - у',1 ^) -(гъ уП- гП УъУ!^?3 (і;г))
ТУ ъ У п У су 0 - хП/3& 1 х%13 (^г) -(гП- хП 3 ^))
Vх; ъ ; п ; СП 0 0 0
Это, в свою очередь, дает следующие выражения для производных скаляра:
/=/л++лч+/—с;=Г - 4 £,
Ґх - /Л: + + /ПпХ + /і; С'х -
(ДуП - Л'у) Ґ (;ъ уП - ;П у'\)
(&)
/■' - /у і /і /'П' і /'Г = (х^/п Хп/ъ) Д (х^;п хп^
!у - /ТТ У +3ъ ъ У + ЗпП У + /^ У -------------уст)-------^
/;
Ґ; - /ТТ'; + + ./пПг + ДСХ - _Т .
(5)
1- № П г"
Гладкость 2—^ = обеспечивает
4 =((4)т -(Т4\ )•
т-> » ГГ »»»» /Гладкость ^ ^п, У^ = У^ обеспечивает
лМуу&)))(пУ4\-(У4)п + у'\(/\-уП(2!
т- п гг п п п п г-
Гладкость ^ ^ ^, х^ = х^ обеспечивает
Гу = ( у ^))) -( 2\ 4 ^ + (у 2\4 )п - Х( ^
Вертикальность координатных линий направления £ у£ = 0 и гладкость
п п п п г-
У^ = У?\, Уп? = У^п обеспечивают
Гх = (у)) 1 ((ПУ4)) - (ч4)п + ((- уП4)/\) • (9)
(6)
(7)
(8)
г
Вертикальность координатных линий направления Z Х1 = 0 и гладкость
п // // ft f-
х^ ^ ^, Xnz = n обеспечивают
/; = (j(Zz)))-(^ (zZ/)n -((-xn)zj. (io)
Поскольку система уравнений хотя и преобразуется в криволинейные координаты, но векторные величины разлагаются не по локальному, а по глобальному декартову базису [2, 3], то компоненты градиента и дивергенций тензорной диады и смешанного тензора 2-го ранга (вязких напряжений) преобразуются независимо друг от друга. Следовательно, уравнение баланса для каждой компоненты импульса преобразуется как уравнение для скалярной функции (дая каждой из которых две другие компоненты преобразуются как такие же независимые скалярные функции, хотя вместе они составляют единый объект - вектор). Это означает, что каж-
(2) -
ная скалярной функции (5)-(10).
Мелкое море. В мелком море R составляет 104 ...105, то есть R >> 1. В этом случае наличие w в уравнениях для и и v, вносимое слагаемыми R 2i а div(xdv/dx“), 2pR Q R- 'V, не существенно из-за коэффициентов R 2 и R 1.
В то же время наличие и и v с коэффициентами R2 и R1 в уравнении для w ,
вносимое этими же слагаемыми, ведет к численной неустойчивости. Причем сгущение сетки в горизонтальных направлениях (уменьшение R ) проблему не устраняет. Дело в том, что сами процессы таковы, что различие в поле скоростей на 1 м
10 .
Наличие компоненты v в уравнении для и и компоненты и в уравнении
для v, вносимого слагаемым R2iа div(xdv/дха), оставляет возможность «перекачивания» энергии между компонентами и и v (строго диссипативным является не
, ). подсеточных процессов развивается еще один тип неустойчивости (или, по край, ).
В то же время при X = const (2) может быть записано в виде
(pv) + iа div(pvаv) = -R2(grad(p - ( + X')div v)- pg) + + )а div(XR2 gradva)+ 2pR Я R-1v .
В гидростатическом приближении и в приближении возмущенной плотности div v = 0 . Это (вместе с модификацией £1) снимает перечисленные проблемы. Но для преобразования слагаемых вязких сил сжимаемой жидкости в слагаемые вязких сил несжимаемой жидкости в силу уравнения неразрывности несжимаемой жидкости на сеточном уровне требуется перестановочность операторов простран, ( ), -бальном смысле интегральной теоремы Гаусса-Остроградского. Этого формат пространственных производных обеспечить не может даже в двумерном случае [2], т.е., несмотря на заявленный принцип [2], упрощаются сами дифференциальные уравнения (2):
(pv) + )аdiv^v^^-R2(gradp-pg) + )аdiv(XR2gradv^+ 2p£2'v, (11) 0 sin Ф 0^
где £2' = Q
J
.
e
Р = a + jpgds, pZ + pg = 0, (12)
г
где a(t, x,;) - атмосферное давление; e(, x,;) - функция возвышения свободной поверхности, производная p г в уравнениях движения (11) отсутствует. Следовательно, уравнение Пуассона для давления MAC-техники [1] pp = pv + x(iа div(XR2 gradvа - pvаv)+ 2pQ' v),
т2 div
/
R grad(2
a +
V. z
e Y\ I pgds
/ I e
ir^p ’V - z£ pp + ((Z - z )pv)z -xz'z R 2grad (2) a + J pgds
a
z
j
где r = p-(^'?) '(Z?p+(i? - Z)p)z ) + Tdiv((г?) ^Z?Pv + (Z - Z)v)Z )); grad(2) / = /i + f;j + 0k , будет двумерным, но потребуется также уравнение для определения положения свободной поверхности, для получения которого проинтегрируем (1) по вертикали от поверхности дна d(, x,;) до e и используем содержащее e соотношение (12):
j( + div(pv))dz - jpdz
dd
+ div.
(2)
j pvdz
d
- - pd®d - 0 ,
e
- fee- + fdd-;
dІV(2)(axІ + ayj + azk)-(ax )х +(ay )y і
где j/t'dz = j/dz
d \d Jt
e f e
J div a dz = div(2) J a dz d K,d J
ne =-exi - e; j + k, nd = dxi + d'; j - k - внешние (по отношению к столбу жидкости) нормали; rae и rad - интенсивности испарения-осадков (rne > 0 - испарение) и фильтрации (rad > 0 -в грунт) в кинематических условиях:
et + ®e =(e,ne), dt + (vd,nd)=®d . (13)
При p = const получается следующее уравнение для e
Н = Н + т2 div(2)(.tf grad (2)(p-1a + g<5))+x(div(2)(x(veffle + vd ®d )-V )-(e + ®d ))
t;
H = e - d - полная глубина; V = Jp dz .
где н = е - а -
Криволинейная сетка для того и вводилась, чтобы связать границы раздела вязких сред (и в первую очередь - свободную) с координатными поверхностями, т.е. подобно [2], граничные условия и кинематическое условие на свободной поверхности (13) получаются из задания по разные стороны границы (радела вязких сред) различных значений плотности или коэффициента вязкости (и, возможно,
-r
e
1ЗЗ
предельного перехода для одного из них к существенно большему или существенно меньшему значению). Рассуждения велись для системы уравнений вязкой сжимаемой жидкости в том числе и с целью возможности этого (задания различных значений). При этом сохранялась и возможность редукции к составляющим основу приближения возмущенной плотности уравнениям вязкой несжимаемой (постоянная плотность) жидкости [1].
[2], следующим образом скорректировать поток через
подвижную £, -грань
РЧ . 1 , = £(^) . 1 , (ри). 1 . 1 , + £(^) . 1 , (ру). 1 . 1 , +
г +-2, ] +1,к г +-1,] +2,£\Г I, +1 ] +±,к г +2,] +-2,к\к I, +1 ] +2, к
+ £1(^1),] + ±,к( (Р^г + 2, ] + 2,к (р*т) 1+2, ] + 2,к)
то конвективные производные аппроксимируются, например, следующими сеточными представлениями
л))л) = ((р/л)т ) + (й1у(ру/)А/. к ,
'1, ],к ' 11, ], к
^(ру\ + МруУ))аУ ^ = ^(р/А)т ^ ^ к + (div(рv/)А).,] к , (14)
где ^(р/А)т у ] к = ^(Л) .,].,к(р/).,].,к -<А) 1,],Лр/>г,]д) /т , а представление
^^(ру/)А). ] к не отличается от аппроксимации на неподвижной сетке [3].
Однако при попытках разрешения системы сеточных уравнений относительно непосредственно компонент вектора скорости и,],к, V,,] к , Т4-г,]к развивается
сеточная неустойчивость. Решение «р^рушается» на первых же итерациях, хотя на сеточном уровне, подобно [2, 3], выдерживается консервативность, сопряженность аппроксимаций слагаемых дивергенции и компонент градиента, строгая диссипа-тивность и энергетическая нейтральность.
Анализ показывает, что в дивергентных формах производных по т, £,, п (6)-(10) функция / входит только в виде произведения / . Действительно, разрешение системы относительно ^и^У к, (г'^.к, ^, где
(Ук)г]к = 2 (у,],к+1 - У,],к-1 Ь,],к , С последующим Выделением иг,м , V,,],к , wi, ]к приводит к положительному результату. Разница в том, что на новом временном слое определяются не г и /, а /. Затем положение г назначается ( , ), значение на новом временном слое / отрабатывается с учетом назначенного положения г .
Характерные длины ребер сетки по горизонтали составляют 1...10км, по вертикали - 0,1.1,0 м. Характерные значения скорости совпадают с реально измеряемыми. Шаг по времени т приходится выбирать в диапазоне 15.30 с (мо). -ранта. При таком значении т в матрице оператора сеточного уравнения для возвышения уровня е образуется настолько сильное диагональное преобладание, что за одну итерацию Г аусса-Зейделя (ни разу не была зафиксирована потребность в более чем одной итерации) относительная невязка снижалась до нуля (то есть при
расчетах с одинарной точностью снижалась более чем на 7 десятичных порядков).
, .
Кроме этого, были исследованы следующие вопросы. В [2] заявлен принцип: не сразу использовать упрощение в соответствующем уравнении, а вначале получить соответствующую сеточную аппроксимацию, наследующую отмеченное ,
упрощенным в силу этой сеточной аппроксимации, а не аппроксиммировать уже упрощенное дифференциальное уравнение. Тем не менее, кинематическое условие на свободной поверхности должно выражаться через испарение, а не через компоненты поверхностной скорости, т.к. в последнем случае возникает сеточная неустойчивость даже при незначительном изменении скорости ветра (вследствие силь-
).
Придонный пограничный слой настолько тонок (и это подтверждают натур),
( ).
Вертикальная компонента вектора скорости w и ] к может определяться из третьего уравнения количества движения, а может из уравнения неразрывности
но при этом возникает проблема проецирования сеточных функций, заданных на , .
Выражение кинематического условия на свободной поверхности через ис,
( ), производную по времени. Но если сеточное уравнение неразрывности не содержит производной по времени, то аппроксимации производных по времени в уравнениях движения (14) имеет смысл согласовывать только между собой (без согласования с производной по времени в уравнении неразрывности).
В случае расчета потоков через £ -грани ячеек не через компоненты вектора скорости, а из уравнений неразрывности для столба (15), нет необходимости в проецировании сеточных функций, заданных на гранях, на пространство сеточных функций, заданных в узлах. Следовательно, полностью отпадает необходимость расчета вертикальной компоненты вектора скорости из уравнений движения.
1. Флетчер К. Вычислительные методы в динамике жидкостей: В 2-х томах: Пер. с англ.
- М.: Мир, 1991.
2. Васильев В.С. Энергетически нейтральная строго диссипативная аппроксимация двумерной гидродинамической системы на подвижной криволинейной сетке в вертикальном поле // Известия ТРТУ. - Таганрог: Изд-во ТРТУ, 2003, № 5. - С. 184-189.
3. Васильев В.С. Энергетически нейтральная строго диссипативная аппроксимация дву-
// , 2002.
2
(15)
й
БИБЛИОГРАФИЧЕСКИМ СПИСОК
№ 1. - С. 171-177.
