CC BY
13
УДК 004.942::524.3-6
УЧЕТ ФАКТОРА ИЗЛУЧЕНИЯ ПРИ МОДЕЛИРОВАНИИ ПРОЦЕССА ОБРАЗОВАНИЯ И РАСПРОСТРАНЕНИЯ ЗАГРЯЗНЕНИЙ В ВОЗДУШНОЙ СРЕДЕ
ПЕКУНОВ В.В., канд. техн. наук
Рассматривается вопрос о включении фактора излучения в математическую модель, описывающую процессы образования и распространения загрязнений в воздушной среде. Приводятся базовые уравнения динамики загрязнений, уравнения переноса прямого солнечного и диффузного (солнечного и теплового инфракрасного) излучения. Даны выражения для вычисления оптических характеристик среды, содержащей капельную фазу.
Ключевые слова: загрязнение в воздушной среде, математическая модель, солнечное и тепловое излучение.
RADIATION FACTOR ACCOUNTING IN SIMULATING THE PROCESS OFAIR POLLUTION FORMING AND SPREADING
V.V. PEKUNOV, Ph.D.
This work is devoted to the problem of radiation factor inclusion into mathematical model, which describes the processes of air pollution forming and spreading. The author gives base equations of pollution dynamics, the equations of direct solar and diffuse (solar and heat infrared) radiation. Formulas for calculating visual characteristics of the environment, containing spot phase, are also given here.
Key words: air pollution, mathematical model, solar and heat radiation.
Моделирование процессов образования и распространения загрязнений требует учета множества факторов, среди которых необходимо особо выделить излучение. Прежде всего излучение влияет на процессы образования вторичных загрязнений в результате фотохимических реакций, что характерно, например, для сухого лос-анджелесского смога. Кроме того, излучение оказывает существенное воздействие на процессы теплообмена, которые, в свою очередь, влияют как на образование, так и на распространение загрязнений. Однако в известных автору математических моделях распространения загрязнений излучение либо просто не принимается во внимание, либо учитывается со значительными упрощениями (например, интенсивность излучения считается постоянной и не учитывается сложная форма расчетной области). Поэтому актуальна задача разработки такой модели образования и распространения загрязнений, которая в достаточно полной мере учитывала бы фактор излучения.
Сформулируем фрагмент трехмерной модели, включающий базовые уравнения переноса и образования загрязнителей, уравнения переноса прямого солнечного и диффузного излучения, а также соотношения для расчета оптических характеристик среды при наличии в ней капельной фазы.
В расчетной области введем прямоугольные координаты (Х1, Х2, хз) таким образом, чтобы ось Охз была вертикальной. Воспользуемся моделью течения несжимаемой вязкой жидкости в системе «вектор вихря - векторный потенциал», которая обладает высокой стабильностью в расчетах. Уравнения для трех компонентов вектора вихря ю запишем с использованием эффективной вязкости Vэфф = vмол + V-, где vмол -молекулярная вязкость, а V - турбулентная вязкость: дю у ^ дю у ^ д ( дю j
~дТ + ^ и1 — = ^ — \уэфф —
i=1
ЭХ;
i=1 dx
ЭХ;
V 8UJ
Ibt i=1 oxi
F, ; j = 1, 2, 3;
Fi = вд ^; F2 = -вд дТ; F = 0,
dx2
дх1
2
где р - термический коэффициент расширения воздуха; д = 9,81 м/с2.
Для расчета турбулентной вязкости может быть применена любая адекватная модель турбулентности, например, были успешно испытаны модели К-Е и Аб-рамовича-Секундова.
Вектор вихря равен ротору скорости и:
dU3 3U,,
dU1 dU
3
дЩ DU,
'2
И>1 — - -~~ , \Х>2 — - - 1 и-»3 — - - .
1 дх2 дх3 дх3 дх1 3 дх1 дх2
Вектор скорости равен ротору векторного потенциала у:
д1з д12 . и -М-д1з . и. -^12 д11
U1 =
дх2 дх3
дх3 дх1
U3
дх1 дх2
Для векторного потенциала запишем три уравнения Пуассона:
I
дЧ j
-- = -ю j
i=1 д х/
j = 1, 2, 3.
При этом граничные условия ставятся согласно [2]. Уравнение для температуры среды Т имеет вид:
-+ I Ui дТ = I JL({Dt + ат Vm)) 1 +i^P, 8t ¿i 1 дх дх, ^ T T mJ дхi J cp
i=1
дxi
где йт - коэффициент температуропроводности; ат -вспомогательный коэффициент; - приток энергии излучения; р и с - плотность и теплоемкость воздуха.
Особо следует сказать о граничных условиях для температуры на стенках. Фактически, в верхнем слое любой стенки, перпендикулярной некоторой оси ху, решается двумерное уравнение теплопроводности
дТ> = а2 удТ + Щ +*и (Т ~тз)
д* "¿у дх} с3Рз •
где Те - температура стенки; а|, с5, р5 - температуропроводность, теплоемкость и плотность материала стенки соответственно; ^ - коэффициент, Вт-м"3-К-1,
характеризующим теплопередачу на границе стенки и окружающей среды; Щь - член, выражающий обмен лучистой энергией с окружающей средой. Выражения для Щ/ и Щь приведены ниже.
Уравнения диффузии веществ с концентрациями С имеют следующий вид:
ВО,
И
L + V uj ° = t-LI(
ti ' dx, fe дх, II
Dr + a^ v
0 mi ÜX;
дО ,
6С; -
-АС, +—-;- = 1,Ы ,
! а ^
где Ц = Ц|, и2 = и2, и3 = Ц + Щ-; Щ - скорость витания --го вещества; йС. - коэффициент диффузии --го вещества; а С , - вспомогательный коэффициент
для --го вещества; N - число веществ; член АС- выражает изменение концентрации --го вещества при
дС,
взаимодействии с водяными каплями; выражает
изменение концентрации в результате химических реакций.
Для фотохимических реакций может быть определен стартовый порог интенсивности излучения. В этом случае реакция начнется лишь тогда, когда уровень солнечного излучения достигнет этого порога. При этом можно ограничиться двумя диапазонами излучения: солнечным видимым (длины волн 0,4-0,75 мкм) и солнечным инфракрасным (0,75-4 мкм), выделив в каждом диапазоне две составляющие: прямое и диффузное (отраженное и рассеянное) излучение. Для корректного рассмотрения тепловых процессов необходимо также принять во внимание тепловое инфракрасное (4-100 мкм) излучение, которое будем считать диффузным.
Интегральная освещенность р0, Вт-м-2, создаваемая прямым солнечным излучением в точке (х,у,г), рассчитывается по следующей формуле:
^ (х,у, г) = е-'(ху,г),
где Г0 - освещенность, создаваемая солнечным излучением на границе расчетной области; '(х,у,г) - безразмерная оптическая толщина среды в направлении падающего солнечного луча, определяемая интегрированием по отрезку I (от точки вхождения луча в расчетную область до точки (х,у,г)):
'(х,у,г) = |ре(з)дЗ ,
причем имеют место соотношения: at (з) = Ре (х, у, г)-(1 -и(х,у,г));
Рз (х, у, 2 ) = и(х, у, г )-Ре (х, у, г),
здесь ре, рз, аt - коэффициенты, м-1, ослабления (экс-тинкции), рассеивания и поглощения соответственно; и - альбедо однократного рассеивания.
Для расчета оптической толщины среды в направлении падающего луча используется упрощенный вариант обратной трассировки луча (метода, широко применяемого в машинной графике). Из каждой ячейки (К,У,Т) расчетной сетки в направлении, обратном направлению падения солнечного луча (Ой), «выпускается» Ыь лучей и отслеживаются ячейки, через которые они пройдут. Для каждого /-го луча вычисляется оптическая толщина
l, (X,Y, Z) =
если луч попал в стенку;
^ ße X,m,n)• ij (k,m,n), иначе,
Pi (xyy ,z )
где P(X,Y,Z) - множество ячеек, через которые прошел i-й луч; i(k,m,n) - длина отрезка луча, ограниченного ячейкой (k,m,n).
Тогда
- (X,Y ,Z ) ye-1; (X,Y ,Z ) Nbü
В рамках диффузионного приближения [1], используя метод сферических гармоник (Р1-приближение), интегральную диффузную интенсивность излучения I, Вт-м"2-ср-1, в направлении Q для любого диапазона можно определить [3, 5] выражением
I (X, у, z, Q) = /0-± ]Т
д/0
■Q,
(1)
Р^ дх/
где pt = ре (1 -ида); да - фактор асимметрии.
Коэффициенты да, и и Ре вычисляются для каждого диапазона излучения по отдельности с учетом наличия в среде капель воды, причем
Ре = Ре0 + Ред ,
где ре0 = ре0(г) - собственный коэффициент ослабления воздуха на заданной высоте г; ред - коэффициент ослабления за счет наличия дисперсной (капельной) фазы.
Выражения для Ре^ и и да будут приведены далее.
Первый компонент разложения интегральной диффузной интенсивности излучения определяется решением следующего уравнения Пуассона для анизотропной среды [3, 5]: з
^ дх-
i=1 дХ1
ßi дх, '3ßsFo
= 3«i/° -
4п
для солнечного излучения,
3at • fT (T), для инфракрасного;
fT (T)= J B(T)dX
где В(Т) - функция Планка.
При типичных расчетах экологии воздушной среды (в диапазоне температур 243-423 К при длинах волн 4-100 мкм) можно упростить вычисления, интерполировав интеграл функции Планка кубическим полиномом:
Т (Т)«1,86 -10-5 - Т3 - 6,932 -10-3 - Т2 + +1,136 - Т - 69,041.
Потоки диффузной освещенности (для любого диапазона) в прямом (+) и обратном (-) направлениях координатных осей задаются соотношением
2 1 д/0
F± х =
JI (х, у, z; q)q x.d Q = п/0 + - п
2п
__0_
ßi дх,
а компоненты вектора Fd потока диффузной освещенности -
X
1 5/0
(Fd) = F+ - F-=-4п-^.
4
3'" pf dXj
Гоаничные условия (в любом диапазоне) для уравнения диффузной интенсивности излучения на открытых границах задаются из условия равенства диффузного потока извне некоторой известной величине Fext:
F+ «(X = 0) = п/00 - !^f = Fext ; F-Xi (X = a) = П/00 + f ^ Щ = Fexf
На закрытых границах (стенках) для солнечного излучения задается условие равенства выходящего диффузного потока Fout полному (включающему прямую солнечную радиацию) входящему потоку освещенности Fin, умноженному на альбедо стенки -s. Пусть n =(n1;n2;n3) - вектор внешней нормали стенки. В зависимости от направления нормали
jn > 0 In, < 0
п/0 - 2 п 1 д/0
= п/0 — — п--
0 3 pf dX;
п/0 + -
1 д/° _
тогда
1 — 2 (1-
3 Pf ÔX; 5) 1 5/0
— F0 -(«0 )x
3 (1 - —)pf Bxi
(«0 )x,. - F0
_ <x- 0
, = —1as (Л '-Г",
П >0) П(1 - — ) П <0J
где Оо - вектор падения прямых солнечных лучей.
На закрытых границах (стенках) для теплового инфракрасного излучения в зависимости от направления нормали
Jn, >0
[n, < 0
(
= п/.
0,/R
0,/R -
2 1
3 neR
d/0°,/R
Pf Ох,
1 д/.
0,/R Л
3 pR dXj
= ("1 -6s ) п/ тогда
0JR — 2 (2-6s ) aÇR 0 — 36s PR dXj
+ 6sfT (T )
= fT (T )
,n, >0)
'n: <0
где в - коэффициент черноты поверхности стенки; верхний индекс ¡И относит величины к диапазону теплового инфракрасного излучения.
Подставив выражения (1) в формулу для вектора потока диффузной освещенности, с учетом вышеуказанного уравнения Пуассона, а также с учетом теплового инфракрасного излучения получим соотношение для локального объемного потока энергии Втм-3, поступающей в среду при поглощении диффузного и прямого солнечного излучения:
Щ (X ^ +Х Р801 + ^ )-
- Х(4(° + агГо) + 4паИ (^ - Т (Т)),
где Г301 - вектор потока прямой солнечной освещенности в заданном диапазоне, причем Ьм (0/) - -Ре£"0 . Суммы берутся по диапазонам солнечного излучения.
Локальный объемный поток энергии излучения Вт-м-3, для элемента стенки с внешней нормалью
n =(n1 n2; n3 )
описывается соотношением
In, >0
Ws j ; !> = 6s s[n, <0
п/.
0,/R
1 д/
0,/R Л
3 pR дх,
-6sT (T ) +
-(1 --s )X
2 1 д/0
П/00 ± 3^ д/^ - F0 («0)x
Рассмотрим оптические характеристики среды, содержащей капельную фазу. Пусть имеется Z компонентов фазы, причем каждому ,-му компоненту соответствуют капли с диаметрами от X, до y, и известна функция распределения капель по диаметрам n,(d).
Фактор асимметрии ga и альбедо однократного рассеивания среды и определим формулами: Z Z
Ре0и0 + X^eZ
и =-=-;
Pe0 + Ped
ga =
Pe0ra0g0 + XPe' Ид
_1=1_
Z
Pe0ra0 + XPe,И
=1
Ped =XPe
I-1
где ио и до - усредненные (по диапазону излучения) альбедо однократного рассеивания и фактор асимметрии незамутненного воздуха; и д, в - альбедо однократного рассеивания, фактор асимметрии и коэффициент ослабления 1-го компонента. Наиболее точно оптические параметры компонентов капельных фаз могут быть рассчитаны по теории Ми, но для прямых оперативных расчетов эта теория малопригодна ввиду ее высокой вычислительной трудоемкости. Существует более простая (но менее точная) теория аномальной дифракции, но при необходимости осреднения оптических параметров по диапазону длин волн (численного нахождения интегралов от таблично заданного комплексного показателя преломления [6]) расчет также будет весьма трудоемким. Поэтому целесообразно применить соотношения, приближающие результаты точных вычислений по теории Ми. Применим интерполяцию непосредственно по диаметру капель, избегая вычисления эффективного диаметра капель компонента (то есть дополнительного осреднения, примененного, например, в [4]).
Определим:
УI
| п (й)- в, (й)Ьй
Pe, -Pa
g, =—
Pe,
У,
Pe,= 4 j n (D)- Qe, (D)
N'k
- D2dD ;
Pa,= 4 j n (D)- Qa, (D)- D2dD,
где в,, Ое, и Оа, - усредненные (по диапазону излучения) безразмерные коэффициенты асимметрии, ослабления и поглощения одной капли 1-го компонента; Ра, - коэффициент поглощения 1-го компонента. Функции в(й), Ое(й) и Оа/Р) являются результатом полиномиальной интерполяции по методу наименьших квадратов. Исходные точки для интерполяции определяются функциями:
s
1 ^
QQei(D ) = р J У (a.)- qe (X,D )d X ;
А-1 1 X2
Qai(D) =- J 3?(X)- qa (X,D)dX ;
Xl
1 X2
G (D) =1 J У (X)- gm (X,D)dX ;
"2
P = J У(X)dX,
где X и ^2 - начальная и конечная длины волн рассчитываемого диапазона излучения; у (X) - нормированная функция Планка для солнечного излучения; Ще, Ща, дт - коэффициенты для одной капли, рассчитанные по теории Ми.
Таким образом, предложенный фрагмент математической модели образования и распространения загрязнений с учетом воздействия солнечного и теплового инфракрасного излучений позволяет принять в
расчет лучистый теплообмен и более корректно учитывать фотохимические реакции, что особенно важно при моделировании вторичного загрязнения.
Список литературы
1. Алексеев Б.В., Гришин А.М. Физическая газодинамика реагирующих сред. - М.: Высш. шк., 1985.
2. Роуч П. Вычислительная гидродинамика. - М.: Мир, 1980.
3. Chen Y., Liou K.N., Gu Y. An efficient diffusion approximation for 3D radiative transfer parameterization: application to cloudy atmospheres // J. Quant. Spectrosc. Radiat. Transfer. -2005. - № 92. - P.189-200.
4. Hu Y. A Study of the Link between Cloud Microphysics and Climate Change: PhD Thesis. - University of Alaska, Fairbanks, 1994.
5. Ou S.C., Liou K.N. Numerical experiments on the Helmholtz equation derived from the solar radiation transfer equation in three-dimensional space // Appl. Math. Comput. -1980. - № 6. - P. 155-175.
6. Segelstein, D.J. The complex refractive index of water: M.S.Thesis. - University of Missouri, Kansas City, 1981.
Пекунов Владимир Викторович,
ГОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина», кандидат технических наук, доцент кафедры высокопроизводительных вычислительных систем, телефон (4932) 26-98-29.
