CC BY
12
Владикавказский математический журнал Апрель-июнь, 2003, Том 5, Выпуск 2
УДК 001.891.57
РАСПРОСТРАНЕНИЕ ЗАГРЯЗНЯЮЩИХ ВЕЩЕСТВ В ГОРНЫХ УЩЕЛЬЯХ
Е. С. Каменецкий, А. А. Радионов
В работе рассматриваются результаты численного моделирования суточных изменений распространения загрязняющих веществ в горных ущельях. Показано влияние режимов течения воздуха в ущелье на приземную концентрацию загрязняющих веществ, источник которых расположен на дне ущелья.
Загрязняющие вещества, выбрасываемые в атмосферу в горных ущельях, рассеиваются значительно хуже, чем на равнине. В то же время, в связи с освоением горных территорий и увеличивающимися транспортными потоками по автомагистралям, расположенным в горах, выброс загрязняющих веществ в ущельях непрерывно увеличивается. Измерение полей концентрации загрязняющих веществ в горах связано с большими трудностями, чем на равнине. В связи с этим необходимо разрабатывать математические модели, которые позволяют оценивать концентрацию загрязняющих веществ при различном положении источников выбросов, различной форме ущелий и учитывать влияние различных термогидродинамических факторов.
Обычно такие задачи решаются на основе системы уравнений Рейнольдса, с некоторыми упрощающими допущениями [1-5]. В данной работе также использовалась эта система уравнений, в которой атмосфера предполагалась несжимаемой и учитывалась только составляющая сил Кориолиса, перпендикулярная к поверхности земли. При учете сил плавучести использовалось приближение Буссинеска.
Для того чтобы расчетная область имела прямоугольный вид, применялось преобразование координат из прямоугольной системы координат {X, У, к новой криволинейной системе координат {х,у,г}:
_ 2 ~ г9 (х-, у)
X Л, у 1, 2 . ,,
1 (ж, у)
где функция гд = гд(Х,У) определяет превышение рельефа подстилающей поверхности относительно некоторой нулевой высоты и не меняется при преобразовании координат, кроме того, полагается, что высота расчетной области равна единице, а ее верхняя граница является прямой линией.
Уравнения Рейнольдса приводились к безразмерному виду, причем в качестве масштаба скорости ветра Ми использовалась максимальная составляющая геострофического ветра на верхней границе расчетной области, масштаба длины Мг — высота верхней границы расчетной области над дном ущелья, которая принималась равной 3.33 -т- 4.0
© 2003 Каменецкий Е. С., Радионов А. А.
высоты ущелья. При меньших значениях отношения этих высот возможно существенное влияние граничных условий на картину течения внутри ущелья, а при больших значениях либо необходимо слишком большое количество расчетных точек, либо течение внутри ущелья рассчитывается с малым пространственным разрешением. Масштаб давления принимался равным Мр = рд(ЗМгМ$. Число Эйлера Ей = Мр/рММасштаб времени М/ = Мг/Ми выбирался таким образом, чтобы он равнялся единице. Масштабом температуры является разность максимальной температуры на дне ущелья и на верхней границе расчетной области. Принятый выбор масштабов гидротермодинамических величин обеспечивал максимальные значения соответствующих безразмерных величин порядка единицы.
В качестве начального приближения использовалось решение одномерных уравнений для входной границы расчетной области. В преобразованных координатах это решение считалось справедливым для всех точек расчетной области.
Уравнение для определения поля потенциальной температуры имело вид
^ иНО. ^-^„А'е 1 1 9и 90
дЬ дх дх Рг Рг (1 — х9)2 дх дх
{ЁО. _ , }-(( дидв
уравнение для проекции скорости ветра на ось х
с)и ди (дт) 1_% дт) \
Ш+"Тх+™Тг = "Еи Ц - Ы*Т=7аТг) + - ^ +
1 ди ди {ди , , 1 — г; дь> \ / ди 1 ^ х ди
(1 — ха)2 дх дх \дх 9 х 1 — хд дх) \дх 9 х 1 — хд дх
Л , 1 — г; \ 2 ди ди + \}г9)УТ^~д) уравнение для проекции скорости ветра на ось у
dv.dv.dv „ , ч . < 1 дуду
——I- и——|- го— = — 1 {< > (и — 1/д) + ^Д'г> +
дЬ дх дх (1 — х9) дх дх
(ди . 1 — г; дь>\ /ду . 1 — г; ду\ Л , 1 — г; \ 2 ди ди
и уравнение неразрывности
ди дт 1 1
^ + Тг ~ {гв)*Т=7ди~ =
Здесь и, V и и) — проекции скорости ветра на оси х, у их соответственно, 11д и Уд — проекции скорости геострофического ветра, 0 — потенциальная температура.
Уравнение Пуассона для определения отклонения поля давления от учитываемого при помощи геострофического перепада давления записывалось без учета вязких членов, поскольку их роль мала [3]
^ л / пд(и,ьи) 1 /ду , , 1 — г ду , , 1 — г; ди\ 1 1 дб ЕйА'р = 2 а) ' + — — - (ха)х--— + {х9)у--— +
д(х,х) Ко \ дх 9 х 1 — хд дх 9 у 1 — хд дх ) ¥х 1 — г9 дг
Zg)x ( du du\ (zg)y ( dv du\ 1 — z du\
-u-
1 — zg \ dx dz J 1 — zg \ dz dx J 1 — zg dz J
Для определения энергии турбулентности использовалось уравнение db db db „ „ ^ 1 . ,, 1 1 du db
1 f du 1 — z du\ f db l ^ z db\ 1 f 1 — 2; \2 du db + 7Г 77- - (Zg)x--- (Zg)x---TJ- + — (Zg^
Б к \ dx 9 х 1 — дг) \дх 9 х 1 — dz) Б к \ 9 у 1 — Zg ) дг дг''
где Б к = 0.7 -т- 1 — константа, И = С — скорость диссипации энергии турбулентности, С = 0.046 — константа. В этих уравнениях коэффициент турбулентной вязкости определяется по формуле Колмогорова и = сЩь1'2, где Ь — энергия турбулентности, I — масштаб турбулентности, ^9)х = дzg(x,y)/дх, (гд)у = dzg(x,y)/ду, Ко = Ми//Мг — число Россби, / - параметр Кориолиса, Рг — число Прандтля. Оператор Лапласа имел вид
^ ~ Я„.2 + ,Л я „2 ¿уЯд/з
dx2 (1 — Zg)2 dz2 9 ж (1 — Zg) dxdz
— ^2 {{Zg)XX + (Zg)yy) ^ ^ + 2 ((Zg)2 + (Zg)y) ^ _ ^ ^ 2 ^ '
Генерация энергии турбулентности потоком воздуха записывалась в упрощенном виде
г< 1
G = и-
/du\2 idvs 2 \dz J \dz
Генерация или диссипация энергии турбулентных пульсаций силами плавучести записывалась в виде
Р = 1.
Гг Рг dz
Для определения масштаба турбулентности использовалась формула [5, 10]
- - - + 1
1 lo I'min
где ¿min — минимальное расстояние до ближайшей твердой поверхности, умноженное на постоянную Кармана, максимальное значение масштаба Iq связано с распределением энергии турбулентности по высоте
оо / 00
10 = 0.1 J Vbz dz / J Vbdz. 0 ' 0
В качестве граничных условий на верхней границе расчетной области принимались следующие:
dp db
u = Un, v = VQ, 0 = Const, — = 0, — = 0.
dz dz
Как правило, поток воздуха втекает в расчетную область через левую границу, а вытекает из нее через правую границу. Но, в некоторых вариантах расчетов, рассматривалось обратное направление воздушного потока, а также возникали разнонаправленные потоки на одной и той же границе. Граничные условия на этих границах имеют вид:
(&=»Г /&=»Г |0=»Г
I ^м — О I ' I — — о I ' I -^ — 0 I
V дх - и / ц<0 ^ дх - и / ц<0 ^ дх - и ) и<0
^ др д'ш ди дт ди
Ей— =и—--т—, — - — . (г)х = 0.
дх дх дх дх дх
В этих выражениях положительное значение горизонтальной скорости соответствует вытеканию из расчетной области, а отрицательное — втеканию.
Эти граничные условия использовались только в расчетах с относительно маленькими шагами по времени (равными 1 секунде, при масштабе времени 200 секунд). При использовании больших шагов по времени (5-10 секунд, при масштабе времени 200 секунд) эти граничные условия уже не могут быть использованы, поскольку в пульсирующих вариантах решения, вихревые структуры, сносимые потоком вниз по течению, приводят к расходимости решения при прохождении границы. Для больших шагов по времени использовались другие граничные условия, получаемые из соответствующих уравнений. В уравнениях вторая производная по горизонтальной координате принималась равной нулю, а в случае, если поток втекает через рассматриваемую границу, также принималась равной нулю и первая производная.
Для того чтобы не ставить граничное условие на верхней границе на вертикальную составляющую скорости ветра, проекция уравнения количества движения на вертикальную ось отбрасывалась. Эта составляющая скорости ветра вычислялась из уравнения неразрывности, а для определения поля давления решалось уравнение Пуассона. Здесь для горизонтальных составляющих скорости ветра выбирались условия и = 11д, V = Уд, для потенциальной температуры, давления и энергии турбулентности первая производная по нормали к этой границе принималась равной нулю: дв/дх = 0, др/дх = 0, дЬ/дх = 0.
В качестве граничных условий на поверхности земли принималось равенство нулю скорости ветра, а для определения давления использовалось значение производной др/дх, полученное из проекции уравнения количества движения на вертикальную ось, записанного на поверхности земли
и = 0, V = 0, ю = 0;
Ей
1-х др 1 + {{г9)1 + {х9)2у){1 - х)2 др
\%д)х „ г. +
1 — Хд дх
(1 -
дх
1__1_
Ег 1 — х„
1 + + ( д2и , (ха)х ди , {х9)у ду
(1 "
+
+
дхдх 1 — х9 дх 1 — х9 дх
2(т\(т\ (1 ~ Х)2 9и 2(т ^ (1 ~ Х)2 9и 2(т ^ (1 ~ Х)2 9У '¿\*д)х\*д)хх ,Л у> д)х ¡л ..я а.. ^Х^д/г
(1 ^Хд)2дХ
(1 - дх К9>* (1 - Хд)3 дх
'(1 дх
(1 ^Хд)*дХ
Поскольку энергия турбулентности в приземном слое, имеющем малую толщину, приблизительно постоянна, ее производная по нормали к поверхности земли принималась равной нулю: дЬ/дх = 0. Температура на поверхности земли вычислялась с помощью уравнения теплового баланса
дТ$ Н4 / 2тг \ , . дТт На
= (ТЗ - Тт)
рзсзаг \ П )
дЬ рзсзЗ 1 \ П / ' дt рзсзЗ-!'
где Тз — температура подстилающей поверхности, Тт — температура слоев почвы, расположенных на глубине проникновения годовых колебаний температуры; Н4 — сумма всех потоков тепла из атмосферы к подстилающей поверхности; рд — плотность почвы; С3 — удельная теплоемкость почвы; й\ = л/К^гг — толщина эффективного слоя почвы, пропорциональная глубине проникновения суточных колебаний температуры; — удельная теплопроводность почвы; т\ — суточный период времени; (¿2 = л/Зб!ГЖЩргг — толщина слоя почвы, пропорциональная глубине проникновения годовых колебаний температуры [2, 7, 8].
При рассмотрении атмосферы без учета влажности, уравнение баланса тепла на поверхности земли можно записать [6, 7]:
НА = (1 — 05) Яя — Яьи + Яьб — ЯЕ,
где аз — альбедо подстилающей поверхности; Яз — поток солнечной коротковолновой радиации на поверхность; Яьи — поток длинноволновой радиации, излученный склоном; Яьв — поток длинноволновой радиации, поглощенный склоном (отраженный небом, либо другими поверхностями); Яе = рсри*Т* — поток явного тепла, р — плотность воздуха, ср — теплоемкость воздуха, и* — динамическая скорость, Т* — размерная температура воздуха (Т, « Мв9*) [6, 7].
При штилевых условиях перепад температуры между уровнями г = 0 и уровнем вязкого подслоя х = хо может достигнуть довольно больших значений. Поэтому для расчета температуры воздуха на нижнем вычислительном уровне использовалась полуэмпирическая параметризационная формула для вязкого подслоя [10]: 9 — 0$ = 0.09620* (и^хо/г^)°'45, где V — кинематический коэффициент вязкости, хо — высота слоя шероховатости, и* — динамическая скорость.
Концентрация загрязняющих веществ определялась путем решения уравнения диффузии в преобразованных координатах.
дс ^ дс ^ дс ^ д/ ^ дь> дс
дЬ дх дх (1 — хд)2 дх дх
д (дь> , , 1 — 2 дь> \ / дс , ч 1 — 2 дс
На входной и выходной границах концентрация определялась из этого же уравнения, в котором вторая производная концентрации по нормали к границе равнялась нулю. В тех точках, в которых наблюдалось вытекание через эту границу, решалось уравнение
дс дс дс . 1 дь> дс
+ и——Ь ь)— = у А с +
ей дх дх (1 — ) дх дх
"Чз* -*ддг)[дх ^я)х1-гддг) + {^я)у1-гд; дхдх
(ду , . \ — г ду\ ( дс , . \ — г дс\ (, . 1 ^ х\'2 ди дс
а в тех точках границы, в которых наблюдается втекание через границу в расчетную область, в этом уравнении принималась равной нулю также и первая производная концентрации по нормали к границе, и уравнение в них имело вид:
дс дс , 1 дь> дс
ей дх (1 — хд)2 9х дх
/ди .. 1 — г ди\ ( . \ — г дс\ ( 1 — г; \2 ди дс
Оператор Лапласа при условии (гд)х = 0 и д2/дх2 = 0 имеет вид:
^ — ч2 я ..О. ¿\^9>УУ /1 .. \ ^¿(¿9)1
(1^9)2 дг* V "а-^)2.
Поверхность земли и верхняя граница предполагались непроницаемыми, т. е. первые производные по нормали к этим границам равнялись нулю.
При расчетах распространения загрязняющих веществ их источник в большинстве расчетов предполагался линейным, что позволяло проводить расчеты в двумерном приближении. Таким источником может быть автодорога, если поток автотранспорта достаточно интенсивен. Другим возможным источником загрязнения, который также можно считать двумерным, является хвостохранилище. Предполагалось, что автодорога расположена в середине ущелья, а выбрасываемые вещества равномерно рассеиваются в следе за автомобилями, высота которого порядка 2-3 метров [4]. Значения концентрации нормировались на их величину над дорогой.
Конечно-разностная аппроксимация
Для решения системы уравнений применялся конечно-разностный метод. Уравнение Пуассона для давления интегрировалось при помощи схемы «классики». Для прогностических уравнений, для проекций скорости ветра, температуры, энергии турбулентности и концентрации загрязняющего вещества использовалась неявная итерационная схема.
Применялась неравномерная сетка с более мелкими шагами вблизи поверхности земли. Предполагалось, что зависимость преобразованной вертикальной координаты от номера точки степенная. Показатель степени принимался равным трем вторым Zj = J3//2. Расчеты осуществлялись на совмещенной сетке с числом узлов М х N. Как правило, по вертикали принималось N = 60 узлов сетки, а по горизонтали М = 180. На входном участке расчетной области нижняя граница полагалась горизонтальной. Этот участок занимал 35-60 слоев сетки, что обеспечивало малую зависимость скорости ветра на входной границе от течения в ущелье. Расстояние от выходной границы до края ущелья также обычно составляло от 35 до 60 точек, что связано с более сильным влиянием выходного граничного условия на картину течения внутри ущелья. Это влияние особенно заметно проявляется при прохождении вихревых структур через выходную границу, что наблюдается во время установления картины течения, а также при расчете суточного изменения картины течения во время смены нагрева подстилающей поверхности на охлаждение и обратно. Если выходной участок занимал 15 слоев по горизонтали и меньше, то картина течения в ущелье существенно отличалась от той, которая наблюдается при большем числе слоев.
Для устранения схемной диффузии конвективные слагаемые во всех уравнениях аппроксимировались конечными разностями первого порядка в направлении против потока. В конечно-разностной аппроксимации производных в диффузионных слагаемых использовались только центральные разности, т. е. производные выбирались второго порядка точности.
Для вычисления инеобходимо пользоваться давлением на том же временном ша-к+1
ге ^ , поскольку определяющим при вычислении скорости является именно давление. Перепад давления распространяется по расчетной области очень быстро, поэтому необходимо учитывать давление с (к + 1)-ого временного слоя.
Использовалась вычислительная схема с весом, то есть для вычисления значения функции на следующем временном шаге использовалась запись
= (1 ^ а) Ф^-+ (1)
где Ф^ • — правая часть уравнения для определения этой функции, вычисленная на предыдущем временном слое, Ф^1 — правая часть того же уравнения, вычисленная на последующем временном слое, а — весовой коэффициент. Обычно весовой коэффициент равнялся 1/2. В некоторых случаях весовой коэффициент равнялся 3/4 для диффузионных слагаемых. В уравнениях, записанных для граничных точек расчетной схемы, весовой коэффициент обычно равнялся единице.
Неявное слагаемое Ф^1 в (1) предполагает, что все члены уравнения записываются на (к + 1)-ом временном слое. При этом нелинейная часть уравнения записывается на промежуточном временном слое /*•, превращаясь, таким образом, в расчетный параметр. В пределах одного шага по времени, далее, при помощи итераций, производится
< в/, где £f -
_ £к +1
вычисление /* • до тех пор, пока не выполняется условие тах
задаваемая точность решения уравнения для функции /. При этом после каждой итерации, в случае не выполнения указанного условия принимается Д^1 = после чего расчет повторяется вновь. Первая итерация начинается при условии = Однако решение, полученное с помощью такой схемы, требует дополнительного согласования полей расчетных функций. Для этого дважды применялась процедура коррекции решения.
Результаты расчетов
В численных экспериментах рассматривалось два ущелья: симметричное, профиль которого приведен на рисунке 1 (а), и ущелья, расположенного на Северо-Кавказском хребте в ущелье реки Ардон в районе унальского хвостохранилища, поперечное сечение которого приведено на рисунке 1 (б).
Рассмотрим суточный ход концентрации загрязняющего вещества, выбрасываемого линейным придонным источником, расположенным приблизительно в центре ущелья.
На рисунке 2 приведены суточные изменения температуры подстилающей поверхности, где по оси абсцисс отложено время, а по оси ординат — температура в градусах Кельвина.
На рисунке 3 представлены временные зависимости приземной концентраций выброса источника загрязняющего вещества для двух равноудаленных от источника точек поперечного сечения симметричного ущелья.
0,12 0,1 0,08 0,06 0,04 0,02
10
20
30
40
50
60
70
0,12 0,1 0,08 0,06 0,04 0,02 0
0 20 40 60 80 100
Рис. 1 (а) Профиль симметричного ущелья. В расчетах поперечное сечение ущелья дополнялось слева и справа горизонтальными участками, занимавшими по 55 точек расчетной сетки.
Рис. 1 (б) Профиль «унальского» ущелья. В расчетах поперечное сечение ущелья дополнялось слева и справа горизонтальными участками, занимавшими по 35 точек расчетной сетки.
23:00
3:00
7:00
11:00
15:00
19:00 час
Рис 2. Суточный ход температуры подстилающей поверхности в проведенных экспериментах
0,14 0,12 0,1 0,08 0,06 0,04 0,02 0
23:00
3:00
7:00
11:00
15:00
19:00 ^ час
х/И=1 с наветренной стороны ущелья
х/И=1 с подветренной стороны ущелья
Рис 3. Приземная концентрация загрязняющего вещества в точках равноудаленных по обе стороны от источника на расстояние х/Ъ. = 1,00.
Вечером приземная концентрация загрязняющего вещества в ущелье низка, поскольку нри рассматриваемом расположении источника загрязнения оно выносится из ущелья вверх. Ночью возникают колебания приземной концентрации большой амплитуды. Эти колебания не имеют определенного периода и при умеренных средних значениях
о
о
приземной концентрации именно ночью достигаются ее максимальные величины. С наступлением утра колебания практически прекращаются и в течение всего дня приземная концентрация загрязняющего вещества достаточно велика, причем оно сносится то к одной, то к другой стенке ущелья.
x/h=0,33 с наветренной стороны ущелья —■— x/h=0,33 с подветренной стороны ущелья
Рис. 4. Суточный ход приземной концентрации загрязняющего вещества в эксперименте с более высокой наветренной стороной ущелья, для точек равноудаленных по обе стороны от источника на расстояние x/h = 0,33.
В экспериментах для реального профиля ущелья, изображенного на рисунке 1.6, источник загрязняющего вещества располагался в середине участка поверхности, которая имеет наименьшую высоту.
х/И=0,33 с наветренной стороны ущелья —■— х/И=0,33 с подветренной стороны ущелья
Рис. 5. Суточный ход приземной концентрации загрязняющего вещества в эксперименте с более высокой подветренной стороной ущелья, в точках равноудаленных по обе стороны от источника на расстояние х/Ъ. = 0,33.
Распространение загрязняющего вещества в несимметричном реальном ущелье существенно другое (рисунки 4, 5). Здесь ночью приземные концентрации загрязняющего вещества стабильно высоки, а утром они уменьшаются. Интенсивные пульсации концентрации загрязняющего вещества начинаются здесь днем и продолжаются до вечера, если ветер над ущельем направлен в сторону более низкой стенки. Если же направление вот-
ра противоположно, то вечером пульсации концентрации возрастают, и пульсирующий режим течения сохраняется до глубокой ночи.
Кратковременное достижение максимальных суточных значений приземной концентрации возникает либо в начале дня, либо поздним вечером в зависимости от направления внешнего ветра. Отметим, что в реальном ущелье расчетные приземные концентрации загрязняющего вещества больше, чем в симметричном ущелье при одинаковой интенсивности источника. Среднее дневное значение приземной концентрации может быть несколько ниже, чем ночью или примерно таким же, что также зависит от направления внешнего ветра.
Полученные результаты показывают, что в условиях горного ущелья измерение концентрации загрязняющего вещества должно производится гораздо чаще и длительность замера должна быть достаточно длительной — примерно от 30 минут до часа.
Очевидно также, что перенос результатов исследования одного ущелья на другое практически невозможен.
Литература
1. Beniston М., Wolf J. P., Beaistoa-Rebetez M., Kolsch H. J., Rairoux P., Woste L. Use of Lidar Measurements and Numerical Models in Air Pollution Research // J. Geophysical Research.—1990.— V. 95, № D7.-P. 9879-9894.
2. Yamada Т., Bunker S. Development of a Nested Grid, Second Moment Turbulence Closure Model and Application to the 1982 ASCOT // Brush Creek Data Simulation. J. Appl. Meteorology—1988 —V. 27, № 5.-P. 562-578.
3. Eager L., Korasia D., Yang X. A Numerical Method of Boundary-Layer Dynamics a Mountain Valley // Boundary-Layer Meteorol.-1993.-V. 66.-P. 357-394.
4. Yamartino R. J., Wiegand G. Development and Evaluation of Simple Models for the Flow, Turbulence and Pollutant Concentration Fields within an Urban Street Canyon // Atmospheric Environment.— 1986.—V. 20, № 11—P. 2137-2156.
5. Svoboda J. Numerical Modeling of the Atmospheric Boundary Layer over a Hilly Landscape // Studia geoph. et geod.—1990.—V. 34.-P. 167-184.
6. Барри Роджер Г. Погода и климат в горах.—Л.: Гидрометеоиздат, 1984.
7. Danard М. A Prognostic Model for the Surface Temperature, Height of the Atmospheric Boundary Layer, and Surface Wind. 1988.
8. Владимиров С. А. Численное моделирование распространения пассивной примеси в атмосфере // Метеорология и гидрология.— 1999.—№ 7.—С. 22-35.
9. Шелковников М. С. Мезометеорологические процессы в горных районах и их влияние на полеты воздушных судов.—Л.: Гидрометеоиздат, 1985.
10. Пененко В. В., АлоянА. Е. Модели и методы для задач охраны окружающей среды.—Новосибирск: Наука, 1985.
Статья поступила 30 мая 2002 г.
Каменецкий Евгений Самойлович, к.ф.-м.н. г. Владикавказ, Институт прикладной математики и информатики ВНЦ РАН, E-mail: eskOalanianet. ru
Радионов Анатолий Анатольевич г. Владикавказ, Институт прикладной математики и информатики ВНЦ РАН
