CC BY
12
УДК 519.5
ВЫЧИСЛИТЕЛЬНАЯ МОДЕЛЬ НЕСТАЦИОНАРНОГО УРАВНЕНИЯ ПЕРЕНОСА ПРИМЕСЕЙ НА ОСНОВЕ МЕТОДА НАИМЕНЬШИХ КВАДРАТОВ И АЛГОРИТМОВ АППРОКСИМАЦИИ
© 2004 г. Н.И. Каргин, В.И. Наац
The construction of calculation model for the unsteadily-state equation of source transferring in the frontier layer of the atmosphere, which is based on method of approximation desired solution using Bernstein polynomial, minimization technique of discrepancy function, least-squares method, to carry out. The analytical model and results of numerical research calculation algorithm are receiving.
Выполняется построение вычислительной модели одномерного нестационарного уравнения переноса примесей, распространяющихся в приземном слое атмосферы в условиях турбулентности. Уравнение представляет собой параболическое дифференциальное уравнение в частных производных, известное под названием уравнения теплопередачи, или уравнения диффузии, или уравнения массопереноса. Особенность задачи, рассматриваемой в данной работе, - исходное уравнение переноса вместо постоянных коэффициентов содержит функции исходных данных, значения которых меняются в каждый момент времени в каяедой точке пространства. Известно, что не существует аналитических решений подобных задач, и поэтому требуется построение вычислительной модели. Известны различные подходы к построению вычислительных моделей данного уравнения. Один из них предлагается и исследуется в данной работе.
1. Основные подходы к построению вычислительной модели. Рассмотрим нормированный вариант нестационарного уравнения переноса примесей с начальными и граничными условиями [1]:
АДА
~У— K(t,x)
ОХ
А
дх у
(1)
V
АА А А А А А А __ А А А А _А
А
q(t = 0,х) = q0(x), q(t,x = 0) = q0(t),q(t,х = 1) = q^(t),
A
A
t e [0,1],x e [0,1],
где р, у, с- нормировочные коэффициенты. В уравнении (1)- (2) нормированные значения функций: q(t,'х) - концентрации примесей; V(t,x) - скорости ветра; K(t,x) - коэффициента турбулентной диффузии; S(t,x) - источника примесей; a(i) - коэффициента, определяющего степень привнесения или вывода примесей из атмосферы; q, V, К, S. а принимают значения из интервала [0,1]. Полагаем, что функции q(t ,'х), V(t,x). K(i,x), S(i,x), a(t) представлены дискретными значениями: q(tj ,х;.) = § .(; V(ij,xj) = Vji;
л Л л л л лАЛ Л ___ ______ J
K(t j ,xt) = Kjt\ S(tj ,xt) = Sji\ a(tj) = aj\ j =0,n, i = 0,m + l; Ax =---;
m+1
A/=—; {? ,x,}- узлы равномерной сетки на интервалах I, е |0.1| и
п
х е [0,1].
В основе вычислительной модели, строящейся в данной работе, лежит
подход, согласно которому искомое решение уравнения (1) q(t,x) = q(t,x) представляется функцией вида:
п
qn{t,x) = C0{t)-u0{x) + '^jCk (0'11 к (х) + си+1 (0-ми+1(х), (3)
к=\
где {ик(х)}„ = {uk(x)}n - совокупность базисных функций, определенная на
отрезке хе [х0,Х] = [0,1], х = х; {Ck(t)}, Л' =()./? +1 - коэффициенты представления функции q(t,x) в базисе {ик(х)}„ в момент времени t =
= t e [0,1]. Представление (3) соответствует концепции конечных элементов и задает аналитическое представление искомой функции q(t,x), которое может осуществляться различными способами. В качестве базисных функций {ик(х)}„ в выражении (3) обычно выбираются положительные многочлены, т.е. искомая функция q(t,x) будет представлена полиномом конечной степени q„(t,x) и, следовательно, будет непрерывной и дифференцируемой. При этом базисные функции обладают свойствами:
и0(х = 0) = 1, м0(х = 1) = 0,
им(х = 0) = 0, ип+1(х = 1) = 1, (4)
ик(х = 0) = ик(х = 1) = 0; к = \,п.
31
Функции и0(х) и и„ | (х) принято называть граничными базисными функциями. В силу свойств (4) будут справедливы выражения:
Со(0 = q(t,x = 0) = goCO; CM(t) = q(t,x = 1) = q„+\(ty, 0 < x < 1.
Тогда соотношение (3) можно переписать в виде
п
<7и(^х) = <7о(0 + £ С, (0 •«* (*) + ?„+!«• (5)
к=\
Для полной определенности значений функции q„(t,x) в (5) требуется найти неизвестные коэффициенты {( '/(0!. Решение этой задачи основано на применении в вычислительной модели уравнения (2), (3) вариационного подхода, согласно которому строится функция невязки, после чего выполняется ее минимизация, позволяющая определить неизвестные коэффициенты
{.cm-
Построение вычислительной модели включает в себя следующие этапы. Выполняется процедура приведения уравнения (1), (2) к системе линейных алгебраических уравнений на основе представления (3), в качестве которого выбирается многочлен Бернштейна. Для определения значений неизвестных коэффициентов {Ck(t)}, к = 0. п + \. входящих в данную систему, строится функция невязки по методу наименьших квадратов, которая затем сокращается с помощью метода минимизации функции многих переменных.
Ниже подробно излагаются этапы построения вычислительной модели, приводятся результаты работы ее программной реализации, выполняются численные исследования алгоритма модели.
2. Построение вычислительной модели. Многочлены Бернштейна, применяемые для аппроксимации функций, представленных дискретными значениями и соответствующих им производных, имеют вид [2]:
вп (*) = Z /(** )Л* W = X )си ■хк (! - хТк,
к=0 к=0
В'„(х) = «Х[/(х,+1)-/(х,)]С*_1х*(1-хГм, (6)
к= О
_В"(х) = п(п- 1)2 [Ж+2 )■- 2Ж+1) + /(X* )]С*2х* (1- хТк-2, X е [0,1].
к=0
При этом полагается, что fix) = {/(xt) j. к = 0,п, fix) ~ Вп(х), т.е. Вп(х) => f{pc), fix) ~ B'nix), fix) и ВЦх). В качестве базисной функции
со
выступает функция p„kix) = Скхк i\-х)п~к для которой справедливы свойства (4): рП: 0(х = 0) = 1, рП: 0(х = 1) = 0, p„t„ix = 0) = 0, рп^х = 1) = 1, p„fix = 0) = Pn,kiX = 1) = 0.
Выбор многочлена Бернштейна обусловлен тем, что он достаточно хорошо исследован [2] и обладает рядом замечательных свойств. Отметим некоторые из них: ИтВп(х) = /(х) равномерно по х для всеххе[0,1]; конструиро-
со
вание многочленов Бернштейна Вп(х) по некоторым специальным значениям Дх/.) в дискретных точках промежутка [0,1], в котором осуществляется приближение функции Дх) многочленами В„(х), выполняется достаточно просто; если у функции Дх), где хе[0,1], существует непрерывная производная /®(х), то Ит ///;) (х) = [<к> (х) равномерно относительно х; порядок приближения
п—> со
Е„(/) = |Дх) - Вп(х) | при п —> да многочленами Вп(х) не зависит от природы (дифференциальных свойств) Дх); если функция Д.х) выпукла (вогнута) на [0,1], то и многочлен Вп(х) также является выпуклой (вогнутой) в указанном интервале. Запишем далее представления (6) для аппроксимации искомого решения уравнения (1) и его производных:
Л Л ОТ + 1 Л л
= )А.+и(*) = вш+1 (Ь,чХ рт+и(*) = Си*С1 -*Г+1~',
(7)
ах /=о
д2 л. л. л Ш~1 л л л
— q^tJ,x) = ^m +1 )т£ [<?(0, */+2) -дх2 /=о
-2 , х,.+1) + д(^, х,. )]/>т_и (х) = 5"+1 (О, х, д).
Аппроксимация исходных данных и их производных выполняется так же, как и (7). Соответствующие аппроксимационные формы обозначим
Вт+1«/х,П Вт+1{1рх,к), /Г. и .х.Г). В'т+1(1рх,к).
Исходное уравнение (1) с начальными и граничными условиями (2) с использованием выражений (7), (8) может быть преобразовано следующим образом. Пусть
5 д(1 х) ^ Вт+1 (?у ,х,д)~ Вт+1 (р-1, х, д)
3? А?
где Д / —>■ 0. Подставляя (7) - (9) в (1), получим
С/ЛЛЛ т л л Л л л Л л л л л
—, х) = (да + 1)У [<ДО,х,.+1 ) - ,х,.)]/?(х) = В’т+1 (О,х,<?),
<*j (*)£ qj,Pm+u (*) + ь, (x)(m +1)£ [q. ,.+1 - q ]/?„, (х) -
,=о ,=о (Ш)
m-l т+1 v 7
-с. (х)(т +1 )да£ [?. ,.+2 - 2?. ,.+1 + qp \рт_и (х) - £ q}_up^ (х) = dj (х),
/=0 /=0
где обозначено:
a(tj,x) = \ + a(t j)At+ р ■ B'm+l(t j ,x,V) = a fix),
b(ij,x) = P- At Bm+l (tj ,x,V)-y-At B'm+l (tj ,x,K)= bj (x), c(tj ,x) = у-At Bm+l (tj ,x,K) = с. (x),
d(t j ,x) = % ■ At Bm+l (tj ,x,S) = dj (x), где x = x, qjt = q(tj,Xi), j =0,n. Вьфажение (10) представляет собой систему линейных алгебраических уравнений относительно неизвестных значений {qjfi, j = 1,п, i=\,m при
заданных начальных и граничных условиях: q0i, /=0,да+1; и ^>т+ь ; = 0 ,п.
Систему (10) можно далее преобразовать, если ввести обозначения:
gj(x) = qj,o{aj(x)pm^ho(x) - Ъ3(х)(т + 1 )рщ0(х) - с3(х)(т + 1)-да-^1>0(х)} +
+ qj,mt\{aj(x)pm+l:m+\(x) + Ь3(х)(т + 1 )рщт(х) - с3(х)-(т + ^-т-р^^х)} -
qj-l,0 Рт+\,0(Х) qj—\rm^\Pm+\,m+\(x),
е3(х) = а3(х)рт+и(х) - Ь3(х)(т + 1 )рщг{х) - с3{х)(т + 1 )-т-рт_и(х);
/3(х) = Ъ3(х)(т + 1) ртз(х) + 2 (да +1 )да-с,(х)^1>1(х); г/х) = -с3(х)-(т + 1 )-т-рт_^(х); г,(х) = Ь3(х)-(т + 1)рт,0(х) + 2 (да + 1 )-т-с3(х)рт_ио(х); м^х) = -с3(х)-(т + 1)-да >0(х);
и/х) = а3(х)рт^\ (х) - Ъ3{х)(т + 1)^т>тч(х) - с3{х)(т + 1)-да ,т-1(х);
.у/х) = я,(х)/>т+1 >т(х) + £/х)(да + 1)/?т,т-1(*) - Ъ3{х)(т + 1)рт>т(х) +
+ 2 (да + 1) -да-с/х) ■ргг&\>т_1 (х).
В этом случае система (10) получит представление:
т-2
X (х) + ?у,/+1/у (х) + 9и+■,?} (*)] + Я,Л (х) + (х) + qJm_luJ (х) +
/=1 А)
+Ч^У] (*) - (Х) - <1м,тРт+1,т (*) + 8у (*) = ^ (*)•
В (11) введем обозначения Ci(j) = qji и С/у-1) = qj_^i, 7=1,п, /'= 1,да. Тогда (11) превращается в систему относительно неизвестных коэффициентов С = {С, О)}:
A (j)fj (х) + С2 (j)Wj (х) + £ [С, (;>_, (х) + См (j)fi (х) + См (j)Zj (х)] +
/= 1
Ст-1 (У')му. (х) + Ст (;)^- (*) - ст_х (у - l)/?m+1>m+1 (х) -“О, (У - 1)Рт+1,т (*) + gj (*) = dj (х).
С = {Ct(j)} находятся из условия минимума функционала (метод наи-
1 _________________________________________
меньших квадратов): f p2(x)dx => min, где невязкар/х), j = 1,п определяется
о CU)
из соотношения
Pj (х) = сх (у Yj (х) + С2 (j)Wj (х) +
т-2
+Yj й (j)ej (х) + А+1 U)fi (х) + См (j)Zj (х)] + Ст_х (j)Uj (х) +
/= 1
+ст (ЛУ; (X) - Ст_х (У - 1)Рт+1т+1 (х) - Ст (У -1 )pm+l m (х) + g. (х) - dj (х). Определив значения коэффициентов ( = {()(,/)\. получим
q(t],x) = Bm+l(t],xAj)\ j=\n, Vxe [ОД]. (12)
Выражение (12) определяет искомое решение уравнения (1).
3. Результаты программной реализации и численные исследования алгоритма модели. Описанная выше модель решения уравнения переноса была детально алгоритмизирована, программно реализована с использованием тестовой задачи и численно исследована. Тестовая задача, разработанная и описанная в [1], условно называется «Блоком исходных данных». В «Блоке исходных данных» по определенному алгоритму формируются массивы
q(tj,xi) = qJr, V(tj ,xi) = Vji; k(ij,xi)=kji; S(tJ,xi) = SJi; a(tj) = aji;
j = 0,n, i = 0.in +1: Дх = —-—; A/= —; вычисляются значения нормиро-
m +1 n
ал л V*T K*T S*T л t
вочных коэффициентов: a(t) = T-a(t), /3 =----------, у-——, Е, -—t-—,
X X с[ Т
Л X Л Л
X =---, t Е [од] , х е [0,1]; выполняется процедура нормирования переменных
X
и функций уравнения переноса. Кроме исходных данных, представленных массивами, в блоке генерируются значения массива qT ={?уг}, j = 0,п,
/ = 0.111 + 1. представляющие собой точное решение. Для генерации значений
указанных массивов требуется задать исходные данные, которые для проведения вычислений в данном примере задавались следующими: Т = 5 с, У0 = 5 м/с, К0 =100 м2/с,Х= 50 м, д0 = 0,75 кг/м-с, а0 =0,75 1/с, п = 10, да + 1 = 20.
Аппроксимация функций исходных данных У(1, х), К(1. х). 8(1. х) и первых производных I "(1.x). к'([.х) по дискретным измерениям {( х )|.
!/\’(/;.х )|. ^(ґ^Хі)}, являясь одним из шагов алгоритма вычислительной модели уравнения переноса, проведена на сетке узлов {£,х,}, / = 0. п.
/ = 0,да +1, ;(є|0.1|. х,є[0,1], п = 10, да + 1 = 20. Точность аппроксимации при этом составила: ст(У ,Вт+1(У)) = 9,54£-03, а(К, Вт+1 (К)) = 1,03 Е - 02,
Ст(^,5)л+1(^)) = 3,71£-03, сг(К’,В’т+1(к)) = 1,04Е-02, где
°ГЛ+1(П) =
= 4,62£-02,
1
а(Е,Вт+1(Ю) = -
(и+1)(да+2)у=0,=0
Результаты работы вычислительного алгоритма по расчету концентрации примесей приведены в виде массивов. Точное решение (/Т(1,х). значения которого генерируются в «Блоке исходных данных», представлено массивом г/7 = \ц,] 1|. / = \.п. /= 0,да+1,3, решение ц = (1.x). получаемое в ходе выполнения вычислительного алгоритма модели - массивом ц = {д„}, ./ = 1 ,п,
і = 0, да +1,3 :
т
ч =
'0,3785 0,4504 0,5657 0,6626 0,6891 0,6310 0,5194 0,4447'
0,3270 0,3891 0,4887 0,5724 0,5953 0,5451 0,4487 0,3841
0,2990 0,3557 0,4468 0,5234 0,5443 0,4984 0,4103 0,3512
0,3014 0,3585 0;4504 0,5275 0,5486 0,5023 0,4135 0,3540
0,3335 0,3968 0,4984 0,5838 0,6072 0,5560 0,4576 0,3918
0,3876 0,4612 0,5792 0,6785 0,7056 0,6461 0,5319 0,4553
0,4504 0,5358 0,6730 0,7884 0,8199 0,7507 0,6180 0,5291
0,5064 0,6026 0,7569 0,8865 0,9220 0,8442 0,6949 0,5949
0,5421 0,6450 0,8102 0,9490 0,9869 0,9037 0,7439 0,6368
0,5486 0,6527 0,8199 0,9604 0,9988 0,9145 0,7528 0,6445у
"0,3785 0,4661 0,5416 0,6379 0,7083 0,6996 0,5862 0,4447^
0,3270 0,3970 0,4474 0,5196 0,5693 0,5801 0,5138 0,3841 0,2990 0,3964 0,4618 0,4941 0,5196 0,5183 0,4586 0,3512 0,3014 0,4123 0,4692 0,5063 0,5380 0,5268 0,4508 0,3540 0,3335 0,4507 0,5131 0,5318 0,5495 0,5493 0,4862 0,3918 4 ~ 0,3876 0,5097 0,5770 0,6159 0,6271 0,6106 0,5445 0,4553 0,4504 0,5779 0,6554 0,6910 0,6946 0,6813 0,6148 0,5291 0,5064 0,6164 0,6949 0,7463 0,7497 0,7267 0,6538 0,5949
0,5421 0,6544 0,7112 0,7643 0,7663 0,7275 0,6753 0,6368 ч0,5486 0,6521 0,7087 0,7600 0,7548 0,7172 0,6764 0,6445у Отклонение а{с[ ,ц) расчетного решения д(1.х) от точного (/г(1,х). вы-
^ пт
числяемое по формуле сг =------У У К//, - д I. составило 5.68/-.’ 02. Сходи-
п-т ]=1 ;=11 ’ •'
мость вычислительного алгоритма проверялась в ходе выполнения расчетов с различными значениями параметров пят, определяющими размерность сетки узлов . / = 0. п. / = 0. т +1. Расчеты показывают, что значительное увеличение размерности сетки узлов приводит к снижению точности получаемых решений за счет накопления вычислительных ошибок. Было установлено, что если задавать значения параметров п и т из интервалов 15 < т < 25, 7 < п < 12, то можно получать решения с приемлемой точностью.
4. Заключительные замечания. Предложенная в данной работе вычислительная модель и основные подходы к ее построению обладают особенностями, достоинствами и недостатками. Особенность данной вычислительной модели - неизвестному искомому решению ставится в соответствие некоторая аппроксимационная форма Ит1 (( (I).х). На ее основе выполняется редукция исходного уравнения переноса к системе линейных алгебраических уравнений, включающей в себя неизвестные коэффициенты ( (I), относительно которых она и решается. Основная идея - для каждого фиксированного момента времени искомому решению ставится в соответствие вектор
коэффициентовС(/) = \('1;(I)|. к = \,т, являющийся представлением решения
в соответствующем базисе {рпк(х)} в узловых точках х/.е |0.11. к = 0. т +1. Поиск значений неизвестных коэффициентов реализуется в модели процедурой минимизации соответствующей функции невязки. По коэффициентам ( ([) с помощью аппроксимационной формы , (( (I). х) получаем непре-
рывный аналог искомой функции </(1,х) \/хе[0,1] для каждого фиксированного момента времени 1е [ОД].
Применение многочленов Бернштейна для аппроксимации искомого решения в вычислительной модели позволяет реализовать их положительные свойства, изложенные выше. В частности, многочленами Бернштейна в вычислительном алгоритме аппроксимируются д'(/,х), (/"(1,х) искомого решения, а также исходные данные I '(1,х). I ’’(1,х). К(1,х). К'(1,х). ,4(1,х). Проведено численное исследование аппроксимации производных. Показано, что аппроксимация многочленами Бернштейна при достаточно большом числе узлов [ОД], к = 0, т +1 позволяет получить результат с точностью порядка 1,0Е- 03- 1,0Е- 02. Кроме того, вычислительная модель, использующая многочлены Бернштейна, позволяет «восстанавливать» значения не только искомого решения (1(1,х). но его первой производной (~/'(1. х). В качестве недостатков можно отметить следующее. Во-первых, необходимо аппроксимировать производные исходных данных I '(/,х). К'(1,х). которые, согласно постановке задачи, представлены дискретными измерениями и сами нуждаются в аппроксимации. Поэтому использование в вычислительной модели аппрок-симационных форм В'т+1(/,х,К(/,х)), В'^^^К^х)) в итоге снижает точность работы алгоритма за счет накопления вычислительных ошибок. Во-вторых, использование многочленов Бернштейна предполагает равномерное распределение узлов хк на всем интервале [0,1]. Эго в свою очередь означает, что измерения исходных данных также должны быть проведены в этих точках, что не всегда соответствует реальным прикладным задачам.
В целом вычислительная модель уравнения переноса «работает» вполне удовлетворительно и позволяет получить решение с приемлемой точностью порядка 1,0Е - 02.
Литература
1. Семенчин Е.А., Наац В.И., Наац И.Э. Математическое моделирование нестационарного переноса примеси в пограничном слое атмосферы: Монография. М., 2003.
2. Микеладзе Ш.Е. Численные методы математического анализа. М., 1953.
Северо-Кавказский государственный технический университет 1 апреля 2004 г.
