Численное моделирование пространственных задач переноса субстанции в пограничном слое атмосферы Текст научной статьи по специальности «Математика»

ПРОСТРАНСТВЕННЫХ ЗАДАЧ ПЕРЕНОСА СУБСТАНЦИИ В ПОГРАНИЧНОМ СЛОЕ АТМОСФЕРЫ

© 2005 г. В.И. Наац

The construction of double-measured model of source transferring in the frontier layer of the atmosphere, which is based on coordinate wise splitting method, to carry out. The analytical model and results of numerical research calculation algorithm are receiving.

К выбору двумерных моделей теории переноса субстанции в пограничном слое атмосферы. Прогнозирование динамики распространения загрязняющих веществ в пограничном слое атмосферы требует решения нестационарных трехмерных задач теории массопереноса в турбулентных средах. Что же касается двумерных моделей переноса, то они являются приближениями к реальным ситуациям. Стандартным приемом, связанным с переходом от трехмерных пространственных к двумерным задачам, является фиксирование одной из пространственных координат как некоторого параметра задачи с последующим решением соответствующей двумерной задачи. Подобный прием можно условно назвать методом сечений. Вместе с тем определенный интерес представляют двумерные модели, при построении которых в той или иной степени учитывается трехмерный оператор рассматриваемой задачи. К подобным моделям в теории переноса можно отнести те, которые получаются из исходных трехмерных моделей усреднением по выбранной пространственной координате. Рассмотрим кратко технику построения двумерной модели переноса субстанции в пограничном слое атмосферы и дадим анализ ее содержательного смысла [1].

Уравнение переноса запишем в виде:

МАТЕМАТИКА, МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

УДК 519.5

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ

dq( x, y, z, t)

+ a(t) q( x, y, Z, t) +—(Vx (X, y, z, t )q( x, y, z, t)) +

d_

dx

dt

dt dx

В этом уравнении величина ш& характеризует скорость оседания частиц переносимого вещества (речь идет о пылевых частицах, переносимых потоками воздуха) в поле силы тяжести. Значение является осреднен-ным значением скорости по всему спектру размеров частиц, взвешенных в воздухе. Считаем далее, что координата г меняется в пределах интервала [21, 22]. В подобных задачах, когда процесс распространения вещества сопровождается непрерывным оседанием частиц, крупных, прежде всего, в поле силы тяжести, вводится так называемое понятие «подстилающая поверхность». Будем ниже считать ее уравнением г = 21, и для функции д(х, у, г, /) принимать в соответствии с этим граничное условие дд( х, у, г, /) |

dz

=z =uq(x yz t )| z=^ (2)

