CC BY
16
УДК 519.5, 519.6
КАЧЕСТВЕННАЯ МОДЕЛЬ ОЦЕНКИ КОНЦЕНТРАЦИИ АЭРОЗОЛЬНЫХ ПРИМЕСЕЙ В АТМОСФЕРЕ, ОСНОВАННАЯ НА ИНТЕГРАЛЬНОМ ПРЕДСТАВЛЕНИИ РЕШЕНИЯ УРАВНЕНИЯ ТУРБУЛЕНТНОЙ ДИФФУЗИИ
© 2012 г. В.И. Наац, Е.П. Ярцева
Ставропольский государственный университет, Stavropol State University,
ул. Пушкина, д.1, Ставрополь, 355009 Pushkin St., 1, Stavropol, 355009
Выполняется построение качественной модели оценки концентрации аэрозольных примесей в атмосфере и соответствующей параметризованной интегральной модели на основе фундаментальной системы решений обыкновенного дифференциального уравнения второго порядка по пространственной переменной. Разрабатывается рекурсивный вычислительный алгоритм, приводятся результаты расчетов и численных исследований.
Ключевые слова: математические модели теории переноса, вычислительный алгоритм, численные исследования.
Construction of qualitative model of the estimation of aerosol impurity concentration in atmosphere and the corresponding parametrized integrated model on the basis of fundamental system of decisions of the ordinary differential equation of the second order on a spatial variable is carried out. The recursive computing algorithm is developed, results of calculations and numerical researches are given.
Keywords: mathematical models of the theory of carrying over, computing algorithm, numerical researches.
Постановка задачи
Оценка концентрации загрязняющих аэрозольных примесей, распространяющихся в атмосфере в пределах ее пограничного слоя, является актуальной задачей в проблеме мониторинга экологического состояния воздушного бассейна и осуществляется в настоящей работе путем численного решения уравнения переноса субстанции, основанного на полуэмпирическом уравнении турбулентной диффузии, которое для случая одной пространственной переменной имеет вид [1]:
д д — q( x, t) + a(t)q( x, t) + — (V (x, t)q(x, t))-
dt dx
д ( д Л
--1 K(x, t)—q(x, t) 1 = S(x, t),
dx ^ dx J
q(x,t ) = q (x) , q(x ,t) = q (t) ,
0 0 0 0
(1)
q(X, t) = q (t) , t e
X
t ,t 0
x , X 0
(2)
В уравнении (1) q(x, I) - концентрация примесей, имеющихся в точке пространства х в момент времени /, еа/1 3 ; а({) - коэффициент, характеризующий степень вывода или привнесения примесей в данный объем за счет химических или других процессов, протекающих в пограничном слое атмосферы, [1/с]; V(х, /) - скорость ветра, 1 /п; К(х,1) - коэффициент
турбулентной диффузии, 1 2 /п, перенос осуществляется вдоль координатной оси Ох; 5(х, /) - источник
примесей, имеет размерность еа/1 3п. Уравнение (1) с начальными и граничными условиями (2) решается относительно распределения q(x, /), остальные функции являются исходными измерительными данными.
Построение параметризованной модели (1), (2)
Важной проблемой при решении задачи (1), (2) является ее обеспечение исходными данными. При этом возникают определенные трудности. Первое, что обращает на себя внимание, состоит в том, что если поле скорости ветра V(х, /) в принципе может быть определено (измерено) в некоторых точках
х , ^ |!> , то поле коэффициента турбулентной
^ ' 1 ШУ.И
диффузии К(х,€) следует считать неизвестным, поскольку оно не поддается измерениям. Значения функции К(х,€) в лучшем случае могут быть определены с использованием полуэмпирических формул, которые связывают в среднем наблюдаемые пространственно-временные вариации компонент скорости ветра со средними значениями К(х,€) в пределах
локальных объемов исследуемой среды. Вторая проблема состоит в том, что в уравнение (1) помимо V(x, t) и K(x,t) входят также и их частные производные (V(x, t)) x и (K(x, t)) x . Что касается определения (V (x, t)) x, то в принципе эта задача решается при использовании соответствующих расчетных методик, если V(x, t) представлена множеством приближенных данных ¡V(x , t )|. Оценка (K(x, t)) x
заметно сложнее, поскольку операция дифференцирования полуэмпирических формул почти бессодержательна в соответствии с теорией аппроксимации функций [2]. Указанные выше трудности в задачах моделирования нестационарного процесса диффузного переноса часто преодолеваются в рамках качественного подхода, основанного на введении так называемых турбулентных состояний пограничного слоя [3]. С формальной точки зрения речь идет об ограничениях на распределения V (x, t) и K (x,t) типа
V(l) < V(x, t) < V(l), K(l) < K(x, t) < K(/), l = 1,2, (3) 1 2 1 2
где V^), v(), k(), K/) - некоторые числа, определяющие границы изменения полей V(x, t), K(x,t) и зависящие от типа состояния (индекса l ). В пределах указанных границ можно говорить о существовании определенной функциональной зависимости между V(x, t) и K (x,t) , представленной той или иной полуэмпирической формулой, в которую, естественно, входят и физические характеристики пограничного слоя атмосферы. Ясно, что говорить о неких средних значениях V и K можно в пределах определенного типа состояния пограничного слоя атмосферы.
Неопределенность исходной задачи требует построения соответствующей адекватной параметризованной модели данного физического явления. Обратимся непосредственно к технике ее построения на основе уравнения переноса (1). Первый шаг преобразования уравнения (1) должен начинаться с нормирования независимых переменных x и t . С этой целью вводим новые переменные x и t :
x(x) = x + (X - x )x , t(t) = t + (T -1 )t , (4)
v ' 0 0 w 0 0
x e [0,l], t e [0,l].
В соответствии с изложенным выше подходом, определяющим ограничения типа (3) для V(x, t) и K (x,t), записываем неравенства V < V(x,t) < V
min max
и K < K(x,t) < K , которые далее позволяют
min max
ввести нормированные значения указанных функций,
V(x, t) = V(x, t)/V , K(x, t) = K(x, t)lK . (5) / max / max
x e
Ясно, что V (x, t) и K (x,t) удовлетворяют условиям
а < V (x, t) < 1 и
aV = Vminl Vm
aK = Kmin/ Km
а < К(x, t) < 1, где К
0<aV<1,
0 <а < 1. Если V(x, t) в области Q = [0,X]x[0,T]
к
(xo = 0 и t^ = 0) является знакопеременной, то ограничения на V (x, t) указанного выше типа должны касаться модуля этой функции. Для определенности решения (1) требуется задание начального распределения q(x,0),
которое выше обозначено как q (x). Определим точку
* * *
0 < x < X условием max q (x) = q(x ) = q . Ис-
q 0<x<X 0 q
*
пользуя найденную константу q , перейдем к функции
*
q(t, x) = q(x(x), t(t)) / q . (6)
*
Аналогично определим точку 0 < x < X условием max S(x,0) = S(x*) = S*. Используя найденную
0<x< X S
константу S , перейдем к функции
~ *
S(t, x) = S(x(x), t(t))/ S . (7)
Подставляя нормированные значения полей (4)-(7) в (1), приходим к параметризованному уравнению вида
x) +a(t) ■ q(t, x) + (V(t, x) ■ q(t, x)) -
,-el\ К(t,x)W-s(t,x),
cx ^ cx J
в котором нормировочные коэффициенты
(8)
^ п V T п KT „ s T ...
a(t) = T - a(t) , ß = —, e= —, S = —, (9)
X
X q
где V = V , K = K . Все коэффициенты, вхо-
max max
дящие в (8), - безразмерные величины. Функции q(t, x), V (t, x), K (t, x), S (t, x) и переменные x, t принимают значения из интервала [0,l].
Относительно параметризованного уравнения (8) сделаем несколько замечаний. Выражения (9) формально вводят 4 безразмерных параметра: а , р, в, £. Если заданы начальные и граничные условия для q(t, x) на Q = [0,l]x [0,l] и фиксированы значения а , р, в, £, то решение (8) следует рассматривать как функцию этих параметров, т.е. q(x,t ,а,р,в,£). Если
V = V и K = K определены однозначно, то
max max
параметры р ив определяют X и T с помощью выражений
* / * * 2 If *Л2
X = к ■р V ■в , T = к ■р \V I ■в. (10)
В результате можно утверждать, что решение , а,р, 9, £,) соответствует множеству ситуаций, ка-
ждая из которых определяется набором ^X ,Т,¥ , К J .
*
Вопрос об оценке вероятных значений параметров q ,
* * *
V , К и £ подробно обсуждается в работе [1].
Построение вычислительной схемы для параметризованного уравнения (8) на основе интегрального представления его решения
В [1] подробно излагается и исследуется так называемый метод интегральных уравнений, применяемый для построения вычислительной схемы параметризованного уравнения (8). Он основан на предварительном интегрировании дифференциального уравнения в частных производных по одной из переменных х или t и сведению его в конечном итоге к интегро-дифференциальному уравнению. При таком преобразовании математической модели явления переноса определяющую роль в ней играет интегральный оператор Вольтерра 2-го рода. Для определения приближенного решения выполняется построение итерационных алгоритмов на основе метода последовательных приближений. Помимо этих интегральных схем решения уравнения переноса могут быть построены и другие формы интегральных уравнений и соответствующие им параметризованные интегральные модели.
В данной работе осуществляется дальнейшее развитие метода инте гральных уравнений и строится параметризованная интегральная модель на основе фундаментальной системы решений обыкновенного дифференциального уравнения 2-го порядка по пространственной переменной. Предлагаемый ниже вариант будет носить сугубо качественный характер. Его основное назначение состоит в разработке метода качественной оценки значений параметров интегральных моделей в задачах математического моделирования явления переноса загрязняющих примесей в турбулентной атмосфере.
Рассмотрим основные этапы построения вычислительной модели для уравнения (8). Для этого перепишем его в виде
t, e-K'-ß-V t a + ß-V
1 + - "
e-K
e-K
q-S- —
e-K
(11)
где q = (q(x, / )) x, q" = (q(x, t)) xx , q = (q(x, / )) t ,
V' = V{x, t)) x, K" = (к(x, / )) x. Для дальнейшего анализа удобно ввести следующие обозначения:
p (x, t )=($■ K'-р■V Ув^ K, p (x, / ) = - (а + р ■ V'^O^ K, f(x, f)=(q-£■ S)lo■ K, в соответствии с которыми уравнение (11) примет вид
q" + p ■ q' + p2 ■ q = f(x,t). (12)
Если считать, что может быть указан способ приближенной оценки величин^1 (g(x, t)) t , а переменная t играет роль фиксированного параметра, то (12) можно рассматривать с определенной долей приближения как обыкновенное дифференциальное уравнение 2-го порядка относительно распределения §(x|?). Последнее
v
обстоятельство позволяет, используя стандартную процедуру построения общего решения подобных уравнений, строить интегральные представления для решений
%(х|?), в которые войдут указанные выше параметры. В
соответствии с теорией дифференциальных уравнений [2], решение (12) представляется в виде совокупности фундаментального решения однородного дифференциального уравнения, соответствующего (12), и частного решения неоднородного уравнения (12).
Ограничимся далее частным случаем коэффициентов р и р , считая их константами. Для этого до-1 2
пустим, что V(х|/)= 1, К(х|/ )= 1 Ух е[0;1] и УТ е[0;1]. В результате получим р^ = -р/в, р^ = -а/в. Фундаментальные решения однородного уравнения, соответствующего (12), примут вид
и (x|F) = e 1 и u (x|F) =
r x „ 2
где гиг- корни ха-1 2
рактеристического уравнения r + pr + p^ = 0:
r =P+ 1 26
pL+a, r =ß— P+a. 2 6 2 26 у 46 6
46
r (X—X') r (X—X')
K(X,X') = e 2 — e 1 , (X'<X),
: e 1
(14)
6 r — r 2 1
затем интеграл (13) можно и в тех случаях, когда р^ и р не являются константами. В частности, эти коэффициенты могут быть полиномами либо другими элементарными функциями. Важно при этом отметить, что если
функции V (Т, х) и К (Т, х) определены табличными зна-
В силу положительности параметров модели а , Р, в, Е корни г иг - вещественные. Частное решение неоднородного уравнения (12) определяется с помощью метода вариации постоянных и после соответствующих выкладок получает представление
~(х|?) = 2 •х ~(х, х') • %(х'|!}&' + <р(х'\7), (13)
0
где обозначено
(р(х'\7)=q(o| 7)- л • е •х ~(х, х') • 5(х'| 7 ,
чениями в узлах |х , Т |, то, используя методы численного дифференцирования, можно рассчитать значения выражений р(х,7) = (в■ К'-р^/в^ К и
р (х,7)=-(а + Р^')/в■ К в указанных узлах, а затем для
них построить аппроксимационные многочлены. Подобное направление расчетно-аналитических исследований также следует считать качественным подходом в теории турбулентного переноса субстанции применительно к задачам экологии.
Построение вычислительного алгоритма. Результаты вычислений
Построение вычислительного алгоритма начинается с дискретизации задачи (13), (14). Зададим сетку
узлов х ( на интервале [0,1]: х = I •Ах , Ах = 1/ т,
О I
I = 0, т . Определим сетку |7 | на интервале [0,1]: Т = j •АТ , АТ = 1/ п , у = 0, п . Введем в рассмотрение
сеточные функции q = х , Т |, V = V\х , Т
У \ г у) у \ г У
К = К\ х , Т I, 5 = 51 х , Т I. Допустим, что
У \ г у) у \ г У
q (x|? )s
q\ x ,t I — q\ x ,t
1 i j) V i j—1
/ —7 .
(15)
Функция ~(х|7), являясь частным решением неоднородного уравнения (12), может быть принята в качестве искомого распределения, связанного с решением исходной задачи (8). Уравнение (13) является интегральным, ядро которого К~(х, х') не зависит от временной переменной Т . Это определяется предположением о том, что V(х|7) = 1, К(х|7) = 1, временная
изменчивость распределения ~(х|7) - лишь временной изменчивостью источника. В этом случае интегральная параметризованная модель (13), (14) носит сугубо качественный характер.
Интегральная модель (13), (14) получена в предположении, что коэффициенты р и р являются константами по крайней мере по переменной х . Вместе с тем подобное предположение не является строго обязательным. Определить фундаментальную систему решений однородного уравнения, соответствующего (12), и построить
Поставим в соответствие интегралам в выражениях (13) квадратурные формулы с коэффициентами ](о |, г = 0, т . Тогда вычислительный алгоритм с учетом (15) получит представление:
q = ßE а К ■ I ~ — ~ \ + p
i,j k=0 k i,k V k'j k'j—1) i'j i ~
(p = q — 7 E а KS ,
i,j 0,j k=0 k i,k k,j
z ч r IX — X I r IX — X T> " " I 2 V i k) 1V i k
К = K\ x , x I = e — e
i,k V ^ k)
ß = k—X , 77 = E ■ AX , X = —--1-,
At 6 r —r
(16)
2 1
r =Z + , 1 26 \
P2 a ß
+ — , r =■
46
2 6 ' 2 26 \
P +a, )
46
26
г = 1, т , у = 1, п .
Данный алгоритм реализован программно при следующих значениях исходных данных: Vo = 101 /п;
К = 100 м/с2 ; % = 0,001 кг/м3; 5 = 0,0001 кг/м3 • с .
Значение V определяется как среднее значение на
интервале
V ,У
тт тах
аналогично задаются K , q
0 0
и £ на соответствующих интервалах. С помощью определенного алгоритма, подробно описанного в [3], генерируются массивы дискретных данных - д
V > , ^г* и - £ !■ . Далее выполняется
У1 тхп ^ У * тхп ^ У * тхп
нормирование заданных распределений в соответствии с процедурой, описанной выше, в результате чего
получаем массивы - д
у
V
K
и
S
элементы которых принимают значения из
* * K и S
интервала [0,1], а также значения V , д Затем задаются коэффициенты а= 0,1, Р = 1, в = 1, т = 15, п = 20, вычисляются Ат = 0,067, А/ = 0,05 , £ = 1,931, Я = 0,845 , / = 1,127, ] = 0,109 , г = -0,092 ,
г = 1,092 2
X = 10 i
и определяются значения
Т = 5п в соответствии с (9), (10). Ниже приведены результаты расчетов для рекурсивной вычислитель-
ной схемы (16) в виде массивов - д =-д
у
i = 1, m
У = 1, п с соответствующим шагом, который представляет собой точное решение уравнения переноса (8) в нормированном виде:
q =
' 0,469 0,524 0,562 0,571 0,545 0,4964
0,556 0,621 0,666 0,675 0,646 0,588
0,637 0,712 0,763 0,774 0,741 0,674
0,711 0,794 0,852 0,864 0,826 0,752
0,774 0,864 0,927 0,939 0,899 0,818
v 0,823 0,919 0,986 1,000 0,956 0,871,
и массив q = <q
i, у
i = 1, m, j = 1, n, полученный в
результате расчетов, проводимых по рекурсивному
алгоритму (16):
' 0,469 0,524 0,562 0,571 0,545 0,4964
0,556 0,565 0,598 0,649 0,761 0,588
_ 0,637 0,646 0,677 0,737 0,831 0,674
q = 0,711 0,719 0,747 0,801 0,886 0,752
0,774 0,781 0,805 0,846 0,907 0,818
v 0,823 0,829 0,846 0,876 0,926 0,871,
Графическое представление данных распределений представлено на рисунке.
Отклонение приближенного решения от точного оп-
1
m n
■sz
toch
i, j
i, j
ределялось по формуле а = -
п • т1=1у=11
=0,082. Полученный результат можно считать вполне приемлемым для данного качественного метода и со-
ответствующего рекурсивного алгоритма. Приближенное качественное решение ~ можно далее использовать как исходные данные для других более точных алгоритмов, описанных в [1]. Ниже приводятся результаты численного исследования алгоритма (16) (табл. 1).
Пространственно-временное распределение в нормированном виде: а - точного решения д ; б - приближенного
toch
решения д
Таблица 1
Влияние размерности задачи на сходимость вычислительного процесса ( Я = 0,845 )
m n M Л а
5 10 1,6903 0,3251 0,0464
10 10 0,8452 0,1625 0,0593
10 20 1,6903 0,1625 0,0726
15 20 1,1269 0,1088 0,0818
Анализируя данные табл. 1, можно сказать, что увеличение размерности задачи приводит к увеличению объема вычислений и, как следствие, накоплению вычислительной погрешности. Приемлемой можно считать размерность 10 х 10 .
у
mxn
j
j
mxn
j
б
Кроме сходимости, исследовалась устойчивость алгоритма к возможным погрешностям в исходных данных уравнения (8), в вычислительном алгоритме -по правой части уравнения (8) [2]. Методика исследования подробно описана в [3], и ее суть состоит в следующем. Предполагается, что вместо точного значения правой части исходного уравнения (8) £^ известно его приближенное значение с погрешностью
8, 8 =
т.е. имеем £
(8)
Величина погрешности
s(°) - S(8)
где S8 = S(0\l + ©е), © - слу-
чайные числа, равномерно распределенные на интервале [-1,1], 0<е< 1. Здесь е имеет смысл процента вносимой погрешности. Вводится коэффициент усиления ошибки 7~, характеризующий степень отклонения получаемого приближенного решения
~(8) [ S(8)
~ (0)
от точного значения q , который опре
делим следующим образом:
с =
8 м
ным погрешностям в исходных данных. При этом максимально возможной погрешностью, вносимой в правую часть (8), можно считать погрешность порядка 10 % и не выше.
Таблица 2
Результаты исследования устойчивости вычислительного алгоритма при т = 15 , п = 20
£, % 8 Л
1,0 0,0016 0,015402 0,01055
2,5 0,00041 0,015405 0,02637
5,0 0,00081 0,015414 0,05271
7,5 0,00122 0,015429 0,07899
10,0 0,00163 0,015450 0,10518
12,5 0,00203 0,015476 0,13125
7 = а / 8. Результаты вычислительного эксперимента представлены в табл. 2.
Анализируя данные табл. 2, можно сделать вывод об устойчивости исследуемого алгоритма к возмож-
Литература
Наац В.И., Наац И.Э. Математические модели и численные методы в задачах экологического мониторинга атмосферы. М., 2010. 328 с. Марчук Г.И. Методы вычислительной математики. М., 1980. 352 с.
Лайтхман Д.Л. Физика пограничного слоя атмосферы. Л., 1970. 366 с.
Поступила в редакцию
11 мая 2011 г.
2
2
