CC BY
24
УДК 519.63
ОБРАТНЫЕ ЗАДАЧИ О МОЩНОСТИ ТОЧЕЧНОГО ИСТОЧНИКА В МАТЕМАТИЧЕСКОЙ МОДЕЛИ РАССЕЯНИЯ ПРИМЕСИ В АТМОСФЕРЕ
© 2010 г. Е.А. Семенчин1, М.В. Кузякина1, Е.О. Лоскутова2
1Kuban State University, Stavropolskaya St., 149, Krasnodar, 355040,
1 Кубанский государственный университет, ул. Ставропольская, 149, г. Краснодар, 355040, rector@kubsu.ru
2Черноморская гуманитарная академия, ул. Орджоникидзе, 10А, г. Сочи, 354000, shga@sochi.ru
rector@kubsu.ru
2Black Sea Humanitarian Academy, Ordzhonikidze St., 10A, Sochi, 354000, shga@sochi.ru
Предлагаются 2 новых алгоритма расчета мощности точечного источника непрерывного действия в рамках математической модели рассеяния примеси в атмосфере. Один из них позволяет производить расчеты мощности источника по экспериментально заданным значениям параметров модели и концентрации примеси; другой — производить расчеты мощности источника по тем же значениям параметров модели, но вместо значений концентрации требует использования значений плотности осадка на подстилающей поверхности и скорости сухого осаждения ее частиц. Оба алгоритма реализованы в программных продуктах (соответственно в средах Maple и Delphi). Приведен сравнительный анализ этих алгоритмов.
Ключевые слова: мощность источника, концентрация примеси, плотность осадка, подстилающая поверхность, фильтрация, регуляризация.
In this paper two new algorithms for calculating the power of a continuous action point source are presented in mathematical model of admixture dispersion in the atmosphere. One of the algorithms in question allows performing source power calculations using experimentally given model parameters' values and concentrations of admixture in the atmosphere; the other one allows making source power calculations using the same model parameters' values, but instead of concentration values it requires the values of sediment density on an underlying surface and the speed of dry deposition of admixture particles. Both approaches have been realized in computer programs (using MAPLE and Delphi respectively). Comparative analysis of the algorithms is given.
Keywords: source power, admixture concentration, sediment density, underlying surface, filtration, regularization.
1. Постановка задачи. Математическая модель, описывающая процесс рассеяния примеси в атмосфере, имеет вид [1]
*+и * - w *=± кх *к * +
dz дх dx dy dy
dt
дх
д dq
+—Kz— + f, t e [t o, T ], q(t о, x, y, z) = p(x, y, z) dz dz
{Kzdq + Wq} z
=z0
={^q}| z
=z0
(1)
q (t, x, y, z) ^ 0, x 2 + y 2 + z 2 ^ да
z > z 0 , где
q(t, x, y, z) - средняя концентрация примеси в точке
(x, y, z) e E+, E+ = {(x, y, z): x, y e (-да; да), z > z0 > 0}, в момент времени t; Kx, Ky, Kz - коэффициенты
турбулентной диффузии соответственно вдоль осей Ox, Oy , Oz ; U - компонента средней скорости ветра вдоль оси Ox; W - скорость осаждения частиц примеси вдоль оси Oz ; <р(x, y, z), f, Vs - соответственно фоновая концентрация, функция источника, скорость сухого осаждения этой примеси; z0 = const > 0 - уровень шероховатости подстилающей поверхности.
Для расчета количества примеси, выпадающей из атмосферы на подстилающую поверхность, можно использовать соотношение [1]
p(t,x,y,z0) = JP(t,x,y,z0)dt .
(2)
где Р(/,х,у,г0) - плотность осадка в точке (х,у), расположенной на плоскости г = г0 за время действия источника /; Р(/, х, у, г0) = ^д^, х,у, г0) - вертикальный поток частиц в точке (х, у, г0) в момент
времени t; Vs = const > 0 - скорость сухого осаждения частиц этой примеси.
Соотношения (1), (2) определяют математическую модель процесса осаждения примеси на плоскость
z = z0, z0 > 0 [2].
Для вычисления P (t, x, y, z 0) авторами данной работы разработан алгоритм в интегрированной среде Delphi, реализованный в программном продукте «DODS». Программа «DODS» зарегистрирована в отраслевом фонде алгоритмов и программ под номером 9293 [3].
Кроме того, в рамках модели (1), (2) подробно исследованы [4, 5] следующие важные (с прикладной точки зрения) задачи для точечного источника мгновенного действия: задача 1 - определить точку (x, y), расположенную на плоскости z = z 0, в которой плотность осадка P (t, x, y, z0) в момент времени t, t e[t0, T] принимает максимальное значение; задача 2 - определить интервал времени (0, t) максимальной длины, в течение которого плотность осадка P (t, x, y, z 0) в точке (x, y, z0) не превысит заданной величины P (т.е. P (t, x, y, z0) < P ), где P - предельно допустимое значение плотности P (t, x, y, z0) на подстилающей поверхности, превышение которого становится небезопасным (согласно медицинским нормам) для здоровья людей. В данной работе будут представлены результаты исследований следующих задач, которые являются обратными в рамках модели (1), (2); задача 3 -определить мощность Q(t) точечного источника непрерывного действия по экспериментально заданным плотности осадка на подстилающей поверхности
z
P (t, x, y, zо), скорости сухого осаждения частиц примеси Vs, основным параметрам модели (1): U, W , Kx, Ky, Kz, и координатам источника (0,0, H); задача 4 - определить мощность Q(t) точечного источника непрерывного действия по экспериментально заданным значениям концентрации примеси q(t, x, y, z), основным параметрам модели (1): U, W , Kx, Ky, Kz, и координатами источника (0,0, H).
2. Методика решения задачи 3. При решении обратной задачи 3 из п. 1 будем предполагать, что коэффициенты турбулентной диффузии и функция источника f (t, x, y, z) имеют вид U = const, Kx = Ky = K0U,
K0 = const Kz = K2zn, K2 = const, n = const, 0 < n < 2, f (t, x, y, z) = Q(t)g(x)S(y)S(z - H), где S(s - s,) -дельта-функция Дирака; Q(t) - количество примеси, выбрасываемой источником в момент времени t .
Тогда, согласно [1], P (t, x, y, z0) можно представить в виде
_ t т
P (t, x, y, z0) = Vs U q(s, x0, y0, H; t, x, y, z0 )Q(s)dsdT, (3)
00
где q(s, x0, y0, H;t, x, y, z0) - функция Грина задачи (1).
Таким образом, задача 3 сводится к решению интегрального уравнения относительно Q(s) при заданных значениях P(t,x,y,z0) и q(s, x0, y0, H;т, x, y, z0)
(определяемой из [1]). Соотношение (3) представляет собой интегральное уравнение 1-го рода относительно неизвестной функции Q(t). Задача построения решения таких уравнений является некорректно поставленной [6]. Поэтому для ее численного решения относительно Q(t) воспользуемся методом регуляризации (по А.Н. Тихонову), согласно которому задачу построения решения Q(s) уравнения (3) можно свести к минимизации функционала вида
sir^ML -
P p(Q) = |\PQ - P "2
L>[0,f ]
HL2[0,i ]
= J (V J J q(s,0,0, H; T, x, y, z)Q(s)dsdT -Ps (т, x, y, z))2 dt +
0 00 t
I _ _ |./=1,.,N
(P T , PT ) +p
I t=1,.., N
t /=1,..,N Л
JT+jdt X = (4)
0 i=1,..,N j ICN j
= | J PTPsdT
V0 ji=1,..,N
Систему (4) можно переписать в матричном виде (Г + /1)С = V, (5)
где Г = ||ггу || - матрица Грамма с элементами = (Р Т,РТ); С - вектор-столбец неизвестных;
V- вектор с элементами вида (P Т,Pg );
R = r^J -
+ Q2(т)аТ , где P(^)=|ds\q(s,x„,y0,H;t,x,y,z0)(«)dt;
0 0 0
Pg - экспериментальные (приближенные) значения плотности осадка P ; ц- параметр регуляризации. Пусть Q(t) представляет собой полином вида
N
Q = ZСТ , С = const, i = 1,...,N .
i=0
Тогда задача определения функции Q(s) сводится к решению системы алгебраических уравнений относительно Ci, i = 1,...,N :
матрица размера N х N, Ту = (Т, Т ).
Так как матрица (Г + /) симметрична и положительно определена, то систему линейных алгебраических уравнений (5) можно решить методом Холецкого [7].
Параметр регуляризации / в (5) определяется методом подбора или с помощью соотношения [8]: т г т
р(/) = | (V Ц д(s,0,0, Н;т, х, у, zо)0/(s)dsdт-0 00
-Р5(т,х,у,z0))2dt = 52 , где Q/(s) - решение интегрального уравнения (3); 5 - заданная погрешность, 5 ^ 0.
Для численного решения уравнения (3) методом регуляризации авторами разработан алгоритм, реализованный в виде программного продукта «1Р18».
3. Пример решения задачи 3. Для оценки качества работы алгоритма из п. 2 воспользуемся экспериментальными данными, взятыми из отчетов Центра лабораторного анализа и технических измерений по Южному федеральному округу (ЦЛАТИ по ЮФО) и содержащими информацию о выбросах в атмосферу загрязняющих веществ (диоксид азота). Согласно этим данным: Н = 20 м, и = 2,5 м/с, К0 = 0,1 м, К2 = 0,1 м, п = 0,15, г0 = 0 с, г = 55 с, V, = 0,1 м/с,
) = 2,71 кг/с.
В начале рассчитаем значения плотности осадка Р численно, решив прямую задачу с помощью программы «ВОВ8». Затем восстановим значения Q(t) (т.е. найдем решение обратной задачи) с помощью программы «¡РК». Сравним значения восстановленной и заданной мощностей (графическая визуализация результатов проведенных расчетов приведена на рис. 1).
160,0 140,0 120,0 100,0 80,0 60,0 40,0 20,0 0,0
О 5 10 15 20 25 30 35 40 45 50 55
Рис. 1. Графическое изображение значений экспериментальной (гладкая линия) и расчетной (ломаная линия) мощностей источника примеси
Согласно рис. 1, расчетные значения мощности источника Q(t) с высокой степенью точности согласуются с экспериментальными значениями Q(t) В
—8
данном примере параметр регуляризации / = 10 .
4. Методика решения задачи 4. Пусть при решении задачи 4 из п. 1 в математической модели (1) фоновая концентрация р(х, у, г) не учитывается, т.е.
р(х, у, г) = 0 , начальный момент времени /0 = 0 .
Из первого уравнения (1) следует, что
/ (/, х, у, г) = ^ + и ^ - W * + ад -д/ дх дг
д к дд д к дд д ^ дд дх х дх ду у ду дг г дг Если источник / является точечным с координатами (х0,у0,И), то, учитывая соотношение /(/, х, у, г) = =ШЖх - х0)8(у -у0)8(г - И), легко убедиться, что если / заменить на $/), то решение задачи (1) не изменится:
Qt) J-LL + и dJL - w дЛ- + aq -
dt dx dz
-4- k.
dL -d K
L
-д к,
dx dx dy dy dz
L
~dz
(6)
Согласно (6), для вычисления значений Q(t) требуется предварительно вычислить значения производных функции д (/, х, у, г) по /, х, у, г (по х, у, г -
до 2-го порядка включительно).
Задача нахождения производной п -го порядка
г(0 функции и(0, т.е. ) = и(п)(/), сводится к решению (относительно )) интегрального уравнения 1-го рода [6]:
1
J——(t -т)n-1 z(T)dT = u(t), u(0)=0. o(n -1)!
(7)
Обозначим Rzz (t, x, y, z ) =
d q(t, x, y, z)
dz
n / Ч д q(t, x, y, z) / \ d q x,y, z)
Ryytt,x,y,z) = ^j-^, R„tt, x, y, z)=-—-;
dy
c,2
n / \ dq(t, x,y, z) , \ dq(t, x,y, z)
Rz (t, x,y, z )= J , Rx (t, x,y, z)= J
z
x
т, ! \ дай, х, у, г) „
Я/ ((, х, у, = ^—-—1. Согласно (7), для определе-
д/
ния(/,х,у,г) будем иметь интегральное уравнение:
г
д(/, х, у, г) = $ (г - т) ■ Ягг (/, х, у, т)йт . (8)
0
Соотношение (8) представляет собой интегральное уравнение 1-го рода относительно неизвестной функции Ягг. Задача построения решения таких уравнений является некорректно поставленной [6]. Перейдем от уравнения (8) к его дискретному аналогу:
q (t, x, y, z) = Е |(zp - zk)Rzz (t, x, y, zk) • rk |, (9)
k=1
rk =
zk + 1 - z
k = 1,
k-1
zp z P-1
2
k = 2,.., (p - 1),
k = p,
zj,..., z p - точки деления интервала [0, z].
Согласно (9), по значениям д(/1, х, у, г),..., д(tN,х,у,г), заданным в точке (х,у,г) в различные моменты времени /1,...,, с ошибками измерения V = ), 1^2 =2), ., Ум =?(/ы), где V (/) - белый гауссов шум, / = /1,...,, требуется найти значения хy,гк,x,У,гк), к = р . Введем в рассмотрение матрицу А = (Ак):
А,к = (гр-гк) ■ гк, к = р, ' = 1,...,м.
Тогда с учетом введенных обозначений и замечаний, (9) можно представить в виде системы
Е [Ak • Rzz (t1, x, y, zk
k=1
)]+ ^ = L(t1
x, y, z),
Е [A2k • Rzz (t2 , x, У, zk Ц + V2 = L^ x, У, z),
k=1
i
X ^Ык ■ (/М , x, У, гк )] + V N = д(/М , x, У, г) ,
к=1
или в матричном виде
АЯгг + ~ = д . (10)
Для подавления влияния белого шума ~(/) на (/к, х, у, г р), к = 1,..., N, можно использовать фильтр Калмана-Бьюси [9]. Тогда оптимальная в среднем квадратическом смысле апостериорная оценка Кгг
решения системы (10) может быть найдена из соотношений: Кгг =ф+ РАТК-1 (д^ - Аф), ф= М [Кгг ], Р = (N-1 + АТЯ-1 А)-1, N = М[(Ягг - ф)(Ягг -ф)Т ], Я = М Т ].
Аналогичным образом определяются Яхх, Яуу , Я/ , Ях, Яг соответственно для Яхх, Яуу , Я/ , Ях,
Яг. Подставляя найденные оценки в (6), получим наилучшую в среднем квадратическом смысле оценку £>0к) значения ^к):
Q(tk) = Я/ + и ■ Ях - W ■ Яг+а ■ д -
-Кх ■ Яхх - Ку ■ Яуу -кг ■ Ягг, к = 1,..., N. (11)
5. Пример решения задачи 4. Для оценки качества работы алгоритма, описанного в п. 4, воспользуемся экспериментальными данными, приведенными в п. 3 настоящей работы. В начале рассчитаем значения концентрации примеси , / = 1,...,N, численно. Для этого решим прямую задачу (1), из которой определим значения концентрации примеси д(/, х, у, г). Затем, используя соотношение (11), найдем Q(tk), к = 1,..., N (т.е. найдем решение обратной задачи).
Результаты проведенных расчетов (их графическая визуализация) приведена на рис. 2.
Сравнительный анализ проведенных расчетов с экспериментальными данными (рис. 2) показывает, что с помощью соотношения (18) мо жно восстановить значения мощности источника Q(tk) с достаточно высокой степенью точности.
zz
2
2
2
Рис. 2. Графическое изображение значений экспериментальной (гладкая) и расчетной (ломаная) мощностей
Основное достоинство алгоритма из п. 2 - совпадение (с малой погрешностью) расчетных и восстановленных значений мощности точечного источника почти на всем рассматриваемом интервале времени. Только в конце интервала времени наблюдения эти значения начинают заметно отличаться друг от друга. При этом этот эффект не зависит от длины интервала (0,Т). Основной недостаток: в случае, когда ядро интегрального уравнения (3) принимает значения, меньшие 10-4, восстановить значения мощности осадка практически невозможно. Этот недостаток алгоритма можно устранить двумя способами: или путем масштабирования, искусственно увеличив значение ядра уравнения (3), умножая обе его части на одно и то же достаточно большое (порядка 104) положительное число, или проводить замеры плотности осадка на подстилающей поверхности в точке, где эта плотность достигает максимального значения.
Достоинство алгоритма, описанного в п. 4, - незначительные расхождения значений экспериментально заданной и расчетной мощностей источника выбросов примеси на всем рассматриваемом интервале времени независимо от его длины. Основной недостаток - необходимость в априорной информации о
Поступила в редакцию_
значениях ковариации ошибок решения N и вектора р. Добыть такую информацию на практике часто оказывается затруднительным.
Оба алгоритма (по отдельности) можно использовать в прикладных исследованиях (оба дают удовлетворительные результаты) на практике. Однако целесообразнее использовать эти алгоритмы одновременно, чтобы уменьшить вычислительные ошибки.
Литература
1. Семенчин Е.А. Аналитические решения краевых задач в математической модели атмосферной диффузии. Ставрополь, 1993. 142 с.
2. Бызова Н.Л., Гаргер Е.К., Иванов В.Н. Экспериментальные исследования атмосферной диффузии и расчеты рассеяния примеси. Л., 1991. 278 с.
3. Лоскутова Е.О. Регистрация в ФГНУ «Государственном координационном центре информационных технологий» разработки, представленной в отраслевом фонде алгоритмов и программ: Электронный комплекс программ определения количества легкой и тяжелой примеси, выпадающей на подстилающую поверхность (БОБ8). Номер гос. регистрации: 9293. Дата регистрации: 14.11.2007 г.
4. Лоскутова Е.О. Прямые и обратные задачи расчета количества примеси, выпадающей на подстилающую поверхность : автореф. дис. ... канд. физ.-мат. наук. Краснодар, 2009. 26 с.
5. Годяева Е.О. (Лоскутова Е.О.), Семенчин Е.А. Обратная задача для плотности осадка, выпадающего на подстилающую поверхность от непрерывного точечного источника // Изв. вузов. Сев.-Кавк. регион. Естеств. науки. Приложение. 2006. № 12. С. 31-34.
6. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М., 1979. 142 с.
7. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М., 2002. 632 с.
8. Денисов А.М. Введение в теорию обратных задач. М., 1994. 208 с.
9. Верлань А.Ф., Сизиков В.С. Интегральные уравнения: методы, алгоритмы, программы. Киев, 1986. 543 с.
20 апреля 2009 г.