где и - некоторая константа, характеризующая взаимодеиствие переносимой субстанции с подстилающей поверхностью [2]. Считаем Р(х, у, г) точкой в области ^ = [0,Х] х [0, У] х [21, 22]. Для простоты считаем, что д(Р, () = = 0 при Р еО, исключая г = 21. В качестве как нетрудно видеть, выбран брус, и г = является его нижней гранью. Проведём интегрирование выражения (1) по координате г в пределах [21, 22]. Процесс интегрирования рассматриваем почленно. В частности для первого члена слева имеем

22 дд(х, у, г,I) , д( ^)

]-д-йг =Я(x, У, 1X

21 д

где введено обозначение

_ 22

д (х, у, /) =| д(х, у, г, ()йг. (3)

21

Ясно, что осредненное распределение (3) является функцией двух пространственных переменных и следовательно процесс интегрирования необходимо выполнить таким образом, чтобы в результате прийти к дифференциальному уравнению относительно этой функции.

Интегрирование второго члена а(/)д(х, у, г, /) не вызывает затруднений в указанном только что смысле, и поэтому обратимся к интегрированию членов уравнения (1), содержащих компоненты скорости. Для первого из

22 д д этих интегралов примем |—(Ух(х,у,г,/)д(х,у,г,/)\йг =—(Ух(х,у,/) х

2 дх дх

х д(х, у, /)), что справедливо только в том случае, если компонента Ух не зависит от г, т.е. полагается, что в исходной задаче Ух = Ух (х, у, /) в пределах области Делая аналогичное допущение относительно компоненты ¥у, получим

22 д д

| — ((у (х, у, г, () д( х, у, г, ()) с!г = — (Уу (х, у, ()д (х, у, ()).

z

В связи со сделанными выше допущениями относительно Vx и Vy можно допустить предположение, что Vz = 0 в пределах [Zb Z2] и перенос вещества в области Q в основном определяется процессом его оседания.

_ Z2 dq(x, y, z, t) , . ,\z2

Остается рассмотреть интеграл - J шg ——--dz = -rngq(x, y, z, t) z2 =

zi dz

= Шgq( x, y, Zi, t).

При его вычислении учитывались условия = const и q(x, y, Z2, t) = 0 (граничное условие). В результате левая часть (1) принимает вид

q( ^ y, t) + a(t )q ( x, y, t) + dx (Vx ( x, y, t ж x y, 0 ) + ) (Vy ( x, У, t )1( x, У, t) ) + Ш( x, У, Z1, t),

где обозначено q( x, y, Z1, t) = q( x, y, Z1, t).

Остается рассмотреть преобразование правой части (1) при интегрировании его по переменной z в пределах [Z1, Z2]. Здесь принимаются те же допущения, что и выше. Считаем, что Kx и Ky зависят только от x, y и t и не зависят от z. В этом случае

Z2 д( dq(x, y, z, t) V d ( ,dq(x, y, t) 1

J dx Kx (x, y, t) dx I dz = т- Kx (x, y, t) 4 кду> I, Z dx ^ dx J dx ^ dx J

ZJ2 ± ( Ky (x, y, t) 1 dz = ± ( Ky (x, y, t)

z1 dy 1 dy J dy { dy

Z21 (Kz < x, У t) dz = ( K (x. * z-1) ^ )| Z2 = (4)

^ , ^ s dq(x, y,z,t)i , ^ ч dq(x, y,z,t)i

= Kz (x, y, Z2, t) dz \z=Z2 - Kz (x, y, Zi, t) dz z=Z1 .

dq( x, y, Z 2, t)

Поскольку q(x, y, Z2, t) = 0, то резонно считать, что и-2— = 0,

dz

... ^ _ dq(x,y,Zi,t) и первый член в (4) будет отсутствовать. Для величины --— в

dz

соответствии (2) следует принять значение uq(x, y, Z1, t).

_ Z2

Аналогично вводя S (x, y, t) =J S(x, y, z, t)dz, окончательно имеем

Zi

уравнение относительно распределения q (x, y, t):

dq (XдtУ, t) + a(t)q (x, y, t) + dx (Vx (x, y, t)q (x, y, t)) + -)y ( (x, y, t)q (x, y, t)) =

= A( Kx (x, y, t) ^Ml+Af Ky (x, y, t) dx ^ dx J dy ^ dy

-(ш% + кг (х, у, г!, Ои ж х, у, 2Х, г) + £ (х, у, г). (5)

Полученное уравнение решает поставленную выше задачу, а именно, строится нестационарная модель переноса, которая, будучи двумерной, давала бы в среднем представление о поле концентрации веществ в пределах трехмерного объема. Совокупность тех допущений о физических характеристиках среды, которые были сделаны выше, дает представление о содержательности модели (5). Если далее предположить, что выполняется условие для компонент скорости ветра, определяемое уравнением неразрывности

дУх (х, у, г) + дУу(х ^ г) = 0 дх ду

то уравнение (5) может быть переписано таким образом:

^y,о + а(г)_(х,у,г) + 11 Ух(х,у,г) Ад(х,у,г) | -

дг [У дх

кх(х,у,г) №МТ1 + {ГУу(х,у,г)Аж , дх У дх ) ду ) (6)

-А (к, (х, у, г) Я!|М I-«*.,.,)-

-(ш& + кг (х, у, г!, г)и )д( х, у, г!, г).

Алгоритм расщепления для уравнения (6) построим следующим образом (метод покоординатного расщепления подробно рассмотрен в [3, 4]):

1. г< г < г+1, у = 0, п,

х, у),} = 0

2. qx{x,y,tj) = •

; |?2(x y, tjX j > 0

dq, (x, y, t) 1 d

1 ' (t)q1 (x,y,t)+Vx(x,y,t)—q,(x,y,t)-

dt 2 dx

-д- ^ Kx (x, y, t) ) + ^ + Kz (x, y, Z,, t)u )q( x, y, Z,, t) =

= S,( x, y, t), (7)

x, y, t) =®,S(x, y, t), q1 (0, y, t) = q (0, y, t), q, (X, y, t) = q4(X, y, t);

3. q2( x, y, tj) = q,( x, y, tj+,),

dq2 (x, y, t) 1 d

2 ' + — a(t)q2 (x, y,t) + Vy (x, y, t)—q2 (x, y, t) -dt 2 dy

-д-( Ky(х,y,t) dq2{X:y t) | + Urng + Kz(x,y,Z,,t)q(x,y,Z,,1) = dy ^ dy ) 2 (8)

= S2( x, y, t), S2(x, y, t) — ®2 S (x, y, t), q2 (x, 0, t) = q5 (x, 0, t), q2 (x, 7,t) = q6 (x, 7,t); ®1 + ®2 = l,

4. j = j + 1; _ _

5. если j < n, то перейти на шаг 2, иначе q (x, y, tj) = q2 (x, y, tj+1).

Алгоритмизация вычислительной схемы, результаты численных исследований. Для того чтобы определить значения исходных данных для последующих расчетов, необходимо создать некий тестовый пример, условно называемый «Блоком исходных данных». Для этого предположим, что поле концентрации примесей определяется функцией вида q(x, y, z, t) =

= q0 • q(x, y, z, t), в которой t, x, y и z - нормированные величины, значения которых меняются в интервале [0, 1] и функция q (x, y, z, t) задается

Л Л Л Л Л Л Л Л

формулой q(x, y, z,t) = (1 + sinx)(1 + siny)(1 + sin z)(1 + sint).

Тогда осредненное распределение концентрации можно определить следующим образом:

Z2 z2 Л Л Л Л Л ЛЛЛ

q(x,y, t) — J q(x,y, z, t)dz и q0H • J q(x,y, z, t)d z — q0 • H • azq(x,y, t),

z, л

1 zi

íH — Z 2 Z,, z, — ^^^, z2 — ^^^,

2 11 H 2 H

az — (z2 - z,) - (sin z2 - sin z,),

q(x, y, t) — (1 + sin x)(1 + sin y)(1 + sin t). (9)

Аналогичным способом зададим значения остальных распределений:

Vx (t, x, y) и VxV (Л, Л, y) — Vo

Vy (t, x, y) * Vy 0Vy (t, x,y) = Vy 0 Kx (t, x, y) * Kx,0Kx (t, x, y) = Kx,0 Ky (t, x, y) * Ky,0 Ky (t, x, y) = Ky,0

(1 + sin t)(1 + cos x)(1 + sin y) (1 + sin t)(1 + sin x)(1 + cos y) (1 + sin t)(1 + sin x)(1 + cos y) (1 + sin t)(1 + sin x)(1 + sin y)

Kz (t, x, y) » KzoKz (t, x, y) = Kz,

(1 + sin t )(1 + sin x)(1 + cos y )(1 + sin Z1)

qz(t, x y) ~ ?0 qz(t > x y) = ?0

(1 + sin t )(1 + sin x)(1 + sin y )(1 + sin z1)

(10)

S (t, x, y) = f S(t, x, y, z)dz » S0H - f S(t, x, y, z)d z = S0 -H -sz - S(t, x, y),

Z A, (ID

sz = (z2 - z1) + (cos z2 - cos z1).

Подставим записанные выше выражения (9), (10) в уравнения (7) и (8). Получим

qp-н-д2 dq1(рt) +1 л q тт g q (P -t)+V V (P -t) q0-н-дг v

-----:-+ -'a(tГq0 -H'gz t) + V0,x'Vx (P1, t)----v

T d t 2 ^

^P, t) q0 - H 'az ■ Kx,0 dKx (P, t) dql(P[, t) ^-H-a,-Kx

X 2

d x

d x

X2

d x

2~ A

x^^^M + \-rng-q0-qz(P,t) + ±-Kzß U q-Кг(P,t)-qz(P,t) =

d x

= S0-H-Sz-S1( P, A);

q0-ü-az dq2(P2,t) + 1 -^(t)-q0 .H.arq2(P2,t) + ^ VPt)х

■Kx ( P, t )>

Г

51

_q0 -H-az xdq2(P2, t) q0 -H-az'Ky,0 dKy (P , t) t)

1 dy ^

q0 -H-az-Ky,0 x K (P a) d2q2(p2,t) 1

5 y 5 y

(12)

12

K (P2, t )-

5 y

2 + 2- qz(P2, t)+

+--Kz,0 Us q -Kz(P2,t)-qz(P2,t) = S0-H-s2 ^(P,,t),

в которых P = P( x

y); P2 = P(y

x); P = P(x, y); S1 (P1, t) + S2 (P2,t) = S(P, t),

S1(P1,t) =a1 -S(P,t), S2(P2,t) =®2 -S(P,t), m1 + m2 = 1, 0 < m, < 1, i = 1, 2. Введем обозначения для нормировочных коэффициентов a(t) =

1 А ™ V0,xT_ V0,yr. _K0,x-r _K0,y'r.

=~7-a(t)-г; ßx =; ßy ; r =

X

X

; r =■

2 ' ' y

Y 2

/I =

z

2

T • и ; M = ; # = So T •

2 • Н • аг 2 • Н • аг д0 • аг

Дифференцируя (10), получим выражения для производных, входящих в (12), которые обозначим:

дШ,х, у) / Л " дд(г,х, у) / л " да (г,х, у) / л лч » = р((г, х, у); чк\:у= Рх(г, х, у); ЧК,Г= Ру (г, х, у);

д г д х д у

д2Ш, х, у) / л " д2д(г, х, у) / л \

* Л 2 = Рхх(г, х, у); 2 = Руу (', х,у);

д х д у

(л, х, у) « л Л Жу ^ у) л л л

л = Ях(г, х, у); — = £у(г, х,у).

д х д у

Выражения для производных и нормировочных коэффициентов подставим в (12).

р (г, Р) + а (г) • дг(Г, Р) + вх • К (г, Р) • Рх (г, Р) -Гх • Ях(г, Р) • Рх (г, Р) --Гх 6 Р) • Кх 6 Р) • Рхх (л, Р) + ^ (л, Р) +1 • К г (л, Р) • Ъ (/, Р) =#• ^ Р),

Р, (г, Р)+« (/) • д[(г, Р)+ву К (/, Р2) • Ру £ Р2)-Гу • Яу (?, Р2) • Ру £ Р2) -

-/у (г, Р2) • Ку (г, Р2) • Руу (г, Р2) + (г, Р2) + .1 • К г (г, Р2) • дг (г, Р2) = £• §2(г, Р2).

Вычисляем значения правой части каждого из уравнений, т.е. значения функции источника §х(г, Р), §2(г, Р2) и § (г, Р).

В итоге имеем распределения всех исходных данных, причем д (Р, г)

далее считаем точным решением и обозначим дт (Р, г). Дальнейшая реализация вычислительного алгоритма проводится в предположении, что

д (г, Р) неизвестно и должно определяться в ходе вычислений при известных исходных данных. Обозначим его дм (г, Р). Точность вычислений можно оценить значением величины отклонения приближенного решения

дм (г, Р) от дт (Р, г), генерируемого в «Блоке исходных данных», по фор-

муле: с =

дт (t, P) - дм (t, P)

Для проведения расчетов с целью исследования влияния скорости ветра и турбулентности на распределение концентрации примесей необходимо задать функцию источника § (г, Р), начальные и граничные условия

qУ, Р), а также значения распределений Ух (/, Р), Уу (/, Р), Кх (/, Р) и

Ку (Г, Р).

Алгоритмизация вычислительной схемы осуществляется путем замены в соответствующих формулах распределений их сеточными функциями. Простейшим способом решения поставленной задачи может служить применение сеточной модели на основе неявной разностной схемы для двумерного случая с последующим подключением метода прогонки [5].

Описанный выше алгоритм численно реализован при следующих значениях исходных данных: q0 = 0,75 кг/м2, а0 = 0,5 с-1, = 0,05 м/с, и = = 0,01 с-1, Ух,0 = 10 м/с, Уу,0 = 10 м/с, Уг0 = 5 м/с, Кх,0 = 300 м2/с, Ку,0 = 250 м2/с, ' К,0 = 170 м2/с, £0 = 0,25 кг/м3с, X = 10000 м = 10 км, У = 10000 м = 10 км, г! = 50 м, г2 = 100 м, Т = 3600 с = 1 ч, X = 0,5, М = К = 10, N = 200. Ниже приведены результаты расчетов в виде массивов qт (т, к) и qм (т, к), т = 0,М, к = 0, К для фиксированного момента времени t(] = 150) = 0,75, что соответствует реальному значению t = 45 мин = 2700 с. Величина отклонения точного решения qт (т, к) от приближенного qм (т, к) для текущего момента времени составила а(] = 150) = 3,76 Е -03, итоговое значение - а = 2,45 Е -03. Это говорит о вполне удовлетворительной работе вычислительного алгоритма.

На рис. 1, 2 представлено распределение концентрации примесей д (^ Р) для различных значений скорости ветра Ух^, Р), Уу(^ Р) при соответствующих начальных значениях констант Уг,0 и Уу,0, а также поля турбулентности Кх(^ Р), Ку(^ Р) в зависимости от значений констант Кх,0 и Ку,0 при фиксированном значении момента времени t.

На рис. 3-5 представлены линии уровней, соответствующие распределениям концентрации д ^, Р), показанным на рис. 1, 2. С усилением скорости ветра происходит смещение «максимума» к началу системы координат, поскольку процесс переноса ускоряется (рис. 3).

8 б

Ч

4

2

Рис 2. Распределение концентрации q(t,P), P = P(x, y) (t = 45 c; Vx,o = 15 м/с; Vy,o = 15 м/с; V,o = 15 м/с; Kxfi = 150 м2/с; Ky>o = 150 м2/с; Kz_o = 15 м2/с)

8-

G

У

4

2

Рис. 3. Линии уровня, соответствующие распределению концентрации д (г, Р) (рис. 1)

В G

У

42

Рис. 4. Линии уровня, соответствующие распределению концентрации q (t, P)

С усилением турбулентности, кроме того, усиливается процесс размывания профиля Р) (рис. 3-5). Линия уровня - окружность в центре рис. 3 - соответствует значению концентрации примесей дтах (/, Р) = 6,03 . Минимальное значение при этом Р) = 1,32.

