CC BY
18
счетной конфигурации с нормальным вектором е1 = ссксп + действуют усилия
{ =71^+71^2/2,^ { 2=°. {1=Ре1'П~1- (4)
Конформное отображение внешности эллипса на внешность единичного круга дается формулой:
ь = {$+т2 0<т<1. (5)
Решение. Подставляя (4) в (1) и решая систему, найдем константы в выражении для комплексных потенциалов:
Р
Hi
= 0, 2а12 +сх =
AßT]
bn = 0, Ьп = 0, 2dn + с1 = 0, ij j = 0, tl2 = 0.
(6)
Подставляя (6), (5) в (2) и применяя аппарат интегралов типа Коши, получаем выражения для комплексных потенциалов:
т2еР2
Ф,=
in?p
1 4 iiri(m2-Ç2)
, Ф2 = -Í-
Vi
8112Г12(т2-¥У'
+т2^2 +т4^ +т6^2)р
4 Ц71 {т2-%2У используя которые вычисляем (3).
Для не очень вытянутых эллипсов тангенциальные напряжения на контуре достигают максимума в вершине с наибольшей кривизной. Коэффициент концентрации К равен:
т2 +1 1т1 — 4 т2 +1 р
К = 2-
- + —
(7)
1-т2 2 (ш2-і) Ч
Область на плоскости (т,р/рі), в которой верна формула (7), расположена под кривой на рис. 1. В области, расположенной над кривой, происходит расщепление. На рис. 2 показаны графики относительных значений тангенциального напряжения Зцхр IР в зависимости от параметра / при р / ц = 0,2,
т = 0,9 в окрестности вершины эллипса с ? = 0. Жирная линия соответствует квадратичному приближению. В вершине появляется минимум напряжений, а симметрично около нее возникают два максимума.
Рис.1
-0,8 -0,6 -0,4 -0,2
0
Рис. 2
0,2 0,4 0,6 0,!
Таким образом, показано, что расщепление одного максимума напряжений на два характерно для нелинейных решений в области конечной, хоть и достаточно большой кривизны контура, а не только для разреза. Литература
1. Астафьев В.И., Крутов А.Н. II Изв. РАН. МТТ. 2001. № 5. С. 125-133.
2. Жуков Б. А. И Изв. вузов. Сев.-Кавк. регион. Естеств. науки. Спецвыпуск. 2001. С. 70-71.
Волгоградский государственный технический университет_____________________________________26 ноября 2003 г.
УДК 519.5
ОБРАТНАЯ ЗАДАЧА ИСТОЧНИКА ЗАГРЯЗНЯЮЩИХ ПРИМЕСЕЙ, РАСПРОСТРАНЯЮЩИХСЯ В ПРИЗЕМНОМ СЛОЕ АТМОСФЕРЫ
© 2003 г. Н.И. Каргин, В.И. Наац
The paper is concerned with methods for solving inverse pollutants source problem as applied to routine monitoring of the pollution condition of the atmosphere. The authors focus attention on the organization of inverse pollutants source problem, development of method for the construction of the calculation model and realizing calculation experiment with describing its results.
Объектом исследования данной работы являются процессы переноса загрязняющих веществ (ЗВ), распространяющихся в приземном слое атмосферы от источников различной длительности действия. Подобные
задачи решаются, в частности, на основе уравнения турбулентной диффузии, исследование которой связано с изучением закономерностей распространения примесей в атмосфере.
Важная роль в исследовании этих явлений отводится разработке математических моделей, построению вычислительных схем и алгоритмов, созданию на их основе соответствующего программного обеспечения, информационно-измерительных систем оперативного контроля и управления экологическим состоянием приземного слоя атмосферы. Эти системы должны также включать в себя технические средства контроля, такие, как стационарные посты наблюдения, оснащенные датчиками концентрации ЗВ в атмосфере, технические средства контроля основных метеопараметров, определяющих состояние пограничного слоя атмосферы, мобильную экологическую лабораторию, оснащенную соответствующим набором измерительных средств.
Новым направлением в решении задач прогноза экологического состояния атмосферы является создание эффективных вычислительных моделей, основанных на решении так называемых обратных задач. С помощью этих моделей можно получить оценку состояния турбулентности атмосферы, определить интенсивность источника по известным замерам концентрации примесей в пунктах наблюдения и др. Во всех случаях большое значение имеет оперативное поступление со стационарных постов наблюдения или мобильных экологических лабораторий измерительной информации.
Представленное ниже исследование имеет своей целью разработку вычислительной модели решения обратной задачи источника ЗВ. Дана постановка обратной задачи источника, разрабатывается вычислительная схема ее решения, приводятся основные этапы проведения вычислительного эксперимента и результаты его программной реализации.
1. Постановка обратной задачи источника. Предполагается, ЧТО В точке (пункте) мо(хо,уо,2о) расположен источник ЗВ интенсивностью Б. Источник 8(^,Мо), 10 ~ъ кг/с, работает в течение конечного промежутка Времени [1о,Т] и производит выбросы ЗВ в моменты времени где ^ < % < 1о+Т. Под действием ветра \(Ух,Уу,Уг), м/с и турбулентной диффузии
К{Кх,Ку,Кг), м2/с идет процесс переноса ЗВ из
пункта М0(хо,уо,2о) в пункт М(х,у,г), в котором за время ЗВ накапливается в количестве <3(МД 10 _3 кг/м3, для расчета которого ранее была получена формула [1]:
']5(£,М0ЖМ0,М,Г,!)^ = б(М,г) , (1)
10
где 11 < I <12 , ^^о+Г-Ть 12=1о+Т+1*+Т2, Т1 и т2 могут определяться, например, как ш , где а=(0,25; 0,5; 0,75); параметр 1* - момент времени, когда в точке наблюдения М будет максимальное значение концентрации ЗВ, для расчета I* ранее были получены следующие расчетные формулы [1]:
•^1+ —1 1 Гт/2
(2)
¿¡(М0,М,1,^) - концентрация ЗВ, фиксируемая в
момент времени I в пункте наблюдения М(х,у,г), для расчета которой применяется функция [2]:
, *, у, г) = б0<7(*> У. г)-
9(г. х, у, г) = --------------г-х
Ы1-10)У2[КХКУК^2
" [(*-*0)-У,(г-/0)]21„
4 Кх((-{0) }
Ь-УоЬМ'-'о)]2]
хехр
хехр
хехр^-
4 Ky(t-t0) [(z-z0)-Vz(í-í0)]2
(3)
^ 4ЛТг(г-г0) ]
называемая гауссовой функцией распределения концентрации примеси.
Обратная задача источника для уравнения (1) состоит в определении интенсивности источника 5(М0,£) по известной правой части уравнения (1) -(2(А/,0. Данная задача относится к классу некорректно поставленных задач, а уравнение (1) является уравнением Вольтерра 1-го рода с переменным верхним пределом. Для решения этой задачи в статье разрабатывается и реализуется алгоритм, основанный на численных методах решения некорректных задач [3].
2. Построение и обоснование регуляризирующего алгоритма для численного решения обратной задачи источника. Вывод основных вычислительных формул для решения обратной задачи источника и обоснование метода начнем с рассмотрения интегрального уравнения Фредгольма 1-го рода;
| К(у,х)(р(х)с1х = /(у). (4)
X
где хех = [х1,х2] и уеУ = [У,,}'2]. Обозначим через Ф множество возможных решений уравнения (4), соответствующих функциям /(у) из исходного множества Б. При решении практических задач множество Р часто неизвестно и представлено в виде множества эмпирических функций , которые не являются дифференцируемыми, в лучшем случае они ограничены. Функции из множества исходных данных Ра обозначим /ст и назовем а - приближением. Если считать, что / = К(р , где (ре Ф , то примем далее, что ||/-/ст||<сг,/представлено в эксперименте своим а - приближением /с , где величина <т характеризует погрешность эксперимента. Далее полагаем, что Ф входит в класс функций Липшица с константой Ь, т.е. любая функция (ре Ф удовлетворяет условию
\(р{х2)-(р{х])\<ь\х2-хх\. В этих условиях возможно применение вариационных подходов к решению операторных уравнений. Для интегрального уравнения Фредгольма 1-го рода вариационный принцип обычно реализуется в виде оптимизацион-
ной задачи для так называемого сглаживающего функционала Та [<р]:
II II2
ТаЫ= У^(у.л)ф(дс))<йс-/<г(у) +
||^(П (5)
+ар\\<р0 -qf^xy
В этом выражении (р0 (х) - некоторая функция из Ф; р - множитель, введенный в (5) из соображения размерности; а > 0 - параметр регуляризации. В качестве решения принимается функция сра (х), доставляющая минимум сглаживающему (квадратичному)
функционалу Та [(р\, т.е. min Та [<р].
<ре Ф
Значение параметра а в (5) не определено и осуществляется на основе тех или иных рекомендаций, формулируемых в процессе соответствующих численных экспериментов. В проблеме выбора «подходящего» значения параметра а наилучшими могут быть приближенные неравенства вида:
min Та [<р] < 2сг 2, а также IIК(ра - fa II2 < а 2
(ре Ф
и ap\(pa-(pQf <о2 .
Соответствующий способ получил название выбора параметра а по невязке [3].
Рассмотрим случай, когда множество исходных данных Fa является конечномерным. Заметим, что последнее обстоятельство в большей мере соответствует нашей прикладной задаче по распространению ЗВ в атмосфере, которая имеет дело с экспериментальным материалом. Теперь постановка задачи формулируется следующим образом: задано множество значений {fa (у1 )}л и при этом
f(yt) = ¡K(yl,x)<p(x))dx , где l = l,n. х
Требуется найти <р{х) таким образом, чтобы ||/- fa\ бьша минимальной. Для решения подобной задачи используем метод сглаживающего функционала
Tabh 1 JК, (х) (p(x)dx- fal
/=\Х
+ар f (<р0 (х) - (р{х))2 dx, х
(6)
и подставим (8) в (7). Получим
Л
2 Ki (х)(-а рС,) - ар(у0 (х) - (р{х)) = 0, откуда имеем i=i
<p(x) = (po(x) + Y,ClKl(x) . /=і
(9)
Считая <р0(х) = 0, результат (9) можно интерпретировать следующим образом: решение вариационной задачи для функционала (6) можно представить в виде
рм-хедм, (ю)
/=1
т.е. имеет место разложение по системе функций {К/(•*)}„. Ясно, что система функций {К1 (х)}и должна быть в совокупности линейно независимой. Последнее обстоятельство в какой-то мере проясняет проблему наилучшего выбора в эксперименте системы узловых точек {у1 }л. Далее необходимо получить линейную алгебраическую систему для определения вектора коэффициентов С = (С^, С2С„). Для этого достаточно подставить (9) в выражение (8). Имеем
0(х')+ 1Кк(х')Ск
X { к=1
Введем обозначения д1к = \К1(х')Кк(х')с1х' и
х
/о/ = \К!(х')(р0(х')с1х'. х
Выражение (11) примет вид:
^9/*^+«рС( =/^ —/ог , 1 = 1,п. (12)
к=1
Решая систему (12), находим коэффициенты С = (С1,С2,—,С„) в представлении (9) искомой функции (р(х) в базисе {АТ; (х)}л . Вычислительный
алгоритм достаточно прост и его достоинством является естественная редукция вариационной задачи к алгебраической системе. Представление (9) позволяет применять расчетно-аналитические методы в прикладных исследованиях, использующих результаты обращения.
Изложенный численный метод решения интегрального уравнения Фредгольма 1-го рода может быть применен к интегральным уравнениям Вольтер-ра 1-го рода. Сформулируем следующую задачу: задана система уравнений вида:
dx'-fol =-<*рС1 . (11)
где Кі(х) = К(уі,х), 1 = 1,п. Для этого квадратичного функционала уравнение Эйлера 8 Та [<р]=0 за-. писывается в виде:
(
2Я/(х) \К1{х')ср{х)с1х'-/а1 1=1 X
-ар(<р0(х)-(р(х)) = 0. (7)
Введем систему ЧИСЛОВЫХ коэффициентов Сі согласно формуле
¡К1(х')(р(х')(іх'-/а, =-арС,,. 1=йі (8)
х
У\
(13)
\К (У1,х)<р(х)с1х=/о1, 1 = 1, п
Л
и требуется найти <р(х), где хб[х,,х2], у1 е [К[,У2] (X = У), из условия минимума сглаживающего функционала
таЫ=2
У\
і=і п У\
\K{yhx)(p{x)dx-fai
+ар£ j(poW“9>W) dx-
i=ix.
Считая, что представление вида (9) применимо для искомой функции и в данной задаче, приходим к следующей квадратичной форме:
„ / „ \2
г«(с)=2
(=1
Б Рік С к + /о/ ~ fol к=1
+
+ ар Ep¿C,C* ,
l,k=1
(15)
(16)
dina
где рш = ¡Kl(x)Kk(x)dx, 1 = 1,п.
Л
В заключение можно сделать несколько замечаний относительно выбора параметра регуляризации а. Помимо упомянутого выше способа оценки значений а по невязке можно указать возможность оценки так называемого квазиоптимального значения ссКвопт . Речь идет о выборе такого значения параметра а из интервала, скажем, [ctj, сс2 ]» в пределах которого решение (ра (х) наименьшим образом зависит от а. Количественным критерием может служить условие
d {х)\\ -> min .
ае[а1,а2\
Выявление подобного интервала [се1, ск2 ] в каждом конкретном случае требует численного анализа свойств параметрического семейства решений {<ра }.
3. Вычислительная схема решения обратной задачи источника. Сопоставляя уравнения (1) и (13), отмечаем, что в качестве ядра K(yt,x) в уравнении (13) выступает функция q(M0,M,t,^), определяемая выражением (3), в котором t = tt, t Е , ¿2 ];
t\ =tü + t* t2=t0+T + t*+t2-, £e[f0,i-f*J;
= y'i=t[-t*; l = 1,n. Искомой функ-
цией (p(x) будет функция источника S(M0,^), в качестве функции fm (у) выступает Q(M,t).
Вычислительная схема определяется соотношениями (13) - (16), в которых первоначально предполагается, что <р0 = 0 . Обобщенный алгоритм решения включает в себя следующие этапы.
Этап I. Решение прямой задачи.
Решая прямую задачу, вычисляем «экспериментальные» значения функции количества ЗВ в пункте М на сетке {у,}, I = 1,л по формуле (13). Обозначим fm(y) = f(y) и (рт{х) = (р{х). Предположим, что функцию (рт (jt) можно задать аналитической формулой вида q>m(x) = q0(l + a •sin(fot + c)).
Этап II. Получение приближенного значения правой части (13).
Моделируя процесс появления случайных ошибок при изменении «экспериментальных» значений fm (у), заменим fm(y) на приближенное значение fa(y) по формуле fa(y) = fm(y)-(l + Q,cj), где
0,- случайные числа, равномерно распределенные
на интервале (-1, 1); а — погрешность правой части
(17). Обозначим /ст (у,) = /Ст;.
Этап III. Решение обратной задачи. Выбор параметра регуляризации по невязке, реализуемой по указанной ниже схеме.
Берется конечный отрезок монотонной последовательности чисел а0,а1,а2,...,ап, например отрезок
геометрической прогрессии ак=а^к, ' к-0,п,
0<4<1 [3].
1) к=0, задаем начальные значения а0 и q;
2) вычисляем ак - а0цк;
3) для каждого значения ак находятся значения коэффициентов Сак (х), минимизирующих функционал
(15); _
4) коэффициенты С„4 ={С;}, 1 = 1,п определяются
путем минимизации функционала тт Та (С), на-
Сей„
пример, с помощью «эвристического симплексного метода» [4];
5) определив С„4 (х), далее получаем значения искомой функции по формуле (10);
6) на основе полученного решения <ра (л:) вычисляем
приближенное к /с(у) (модельное) значение функции/(у): ](у)=\(ра(х)К(у,х)<1х\
Уо
7) далее, согласно алгоритму нахождения параметра регуляризации а по невязке, вычисляем норму:
Ц/ОО-ЛгОО^
1/2
Если ||/(у) - fa (у)|| ~ о , то выбираем а* = ак и вычисляем соответствующие значения (Ра- (л) и
/(у). Если /(у) - /ст (у)|| > а , то к=к+1 и возвращаемся на шаг 2.
4. Постановка численного эксперимента для обратной задачи источника. Для реализации рассмотренного выше алгоритма необходимо прежде всего задать исходные данные и провести их нормировку, а также основных переменных и полей задачи.
Исходными данными для уравнения (1) будут: а) время работы источника [1о,Т], с; б) координаты источника Мо (хо, уо, го), м; в) расстояние от источника до пункта наблюдения Я, м - расстояние от точки М0 до точки М; г) (2о, 10 ~3 кг/с - количество ЗВ, выброшенное источником в начальный момент времени д) турбулентность атмосферы в приземном слое, определяемая коэффициентом турбулентной диффузии К={КХ, Ку, К,}, м2/с; е) скорость ветра У= {Ух, Уу, У2}, м/с.
Помимо исходных данных, необходимо также вычислить следующие параметры задачи:
а) координаты пункта наблюдения М, где проводятся замеры концентрации ЗВ, поступающих от источника: х = Я/ -\/2, у = Я / л/2, г > г0;
б) время прихода в пункт М «пика» (максимума) ЗВ, поступающих от источника, делающего выброс в момент времени £,, по формуле (2);
в) интервал времени, в течение которого идет накопление ЗВ в пункте наблюдения рь 12]:
г, =г0+г*-т1, г2 = г0 + Г + Г* + т2, где Т| = т2 = 0,1-г*;
г) задать разбиение отрезка [1!, 12] и сформировать
массив {у1 }= {^}, 1 = \,п, 1М=11+Ы,к=———.
п
Далее проводится нормировка переменных задачи. Положим, что д(М0,М,г,^) = <7*<7(М0,М,Г,|,У,ЛТ), где $,М0,МV,К - нормированные величины, принимающие значения в интервале (0,1); х0 = х01 Я, Уо=Уо/К’ г0=г0/Я; х = х/Я, у = у/Я ;
К = КХ=КХ/К\ Ку=Ку1К*,
Кг =Кг1К\где К* = таях{кх,Ку,К1}‘,
У = $х>Уу,У'}, УХ=УХ/У\ Уу=Уу/У\ Уг=У1!У\ где V* = тах^х,Уу,Уг};
Г, = Г, /Г, Г =г*/Т , Iо =10 П, Г2 =^!Т,
£0<|<г-Г, |=|/г, <*§ = га|.
Тогда можно записать следующее выражение:
1
q(M0,M,t,Ç,V,K) =
bn(t-i)Y2[KxKykzY2
xexp
xexp -
xexp-!—
4Кг(Г-|)
[(x-x0 )-Vx{t-$ 4Kx(t-b
X
X
4*y(i-S)
1
\3/2
(Ж )J
Далее рассмотрим функцию источника. Допустим, что
S (м о. £) = So о, |) = qa (1 + а sin(bx + с)). Тогда S0=<70,a 5(М0,|) = g0(l + asin(£x + c)). В итоге имеем:
где /5(М0,£)д(М0,М,Г,£)Л; = Q(M,t),
| 5(М0,|)^0,М,Г,|,У,ВД = б(М,Г),
*0
О* =т q0 <?*.
При этом /т(>’) = 2(М,£), <рт(л:) = 5(М0,|),
К(у,х) = !}(М0,М,г,1;,У,к), у = !, х = £.
В завершение постановки численного эксперимента определим аналитический вид параметра р в
сглаживающем функционале Та\р] исходя из соображений его размерности следующим образом:
Г-Гп
Р =
(лг-^оХР — Ро)Сг —г0)
5. Результаты численного эксперимента. Вычислительная схема была реализована программно при следующих исходных данных: г» = 0; Г=30мин; хо=Уо=0 м; 2*7=50 м; Я=1 км; £2о=50-10 “3 кг/с; Кх= Ку = Кг— 5 м2/с; Ух = Уу= Уг =2 м/с; г = 52 м; п = 10; а = 0,5; ¿> = 8; с = 0,5.
На рис. 1 представлены результаты расчетов значений искомых функций (р (х) и (ра> (х). Относительные ошибки от измеряемой величины /ст (у) в узловых точках {у,} составили 5 % (ст = 0,05). Оптимальное значение параметра регуляризации а*=9,78Е-08. При этом норма \фт-<ра\-0,12,
Цфт^-фа* (*)Ц
что составляет
IM*) I
!-х100 = 11'
фа*
Рис.1. Графики функций (р (х) и фа‘(х) при
а* =9,78Е-08 и сг = 0,05
На рис. 2 приводятся графики функций /а (у) и q = /(у), где значения ^ = /(у) рассчитываются на основе (ра-(х). Норма ||/ст (з1) — /(у)]) — 7,01£Г — 04, а
оо-7оо!
" >
'-X100 = 1,9%.
6.10Е-09 2.44Е-08 9.76Е-08 3.90Е-07 1.56Е-06 6.24Е-06 2.49Е-05 - • в- - -г1 — -•— ■ г2 - - -А- - -гЗ
Рис. 3. Зависимость Г = |фт ~(Ра\ от величины погрешности измеряемых данных сг = 0,05, (7 = 0,10 и
сг = 0,15 соответственно при различных значениях а
В этой ситуации повысить точность получаемого решения можно дальнейшим подбором параметра регуляризации а. Так, при а = 0,10 оптимальным
значением а будет значение а = 1,52Е-09, а при <т = 0,15 дальнейшее уменьшение параметра а уже не приводит к снижению нормы г = ||(рт - (ра |.
/<,00 я
Рис. 2. Графики функций fa(y) и q = /(у) при а* = 9,78Е-08 и сг = 0,05 На рис. 3 показано изменение точного решения от расчетного г = ||<рт —<ра || при различных а для трех значений погрешности измеряемых данных а = 0,05, (7=0,10 и (7 = 0,15 соответственно. Расчетные данные показывают, что для каждого значения а существует свое оптимальное значение а, при котором достигается минимум г = ||<рт - (ра ||. Однако величина а оказывает значительное влияние на точность получаемого решения. Так, при а = 0,05 отклонение получаемого решения от точного составляет 11 %, для сг = 0,10 оно составляет 25 %, а для а = 0,15 - 28 % •
при одном и том же значении а* = 9,78Е-08 .
Г 0,8,
На рис. 4 показано изменение g = ||/сг(>’)-/(у)|
в зависимости от выбора а при различных значениях (7 = 0,05, (7 = 0,10 и (7 = 0,15. Для значений
а*=2,44£-08 и а* =9,78Е-08 точность приближения ]{у) к /ст(у) - !/ст (у) -1 (у)! = 5,88£ - 04 для (7 = 0,05 и <7 = 0,10 и \\/АУ)-?(У)\ = 5Л4Е-04
для а =0,15. Это позволяет говорить о вполне удовлетворительных результатах работы вычислительного алгоритма.
1.6Е-02
1.4Е-02
1.2Е-02
1.0Е-02
8.0Е-03
6.0Е-03
4.0Е-03
2.0Е-03
щ
/1
НА-**-*
О.ОЕ+ОО 6,1 Е-
09 2.4Е-08 9.8Е-08 3.9Е-07 1.6Е-06 6.2Е-06 2.5Е 05
■ -А" “ д1 - - - 'д2 • ■ •Я- ■ дЗ
Рис. 4. Зависимость g = ||/ст (у) — /(зОЦ от величины погрешности измеряемых данных (7 = 0,05, (7 = 0,10 и
(7 = 0,15 соответственно при различных значениях <Х
Литература
1. Паац В.И. II Материалы междунар. науч.-техн. и Российской науч. школы. Ч. 3. М., 1998.
2. Берлянд М.Е., Генихович Е.Л. Метеорологические аспекты загрязнения атмосферы. Л., 1971.
3. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М., 1979.
4. Доннис Дж. мл., Шнабель Р. Численные методы безусловной оптимизации и решения нелинейных уравнений. М., 1988. .
Северо-Кавказский государственный технический университет
26 ноября 2002 г.
