CC BY
13
УДК 551.465
С.Г. Демышев
Исследование чувствительности параметризации Меллора — Ямады к выбору конечно-разностных аналогов в численной трехмерной модели оперативного прогноза течений в Черном море
В рамках подхода Меллора - Ямады реализована численная схема расчета коэффициентов турбулентной вязкости и диффузии в z-системе координат для трехмерной модели оперативного прогноза течений в Черном море. Исследованы дискретные аналоги уравнений для кинетической энергии турбулентности и макромасштаба турбулентности. Показана их высокая чувствительность к выбору конечно-разностных аппроксимаций. На основе сопоставления результатов прогностических экспериментов с данными наблюдений выбрана лучшая аппроксимация слагаемого, описывающего генерацию энергии турбулентности.
Введение. Правильное описание процессов, происходящих в верхнем перемешанном слое Черного моря, принципиально важно для адекватного воспроизведения морской термодинамики и, следовательно, для прогноза его состояния. Формирование и эволюция верхнего слоя в модели динамики [1, 2] ранее описывались на основе аппроксимации Филандера - Пакановского [3].
В серии численных расчетов (например, в работе [2]) было показано, что в случае гладкой структуры атмосферных полей использование этой аппроксимации оправданно. В то же время при резких изменениях атмосферной ситуации приближение Филандера - Пакановского приводит к неадекватным результатам. В первую очередь это связано с тем, что в аппроксимации Фи-ландера - Пакановского в отличие от параметризации Меллора - Ямады [4] при расчете коэффициентов турбулентности влияние атмосферного воздействия учитывается опосредованно, через число Ричардсона. Для реализации численной модели оперативного прогноза течений в море, когда необходимо воспроизводить влияние резких изменений атмосферной ситуации, такой недостаток должен быть преодолен.
Параметризация Меллора - Ямады используется в численной модели динамики океана, разработанной в Принстонском университете (POM) [5]. Эта модель выписана в ст-системе координат и применяется для решения задач диагноза и прогноза состояния морской среды. В отличие от POM в данной работе на основе подхода Меллора - Ямады разработана и реализована численная схема расчета коэффициентов турбулентной вязкости и диффузии для трехмерной модели оперативного прогноза течений в Черном море в z-системе координат. Рассматривался одномерный по пространству вариант параметризации Меллора - Ямады, при котором в уравнениях для кинетической энергии турбулентности и макромасштаба турбулентности адвективные и диффузионные слагаемые не учитывались. В качестве первого шага такое приближение оправданно, так как основные проблемы обусловлены слагаемыми, описывающими генерацию и диссипацию турбулентной энергии. Про© С.Г. Демышев, 2010
ISSN 0233-7584. Мор. гидрофиз. журн., 2010, № 3
29
веден анализ конечно-разностных аналогов уравнений для кинетической энергии турбулентности и макромасштаба турбулентности и предложены различные аппроксимации для слагаемых, обеспечивающих генерацию турбулентной энергии.
1. Постановка задачи. В приближении Буссинеска, гидростатики и несжимаемости морской воды рассматривалась следующая система уравнений:
и, - (X + f)v + wuz = -ggx - — (Р'+Е)х + (ууиг)г - ун У4и , (1)
Р0
V + (X + /)и + wvz = -gдy - — (Р'+Е)у + (кЧ)z - Ун, (2)
Р0
их + Vy + Wz = 0, (3)
z
Р = gРoV + g I Рт = gРoV + P', (4)
0
Т + (иТ)х + (VI)у + (wT)z = -к У4Т + (кТг)z, (5)
5,+(и5)х + (^)у + (wS)z = -к1 У4^ + (кХ )z , (6)
р = р0 + аТт + а^ + аТт2 + а5ТБТ . (7)
Обозначения в уравнениях (1) - (7) - общепринятые [1]. Краевые условия имели следующий вид.
На поверхности при г = 0
ууи2 = -тх, vvvг =-^у, кТ2 = дТ, кХ = (еу-Рг)^ + /(5с1 -5о), (8)
где О - потоки тепла, Еу - испарение, Рг- осадки, 50 и 5с1- рассчитываемая и климатическая соленость, / - заданный релаксационный параметр.
На дне при 2 = Н(х, у)
и = V = w = 0, Тг = 5 = 0. (9)
На твердых боковых стенках: для меридиональных участков границы
и = У 2и = ^х = У\ = 0, Тх = (У2Т)х = ^ = (У25)х = 0, (10)
для зональных участков границы
V = У^ = иу = У 2иу = 0, Ту = (У2Т)у = £у = (У25)у = 0. (11)
На части границы, где вода втекает, используются следующие условия: для меридиональных участков
и = ир, У2и = vx =У\ = 0, Т = Тр, 5 = 5р, (У2Т)х = (У25)х = 0, (12) 30 0233-7584. Мор. гидрофиз. журн., 2010, № 3
для зональных участков
V = ур, У2У = иу =у2иу = 0, Т = Тр, 5 = 5Р, (У2Т)у = (V25)у = 0. (13)
Для верхнебосфорского течения и Керченского пролива, когда течение направлено из Черного моря в Азовское:
V = Vх, V^ = иу ^2иу = 0, Ту = 0,= 0, ^2Т)у = (V2= 0. (14)
При 1 = 10 потребуем выполнения следующих условий: и = и°(х,у,г), V = v0(x,у,г), д = д0(х,у), Т = Т0(х,у,г), 5 = 50(х,у,г), (15)
где и0,V0,д°,Т0,50- известные поля на начальный момент времени счета 10.
Полагаем, что У - неизвестные функции. В соответствии с теорией Меллора - Ямады [4] для их определения необходимо знать кинетическую энергию турбулентности (е2/2) и макромасштаб турбулентности (/). Соответствующие соотношения для расчета коэффициентов имеют вид
У = 1е5н, к = 1е5м, (16)
где 5н, 5м - функции устойчивости, которые определяются из эмпирических соотношений [4]:
5м = Л1
( 6Л ^ 1 _ ^А - 3С
В1
+ 9(2 А + Л^н
(1 - 9 А Л2Он )-1,
(17)
5Н = Л2
1 -
6 А
В
[1 - 3 Л2вн (6 Л + В2)].
1 /
N1 е
Здесь Он =-\ — | =N =
,11 А. дР
е2 р0 дг
Р0 дг
(18)
- частота Вяйсяля -
0
Брента. Эмпирические константы Ль Л2, £ь В2, Су определялись из эксперимента [4].
Для нахождения е2 /2 и I решаются уравнения:
де!=_д_
д1 дг
(
т
де
2
д(е2!) =д_ д1 дг
т
дг
V д(е21)
+ 2п
\ J
2;
ди | ( дv
Ъг I +1 дг
+22к к др-
Р0 " '
дг В1 '
дг
+1ЕУ
дг
+
+ кдр-н. р0 дг В
(19)
(20)
где н - эмпирическая функция и Е1, Е3 - эмпирические константы. В общем случае [4] в уравнениях (19), (20) необходимо учитывать адвективные слагаемые и турбулентную диффузию по горизонтали. На данном этапе мы ограничились рассмотрением упрощенного варианта (19), (20), так как основные вычислительные проблемы связаны именно с правыми частями этих уравнений.
ISSN 0233-7584. Мор. гидрофиз. журн., 2010, № 3
31
2
2
2
2
Эта система уравнений должна быть дополнена соотношением для расче-
та ¡1 , которое имеет вид
(21)
т = 1е5е,
где эмпирическая константа 5е =0,2.
Систему уравнений (16) - (21) необходимо дополнить краевыми и начальными условиями:
е2 = 512/3 \(тх +ту У Го2 j1'2, e2l = 0 е2 = 0, е2! = 0 е = е0, ! = !0
при z = 0, (22)
при z = H (x, y), (23)
при t = t .
(24)
Таким образом, решается система уравнений (1) - (7), (16) - (21) с краевыми (8) - (14), (22), (23) и начальными (15), (24) условиями.
2. Конечно-разностная схема параметризации Меллора - Ямады.
Конечно-разностная схема модели (1) - (15) подробно описана в работах [1, 6]. Поэтому остановимся на особенностях разностных аналогов уравнений
(19), (20).
В правых частях уравнений (19), (20) содержатся нелинейные слагаемые. В силу соотношения (21) наряду с последним членом в уравнении (19) нелинейным является первое слагаемое в правой его части. С учетом (16) в уравнении (20) все четыре члена в правой части являются нелинейными. Причем последнее слагаемое в (19) имеет особенность в нуле, второе и третье в (20) обладают квадратичной нелинейностью. Именно поэтому использование явных схем для аппроксимации по времени этих уравнений, как правило, приводит к неустойчивости решения.
По аналогии с работой [5] распишем разностные аналоги уравнений (19),
(20) с учетом неявного представления по времени диффузионного члена. Имеем
^ „ е7;'/2 = °z [(¡Г+1/2) " "Ч 0^2)]+ 2(^+1/2)
2т
о n-1 °z I uk+1/2
+
(25)
+
°z I vk+1/2
2g Po
2p n+1 (pn-1 )3/2 k+1/2 V k+1 / 2 /
i ^S i~V \n-1 S l r,n-1 \ ^k+1/2\^k+1/2
+-(kk+1 / 2 ) °z (Pk+1 / 2 )
n -1 k+1/2
Lk+1 /2 Lk+1 /2 =°\(nV \n 1 ° i ДП+1 \] + /n-1 77 (~)V \n-1 -(- z L(mk+1/2) °z (Lk+1/2^ + !k+1/2E1 nk+1 /2 ) <
n—1
°z I uk-+1/2
+
+
°z I vn+1/2
+
jn-1 ^ ДП+1 n-1 \3/2
Lk+1/2E3g / ~V \n-1 с / .n-1 ч Lk+1/2 (ek+1/2)
(26)
Po
(kk+1/2)n °z (Pn+1/2) -
n-1 k+1/2
H
k+1/ 2'
где n - временной уровень и т - шаг по времени.
32 ISSN 0233-7584. Мор. гидрофиз. журн., 2010, № 3
2
x
n-1
2
2
x
2
Разностные операторы имеют вид
d f — fk +1/2 — fk—1/2 d f uzJk — ■> uzJk+1/2
fk +1/2 — fk- 1/2
zk +1/2 — zk - 1/2
fk+1/2 + fk -1/2
2
_ fk+1 fk Zk+1 — Zk fk+1 + fk
х _ J к+1/2 ' J к-1/2 х _
Jk 2 ' ^к+1/2 2 '
Здесь введены следующие обозначения:
К"
_ (е2)" К" _ (|е2)" I" _ +1/2 (27)
ек+1/2 -Vе ^к+1/2> к+1/2 — ^к+1/2> к+1/2 _ „" •
ек+1/2
Коэффициенты У, х^и т могут аппроксимироваться следующим образом:
(~У \" _ 7" ^" / Г» \ " \" _ 7" ^" ^ Ч*
(пк+1/2) _ 1к+1/2ек+1/2 (8И )к+112, (Кк+1/2) _ 1к+1/2ек+1/2(8М )к+112,
(28а)
(Ж) _ I" е" ;
-z-z
w)" — in en (SB)n, (k )n — k en (SM)n, ,28_.
(28б)
-z -z
) — lnk en St.
В известной схеме POM [5] расчет этих коэффициентов проводится для одних и тех же точек (i, j , k+1/2) (формула (28а)). В настоящей работе используется аппроксимация (28б), которая согласуется с разностной формулировкой уравнений (1), (2), (5), (6) [6].
В уравнениях (24), (25) преобразовано последнее слагаемое. Цель такого преобразования заключается в следующем. Нетрудно видеть, что эти разностные уравнения сводятся к уравнениям прогонки. Условием их разрешимости является свойство диагонального преобладания. Преобразование последних членов этих уравнений приводит к усилению этого свойства и, следовательно, к повышению устойчивости решения конечно-разностной задачи.
Для анализа численных расчетов необходимо указать на особенность аппроксимаций коэффициентов турбулентности в уравнениях (1), (2), (5), (6). В соответствии с распределением переменных на сетке С [1] компоненты горизонтальной скорости u и v рассчитываются для разных относительно друг друга узлов. Поэтому, строго говоря, коэффициенты вертикального турбулентного обмена в конечно-разностных аналогах (1), (2) должны быть определены для различных точек сеточной области. В уравнении (1) nV аппроксимируется в точках (/+1/2, j, k+1/2), в уравнении (2) - в точках (i, j+1/2, k+1/2). В свою очередь коэффициент вертикальной турбулентной диффузии определяется в точках (i, j, k+1/2). Поэтому возникает неоднозначность в аппроксимации второго слагаемого в правых частях уравнений (24), (25).
3. Численные эксперименты. Предварительный анализ показал, что результаты расчетов чрезвычайно чувствительны к некоторым параметрам уравнений Меллора - Ямады. Поэтому были проведены специализированные численные эксперименты. В качестве иллюстрации чувствительности турбулентной модели Меллора - Ямады к выбору аппроксимаций приведем в со-
ISSN 0233-7584. Мор. гидрофиз. журн., 2010, № 3 33
поставлении результаты двух расчетов. Их отличие состояло в следующем. В I эксперименте второе слагаемое аппроксимировалось в соответствии с разностными аналогами уравнений движения и имело вид
2/
Эи
+
■2К1/2 Г
п—1
dz | ик+1/2
+
((Г+1/2 )П-1
п-1 ук+1/2
Вторая аппроксимация (II эксперимент) использовалась в согласовании с уравнениями (25), (26) и соотношениями (28б):
2(
Эи ^ ( ^
I +I &
2(+шГ1<
1к+1/2
+
к+1/2
3.1. Схема проведения расчетов и параметры модели. Для проверки качества модели оперативного прогноза течений в Черном море необходимо сравнить результаты эксперимента с данными наблюдений. Такое сопоставление целесообразно, когда модель воспроизводит близкую к наблюдаемой гидродинамическую ситуацию. Поэтому необходимо задание реальных краевых и начальных полей.
В качестве краевых условий использовались поля тангенциального напряжения трения ветра, полученные на каждые сутки прогноза [7]. Имелись данные о ежесуточных потоках тепла, осадках и испарении, которые были получены по натурным измерениям [7 - 8].
Расчеты проводились на горизонтальной сетке 5 X 5 км, по вертикали использовались 38 горизонтов с глубинами 2.5, 5, 10, 15, 20, 25, 30, 40, 50, 62.5, 75, 87.5, 100, 112.5, 125, 150, 175, 200, 250, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000, 2100 м. Шаг по времени равнялся 5 мин.
Коэффициенты горизонтального турбулентного обмена импульсом и горизонтальной турбулентной диффузии брались следующими:
(н = 1017пх 4 • п-1,
к = 1016ш 4 -1
•п
Коэффициенты турбулентности по вертикали вычислялись в соответствии с параметризацией Меллора - Ямады.
Схема проведения численных экспериментов была следующей. Задавались начальные поля горизонтальных компонент скорости течения, уровня, температуры и солености (условие (15)), которые соответствовали 1 января 2006 г. Эти поля были получены в рамках проекта «Оперативная океанография» [9]. Начальные условия (24) для уравнений (19), (20) задавались следующим образом:
в° = 10ш 2
10 = 10ш
3.2. Результаты численных экспериментов. В обоих расчетах наблюдается качественное соответствие временной изменчивости полей течений, температуры и солености. Как пример рассмотрим эволюцию уровня моря и холодного промежуточного слоя (ХПС).
2
2
2
2
£
2
2
2
z
z
п—1
п-1
2
34
0233-7584. Мор. гидрофиз. журн., 2010, № 3
Приведенный уровень отражает структуру геострофических течений в верхнем слое моря. Глубоководную область моря охватывает циклонический круговорот с центрами в западной и восточной частях бассейна. По его периферии формируются и эволюционируют антициклонические вихри, среди которых наиболее известны Севастопольский и Батумский антициклоны.
В обоих вариантах воспроизводятся известные особенности в структуре полей уровня, течений, температуры и солености. На рис. 1 приведено поле уровня для I эксперимента в различные моменты времени.
Р и с. 1. Рассчитанная уровенная поверхность (см) в 2006 г. для различных моментов времени
По периферии основного циклонического круговорота располагаются Севастопольский и Батумский круговороты. Около Анатолийского и Кавказского побережий квазипериодически возникают и эволюционируют антициклонические вихри.
В зимний сезон происходит интенсивная зимняя конвекция, которая проникает до глубины 100 м (рис. 2, а). К весне ХПС уже сформировался в виде непрерывного слоя холодной воды толщиной примерно 80 м (рис. 2, б). В нем четко выражено ядро, температура которого меньше 7°С. В период весны -лета формируется приповерхностный теплый слой. Ниже его сохраняется ХПС, температура которого постепенно повышается (рис.2, в). К осени ХПС сохраняет качественную структуру, но температура в нем выше 7,5°С. Приповерхностный теплый слой увеличивается по толщине, что приводит к опусканию ХПС (рис. 2, г). Указанные особенности хорошо согласуются с данными наблюдений и современными представлениями об изменчивости температуры в верхнем перемешанном слое.
ISSN 0233-7584. Мор. гидрофиз. журн., 2010, № 3 35
Р и с. 2. Вертикальное сечение в поле температуры (°С) по 43,7° с.ш. для различных моментов времени
Анализ полей течений и солености также подтверждает качественное соответствие результатов обоих вариантов. Иная ситуация наблюдается в распределении кинетической энергии турбулентности и макромасштаба турбулентности. На рис. 3 приведены графики изменения со временем средней по горизонту кинетической энергии турбулентности е2 . Видны качественные различия в их изменчивости. В I эксперименте прослеживаются нерегулярные колебания со значительной амплитудой. Существенные различия между экспериментами наблюдаются в величине энергии. Во II эксперименте ее значения в два раза меньше, чем в I.
В пространственном распределении е2 в I и II экспериментах имеют место качественные и количественные различия. Иллю-
Р и с. 3. Изменение со временем в I и II экспериментах кинетической энергии турбулентности (см2-с-2), осредненной по горизонту 81,5 м
36
0233-7584. Мор. гидрофиз. журн., 2010, № 3
страцией служит рис. 4, на котором представлено распределение энергии на горизонте 81,5 м. В I эксперименте (рис. 4, а) большие значения е2 наблюдаются вдоль свала глубин, они достигают 34 см2*с-2. Во II эксперименте
(рис. 4, б) значительные величины е2 отмечаются вдоль периферии области,
2 -2
максимальное значение энергии не превышает 5 см -с .
28° 30° 32° 34° 36° 38° 40° в.д.
Р и с. 4. Пространственное распределение кинетической энергии турбулентности на горизонте 81,5 м для 15 мая в I (а) и II (б) экспериментах (интервал между изолиниями составляет 0,2 см2-с-2)
188М 0233-7584. Мор. гидрофиз. журн., 2010, № 3
37
Аналогичные различия наблюдаются в пространственной и временной изменчивости макромасштаба турбулентности. В свою очередь это привело к тому, что коэффициенты турбулентности по вертикали в обоих вариантах различаются. Поэтому существует разница в воспроизведении перемешивания в ХПС.
Б 8 10 12 14 16 18 20 22 6 8 10 12 14 16 18 20 22
z, м а Т, °С б Т, °С
Р и с. 5. Профили модельной (кривая 1) и измеренной (кривая 2) температур в I (а) и II (б) экспериментах для 7 июля 2006 г.
В качестве примера на рис. 5 приведены профили измеренной [10] и модельной температур для 7 июля 2006 г. Для слоя 40 - 60 м наблюдается большое отличие в результатах двух расчетов. Этот слой соответствует ядру ХПС, которое четко прослеживается по наблюдениям в этот период. Во II эксперименте модель воспроизводит эту особенность, в I - нет.
4. Заключение. Проведенный анализ уравнений для кинетической энергии турбулентности (25) и макромасштаба турбулентности (26) в рамках параметризации Меллора - Ямады показал высокую чувствительность результатов расчетов к выбранным конечно-разностным аппроксимациям этих уравнений. В первую очередь она определяется их существенно нелинейным характером. Правая часть уравнения (19) содержит два таких члена (первый и четвертый), в уравнении (20) все слагаемые в правой части являются нелинейными. Важно отметить, что последнее слагаемое в (19) имеет особенность в нуле, второе и третье в (20) обладают квадратичной нелинейностью. Поэтому решение сформулированной задачи (1) - (24) зависит от используемых аппроксимаций.
На примере двух конечно-разностных аналогов коэффициента переме-
V
шивания n , который участвует при расчете второго слагаемого в правых
частях уравнений (25), (26), показана такая зависимость.
Несмотря на то, что основные характерные черты в термохалинной и динамических структурах воспроизводятся в соответствии с данными наблюде-
38 ISSN 0233-7584. Мор. гидрофиз. журн, 2010, № 3
ний и имеющимися представлениями, более тонкие особенности зависят от вида конечно-разностных аппроксимаций уравнений (25), (26). Анализ двух типов разностных аналогов правых частей этих уравнений позволил выбрать наиболее подходящую из них.
СПИСОК ЛИТЕРАТУРЫ
1. Демышев С.Г., Коротаев Г.К. Численная энергосбалансированная модель бароклинных течений океана на сетке С // Численные модели и результаты калибровочных расчетов течений в Атлантическом океане. -M.: ИВМ РАН, 1992. -С. 163 -231.
2. Демышев С.Г., Кныш В.В., Коротаев Г.К. Численное моделирование сезонной изменчивости гидрофизических полей Черного моря // Морской гидрофизический журнал. -2002. -№3. -С. 12 -26.
3. Pacanowski R.C., Philander S.G.H. Parameterization of vertical mixing in numerical models of tropical oceans // J. Phys. Oceanogr. - 1981. - 11, №11. - P. 1443 - 1451.
4. Mellor G.L., Yamada T. Development of a turbulence close model for geophysical fluid problems // Rev. Geophys. Space Phys. - 1982. - 20. - P. 851 - 875.
5. Mellor G.L., Yamada T. Users Guide for Three-Dimensional Primitive Equation Numerical Ocean Model // Available on the Princeton Ocean Model web site.
6. Демышев С.Г. Численная энергосбалансированная модель бароклинных течений океана // Дис. ... д-ра физ.-мат. наук. -Севастополь: МГИ НАН Украины, 1996. -343 с.
7. Ратнер Ю.Б., Мартынов М.В., Баянкина Т.М. и др. Информационные потоки в системе оперативного мониторинга гидрофизических полей Черного моря и автоматизация процессов их обработки // Системы контроля окружающей среды. - Севастополь: ЭКОСИ-Гидрофизика, 2005. - С. 140 - 149.
8. Ратнер Ю.Б., Иванчик М.В., Баянкина Т.М. и др. Структура системы и управление вычислительным процессом моделирования динамики Черного моря // Системы контроля окружающей среды. Средства и информационные технологии. - Севастополь: ЭКОСИ-Гидрофизика, 2006. - С. 150 - 158.
9. Дорофеев В.Л., Коротаев Г.К., Ратнер Ю.Б. Система мониторинга гидрофизических полей Черного моря в квазиоперативном режиме // Системы контроля окружающей среды. Средства и информационные технологии. - Там же. - 2004. - С. 9 - 23.
10. Толстошеев А.П., Коротаев Г.К., Лунев Е.Г. Термопрофилирующий дрейфующий буй // Экологическая безопасность прибрежной и шельфовой зон и комплексное использование ресурсов шельфа. - Севастополь: МГИ НАН Украины, 2005. - Вып. 11. - С. 143 - 154.
Морской гидрофизический институт НАН Украины, Материал поступил
Севастополь в редакцию 03.02.09
E-mail: sgdem@stel.sebastopol.ua После доработки 27.04.09
АНОТАЦ1Я У рамках тдходу Меллора - Ямади реалiзована чисельна схема розрахунку кое-фщенив турбулентно! в'язкост i дифузп в z-системi координат для тривимiрноI моделi оперативного прогнозу течш в Чорному морь Дослщжеш дискретш аналоги рiвнянь для юнетично! енергп турбулентност i макромасштабу турбулентность Показана !х висока чутливють до вибору юнцево^зницевих апроксимацш. На основi зютавлення результата прогностичних експерименив з даними спостережень вибрана краща апроксимащя доданку, який описуе ге-неращю енергп турбулентность
ABSTRACT Within the framework of the Mellor-Yamada approach, numerical scheme of calculating turbulent viscosity and diffusion coefficients in the z-coordinate system is realized for a three-dimensional model of operative prediction of the Black Sea currents. Discrete analogues of equations for kinetic energy and macro scale of turbulence are investigated. Their high sensitivity to choice of finite-difference approximations is shown. The best approximation of the item describing turbulence energy generation is chosen based on the comparison of the results of prediction experiments and the observation data.
ISSN 0233-7584. Мор. гидрофиз. журн, 2010, № 3
39
