CC BY
29
УДК 517.2
ИТЕРАЦИОННЫЕ МЕТОДЫ ЧИСЛЕННОГО РЕШЕНИЯ ЗАДАЧ ПЕРЕНОСА НА ОСНОВЕ ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ
© 2004 г. Н.И. Каргин, В.И. Наац
The method for the numerical solution unsteadily-state equation of source transferring, which is based on iteration calculation algorithm, is developed and research. The analytical model and results of numerical research calculation algorithms are receiving.
В статье рассматривается и исследуется метод численного решения нестационарного уравнения переноса на основе итерационных вычислительных алгоритмов, построение которых осуществляется с использованием интегральных представлений для решений уравнений параболического типа. Основная идея метода достаточно подробно изложена в [1]. Она заключается в том, что для нестационарного уравнения переноса субстанции в пределах пограничного слоя атмосферы можно построить эквивалентное ему интегральное уравнение Вольтерра II рода. Далее решающий алгоритм для определения поля концентрации (/(I, х) в пределах замкнутой области I) можно строить на основе метода последовательных приближений, сходящегося при выполнении условия ограниченности ядра. В работе [1] приведены условия сходимости метода. Цель данной работы - построение вычислительных алгоритмов метода и их численное исследование. При этом предлагаемый в данной работе подход рассматривается в двух вариантах - первой и второй итерационных схем. Ниже приводятся соответствующие аналитические построения и результаты численных исследований алгоритмов метода.
Метод интегральных уравнений.
Первая и вторая итерационные схемы
Уравнение нестационарного переноса примесей имеет вид [2]:
в котором с/(1, х) - концентрация примесей, имеющихся в атмосфере (в приземном слое атмосферы) в момент времени t в точке с координатой х; I '(/, х), K(t, х) и S(t, х) - соответственно функции скорости ветра, турбулентной диффузии и источник примесей; a{t) - параметр, описывающий процессы вывода или привнесения примесей за счет химических реакций.
Если уравнение (1) переписать в виде
q{t g, х) = tp(x) , q(t, x0) = (f>(t) , q(t, X) = ^(/), / е[о,г], xe[o,j], (2)
дК (t, л) 1 dq(t, х)
дх
дх
(3)
то выражение (3) можно рассматривать как линейное дифференциальное уравнение первого порядка по переменной t при фиксированном значении переменной х. Известно, что для подобного вида уравнений существует аналитическое решение, которое для уравнения (3) записывается так:
\ ‘ 5V(t' х) ]
q(t, х) = ехр < - ¡(а ((')-{-----:--)dt' ^ х
I 'О дх \ (4)
' dV(t" х) ]
q(t0,x)~ J x) - S(t', xj) ■ exp < J (a (/") н---------:—)dt">dt'
>0 [»o 8x )
где w{t,x) = lv(t,x)- ЁМЩMhll_ {к(/>.
{ дх J дх дх2
Обозначим в выражении (4)
K(t,t',x) = exp j- И a (t") + 8V(t8"x’X))dt" j > С5)
(p(t,x)= q(t0,x) ■ exp J — J(a(i') + ЁИИ.liL^-)dt' > + J S (t', x) ■ К (t ,t', x)dt'.
I <o dx \ ,0
Окончательно получаем:
q(t,x) = <p(t, x) - J К (f, f' ,х)цг (f' ,x)dt'.
(7)
г0
Выражение (7) можно рассматривать как интегральное уравнение относительно функции (/(I, х). При построении метода численного решения (7) относительно искомой функции q(t, х), удобно вводить функции q(l х). где х играет роль параметра. Тогда для кавдого фиксированного х интегральное уравнение относительно q(t|х) является интегральным уравнением Вольтерра П рода, численное решение которого можно осуществить методом последовательных приближений. Данный метод реализуется в виде следующей итерационной схемы:
q<v>(t,x) = ip(t,x) - j K(t,t' ,x)i//(t' ,x,q^v~r> (t'\x))dt', (8)
>o
где v - номер итерации. Для вычисления значений выражения i//(i. х, q(t, х)) на множестве функций q(t\x) необходима аппроксимация производных
дх ’
д . С этой целью введем равномерную сетку узлов {(х,)} на интервале
дх2
х е [0, X], / = 0, т +1. Далее будем полагать, что на каждом шаге итерации V вычисляется последовательность 1 (I) = = q('v> (;|х/;), к = 0, т +1. Эго позволяет для каждого у воспользоваться следующими конечно-разностными
аппроксимациями для производных
дх
dq(t\x) и d2q(t\x) ;
дх
3q(j\x)
дх
d2q(t\x)
(у)
Дх
дх1
Чк+і ї)(0~2‘іі'/)(0 +д(к\(0
Ах2
где к = 1, т. При этом ^(О = ?(/,|х0 = 0) и <7^(0 = <?(^т+і = X) определены граничными условиями (2) для исходного уравнения (1) (задача Коши). В итоге имеем следующую формулу для расчета значений функции і//(/, х, х)) В точке X = Хк.
у X,, 9 ^\г\хк ))= (к (ф*) -
V )
9(у)(ф*) - Ч<У>(і\хк_і)
V_______!___________________
Ах
дІ>~1)о|хі+і)-2 д'>)(фі) + ^{у)^\хк-1)
(9)
- К (t
хк)х-
Ax2
Таким образом, вычислительный процесс для итерационной схемы (8), (9) полностью определен и его результатом является последовательность функций
(У)
«},
v = 1,2,..., к = \,т .
Вычислительная схема (5), (6), (8), (9) не единственна и возможен другой вариант. Действительно, если выражение (1) переписать по-другому, а именно:
— \уЦ,х^(!,х)-КЦ,х) д<}{-і’х')\ = £(/,*)- - а(/)?(/,*) (Ю)
дх I дх \ д(
и ввести обозначение
J(J,x,q(J,x)) = Г х) - К (1,х) 3 д ^’ ’ * '> ’
О X
то, интегрируя обе части выражения (10) с учетом обозначения (11), получим уравнение вида:
X £) X
J—J(t,x,q(t,x))dx = J
х0 дх Хо
откуда имеем
J (t, х, q(t, х)) - J (t, Xq ,q(t, x0 )) — J
S{t,x)~ ,X^ - a(t)q(t, x)
dt
S(t,x)~ 5g(f’x) - a(t)q(t,x) dt
dx.
Согласно (11), запишем выражение для определения J(t. х0, q(t, х0)): J(t,x0,q(t,x0)) = V (t, х0 )q(t, х0 ) - К (t, х0 ) Х° -* = J0 (f ).
OX
Перепишем (12) с учетом (11) следующим образом:
дq(t, х) V (?, х)
дх
1
К{1,х)
К{1,х)
У о (О + I
q (?, х) +
5(1,х)- ад(?’х) -а(1)д(1,х) 3?
сЪс
(13)
= 0.
Выражение (13) можно рассматривать как линейное дифференциальное уравнение первого порядка по переменной х при фиксированном значении переменной Для этого уравнения аналитическое решение имеет следующий вид:
V (t х') I * 1
q(t,x) = ехр - | ———с1х' ^ • [д(?, х0) - |
+ \
К (?, х')
.10(1) +
(14)
Введем обозначения:
х0 ) ■ ехр -{ / ^ ' с1х '
с о ^ * )
К(/, х, х') = ^ ~ехР ] | г-<з?х " ¡-Л",
КЦ,х')
¥(1,х') = J0(t)+ |
5 (/, х") - ^^,х ) _ a(t)q(t,x") д(
с1х".
В этих обозначениях соотношение (14) перепишем в виде:
q(t,x) = <рЦ,х)~ \ К Ц, х, х')ц/ Ц, х')йх'.
(15)
(16)
Уравнение (16), также как и (7), есть уравнение Вольтерра II рода, для численного решения которого естественно воспользоваться методом последовательных приближений. Переменная / играет в уравнении (16) роль параметра, поэтому удобно ввести параметрическое семейство функций ц(1 х). где I е [0, 7], х е |0. А’|. Итерационная схема метода последовательных приближений примет вид:
ql'v\x\t) = <р(х^) - ¡К (х, х'|?)(// ^_1^(х',?,^^_1^(х'|?))^/х', (17)
*0
где V- номер итерации, г = 0. 1.2.3... Нулевое приближение (/'" (х. I) = = <р(х, /). Далее примем предположение о том, что параметр / пробегает значения I/. / = 0. п. для которых интервал \/ = ] достаточно мал в том смысле,
что конечная разность
дОМр-дОМу-О
А?
вполне приемлемо оценивает вели-
да(х,1) „ , , (у\, . чч
чину -----------. В этом случае вычисление значении у/(х. I. ц (х. I)) можно
3/
оценивать по формуле:
V (у)(х', *_,■) = J0 а ]) +
х' Ч ^ ^ (х ",! ,) - д ^ -1 (х ", г ■_ 1 ) , ,
+ | 5(х", I¡ )----------------—-----------------а(^ )г( }(х", I ¡)
(18)
с1х''
В выражении (18) считаем, что к моменту времени I,. когда последовательно по V вычисляется значение (¡ " (х,/,). значение (¡ " (х, /,ч) уже известно.
Параметризованные интегральные модели и вычислительные итерационные схемы
Выполняя процедуру нормирования переменных и величин [2], входящих в (5), (6), (8), (9), получим основные расчетные формулы первого итерационного метода:
г _ л^-1)
q (?, х) = <р(Х, х) - I К (1,1', х)ц/(Х:, х, q (?'|х))Л',
~ * л д У(1" х")
К (?, х) = ехр | (а (?") + /?------------------------------)Л" >,
\ ? дх \
а г а д V а' х)
<Р((,х) = q(tо,х)-ехр }(«(?')+/?----------------------)^'\ +
t дх \
(19)
¥
Г (Г
*к)~ К(1\хк-\)
Ах
л(-)
<? (1
:к)~Ч ц\хк-1)
Ах
-у-К(1
хк)х'
л(^-1)
<? (Г хА + 1)-29 (?|хА-1)
Ах"
где а(0 = Та((), р = ^ ^ - к т ? = хе[о,1],
X х2 д* 7 К*
л 1 л 1 ,
, х = —, I = — .В выражениях (19) подразумевается, что х = х . 1 = 1.
X т
Значения функций с{ , V,К ,8 - нормированы и меняются в диапазоне [0,1].
л(У) ~ л(у)
Заменим функции д (/, х), К(/, х), (р (/, х) и (//(/, хк, д (ф&)) сеточными
Л1» Л^)
функциями д (1],хк) = Ч]к, К(1],11,хк) = К]1к, ср(}? ,Хк) = ср,-к И |//(г)(^,хк) = |//д;, где (/,.X/.)- внутренние узлы равномерной сетки, у = 1, и,
к = 1,да , Д/ = —-— , Дх = —-—, X = -—г ; / ,■ = у • Д/, У = 0,п ; хк = к - Дх, и +1 да +1 Дх
------- л(г)
к = 0,да +1. Начальными условиями определяются значения ц = 0, хк).
------- /у(г)
к = 0, да +1; граничными условиями задаются значения ц (/,, х0 = 0) и
Л(У) ----
д = 1), у = 0,/?. Тогда выражения (19) получат следующее представ-
ление:
'=[од]
л(°) _ ___
где v = 0,1,2, q д = (p jk , / = 1. п. к = 1. in , интегралы, входящие в выражения (19), приближенно заменены интегральными суммами с квадратурными коэффициентами {ю,}, i = 0. /. j = \,п .
Преобразуем (20) следующим образом:
К jk, (Л) = ехр |- А ■ £ ®, ■ Дх ■ (а Дх + Р ■ (V ¡к - V 1гк- 1))
Л Г j Л АЛ |
<Pjk W = Уок'ехР 1 - ^ • S ®г • («г • Ах + /в ■ Ах -(V гк ~ V г,к-\)) +
j л _
+ £ • Л • X a>iSik К пк Ах і = 0
ИГ0 = —
^ Дх:
.. (V - 1) „(V-1)
/3 ■ &х -V it - Г ■ (К ik - К /д-l) ■ q ik - quк_
(21)
Ax"
-К1к-
лС1-1-1) a(v_1) л^-1)^
с1і,к + \-2с1ік +с1г,к-\
Ах
ч(*0
где
я =
¿1 X
Выражения (21) составляют параметризованную вычислительную схему первого итерационного метода интегральных уравнений. Роль параметров в ней
играют а (/), Д у, с. X. Т и параметр сетки -Л =
At
Аг2
, определяющий взаимо-
связь между величинами Ах и А/. Должно выполняться условие 0 < Я < 1.
Выполним аналогичные построения для второй итерационной схемы метода. Выполняя процедуру нормирования величин, входящих в выражения (15-18), получим следующее представление вычислительной схемы:
9?(х|г) = q* q(x0 |i) • ехр К (х, x'|f)
____і______expkj L^LdA
К * К (х'|t) [ 7 *'К (х"|t) J
* л * л * л д(а q(t,Xn))
J(t, x0,q(t, х0У) = V V(t,x0)q q(t,x0)-K K(t,x0)----------——-----
X ■ ox
q Jo (O
(22)
л I дq (х"|0 л л , |/)------------—-----'■--a(t)q(x \0
с1х"
д *■ \х|?) = д (р(х^)-д "¡К (?, х, х')цг ^ ^ (х'|?, д ^ ^ (х'|?)) • X ■ с1х' '
<3'^'|<3'(х|?) = ^)(х|?)- \ К х, х')уг 1 -1 (х , «з' (х'У)) ■ X ■ с1х' ’
х0
Запишем сеточную модель итерационной вычислительной схемы (22):
<Рц = с1}, о’ехР
Р * V ]к
— X —^а>кАх
Гк-°к]к
К(1],хг,хк) =-----^---ехрX х^ = К]гк
К*К]к \Г1 = кк]I
■ (у -1) _ 7 (0)
X *
J)'" + ^гY. (-Ял-¡1= о
Л (V - 1)
ч и 1-1.1
Д I
а 1 1 ¡1
■ ю I Дл
(23)
(24)
7(0) тл*т) (V —1) ^
-/) =¥ ¥ ]0 ч)о ~к к ]0
л(г)
л(°)
(1-^2)? (^■^1)~4'1 д(.^,х0)
X -(2 Ах)
* ~ -0-1) . Чр =<Рр- Е кркф)к ®ккх-. к = 0
(25)
(26)
л(°) _ _ лМ ЛМ
где у = 0, 1,2, <7-г- =<Рр, 7=1,и, 1=\т, #-0 и д0г - определяется из
начальных и граничных условий. В выражении (25) для аппроксимации произ-
*Л
„ 5(<7 <7(/,х0))
ВОДНОЙ --------------- применяется ПОДХОД, В котором ¡-] И ¡_2 - числовые коэф-
X-дх
фициенты И с”, = 1 - с 2.
Преобразуем выражения для функций 1//д 1; (24) и (26) таким обра-
где ^(v-1) =
£ ' S ji At —
( rSv-l'i
Ifl -
Vj-U
л aO-D ■a j q jj A t
(27)
'■{V) >■ -( ji ~ Ф ji ~ 2 jik J k = 0
a>k Ax ■
где
Ai'
Ai
В итоге вычислительный алгоритм (23), (25) и (27) представляет собой параметрическую вычислительную модель, в которой роль параметров также
играют а (/), Д у; с Л'. 7: и параметр сетки - „ =
Ах
Ai
определяющий взаимо-
связь мслсду величинами Дх и А/. Заметим, что параметр , где Я - анало-
гичный параметр в первой итерационной схеме. Ясно, что должно выполняться условие О
Численное исследование итерационных алгоритмов
В работе [2] приведено описание тестовых примеров, предназначенных для проведения вычислительного эксперимента. Тестовые примеры позволяют генерировать значения массивов исходных данных и точного решения: <7,, = <7(/„ х,), Vji = х,), X,, = Щ, х,), ^ х,), а, = а(/,), ] = 0,п, / = 0, да +1
при заданных значениях переменных X, ^ ^сь ^Г0, а0. Для первой итерационной вычислительной схемы задаются Дх и 0<1<1, далее вычисляется значение 4 = л- Дх2. Для второй итерационной вычислительной схемы задаются значения \/ и 0 < и < 1. затем вычисляется значение Дх = у[/й~М . После этого опре-
X 1 1
деляются значения Т « —, да + 1 =------ь 1; п =---ь 1. Формируются масси-
¥0 Дх А1
вы
j = О, п , где t j = j • At ;
/ = 0, да +1, где хг = / • Дх. Вычис-
ляются значения нормировочных коэффициентов a(t j) = Т ■ a (t j) ,
j = О , п
г = ■
К т
= s т . Выполняется процедура нормирования пере-
менных и значений массивов. Считаем, что цт (/х,) = ц(1. х,). / = 0.п.
/ = 0, да +1 - точное значение искомого решения нам известно.
Далее программно реализуются первая итерационная схема метода интегральных уравнений (21), а затем вторая - (23-27). При этом теперь полагается искомое решение д (/;) не известным. На кавдом шаге итерации V проверяется условие сходимости:
1
S S
« ■ т к = 1 г= 1
Л<>') л O'-D
jk - ч jk
(28)
если /)<£■, то вычисления заканчиваются, и мы получаем а = а(1 1 - искомое
решение. В противном случае переходим к следующему шагу итерации. По окончании выполняем сравнение полученного решения с точным решением:
£ £ \q jk ~ 'я jk
(29)
п -т к = \ /=11
Для первой итерационной вычислительной схемы расчеты проводились при следующих значениях исходных данных: Л'= 900м, Го= 15м/с, /(,, = 40 м/с2. а0 = 0,5 l/c, Yr = 0.1. л= 1/2. При этих значениях вычислены \/ = 0,005, п = 200, т + 1 = 10. Для второй вычислительной схемы исходные данные задавались такими: Л'= 10 м, Г,, = 1 м/с, К0 = 100 м/с2, а0 = 0,5 l/c, At= 0,02,
/î = 1/800 = 1.25/-.’ - 03 и получены значения переменных Лг = 0,005, п = 50, т + 1 = 200.
На рис. 1 и 2 показано графическое представление точного решения и искомого решения, вычисленного по первой итерационной схеме. Графическое представление искомого решения q = ^ |, полученного по второй итера-
ционной схеме, выглядит аналогично. При этом значение cr-отклонения точного решения от расчетного значения (29) составило 2,0Е-02 для первой итерационной схемы и 8.85Е-03 - для второй вычислительной схемы.
Рис. 1. Графическое представление точного решения qT =
Рис. 2. Графическое представление решения, полученного по первому алгоритму итерационного метода q =
Для метода интегральных уравнений, соответствующего первой и второй итерационным вычислительным схемам, проведено численное исследование сходимости. В ходе проведения вычислительного эксперимента было установлено, что итерационные алгоритмы сходятся достаточно быстро. Процесс сходимости отражен на рис. 3. На рисунке показано как меняется значение величины р (28) для каждой следующей итерации V.
Существенное влияние на вычислительный процесс оказывает выбор значений исходных данных и количество узлов сетки. В ходе численных исследований было установлено, что для метода интегральных уравнений, соответствующего первой итерационной схеме, необходимо при задании исходных данных выполнять условие Т<Л' например задавать значение Л'~ I \,'Г. где Г,, > 10. Значение коэффициента турбулентной диффузии не должно быть слишком большим, например, Ки<50. Таким образом, подразумевается, что перенос примесей осуществляется на большие расстояния в основном за счет ветра. Для второго варианта метода интегральных уравнений должны соблюдаться обратные установки, т.е. Л'< Т. например, Л'~ V,-,Т. где Го < 1. а значение коэффициента турбулентной диффузии должно быть достаточно большим, например, К0 >100. Здесь полагается, что перенос примесей идет на небольших расстояниях, достаточно медленно, практически при отсутствии ветра. Перенос происходит за счет турбулентности (табл. 1).
Рис. 3. Сходимость итерационного процесса, соответствующего первой и второй вычислительным схемам метода интегральных уравнений соответственно —го\ = р\ и го2 = Ръ иш = V — число итераций
Таблица 1
Результаты расчетов искомого решения уравнения переноса при соответствующих исходных данных притом, что общее количество узлов сетки составило 5000 для первого и второго итерационных алгоритмов
Название метода X, м т, с Vo, м/с Ко, м2/с а
Метод интегральных уравнений (первая итерационная схема) 900 60 15 40 1,46Е-02
Метод интегральных уравнений (вторая итерационная схема) 10 20 1 100 1,45Е-02
Для метода интегральных уравнений построены их параметризованные модели и соответствующие им алгоритмы. Каждый вычислительный алгоритм в качестве одного из параметров содержит параметр сетки а г или
Ах2
Л 2
и = . В ходе проведения вычислительного эксперимента установлены
А/
числовые значения параметров, при которых обеспечивается приемлемая точность расчетов (табл. 2). Так для первой итерационной схемы-
з Г1 1
А £
6 4
_1______1_
400 ’ 200
= [о 75- 0 16 ]’ ДЛЯ второй итерационной схемы-= [о,002 ; 0,005 ] • В этом случае общее количество узлов сетки
ц е
будет 4-6 и 5-8 тысяч узлов соответственно для каждого метода.
Таблица 2
Влияние значений параметров X п на размерность сетки узлов и точность расчетов а для первой и второй итерационных вычислительных схем
Первая итерационная схема m = 10 Вторая итерационная схема и = 50
/1 п а ß m а
1/2 200 2,00Е-02 1/200 100 1,45Е-02
1/4 400 1,50Е-02 1/450 150 1,07Е-02
1/6 600 1,43Е-02 1/800 250 8,85Е-03
1/8 800 1,41Е-02 1/1250 250 7,71Е-03
Оба вычислительных алгоритма исследовались на устойчивость к случайным погрешностям правой части исходного уравнения переноса (1)-(2) следующим образом. Пусть вместо точного значения правой части (1)
5<°> исходного уравнения известно его приближенное значение с по-
грешностью 8, т.е. имеем
8($) = |, Величину погрешности д определим
так:
5(<у) - 5(0)
8 =
к(°)
где
(30)
=5,(0)(1 + ©е).
(-) = |(-) /,). ] = 0, п , / = 0. т + \. - случайные числа, равномерно распределен-
ные на интервале [-1,1] (генерируются датчиком случайных чисел), 0 < /;< 1. Введем в рассмотрение коэффициент усиления ошибки ц, характеризующий
На рис. 4 показаны графики роста коэффициента усиления ошибки //(31) с увеличением погрешности правой части исходного уравнения д (30) для каждого метода: щ - для первой итерационной схемы метода интегральных уравнений и ц2 — Для второй итерационной схемы метода интегральных уравнений. Рис. 4 показывает, что более устойчивой оказалась вторая итерационная схема метода интегральных уравнений.
В целом можно утвсрлсштъ. что рассмотренные вычислительные методы и алгоритмы обладают свойствами сходимости и устойчивости и позволяют получать решения с достаточной точностью.
степень отклонения получаемого решения в ходе реализации вычислительного алгоритма от точного значения д(0), который определим следующим образом:
где
(31)
Рис. 4. Зависимость коэффициента усиления ошибки метода eta = ц от величины погрешности правой части исходного уравнения d = S
Литература
1. Насщ Н.Э., Насщ B.II. II Сб. науч. трудов СКГТУ, сер. «Физ.-хим.». Вып. 6. Ставрополь, 2002. С. 99-101.
2. Семенчпн Е.А, Насщ B.II., Насщ Н.Э. Мат. моделирование нестационарного переноса примеси в пограничном слое атмосферы. М., 2003.
Северо-Кавказский государственный технический университет 10 ноября 2003 г.
