УЛУЧШЕННАЯ ВЫЧИСЛИТЕЛЬНАЯ МОДЕЛЬ ДЛЯ ТУРБУЛЕНТНОГО АТМОСФЕРНОГО ПОГРАНИЧНОГО СЛОЯ С ПАРАМЕТРИЗАЦИЕЙ ГОРОДСКОЙ ШЕРОХОВАТОСТИ*
А. Ф. КУРБАЦКИй Институт теоретической и прикладной механики СО РАН,
Новосибирск, Россия e-mail: [email protected]
An improved mesoscale model for the turbulent atmospheric boundary layer is formulated. The momentum and heat turbulent fluxes are calculated from the completely explicit algebraic expressions. The urban roughness is parameterized. A model reproduces the observably transformations of a wind field above the urbanized surface, in particular, increase of wind speed at night.
Введение
Для моделирования процессов переноса импульса, тепла и рассеяния примесей в городском пограничном слое в последнее время использовались модели различной степени полноты описания турбулентности и различные схемы параметризации городской шероховатости. Так, например, в [1] использована K — е-модель турбулентности, а воздействие стратификации на турбулентный перенос импульса и тепла учтено введением поправок на стратификацию в коэффициент турбулентной вязкости. Недостатки такого представления общеизвестны. Помимо двух параметров: кинетической энергии турбулентности (КЭТ) и скорости ее диссипации — турбулентная вязкость зависит еще от градиента средней скорости и вертикального турбулентного потока тепла (потокового числа Ричардсона), и, следовательно, турбулентные потоки импульса (и тепла) явно не выражаются через градиенты средних полей, и требуется итерационная процедура. Для учета воздействия шероховатости на процессы переноса тепла и их влияния на городской климат определяющие уравнения Навье — Стокса и притока тепла усредняются не только по ансамблю, но также и по пространству с введением некоторой функции эффективного объема. В другой схеме параметризации [2] используется приближение “пористой городской шероховатости”, при котором индуцируемые зданиями (различной высоты) силы давления и вязкого трения учитываются в виде дополнительных аддитивных членов в уравнениях движения,
* Исследование выполнено при финансовой поддержке Российского фонда фундаментальных исследований (грант № 03-05-64005), Интеграционного проекта Президиума СО РАН (№ 130) и Научной программы “Университеты России” (№ 01.01.190).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.
притока тепла и влажности способом, предложенным в [3]. Такая схема параметризации реализована в простом двумерном тесте эволюции атмосферного пограничного слоя (АПС) с использованием однопараметрической модели турбулентности, в которой только КЭТ определяется из уравнения переноса. Для всех турбулентных потоков использована градиентная модель теории с линейным турбулентным масштабом, рассматриваемым функцией только вертикальной координаты.
В настоящей работе представлена улучшенная модель турбулентного АПС, в которой использованы более полные, чем в [4], модели для корреляций между пульсациями давления и скорости и пульсациями давления и температуры [5] (подробности см. в [6]). Модель позволяет воспроизвести структурные особенности трансформации поля ветра над урбанизированной поверхностью при суточной эволюции АПС, не воспроизводимые с помощью К — Ь- и К — б-технологий моделирования турбулентности.
1. Уравнения для среднего поля
Для двумерного течения в планетарном пограничном слое определяющая система уравнений в приближении свободной конвекции записывается в виде
их + Wz = 0; (1)
Щ + 1111х + \VlJz =--Рх — (ши) г + /V + Ии] (2)
Ро
V* + иУх + WVz = — — /и + В; (3)
Wt + UWx + WWz = -—Pz-{ww)z + |Зeg^, (4)
Ро
0* + иОж + W0z = — {и@)х — {w@)z + Пв ■ (5)
Зависимые переменные в (1)-(5) — это осредненные по Рейнольдсу (т. е. по времени) скорости и, V и W в направлении осей х, у, и г соответственно; Р — среднее давление; 0 — среднее отклонение потенциальной температуры от стандартного значения Т0; в — коэффициент объемного расширения воздуха (3.53 ■ 10-3 К-1); р0 — средняя плотность воздуха; Р — отклонение от гидростатического давления; строчными буквами обозначены турбулентные флуктуации величин. Слагаемые Пи, представляют дополнительный
источник сил (трения, сопротивления формы), индуцируемых взаимодействиями между твердыми поверхностями (земной поверхностью, зданиями) и потоком воздуха, а слагаемое Пв учитывает воздействие явных потоков тепла от твердых поверхностей (зданий или земной поверхности) в балансе потенциальной температуры. Конкретный вид этих дополнительных членов взят, следуя параметризациям [2], и здесь не приводится.
2. Полностью явные алгебраические аппроксимации для турбулентных потоков импульса и тепла
Полностью явные алгебраические модели для напряжений Рейнольдса т^- и вектора турбулентного потока тепла hj в системе (1)-(5) могут быть получены [7-9] путем упрощения в приближении слабо равновесной турбулентности уравнений переноса для турбулентных
потоков к системе связанных алгебраических уравнений и ее решения для АПС с помо-тттью символьной алгебры [6]. Ниже приведены выражения для напряжений Рейнольдса и вертикального потока тепла (формулы для других моментов второго порядка см. в [6]):
((игу), (игу)) = -Км ) ; (6)
'ои тл
дг' дг )
(юв) = -Кн^ + 7С; (7)
Км = ЕтБм, Кн = ЕтБп; (8)
Бм = — /^0 [1 + 8\Он ($2 — $3Сн)\ + ^4^5 (1 + ЗеСн) (тРд)2^ТТ~\ ; (9)
в | и ь , - V * о ^ „ О ч и V , ^ Е
1(2 1
где
$н — у: ^ о— (1 + ^бСя) ^ , (10)
В [3 Сю
7с = — | 1 + -сх2Ом + $бСн ^ <у.ъ{т(3д) І92') (11)
есть противоградиентный член, который в моделях уровней 2 и 2.5 замыкания [4, 6] отсутствует. Величины Он и Ом определяются как
Он = (тN)2 , Ом = (тв)2 , (12)
2 п дв п2 (ьи\2 (0УХ 2
ы 52-(эг] +(э7
и для уравнений (6)—(11):
В = 1 + й\Ом + (І2Он + Л^ОмОн + й^О^Н + \_(І5,О2н — й^ОмОн] Он; (13)
; 2 2 л 10 аз і 2 а3 / \ і-,
йх = -а2, й2 = —---------, й3 = -а2 — {а2 - а5); (14)
3 3 С19 3 С19
11 /а3\2 , 4 /а3\3 2 /а3\2
4 = — — , 4 = - — , Об = оа2«5 — ; (15)
3 \С1^ 3 VС1®/ 3 \С1^/
_ 2 _ 1 /скз\ _ — [ аз\ _
50 — “«2) 51 — ---- ----- , в2 — (У.2 — СК 5, ^3 — «5 ---- , ^4 — 0!з0!5- (16)
3 а2 VС1® / \С1® /
|2
3. Е — є — (в )-модель турбулентности
Для вычисления параметров Е, є и <92 > в (6)—(11) используются замкнутые уравнения баланса [8, 9] для КЭТ, скорости ее спектрального расходования и дисперсии турбулентных флуктуаций температуры соответственно:
ВЕ/Ві = Ви + Ри + Ои — є + Ве ; (17)
Вє/Ві = ВгЦе + Ре + Оє — Г + В; (18)
В<92 >/Ві = Вв + Рв — єв, (19)
где Ве и Ве — источники, учитывающие механические факторы городской шероховатости, вычисляемые, следуя [2]. Подробности вывода модели (6)—(19) можно найти в [6, 8, 9].
4. Вычислительный тест
Горизонтальная протяженность области интегрирования равна 100 км с разрешением 1 км. Вертикальное разрешение равно 10 м в пределах первых 50 м от подстилающей поверхности с последующим растяжением сетки в вертикальном направлении вплоть до высоты 1000 м (выше, до 5000 м, шаг сетки постоянен). Топография поверхности плоская с урбанизированной областью (моделью города) протяженностью 10 км, расположенной в центре вычислительной области с абсциссой от 45 до 55 км. Метеорологические начальные условия определялись заданием геострофического ветра (скорости 3 и 5 м/с) в направлении с запада на восток и атмосферной термической стратификацией, равной 3.5 К/км для потенциальной температуры. Температура поверхности Земли задавалась в виде
<^д(х, 0, £) = 6 вш(п£/43 200), (20)
где £ — текущее время в секундах. Это единственное нестационарное граничное условие задачи, которое моделирует 24-часовой цикл нагревания солнцем земной поверхности. Остров тепла задавался в виде контраста температуры по отношению к температуре поверхности по тому же закону (20), но с амплитудой, увеличенной на 4 град. На поперечных границах для всех искомых функций нормальные производные полагались равными нулю. Такому же граничному условию удовлетворяли искомые функции и на вертикальной границе. Определяющие уравнения модели решены методом переменных направлений в сочетании с методом прогонки на смещенной разностной сетке. Адвективные члены уравнений аппроксимированы второй схемой с разностями против потока [10]. Распределение давления можно вычислить одновременно с вычислением поля скорости из диагностического уравнения Пуассона. В настоящем исследовании при применении модели к АПС с плоской топографией подстилающей поверхности можно предполагать справедливым гидростатическое приближение для вычисления распределения давления. Вертикальная скорость ветра вычисляется квадратурой из уравнения неразрывности, а распределение давления находится в конце каждого цикла вычислений путем интегрирования уравнения для вертикальной скорости. Не зависящее от вычислительной сетки решение получено для сетки 120 х 50. Шаг по времени выбирался из условия сохранения точности; вычисления проведены с шагом, равным 0.625 с.
5. Результаты численного моделирования
Результаты численного моделирования 24-часовой эволюции АПС представлены в виде графиков. Кривые локальной скорости трения на рис. 1 получены осреднением вычисленных значений в течение 24-часового цикла моделирования. Вычисленные профили средней скорости горизонтального ветра для двух значений скорости геострофического ветра (и<з = 3 и 5 м/с) показаны на рис. 2. Вертикальные разрезы отклонений потенциальной температуры (а, в) и среднего горизонтального ветра (б, г) представлены на рис. 3 на 12 часов полуденного времени и для ночного АПС. В согласии с измерениями [18] вычисления (рис. 3,г) указывают на возрастание скорости ветра (как и для дневного АПС, рис. 3,б) на подветренной стороне города (из-за термической циркуляции, перемещенной ветром на подветренную сторону ). В [2] получена несколько иная картина поля ветра, которая не фиксирует повышение скорости ветра для сходных синоптических условий. На рис. 4 приведены контуры КЭТ. На рис. 5 векторное поле среднего ветра и изотахи вертикальной скорости ветра ночного АПС (24 часа полуночного времени) показывают термическую
и* ! и*тах
Рис. 1. Вертикальные профили локальной скорости трения, определенной как (ит) + (ьт))1/4 в центре урбанизированной поверхности, нормализованные на ее максимальное значение. Символами различной конфигурации показаны данные измерений: ♦ — [11—14], I — [15], ▲ — [16]. Линия 1 — скорость геострофического ветра и а = 3 м/с, линия 2 — и а = 5 м/с. Вертикальная координата нормализована на среднюю высоту зданий урбанизированной поверхности.
и*/и
Рис. 2. Вертикальные профили отношения локальной скорости трения к средней скорости горизонтального ветра в центре урбанизированной области. Символы — данные измерений различных авторов, приведенные на рис. 1,бв [17]. Остальные обозначения, как на рис. 1.
20 40 60 80 100 120
X, км
Рис. 3. Вертикальные разрезы отклонения потенциальной температуры (а, в) и скорости горизонтального ветра (б, г) для моделирования со скоростью геострофического ветра и а = 3 м/с на
12 часов полуденного времени (а, б) и 24 часа ночи (в, г).
X, кт
Рис. 4. Контуры кинетической энергии турбулентности ( т2с 2) для ночного АПС (24 часа) при моделировании со скоростью геострофического ветра ис = 3 м/с.
40 50 60 70 80 90 100
X, KM
Рис. 5. Векторное поле скорости ветра и изотахи вертикальной скорости для ночного АПС (24 часа) при моделировании со скоростью геострофического ветра Ug = 3 м/с.
циркуляцию на подветренной стороне города, генерируемую горизонтальным градиентом температуры между городом и его окрестностью.
Выводы
Результаты вычислительного теста воздействия механических факторов городской шероховатости и городского острова тепла на глобальную структуру поля ветра показывают, что улучшенная модель турбулентности для атмосферного пограничного слоя позволяет получить картину трансформации поля ветра в суточной эволюции городского АПС, лучше согласующуюся с данными наблюдений по сравнению с K — L- и K — £-моделями турбулентности.
Список литературы
[1] Vu T.C., Ashie Y., Asaeda T.A Turbulence closure model for the atmospheric boundary layer including urban canopy // Boundary-Layer Meteorology. 2002. Vol. 102. P. 459-490.
[2] Martilli A. An urban exchange parameterization for mesoscale models // Boundary-Layer Meteorology. 2002. Vol. 104. P. 261-304.
[3] Raupach M.R., Antonia R.A., Rajagoplan S. Rough-wall turbulent boundary layers // Appl. Mech. Rev. 1991. Vol. 44. P. 79-90.
[4] Mellor G.L., Yamada T. A hierarchy of turbulence closure models for planetary boundary layer // J. Atmos. Sci. 1974 . Vol. 31, N 10. P. 1791-1806.
[5] Launder B.E. Simulation and Modeling of Turbulent Flows. N.Y.: Oxford Univ. Press, 1996. P. 243-310.
[6] КУРБАЦКИЙ А.Ф. Численное исследование воздействия поверхностного теплового пятна на структуру атмосферного пограничного слоя // Теплофизика и аэромеханика. 2005. Т. 12, № 1. С. 41-60.
[7] Cheng Y., Canuto V.M., Howard A.M. An improved model for the turbulent PBL // J. Atmos. Sci. 2002. Vol. 59. P. 1500-1565.
[8] Kurbatskii A.F. Computational modeling of the turbulent penetrative convection above the urban heat island in a stably stratified environment // J. Appl. Meteor. 2001. Vol. 40, N 10. P. 1748-1761.
[9] Курбацкий А.Ф., Курбацкая Л.И. Проникающая турбулентная конвекция над островом тепла в устойчиво стратифицированной окружающей среде // Изв. РАН. Физика атмосферы и океана. 2001. Т. 37, № 2. С. 149-161.
[10] Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. 516 c.
[11] Rotach M.W. Turbulence Within and Above an Urban Canopy. ETH Diss. 9439, 1991. 249 p.; published as ZGS, Heft 45, Verlag vdf, Zurich.
[12] Rotach M.W. Turbulence closure to a rough urban surface. Pt I: Reynolds Stress // Boundary-Layer Meteorology. 1993. Vol. 65. P. 1-28.
[13] Rotach M.W. Turbulence closure to a rough urban surface. Pt II: Variances and gradients // Boundary-Layer Meteorology. 1993. Vol. 65. P. 1-28.
[14] Rotach M.W. Profiles of turbulence statistics in and above an urban street canyon // Atmos. Environ. 1995. Vol. 29. P. 1473-1486.
[15] Oikawa S., Meng Y. Turbulence characteristics and organized motion in a suburban Roughness Sublayer // Boundary-Layer Meteorology. 1995. Vol. 74. P. 289-312.
[16] Feigenwinter C. The Vertical Structure of Turbulence above an Urban Canopy. Ph. D. Thesis. Univ. Basel, 1999. 76 p.
[17] Roth M. Review of atmospheric turbulence over cites // Q. J. R. Meteorol. Soc. 2000. Vol. 126. P. 941-990.
[18] Bornstein R., Johnson D.S. Urban-rural wind velocity differences // Atmos. Environ. 1977. Vol. 11. P. 597-604.
Поступила в редакцию 2 июня 2005 г.