Рис. 5. Линии уровня, соответствующие распределению концентрации д ^, Р) (рис. 2)

С увеличением значения коэффициента турбулентности при том же значении скорости ветра (рис. 4) уровень концентрации снижается: Ятш(*, Р) = 0,87 - дшах(/, Р) = 5,3.

При дальнейшем усилении скорости ветра (рис. 5) при том же достаточно высоком значении коэффициента турбулентности уровень концентрации несколько повышается за счет ускорения процесса переноса: дтт (/, Р) = 1,48 - дтах (/, Р) = 5,57, но при этом размывание профиля д (/, Р) усиливается.

( 1,888 1,510 1,369 1,572 1,965 2,251 2,213 1,880 1,504 1,370 1,579 ^ 1,414 1,131 1,025 1,178 1,472 1,686 1,657 1,408 1,127 1,026 1,182 1,023 0,818 0,742 0,852 1,065 1,219 1,199 1,018 0,815 0,742 0,855 0,810 0,648 0,587 0,675 0,843 0,966 0,949 0,807 0,645 0,588 0,677 0,828 0,662 0,600 0,690 0,862 0,987 0,970 0,824 0,660 0,601 0,692 дТ (т, к) = 1,072 0,857 0,778 0,893 1,116 1,278 1,257 1,068 0,854 0,778 0,897 1,482 1,186 1,075 1,235 1,543 1,768 1,738 1,476 1,181 1,076 1,240 1,959 1,567 1,421 1,632 2,040 2,336 2,296 1,951 1,561 1,422 1,638 2,385 1,907 1,730 1,986 2,483 2,843 2,795 2,375 1,900 1,731 1,994 2,655 2,124 1,926 2,212 2,765 3,166 3,113 2,645 2,116 1,927 2,221 2,705 2,163 1,962 2,253 2,816 3,225 3,171 2,694 2,155 1,963 2,262,

qM (m k) =

( 1,888 1,414 1,023 0,810 0,828 1,072 1,482 1,959 2,385 2,655 2,705

1,511 1,131 0,818 0,648 0,662 0,858 1,187 1,568 1,910 2,127 2,167

1,375 1,029 0,743 0,589 0,602 0,779 1,079 1,427 1,739 1,937 1,974

1,577 1,180 0,853 0,676 0,691 0,895 1,238 1,637 1,995 2,222 2,264

1,964 1,471 1,064 0,843 0,862 1,116 1,543 2,038 2,481 2,763 2,815

2,244 1,681 1,217 0,965 0,986 1,276 1,764 2,328 2,832 3,153 3,212

2,206 1,653 1,196 0,948 0,969 1,255 1,734 2,289 2,784 3,099 3,157

1,877 1,406 1,017 0,806 0,824 1,067 1,474 1,947 2,369 2,638 2,687

1,506 1,127 0,815 0,646 0,660 0,855 1,182 1,563 1,903 2,119 2,159

1,376 1,029 0,744 0,589 0,602 0,780 1,079 1,428 1,740 1,939 1,975

1,579^ 1,182 0,855 0,677 0,692 0,897 1,240 1,638 1,994 2,221 2,262

Литература

Наац В.И. // Тр. Международной школы-семинара по геометрии и анализу памяти Н.В. Ефимова: Тез. докл. Ростов н/Д, 2004.

Берлянд М.Е. Современные проблемы атмосферной диффузии и загрязнения атмосферы. Л., 1975.

Наац В.И. // Сб. научн. тр. У-го Всероссийского симпозиума «Мат. моделирование и компьютерные технологии»: Тез. докл. Кисловодск, 2002. МарчукГ.И. Математическое моделирование в проблеме окружающей среды. М., 1982. Марчук Г.И. Методы вычислительной математики. М., 1980.

Северо-Кавказский государственный технический университет

14 февраля 2005 г.

УДК 517.9

О РАВНОСИЛЬНОЙ РЕГУЛЯРИЗАЦИИ НЕКОТОРЫХ КЛАССОВ ДВУМЕРНЫХ ОПЕРАТОРОВ ТЕПЛИЦА

© 2005 г. А.Э. Пасенчук, Н.А. Еволенко

A two-dimensional Teplits's operator with a symbol that is analytical with respect to each variable is examined in the space. A construction diagram of an equivalent regularizer is shown. The examples of this diagram application are given.

1. Будем пользоваться следующими обозначениями: С - поле комплексных чисел; Г = {t е С: Ы = 1}, Г2 = Г х Г, Ь2(Г), Ь2(Г2) - гильбертовы пространства измеримых, суммируемых с квадратом функций на Г, Г соответственно; Г), L++ (Г2) - подпространства функций, аналитически продолжимых внутрь единичной окружности по каждой переменной в пространствах Ь2(Г), Ь2(Г2); Р+ Р+ - операторы проектирования Ь2(Г),