CC BY
13
УДК 519.5, 519.6
ПЕРЕНОС И КОАГУЛЯЦИЯ ДИСПЕРСНЫХ ЧАСТИЦ В АТМОСФЕРЕ ПРИМЕНИТЕЛЬНО К ЗАДАЧАМ ЭКОЛОГИЧЕСКОГО МОНИТОРИНГА
ВОЗДУШНОЙ СРЕДЫ
© 2010 г. И.Э. Наац1, В.И. Наац1, Р.А. Рыскаленко
2
1Ставропольский государственный университет, ул. Пушкина, д. 1, г. Ставрополь, 355009, info@stavsu.ru 2Ставропольский военный институт связи РВ, ул. Артема, 2, г. Ставрополь, 355017, svisrv@mail.ru
1Stavropol State University, Pushkin St., 1, Stavropol, 355009, info@stavsu.ru 2Stavropol Military Institute of Communications, Artem St., 2, Stavropol, 355017, svisrv@mail.ru
Рассматривается математическая модель переноса дисперсных частиц в атмосфере, в которой учитывается процесс коагуляции частиц. Для данной модели предложен вычислительный алгоритм, в основу которого положен метод расщепления по физическим факторам. Для уравнения коагуляции выполняется построение рекурсивного алгоритма.
Ключевые слова: экологическое состояние атмосферы, пограничный слой атмосферы, математическая модель, перенос дисперсных частиц, коагуляция частиц, вычислительный алгоритм, уравнение.
Mathematical model disperse particles transfer at the atmosphere was take with coagulation process. For it computational algorithm was devised and used splitting method at physical parameters. Further, recursive algorithm was executed for equation of coagulation.
Keywords: ecological states of atmosphere, border air layer, mathematical model, disperse particles transfer, coagulation, computational algorithm, equation.
В настоящей работе выполняется построение вычислительной схемы модели переноса дисперсных частиц в атмосфере на основе метода расщепления. Особое внимание уделяется разработке рекурсивного алгоритма для нелинейного уравнения коагуляции, входящего в состав исходной модели переноса.
Постановка задачи
Математическая модель переноса многокомпонентной примеси в турбулентной атмосфере с учетом коагуляции частиц может быть представлена в операторной форме:
ф + Оср-Аср = 8, (1)
где В , А - дифференциальный и интегральный операторы вида
1 -3 У а(1)ср(Р. I . 3) +
-X — (yt(P,f)<p(P,t,0))-
OXi
cXj
Kt (P, t )
dcp(P,t, 3) Sx,
(2)
l9
= — J Ф(й>, 3- a>)ç(P, t,a>)(p(P, t,3- a>)da> -
2
3
- | Ф{со, 3)<р{Р, со)(р(Р, и 3)йсо . (3)
Здесь Р - 1'(Х\,х2, х3), Р е с Щ - точка про-
странства; /е - временная переменная;
<р(1'. 1.3) - концентрация частиц массой 3 в точке Р
в момент времени г; начальные условия:
ср{Р^,3)| =(р0{Р,3) для РеО; граничные условия: ср(Р,Г,3) = р(Р,Г,3) для Р<еП, где <р{)(Р.9) и <р(Р, !■ 3) - заданные функции; О - граница области О : а([) - коэффициент, характеризующий степень вывода или привнесения примесей в данный объем за счет химических или других процессов, протекающих в пограничном слое атмосферы; ¥±(Р,г), У2(Р,О, У3(Р, г) - компоненты вектора скорости ветра; К (Р, г), К2 (Р, О, Къ (Р, г) - турбулентность, характеризуемая коэффициентом турбулентной диффузии; перенос осуществляется вдоль координатных осей Ох1, Ох3, Ох3; Б(Р, 1.3) - источник примесей (частиц). Уравнение (2) описывает процесс переноса частиц под действием ветра при наличии турбулентности в атмосфере, а также с учетом осаждения крупных частиц под действием силы тяжести; (3) - явление коагуляции частиц в воздушной среде в процессе их движения (турбулентного перемешивания). 1-й интеграл в (3) описывает образование частиц массой 3 за счет соединения их с частицами меньшей массы со и ^ - о , где со - текущая переменная масса, в общем случае со е ^. 32 2-й интеграл в (3) - распад частиц массой & на более мелкие; Ф 4р,3 - частота столкновения частиц с массами 3 и со . Оператор А является нелинейным. Ядро интегрального оператора явно не определено. В общем случае оно может быть задано с помощью некоторой аналитической ее функции, пример приведен в работе [1]. Набор возможных аналитических функций позволит получить семейство | и исследовать в вычисли-
3
тельном эксперименте влияние выбора той или иной функции на форму и временную изменчивость спектра дисперсных частиц <р(Р,/. $).
Построение вычислительного алгоритма модели (1)
В основу вычислительной схемы, разрабатываемой для модели (1)-(3), положена процедура расщепления (1) по физическим факторам [2]. В соответствии с этим подходом полагается, что /. меняется в пределах дос-
таточно малого интервала
I,
7+1
= В соответствии с методом расщепления формулируются две последовательно решаемые задачи. Задача I:
ф1 СР,/, 3) + С<Р1 ~$РА <9) = РА 3),
7=0, РеП
j = 1,2,...
<pl(P,tJ,3) = ^(P,tJ,3), Ре п.
Задача II:
ср2 (Р, t,S)~ 41 ср2 t, 3) = a>2S(P, t, 3).
((p2(P,tJ,3) = cpl(P,tJ+l,3), PeCl \p2(P,tj,S) = p(P,tj3), Pen '
3^<3<32, щ+а>2=\,
(4)
(5)
(6) (7)
tj<t<tj+1
j=l,n.
В результате (р(Р,(,3) = ср2(Р,г,3). Таким образом, изложенный подход позволяет последовательно решать уравнения (4) и (6) в пределах интервала , / -+1 _ И при ЭТОМ функции (Р\ (Р. I . 3) и <р2 (Р. I . 3)
связаны условиями (5) и (7).
Построение вычислительных моделей для решения первой задачи (4), (5) приведено в [3-5]. Ниже рассмотрим способ построения вычислительного алгоритма для 2-й задачи, а именно для уравнения коагуляции (6), (7).
В основе предлагаемого авторами вычислительного алгоритма лежит идея построения рекурсивной схемы по переменной г (метод последовательных слоев), на каждом шаге которого решается соответствующая система линейных алгебраических уравнений относительно
ср2 = $2(Зк , к = \,т. Построение алгоритма выполняется следующим образом. По переменной t определяем последовательность точек I , в пределах
интервала [). '/ . Для каждого фиксированного момента времени I , имеем функцию <р2 (3 | /7.Р) - <р2 (3). Независимые переменные задачи пронормируем следующим образом: & = + (9^, а> = 31+(32-31)-3,
/=/0 + (Т-(0)-7, где 0<,9<1, 0<®<1, 0<Г<1. Подобным же образом могут быть пронормированы и
значения искомого распределения ср2(Р,1,&) (процедуры нормирования приведены в [4, 5]). В итоге уравнение (6) с учетом (3) запишется в виде
2 ^
3
х | Ф(~, ~ - ~) • (р(3 - г )(р(~, г ~ -о
~ 1 ~ ~ - <р(3,Г)|Ф09,3) ■ ср(3,7)й3 + 8(3,7). (8) о
В (8) опущена точка Р для простоты представления формулы, однако при этом имеется в виду, что
<р(3,1) = <р2(Р,1,3).
Введем Я = (Т — ^)(32 -3±) - масштабный множитель. Уравнение (8) решается относительно искомого распределения <р(3^). Предварительно необходимо провести линеаризацию нелинейного уравнения (8) и включить соответствующую процедуру в общую схему вычислительного алгоритма. Одним из возможных вариантов линеаризации (8) может быть следующий. Выберем на множестве три точки
, tj, 1,+\ и будем рассматривать 3 функции
<р}~1(&), (р3 (3), срУ+1 (3). На основе (8) свяжем эту
совокупность уравнением
<р3
\3)-ср 1(~)
2т
3ф (9-3,3)
о
2
1
• <р(3 - 3, t j)
(pJ+1(3)-ср J~ 1(3)
d ~ -
- ÄtpJ (3)\Ф(3,3) о
<pJ+l (3)-<pJ'1(3) 2
d3 + SJ (3),
T ~ O+i " h ■■
7=1,«. (9)
Введем интегральные операторы, соответствующие интегралам в правой части (9). 1-й интеграл соответствует нелинейному интегральному оператору Вольтерра и может быть записан как А(<р)<р. 2-й оператор - линейный интегральный оператор Фредголь-ма с ядром Ф(,9,ет), где (3,3) <е ¡¡Д^ |;1 . Если последний оператор обозначить через В , то исходная модель (8) примет вид ф — АА(ф)ср - /¿р ■ Вер + 5" или
-ср
7-1
2 т
-кр]В Ядро
(,oJ+l -ср]Л 2
2
7 = 1," •
(10)
интегрального
оператора
A
Ф(<9 - З)^ (3 - 3) . Далее обозначим его как
1
_I
2
А3 =А((р3). В этом случае (10) можно записать в виде
<о3+1 -(¿А3 -тА<р3В=
{¡А ' "
= <pj~x
+ \?AJ -TÄtpJB(ßJ~l + 2tSj
(11)
или, обозначив через К3 = ср3В - А3, на основе (11) получим ( + тЖ] = (- тЖ3 + 2тБ3.
2
Теперь можно записать схему прогноза <р3+1 на основе т}~1 и
, а именно
ср}+Х = {+тЖ< (-тЖ1 +
+ 2г( + гЖ-/^15,-/'. (12)
В соответствии с теорией эволюционных моделей
1 т Ii
2 к=0
(15)
Аналогично формируем следующее выражение: <р3{31)^<р3^1) = ~ 1 ~ = (р] Нф(~, ®)<Р1
о
. ~ т ^
(<%) 1Ф (Зиюк)<р-> (юк)ук1г =
к=О
(16)
=т х фик<Ркук^-
к=О
С учетом (15) и (16) выражение (11) после дискретизации получит представление
- тМ х
оператор { + тЖ' ^ { - тЖ ' в (12) обозначается через Т3 и называется оператором шага, а
¡5-3 = ^ + тЯК-3 ^ - источника. Окончательно имеем рекурсивную вычислительную схему вида
<р]+х =Т](р]~х . (13)
Оператор шага Т определяется на С |,1_. Выполнение данных условий обеспечивает сходимость рекурсивной вычислительной схемы (13).
Завершающим этапом в построении вычислительного алгоритма для уравнения (13) является дискретизация данной задачи по переменной ,9 с использованием квадратурных формул. Поскольку 19 е |,1_ и о>(= |д , на интервале |,1_ определим сетку узлов , / = 0,т и , к -0,т . Тогда имеем
С}<р] ^г У - / ф(Д - ®)<Р3 (Д - аур3 {со)йа> » 2 о
\ г ~ . ~
и- &к)уф ■ (14)
2 к=0
В (14) г{к , к = 0,т - квадратурные коэффициенты, // - сок | - о)к. Введем обозначения: ср! = <р ' (31):
<*>/_* = (Р] (Д ®г-к,к = Ф(Д -Щ,юк)\
80г ,щ)= 1 п ¿г,к = ¿(Д = Щ) ■ Тогда
[О,
последнее выражение в (14) можно переписать в виде 1 *
хф(# -Sk,Sk)<pJ(3j -Sk)<pJ(юк)vkh = 2 £=0
1
X X I ^Ф+<РгФ1,к kk<PJk+l к=0 V 2
= (pj 1 + täh x 1
X 2 I ^Ф 1-к,к(р1А,к -^фа М" • (17)
Обозначим В (17) Кг,к=\ Ф1-к,к<Рг-к^г,к ^
-(р'Ф,к- а = тЖ. Тогда окончательно можно записать
J+1
-а X 1 к=0
ЯП
j-1 ^ = щ + а х
к=о
с
i
г,А
>
ЧкУкП
/ 1 + 2Л,' .
(18)
Выражение (18) представляет собой систему линейных алгебраических уравнений относительно
/ = 0, т, у = 1 ,п .
Вычислительный алгоритм (18) предполагает для каждого значения 1 решение системы алгебраических уравнений, в которой матрицы К] и Я3 также рассчитываются для каждого значения 1 . Наиболее предпочтительным может быть метод последовательных приближений. Однако при этом необходимо определить условия его сходимости. Одним из таких условий является выполнение неравенства
KJ
<1 V/', j =\,п , где К3 = <pJВ - А3. Так как
а = тЖ и при этом тик- достаточно малы, то обеспечить выполнение условия а < 1 можно, подбирая соответствующее значение Я в вычислительном
К]
эксперименте. Выполнение условия
< 1 - зада-
Поступила в редакцию
ча значительно сложнее. Здесь требуются дополнительные исследования, в том числе проводимые в вычислительном эксперименте.
Литература
1. Алоян А.Е. Математическое моделирование взаимодей-
ствия газовых примесей и аэрозолей в атмосферных дисперсных системах // Вычислительная математика и математическое моделирование: тр. междунар. конф. М., 2000. Т. 1. С. 214-230.
2. Марчук Г.И. Математическое моделирование в пробле-
ме охраны окружающей среды. М., 1982. 320 с.
3. Наац В.И. Вычислительная модель нестационарного
уравнения переноса примеси на основе метода взвешенной невязки // Изв. вузов. Сев.-Кавк. регион. Ес-теств. науки. 2004. № 5. С. 3-15.
4. Наац В.И., Рыскаленко Р.А. Разностная аппроксимация в
задаче Коши для нестационарного уравнения переноса // Вестн. СевКавГТУ. Серия Физико-химическая. 2004. № 1. С. 93-99.
5. Семенчин Е.А., Наац В.И., Наац И.Э. Математическое
моделирование нестационарного переноса примеси в пограничном слое атмосферы. М., 2003. 291 с.
_21 января 2009 г.
т
т
