CC BY
29
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
УДК 519.5
ВЫЧИСЛИТЕЛЬНАЯ МОДЕЛЬ ОБРАТНОЙ КОЭФФИЦИЕНТНОЙ ЗАДАЧИ ДЛЯ УРАВНЕНИЯ ПЕРЕНОСА ЗАГРЯЗНЯЮЩИХ ПРИМЕСЕЙ В АТМОСФЕРЕ
© 2005 г. В.И. Наац
The paper «Inverse coefficient problem of equation turbulent diffusion» is concerned with calculation methods for estimating coefficient turbulent diffusion as applied to operative control of the pollutants transfer in the frontier layer of the atmosphere. The new approach for the construction of the calculation models for coefficient turbulent diffusion based on the ideas of the theory inverse problems of atmospheric optics is described in this paper.
Определение коэффициента турбулентной диффузии из уравнения переноса методом обратной задачи
Уравнение переноса загрязняющих примесей в приземном слое атмосферы в условиях турбулентности в направлении координатной оси Ох записывается в виде:
д
q (x, t) + aq( x, t) +—(Vx ( x, t )q( x, t) - K ( x, t )q'( x, t) ) = S (x, t). (1)
dx
В (1) q(x, t) описывает поле концентрации переносимого вещества; Vx(x, t) - скорость ветра в направлении оси Ох; K(x, t) - коэффициент турбулентной диффузии. Функция S(x, t) в правой части уравнения (1) описывает распределенный источник в рассматриваемой задаче диффузионного переноса. В отличие от так называемых прямых задач, заключающихся в определении концентрации q(x, t) [1], обратная коэффициентная задача, связанная с определением поля коэффициента K(x, t), предполагает известными функции q(x, t), Vx(x, t) и S(x, t). Подобные задачи связаны с развитием теории и методов дистанционного зондирования атмосферы на основе радиолокационных наблюдений, а также с использованием средств лазерной локации [2]. Цель настоящей работы - разработка вычислительной модели и её численная реализация, позволяющая оценить значения коэффициента турбулентной диффузии на основе метода обратной задачи применительно к задачам оперативного контроля переноса и содержания загрязняющих аэрозольных примесей в приземном слое атмосферы вблизи промышленных предприятий.
Формальные аспекты решения так называемой обратной задачи для уравнения переноса и построение соответствующего вычислительного метода подробно рассмотрены в [3]. Приведем здесь кратко содержание этого метода. Уравнение (1) преобразуется к виду
J(x, t) - J(xo, t) = Ф(х, t), (2)
где J(x, t) = V(x, t) ■ q(x, t) - K(x, t) ■ q'(x, t), J(x0, t) = V(x0, t) ■ q'(x0, t), x0 = 0,
x
<p(x,t) = j(S(x,t) -q(x',t)-a(t) ■ q(x',t))dx\ V = Vx и K = Kx;
x0
или
K (x, t) ■ q ' ( x, t) - K ( xo, t) ■ q ' ( xo, t) = ^
= V ( x, t ) ■ q( x, t ) - V ( x0, t ) ■ q( x0, t ) - <( x, t). Уравнение (3) позволяет найти значения коэффициентов турбулентной диффузии K(x, t) в точках (x, t) области Q = [0, Х] х [0, Т] при условии, что известны функции q(x, t), q(x,t) = q't(x,t), q'x(x, t), V(x, t), S(x, t) и a(t). Неопределяемыми остаются значения q'x(x, t) в точках (x, t), в которых она обращается в ноль.
Для решения задачи (3) в [3] построен вычислительный метод, основанный на том, что искомая функция K(x, t) в пределах области Q представляется в виде
m
Km (x, t) =Z Ck (t) ■ uk (x), (4)
k=0
где {uk(x)}m - система непрерывных базисных функций на интервале [0, Х]. Считаем далее, что требования к K(x, t) и указанному базису, обеспечивающие сходимость по норме Km(x, t) к K(x, t) во всех точках (x, t), выполнимы. Подставляя (4) в (3), получим
m m
X Ck (t) ■ Uk (x) ■q'(x, t) - 2 Ck (t) ■ Uk (x0) ■ q'(x0, t) =
k=0 k=0 (5)
= V(x, t) ■ q(x, t) - V(x0, t) ■ q(x0, t) - <(x, t). В итоге задача сводится к определению коэффициентов разложения {Ck(t)} из данного функционального уравнения. Применяя далее метод взвешенной невязки и выбирая систему весовых функций {mh Q}m,
m
U Q = [0,X] [1], (5) приводим к системе уравнений относительно коэф-
I=0
фициентов разложения {Ck(t)}, а именно
m
X Ck (t) ■ J (x) ■[Uk (x)q' (x, t) - Uk (x{))q'(x{), t)]dx = k=0 Q',k _ (6) = J a>l (x) ■ [V(x, t)q(x, t) - V(x0, t)q(x0, t)] dx - J соt (x)<(x, t)dx, l = 0, m.
q, Qi
С помощью обозначений
aik (t ) = J Cl ( x) ■[Uk ( x)q ' ( x, t ) - Uk ( x())q ' ( x0, t )] dx,
Qi ,k
b(1) (t) = J (x) ■ [V(x, t)q(x, t) - V(x0, t)q(x0, t)]dx,
o,
(7)
b}2) (t) = J col (x)^>(x, t)dx,
o,
можно представить (6) в виде алгебраической системы уравнений относительно вектора С(() для каждого момента времени t е [0, Т].
А(0 • C(t) = Ьф, (8)
где Ь(t) = Ь(1)(0 - Ь(2)(t).
Таким образом, определение коэффициента турбулентной диффузии методом обратной задачи свелось к решению системы алгебраических уравнений (8) для каждого момента времени t. Формально поставленную выше задачу можно считать решенной. Однако необходимо иметь в виду, что исходное распределение д(х, () определяется в эксперименте по постановке задачи и, следовательно, д(х, () и тем более д(х, t) = д' (х, t), дх(х, Г) известны приближенно. Значит в системе (8) матрица А(() и её правая часть известны с ограниченной точностью. Это лишает возможности увеличивать т в целях уменьшения ошибки аппроксимации функции К(х, Г) рядом Кт(х, (). Другим фактором, ограничивающим эффективность предложенного алгоритма, может являться так называемая слабая обусловленность матрицы А(() для некоторых t е [0, Т]. Эти обстоятельства указывают на то, что в качестве обратной матрицы А~1(() в системе (8) предпочтительно брать обобщенную обратную матрицу и обеспечивающую получение так называемого нормального решения системы (8), определяемую согласно выражению
-1 Т -1 Т
Аа = (А А + а • I) А . Нетрудно видеть, что при а ^ 0 матрица Аа ^ А . При а > 0 нормальное решение Са(() по сравнению с С(^ характеризуется большей устойчивостью к ошибкам как в матрице А((), так и в правой части системы (8). Доказательство этих утверждений можно найти в [4]. В задаче (8) речь идет о последовательности решений {С(Щ, ] = 1,п. Подобная задача сложнее в вычислительном плане, чем задача решения отдельно взятой
алгебраической системы, и является более содержательной. Действитель-
-1
но, появление оператора Аа , зависящего от параметра регуляризации а, связано с вариационными подходами к решению операторных уравнений вида Аф = / а именно с оптимизационной задачей для так называемого сглаживающего функционала [4]
Та(ф) = (Аф - / Аф - / + а(ф, ф). (9)
Уравнение Эйлера для этого функционала имеет вид А*Аф + аф = А*/, откуда имеем
фа = (А*А + а1)-1А */ (10)
где А* - оператор, сопряженный к А. Выражение (10) делает понятным
. -1
появление оператора Аа при решении указанных операторных уравнений. Ясно, что оператор Аа определяет в качестве решения некое параметрическое множество функций {фа}, из которых выбором значения а определяется наиболее приемлемое с практической точки зрения. Квадратичный функционал (ф, ф) в теории некорректных задач принято обозначать 0.(ф) и называть стабилизатором задачи, в то время как функционал (Аф -/ Аф -/) по существу является квадратом невязки исходного операторного уравнения Аф = /. Для сокращения записи последний можно обозначать р2(ф). Это замечание касается функции фа, доставляющей минимум квадратичному функционалу (9). Поскольку (ф, ф)^2 является одновременно нормой в классе функций Ь2, то фа есть решение с минимальной нормой (т.е. фа наименее уклоняется от нуля). Поскольку в решаемой задаче (8) имеем дело с последовательностью {Са(Щ, О = 1,н, то естественно в этом случае потребовать от {Са(у)} наименьшего уклонения не от нуля (нулевого элемента в соответствующем векторном пространстве), а от {Са(О-г)}. Этот способ регуляризации ведет к согласованному функционалу вида
Та (с(О)) = (А(]) ■ С(]) - Ь (]), А(]) • (С(]) - Ь(])) +
+а-р- (С(у) - С (У -1), С (у) - С (у -1)),
где р - некоторый масштабный множитель. Нетрудно показать, что решение Са(у), доставляющее минимум функционалу (11), имеет представление Са(У) = (I - Аа (!) • АО) • Са( - 1) + Аа О • О т.е. искомый вектор СаО может быть определен по рекурсивной вычислительной схеме, в которой
-1 Т -1 Т
Аа (!) = (А (!) • АО + а • р • I) • А (О), и С(0) определяется начальными условиями К(х, 0). Можно указать качественный подход к выбору значения параметра регуляризации а. Допустим, что в уравнении Аф = /правая часть известна с погрешностью а, т.е. определена функция /а, такая, что \/а -/0\\ < а\/0\\, где/0 - это точная правая часть уравнения. При этом естественно ф0, удовлетворяющее уравнению Аф0 = /0, считается точным решением поставленной задачи. При минимизации невязки р2(ф) = (Аф -/а, Аф - /а) на ф е Ф, ее наименьшее значение в лучшем случае достигнет уровня а, в то время как р2(ф0) = 0. В связи с этим естественно ожидать, что величина а • (ф - ф0, ф - ф0) также не должна превышать а2 Откуда получим оценку параметра регуляризации
а2
а <--.
Ь>а-Ф0>\
Если множество исходных функций ф е Ф ограничено константой Мр = вир 1, то последнее неравенство примет вид а < и21М2, т.е. ве-
среФ
личина параметра регуляризации имеет порядок О(о^), где а - относительная ошибка правой части решаемого уравнения. Подобный способ выбора а называется выбором по невязке. Предварительная оценка значения а находится естественно в соответствующих вычислительных экспериментах на моделях.
Последнее замечание следует сделать в связи с использованием в излагаемой здесь теории метода взвешенной невязки. В его рамках можно упростить соответствующие вычисления и одновременно ослабить требование, предъявляемое к исходным данным. Действительно, вычисление элементов главной матрицы А(^, т.е. интегралов (7), можно осуществить с использованием равенства
J с, (x)uk (x)q'(x, t)dx = ( (x)uk (x)q(x, t))
lk - J ((x)uk(x)) q(x,t)dx,
,k O,
где {о1 к = |Х^кг, х\2к Л, 1, к = 0, т - система конечных элементов; 0,к =
= 01 п 0к. Этот прием избавляет от необходимости знать значения д'(х, (), что немаловажно в случае вычислений с экспериментальными данными. Конечно, это ведет к получению так называемых слабых решений, о чем уже говорилось выше.
Вычислительная модель обратной коэффициентной задачи
В соответствии с методом параметризованных моделей, подробно изложенным в [1], для задачи (3) можно ввести нормировочные коэффициенты {а, в, у, £} и нормированные распределения полей и переменных, принимающие значения из интервала [0, 1]. После этого задача (3) может быть преобразована в параметризованную модель:
у • (К(х, t) • д'(х, t) - К?(х0, t) • д'(х0,0) = = р\у (х, 0 • д(х, t) - У(х0, t) • д(х0, t) )-р(х, 0, (12)
р(х,t) = | &(х',^)-дх',t)-а^)• д(х',0)х'.
х0
Далее, в соответствии с (4), искомую функцию К(х, t) в пределах области О = [0,1] представим в виде
л л л m л л
Km (x, tj) =2 Ck (tj) ■ Uk (x), (13)
k=0
где х е О, ] = 1, п, t]■ е [0,1]. Аналогично аппроксимируем остальные поля исходных данных по дискретным измерительным значениям (экспе-
25
риментальным данным). Тогда система (8) получит представление:
2 a,,k (tj )Ck (t.) = ь, (tj), l, k = 0m, j = M,
tfldj) =
k=0
at,k (tj) = ag(tj) - alk (t,),
л (2) л (2) л л (2) л л (1) л (1) л л (1) л
с (Xi,k )uk (Xi,k ) qm (Xi,k, tj ) - с (Xi,k)uk (Xi,k) qm (Xi,k, tj )
- J (x)uk (x) qm (x tj)dx- J cl (x)uk (x) qm (x= tj )dx>
afttij ) =
л (2) л л л л л (1) л л л л
с/ (Xl,k )uk (x0 ) qm (x0 , tj ) - с/ (Xl,k)uk (x0 ) qm (x0, tj )
- J с(x)uk(xo)qm(xo,tj)dx- J с(x)uk(xo)qm(xo,tj)dx,
Alk Alh
b (tj) = b(1)(tj) - bl(2)(tj),
bff>(tj) = £■ J®l(x)■
Al
Vm (x=t] ) qm (x= ^ ) - Vm (x0, ^ ) qm (x0 , tj )
bl(2)(tj) = J с (x)>
d x,
(14)
Sm (x'= tj ) -
qm (x> tj-1) - qm (x> tj ) At
-а (tj ) ■ qm (x = tj )
d x
d x.
Соответствующее выражение для сглаживающего функционала (11) примет вид
та (() 1=2
l=0
k=0
2 (a,,k (? )Ck (tj) - b (tj)
+а p 2
k=0
Ck (t;) - Ck (t;-0
C (tj) = (C (tj)}.
Решая далее задачу минимизации
Та. (c(tj) min,
C (tj)
(15)
(16)
находим искомые коэффициенты Са, ) при соответствующем значении параметра регуляризации а] для каждого фиксированного момента вре-
A
мени tj, j = 1, n. Решение задачи (15), (16) осуществляется при заданных начальных и граничных условиях:
Ck (t0 ) = Ck (0) = К ( xk ,0), k = 0m, j = 0;
Cq (tj ) = К (¿0, tj ) = К (00, tj \ j = 1П; (17)
Cm ( t j ) = К ( Xm , t j ) = К (1, t j X j = ïjl.
В итоге искомое решение обратной коэффициентной задачи (12) по определению коэффициента турбулентной диффузии запишется в соответст-
л л m л л
вии с (13): Km (X, tj ) = 2 Ca k (tj ) ■ uk (x).
k=0 J
Алгоритмизация вычислительной схемы, постановка вычислительного эксперимента, результаты его численной реализации
Для проведения вычислительного эксперимента воспользуемся процедурой генерации исходных данных q(x, t), V(x, t), K(x, t) и a(t), составляющей содержание так называемого тестового примера (блока исходных данных) [1]. Данный блок включает в себя алгоритмы нормирования исходных распределений и переменных задачи, вычисление нормировочных коэффициентов а, в, у и вычисление значений функции источника S(x, t)
с помощью уравнения переноса. В итоге в блоке исходных данных хра-
* л л л
нится «точное решение» KT (x, t) = K K(x, t). Исходные данные представлены массивами дискретных значений {q(i, j)}, {V(xi,tj)}, {K(xi,tj)},
{s(xi,tj)} и {a(j)}, {x(i)}, {t(j)}, i = 0,m, j = 0,n, нормированных на
максимальный элемент каждого из них соответственно.
Исходные данные: Х = 1600 м = 1,6 км; Т = 600 с = 10 мин; q0 = = 0,75 кг/м; V0 = 5 м/с; К0 = 150 м2/с; а0 = 0,025 1/с; m = 25; n = 25. Они используются в блоке исходных данных для генерации значений упомянутых выше массивов. Расчетные значения нормировочных коэффициентов: в = 4,33; у = 6,58 Е - 02; £ = 16,91. Вычисляются диапазоны изменения значений концентрации примесей q(x, t) е [qmm, qmax\ = [0,37; 1,26],
V x е [0, X], V t е [0, T]; поля скорости ветра V(x, t) е [Vmm, Vmax] = [2,64; 11,55], V x е [0, X], V t е [0, T]; коэффициента турбулентной диффузии K(x, t) е [Kmin, Kmax] = [93; 281], V x е [0, X], V t е [0, T]; интенсивности источника Sx, t) е [Smin, Smax] = [0; 3,57 E - 02], кг/м • с, V x е [0, X],
V t е [0, T], a(x, t) е [amin, amax] = [2,26 Е - 02; 2,75 Е - 02], V t е [0, T], a (t) = T ■ait ).
Для перехода в вычислительной модели к методу взвешенной невязки необходимо выполнить аппроксимацию исходных данных по их дискретным измерениям и соответствующих производных:
Vm (X, j) = 2 Vi, j ui (x), Km (X, j) = 2 j u, (x), Sm (X, j) = 2 u, (x),
i=0 i=0 i=0
qm(x j) = 2 qi j u(x), qm(x, j) =
i=0
qm (x, j) - qm (x, j -')
At ;
(18)
+2 q(i + 1j)-^-1j)Uk(x) +
i=1 2Ax
qm(x, j) = q(1, j)A-q(0,j) u0(x) +
Ax
q(m, j) - q(m -1, j)
(19)
Ax
Дx), V x e[0,1].
Выражения (18) и (19), аппроксимирующие соответствующие производные, приводились и исследовались ранее в [1]. В этом случае выражение (2) получит следующее представление:
J (x, j) - J (0, j) = < (x, j ), J(x, j) = ß ■ Vm (x, j') ■ qm (x, j') - Y ■ Km (^ j) ■ q'm (^ Л
J(0, j) = ß ■ Vm (0, j) ■ qm (0, j) - Y ■ Km (0, j) ■ qm (0, j),
i(
<(x, j) = J( ■ ^^ (x', j) - qm (x', j) - а j ■ qm (x', j))d x'.
(20)
(21)
Интеграл в выражении (21) будем вычислять по методу Гаусса. Отклонение левой от правой части первого уравнения (20) обозначим
сг, =
1
2 21 J(x,, j)- J(x,, j) -<p(x,, j)
1
2\/2
п • (т +1) ]=1 г=01
Величина Ст1 будет включать в себя не только погрешность вычисления интеграла в (21), но также ошибки аппроксимации исходных данных и,
особенно, ошибки аппроксимации производных {д'х(г, ])}, (г, ])} рядами д'т (х, ]) и ¿¡т (х, ]) соответственно. Ошибку аппроксимации будем определять по обобщенной формуле вида:
1
п т -у
2 2(f (x, j) - fm (x , j) )2
^ (п + 1) • (т + 1) ]=0 ¿=0 Далее предположим, что экспериментальные значения поля концентрации примесей {д(г, ])} измерены с погрешностью 3, т.е. нам известны
(л (5)
значения <! ч (г, Л) ^ (предполагается, что все прочие исходные данные
измерены точно). Тогда в расчетных формулах (20), (21) необходимо зал л (5)
менить qi ■ на qi ^ . В итоге мы придем к уравнению
„ л л (5) „ л (5) л л (5)
3(X, Л, qг,] ) - 3(0, Л, qг,] ) = Ф(X, Л, qг,] ). (22)
Соответственно для оценки погрешности будем иметь выражение
1/
a2(S)--
1
(S) „ Л Л (S) Л л Л (S) ^ J (Xi, J, qt j ) - J (Xi, J, qi j ) \-ф( Xi, J, q. j )
2
(23)
п ■ (т +1) г=0 ч
\ / Процесс появления случайных ошибок д при измерении «экспериментальных» значений q(г, Л) моделируем следующим образом. Заменим
л л (5)
q(г,Л) на приближенное значение q (г, Л) по формуле:
л (5) л
q (г, Л) = ^(г, Л) 41 + ©¡5),
где 0/ - случайные числа, равномерно распределенные на интервале (-1, 1), генерируемые датчиком случайных чисел.
На рис. 1 показано поведение погрешностей 0\(Л), а2(Л, д) при д = 0,01 с течением времени. Точность аппроксимации полей исходных данных достаточно высокая благодаря свойствам базисной функции {ик(х)}т (табл. 1) [1]. Основной вклад в итоговое значение погрешностей делают ошибки аппроксимации производных. Очевидно, что необходимо по возможности исключать производные из расчетных моделей. Поведение кривой а2(], д) характерно для вычислительного процесса, связанного с решением некорректных задач, каковой и является обратная коэффициентная задача (3). Для количественной характеристики влияния величины д на погрешность
л л (5)
(ч, q )
а2(Л, д) можно ввести «коэффициент усиления ошибки» п =-,
значения которого приведены в табл. 1.
Полагаем далее неизвестными значения коэффициента турбулентной
диффузии
{ г
Л У, кроме значений, определяющих начальные и граничные
условия: К(хк ,0), к = 0, т, Л = 0; К(Х0, г^) = К(0, г^), К(Хт, г}-) =
= К (1, г л ), ] = 1, п. Вычисление приближенных значений {Кг Л} осуществляется, как уже говорилось выше, путем решения обратной коэффициентной задачи (17).
1,80E-01 -
1,60E-01 -
1,40E-01 -
1,20E-01 -
| 1,00E-01
■j? 8,00E-02 -tfl
6,00E-02 -I
4,00E-02 -2,00E-02 -
0,00E+00 -.........................J
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
a1 ..-■-.- a2
Рис. 1. Распределение погрешностей oj и o2(j, ö) и при ö = 0,01 по времени
Таблица 1
Зависимость погрешности вычислений а от ошибки аппроксимации полей исходных данных и погрешности измерений поля концентрации примесей а
Погрешность аппроксимации исходных данных и производных Ошибка в измерениях поля концентрации примесей q
S = 0 (0 %) S = 0,01 (1 %) S = 0,1 (10 %)
(V v ) 1,022Е-08 1,022Е-08 1,022Е-08
^K (K, Kn ) 1,110Е-08 1,110Е-08 1,110Е-08
(S, sm ) 9,280Е-09 9,280Е-09 9,280Е-09
fe qm ) 1,093Е-08 1,093Е-08 1,084Е-08
(.q', q'm ) 7,815Е-02 8,869Е-02 0,535
fe qm ) 0,106 0,133 0,963
Значение о 01 = 7,716Е-02 02 = 8,73 8Е-02 02 = 0,374
л л (ß (q, q ) 0 3,66Е-03 3,66Е-02
Коэффициент усиления ошибки n - 0,042 0,097
Алгоритм решения обратной коэффициентной задачи. Выбор параметра регуляризации по невязке:
1)] = 1;
2) берется конечный отрезок монотонной последовательности чисел а0(]),а1(]),а2(]),...,ап(]), например, отрезок геометрической прогрессии а(]) = «„(]■) • Гк, к = 0 < г < 1 [4];
3) к = 0, задаем начальные значения а0 и г;
4) вычисляем значение параметра регуляризации ak (j) = а0 (j) ■ rk;
5) определяется сглаживающий функционал Tät(j) (C(j)) в соответствии с (15), в котором элементы матрицы A(j) = {alk(j)} и вектора b(j) = {bl (j)} вычисляются по формулам (14) и с учетом начальных и граничных условий (17);
6) для каждого значения äk (j) находятся значения коэффициентов C äk(j)(j) ■> минимизирующих функционал (15). Коэффициенты
Cäk(j)(j) = {Cät(j),i (j)}, l = m-1 определяются путем минимизации
функционала Ta, (C( j)min с помощью «эвристического симплексно-
C (tj )
го метода» [5];
7) вычисляются значения искомой функции Km (x, j) по формуле (18);
8) с помощью (22), (23) вычисляется величина погрешности c3(S);
9) согласно алгоритму нахождения параметра регуляризации ak (j) по невязке, вычисляется значение
(KT , Km ) =
/ \У
I 1 n m i л \2 к 2
-- ZZ(t (X , j) - Km (X, j) )
m +1) j=1 i=o \ !
п ■ (т +1) ] =1 г=о\
л ~ *
10) Если сК (Кт, Кт) < е(8) и/или а3(д) ~ а2(д), то выбираем а (]) =
„ ~ * л
= ак (]) и вычисляем соответствующие значения Кт (х, у). В противном случае к = к + 1 и возвращаемся на шаг 4.
11) ] = ] + 1; если ] < п, то переходим на шаг 2, иначе - конец алгоритма.
Для реализации данного алгоритма зададим значения д = 0,01, ао(]) = 0,1 У] = 1,п, г = 0,75, е = 0,1. Ниже приведены результаты расчетов. На рис. 2 приведены распределения по оси Ох при фиксированном моменте времени нормированных значений «точного решения»
Кт (х, ] = 5) и приближенного, вычисленного по выше приведенному алгоритму Кт (х, ] = 5).
На рис. 3 приведено пространственно-временное распределение поля турбулентности - его «точного» и «модельного» решений, их взаимное
расположение. Профиль распределения Кт (хг, (у) в целом повторяет профиль распределения Кт (хг, (у). В каждом случае приведены значения
отклонения «точного» от приближенного решения. На рис. 4 отражено временное распределение погрешности вычисления коэффициента турбулентной диффузии. Со временем процесс расходится, что характерно для нестационарных процессов. На рис. 5 показаны кривые а3(3) и а2(3). Кривая а3(3) повторяет ход а2(3), незначительно уклоняясь от нее.
ktd - - - - ktdmod
Рис. 2. Графики Kt (x, j = 5) (на рисунке обозначено «ktd») и Km (x, j = 5) (на рисунке обозначено «ktdmod»), ак (Kt , Km, j = 5) = 2,91E — 02
Рис. 3. Взаимное расположение полей Kt (x,, tj) и Km (xt, tj), i = 07m, j = 1"n, nK (Kt , K m) = 8,73E - 02
1 2 3 4 5 6 7
9 10 11 12 13 14 15
Рис. 4. Изменение ак (Kt , Km, j) по времени
1,20E-01
1,00E-01
8,00E-02
6,00E-02 --
4,00E-02
2,00E-02
0,00E+00
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
------- a2 a3
Рис. 5. Изменение а3(3) и а2(3) по времени
В табл. 2 приведены возможные значения параметра регуляризации ак (]) для некоторых фиксированных моментов времени и соответствующие им значения погрешности аК (Кт, Кт). Для каждого момента времени может быть выбрано свое оптимальное значение а*( ]), при котором аК (Кт, Кт) минимально.
Таблица 2
Выбор параметра регуляризации по невязке
ak(j) ®k (Kт, Km )
j = 1 j = 6 j = 8
7,50Е-02 1,95Е-02 2,21Е-02 5,91Е-02
5,63Е-02 1,37Е-02 1,97Е-02 5,78Е-02
4,22Е-02 8,39Е-03 1,83Е-02 4,87Е-02
3,16Е-02 2,85Е-02 2,46Е-02 4,19Е-02
Литература
1. Семенчин Е.А., Наац В.И., Наац И.Э. Математическое моделирование нестационарного переноса примеси в пограничном слое атмосферы. М., 2003.
2. Зуев В.Е., Наац И. Э. Обратные задачи оптики атмосферы. Л1990.
3. Наац В.И. Обратная коэффициентная задача для уравнения турбулентной диффузии // Вестн. СевКавГТУ. 2003. № 1(7). С. 89-94.
4. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М., 1979.
5. Дэннис Дж., Шнабель Р. Численные методы безусловной оптимизации и решения нелинейных уравнений. М., 1988.
Северо-Кавказский государственный технический университет,
г. Ставрополь 14 февраля 2005 г.
