Применение итерационно-интерполяционного метода для решения задач математической физики

В работах [1 - 7] предложен и развит итерационно-интерполяционный метод (ИИМ) для решения различных уравнений математической физики. Так при переводе работы [2] на английский язык [3] авторами сделано дополнение, в котором для одного частного случая утверждается, что погрешность аппроксимации составляет O ( h ^{2 k} ), где h — шаг разностной сетки, а k = 1, 2 — номер итерации. В статье [4] даны алгоритмы ИИМ для решения уравнений параболического, гиперболического и эллиптического типов, допускающие точное решение, и дана оценка точности полученных численных решений, в публикации [5] получили развитие ИИМ для решения трехмерного уравнения теплопроводности. Обобщение ИИМ на трехмерный случай для уравнения параболического типа сделано в работе [6]. В книге [7] излагается современное состояние ИИМ.

Суть алгоритма ИИМ заключается в использовании метода последовательных приближений на каждом элементарном отрезке разностной сетки, соответствующей области определения исследуемой краевой задачи, в результате чего точность приближенного решения повышается как путем уменьшения шага разностной сетки, так и путем увеличения числа итераций [1]. В дальнейшем устанавливается связь ИИМ с теорией сплайнов [3] и даны примеры применения к решению некоторых нелинейных краевых задач [4].

2. Построение разностной схемы для уравнения типа пограничного слоя на основе ИИМ

Пусть требуется решить параболическое уравнение вида

5 L 5X'    ГX г v . , 5X _ 5X f+ f^+ f X + f-+ f-

5j( ¹ 5Z )   ² 5Z ^{3      4    5} 55    ⁶ 5/

в области 0 <  Z <  S , 0 < ф <Ф , 0 <  Z <  L с начальным условием

X I _e — 0 — X ₆ е ( Ф , Z ),    X I ,— 0 — X „е, J )                                              (2)

и, для простоты выкладок, с граничным условием первого рода

X I Z— 0 — X , ( 5 , ф ),    X Z— L - X e ( 5 , ф ).

Алгоритм ИИМ для уравнения параболического типа (1) приводит к следующей системе линейных уравнений:

AX i + 1 - B i X i + C i X i - 1 + D i — 0     ( i — 1, N - 1) .

Коэффициенты A i , B i , C i , D i согласно алгоритму, подробно описанному в [7], с учетом трехмерного характера течения и заменой частных производных конечными разностями д X /55 — ( X - X ₅ )/ А5 , д X /5ф — ( X - X _ф )/ Аф определяются следующим образом:

A i ^— Y 1 ^А5Аф( ^f ,, , + , ^{+ f} „■)+? 1 -Ч-Аф^у-¹ ( ^f 2, , „ + ² f ₂, ⁾⁺

+ у,А5Аф ^ ( f + f 3, )^ /^ ( f _w + f „) - , А ^A J • ( f _w + f J ; 6                      6                      6

B, — Y . А5Аф ( f z , „ + f _U) +Y . А5Аф^А^₃+ 1. ( f ₂,„ + 2 f ₂, , ) -

- _{Y l} А5Аф^А| +. ( f_U1] + 3 f _u) + _{Y ]} A

1.(f_:,i„ + 3f,,.)+ 66

+Y]А5^А|⁺■ (f,,i+1 + 3f,,)₊y,a^A; (f6,-1+3f,,i)+ 66

• ⁺ Y 2 ^А5АФ⁽f1,i--1 ^{+ f} J-Y 2 ^А5АФ^А|¹(f2,i-1 ^{+ 2 f} 2,i-)^-

• ^-Y2 ^А5АФ^А^(^fVi-!^{+ 3f} ..J+Y2^АФ^А|^L(^f;,M ^{+ 3f},.);

Ci^— Y 2 ^А5АФ^(f|,i-1 ^{+ f}u)-r 2 ^А5АФ^А|^{А (f}l.i-1^{+ 2 f}2i⁾⁺

+ Y2А5Аф^|2(f,„ + f J-Y А^/ 'J (f_U.1 + f_;,i)- Y А- 'J (f6,-1 666

;

AZ 2, /               x             AZ2 /               x

Di =-Y i A5Aф -|+L (f_{w +} 2 f„)-, 2 A5Aф -3- (f4.—1 + 2 f'J+

AZM

12 A5

5, i+1

-

+

+ Y \ф 7 [(f5,,-i + 3f,,,)X„ +(f,5,,—1 + f,>_e—1 ]+

■ A1^Az' [(fs.i« + 3 f.,,)Xф, +(f, .1 + f.,,)X.,+1 ]+

+ Y \_ ^A6 [(f_u_, + 3f,.)X_w +(fу,, + f;>.,—1 ].

В соотношениях (5) и ниже обозначено: X^ — значение искомой функции на предыдущем слое по 5; X^ — значение искомой функции на предыдущем слое по ф.

При добавлении к системе (5) граничных условий

X0 = Xw (5, ф), Xn = X. (5, ф)

получается система (N+1) алгебраических уравнений для определения неизвестных величин X₀, ... ,X_N . Можно показать, что приведенная разностная схема аппроксимирует исходную краевую задачу (1) с погрешностью O(д5+Аф + (AZ)²) в случае постоянного шага по переменной Z и с погрешностью O(AZ+Aф + AZ) — в случае переменного шага по Z , AZ = max(AZ,-). Здесь и ниже

Y1 = A£,/(AC,+A£_M), Y 2 =AZ„1 /(AZi+AZ„1), AZ, = Z" Z,-1 (i = 1, N-1). (7)

Так как при численном интегрировании системы уравнений пограничного слоя необходимо знать потоковые величины, связанные с градиентами искомых функций. Ниже приводятся выражения для вычисления указанных градиентов искомой функции:

=( X,- X- )•

Z=Z.

J1,,-1     Ji,i V2,i-1       J 2,i

2 AZ,          6

-

—7Y^[(fi,i-1^{+ 3 f}m )^Xi^{+ (f}i,i-1^{+ f}i,,)^X,-1^]+(f4,,-1^{+ 2 f}4,,⁾⁺ 12                                     6

^{+       [(f} 6,,-1^{+ 3 f} hi)(^X i^{- X}ф, )^{+ (f} hi-1^{+ f} v X^Xi-1^{- X}фi-1^)];

12 Aф

= ( X+1 - Xi)-

Z=Zi

^f,i+1 ^{+ f},i + ^f2,i+1 ^{+ 2 f}2,i

²AZi+1          6

-

^A;-' [(f3,_M+ 3f3,i)Xi +(f,+1 + fч)Xi,1 ]+ ^T"(f>,i,1 + 2f4,i)-12                                      6

7^AT_L [(f , 12 A^ ⁷5’i⁺¹

i +1 - ^XEi+1 ^{)] +}

⁺₁    ^{+1 [(f}6,/' + 1 ^{+ 3 f}, ii X^Xi - ^Xф/ ) ^{+ (f}6,i' + 1 ^{+ f},,i )(^Xi' + 1 - ^Xф/ + 1 )] .

12 Аф

В силу того, что при турбулентном режиме течения в пограничном слое коэффициент f₁ в уравнении сохранения количества движения во внутренней области турбулентного ядра имеет структуру

**

J1 = J1 + J 2 ЧТТ ’ 9z уравнения сохранения количества движения сводятся к виду

д

9Z

дXАдX + f

^{1    2}dZ JdZ

, дX  r v , г XX    XX

^{+ f}2 7"7 ^{+ f}3^X = ^f4 ^{+ f}5 777 ^{+ f}6 777 ,

9Z              9Z     дф

а выражения, содержащие градиенты искомой функции, записываются как:

9X АдX

= ^(X,-1 -

Z=Z.

X f²i-1+²f² > +(X,- X.-1 )x

X

***

I, i-1 ⁺J1, i⁺( J 2, i-1

2 AZ,

i i-1

AZ,

^{-      [(}fi,i-1 ^{+ 3}fi,i)^Xi^{+ (}fi,i-1 ⁺fi,i)^Xi-1^{] +} "77"⁽f4,i-1 ^{+ 2}f4,i) ⁺

12                                     6

⁺|2A^|■ ^[(f5,i-1 + ³f5,i )(^Xi - X^i^{)+ (}f5,i-1 + f5,i )(^Xi-1 - X^i-1 )

^{+        [(}f6,i-1^{+ 3}f6,i)(^Xi^{- X}фi) ^{+ (}f6,i-1 ⁺f6,i ^^X^Xi-1 - ^XФi-1 ^)];

12Аф

+

X

дX АдX

= (Xi - X +1i) x f^,,M+ ²f2i + (X„, - X, )x

Zz.                       ⁶

***

I,/+1 ⁺f 1,i⁺( f 2,/+1 ⁺ .

2А<1

i+1 ~ -AZi+1

-

Tr [(f3,/.1 + 3 f3,/)X/ +(f3,/+ + f3,/)X/+, ".'■-■ (f,,/„ + 2, f,,/)+

12                                      6

₊^А^1+¹[(f

12 А. ^i

+ 3 fJX_- X J+(f5,_;+1

+

-

+

+

^AZ/+1[(f 12Аф ^/6,z+1

^{+ 3}fмДXi^-Xф)^{+ (}f_v+1 ⁺fvДX+1

^{- X}фi+1^)].

С использованием условия сшивки потоков (11)-(12) в i-ом внутреннем узле далее записывается разностный аналог уравнения (1) в виде (4), где коэффициенты Ai,Bi,Ci,Di равняются:

Ai = Y1 ^А.^Аф^(fi*-+1 ⁺fi**⁾⁺ Y1^A.^AФ⁽f^+1 ⁺f J Xi+* Xi⁺

^AZi+1

+ Y1 А^ф^(f22+1 + 2f2.,)+Y1 '-Аф Д/ (f„+1 + f3,)-36

Y АфА; - (fw+fJ-m.^(f„+1+f„.); 66

Bi = Y1 А.Аф(f+ f'J+Y 2 А.Аф(f+ f')+

А^+1(.

• ^{+    Y}1 ^А.^Аф3 ^(f 2,i+1 ^{+ 2 f} 2,i ) ^Y2 ^А.^Аф3 ^(f 2,i-1 ^{+ 2 f} 2,i )

• - Y1 А^Аф^А|^,+1- (f,,_M+ 3 f,,,)-y 2 а;Аф^А)ф (f„-1 + 3 f„)+ 66

• + Y1Aф^A|²⁺¹(f_+ 3f_i.,)+ y Аф ^Л;^- (f,,_1 + 3f_i.,)+ 66

+ Y₂А^ ^i- (f ,_и + 3f_ы)+ . A. '; ' (f„+1 + 3,f„); 66

^Ci = Y2^А.^Аф^(f!*-1 ^{+ f}1*-⁾⁺ Y2^А.^Аф⁽Гтд-х^{+ f}2',i^{) Xi} X^iA ,            ,                          •            ,        ^А^.

-

А^,/                     А^2/

-Y 2 А.Аф —(f₂,_z_1 + 2 f2,J+y 2 А.Аф —(f₃,_z_1 + f.J-

- Y 2 Аф ^А|²(f_;,_1 + f_;,)- Y 2 А. ^ (f;,_и+ f„); 66

АС2.,                          АС2

^Di =^-Y1 ^А.^Аф3 ^(f4,i+1 ^{+ 2 f}4,i^)-Y2 ^А.^{Аф     (f}4,i-1 ^{+ 2 f}4,i^{) +}

-

^- [(f.2.1 + 3 f„JX- Xф,)+ 12Аф

+ Y 2Аф^А|²[(f.2-1 + 3 f^JX.2+(f 52-1

+

-

[(f„+1 + 3 f_UXx,-X_5i)+

+ Y А_ ^А6 Кf_u4 + 3 fJx_<+(f,._A + f,)x_w_i ].

Граничные условия при этом имеют вид (6). Порядок аппроксимации разностной схемы (13) совпадает с порядком аппроксимации схемы (5).

Схема (13) позволяет рассчитывать турбулентные режимы течения. Однако присутствие в коэффициентах искомой функции, которая берется с нижнего итерационного слоя, может приводить к существенному уменьшению скорости сходимости, а в отдельных случаях, и к расходимости итерационного процесса, так как в области турбулентного ядра коэффициент турбулентной вязкости много больше величины коэффициента молекулярной вязкости, то есть dX f1 << Л — .

dZ

С учетом неравенства (14) ниже приводится один из вариантов схемы, свободной от указанного недостатка.

Потоковые величины (11)–(12) переписываются в следующем виде:

^a. (X-x,- )²[ V^ ]²;

^Zi

4- (x_m - x, )² [ /+7; ] ²,

где величины, входящие в (15)–(16), определяются из соотношений:

£i

^a2 ^(X,^{- X}1-1 ^)+a3 ^{AZ a}1 ^(X,^-X,-1 ⁾²

_: P2 (X..1 -X, )+P3AZ_M Pi (X..1 - X, )²

«1 = ¹(f2*,-i+f2*,);                   «2 = ¹(fi'-i + fX^t f_U-1+2 f2,,);

a = — (f, 1 + 2 f /)- ^^ [(f Z 1+3 f z)Xz- + (f Z 1+f3 z)Xz-1 ]+ 3                       4,1 1            4,,                         3,, 1            3,Z                   3,Z 1         3,,,

+      ^[(/5,,-1 ^{+ 3f}5?k^X,^-X Z!'^)+(f5,-1^{+ f}5ДX-1 ^-X y-1^{)] +}

+ -^^[(f„,-1 12 Аф

P1 = ¹(f*,+1 + f2");                    P2 = 1 (/1 *+1+/1^x- ^X; ' (f_u+1+2f„);

2     26

Рз = -^7+1 (fм+1 + 2 fм)+ ^ [(fм+1 + 3 fм)Xi + (f«*1 + f’■ > X+1 ]6                  12

AZ+1

12 A^

[(f5,/+1

+ 3 f,XX,- X_w)*(f,.„1

^-X«+1^)]-

AZ*!

12 Аф

[(f6J.1

На основе условия сшивки потоков

0,5

а также разложения функций (1 + £1)0’5 и (1 + £₂)^0,5 в ряд по степеням £1 и £₂, соответственно, с удержанием первых двух членов ряда с точностью O(£²), где £ = max( £1, £₂), получается следующий вид коэффициентов разностной схемы:

A=Y1 ^1 ,     B, = Y 17^1 *Y 2 7°1 ’     C, = Y 2 Т°Г ’

Di =Y 1^AZ+1 00^^-Y 2 ^ТТ^ *

2 -7P1          2^ a 1

Y 1в з AZ 2+1

2 7p1 (X+1 - X)

Y2^aз^AZ2

2701 (X,- XM )■

Для использования схемы (19) необходимо выполнение условия

£ = max(£1, £2 )<£₀,

где £o — некоторая малая величина, выбираемая исходя из необходимого порядка аппроксимации по Z исходного уравнения в рассматриваемой области интегрирования. Данная схема является монотонной, так как справедливы соотношения А> 0, C, > 0, В, = A, + C,, что делает метод прогонки абсолютно устойчивым.

Тестирование схемы (19) показывает, что эта разностная схема позволяет рассчитывать систему уравнений пограничного слоя для чисел Рейнольдса > 10¹² , при этом скорость сходимости увеличивается в несколько раз.

Область интегрирования по нормальной координате разбивается на подобласти 0< Z- Z* и Z* L, где точка перехода Z* находится из условия равенства коэффициентов турбулентной вязкости pt, вычисленных по формулам для внутренней и внешней областей. Для внутренней области используется разностная схема, для которой коэффициенты A ,B , C , D определяются по формулам (13), а при выполнении условия (15) — по формулам (19). Для внешней области используется схема (5). Таким образом, при расчете турбулентного течения происходит стыковка трех разностных схем.

Отметим, что в силу физического смысла условий «cшивки» (18) разностные схемы, получаемые при помощи ИИМ, обладают свойством слабой консервативности [8], то есть консервативность удовлетворяется с порядком аппроксимации [7].

В работе [7] методом Фурье [9] в приближении замороженных коэффициентов получено, что разностные схемы ИИМ безусловно устойчивы.

3.    Сходимость и пример применения метода

Сходимость разностных схем ИИМ устанавливалась практически на последовательности сгущающихся сеток по т и hm (m = 1, 2 ,3) при решении частного случая трёхмерного уравнения (21) со смешанными условиями (22) в кубе Q [0 < x_m< 1, m = 1, 2, 3] при 0 < t< tk :

g

д X dt

- Z⁽Z - 1)vaxa-2 ] ^ exp⁽t);

X

X

= r, t=0

= (1+ x 1 =0

r=¹+E xa, a=1

x Z + xZ )exp(t) ,

X

= (1 + xZ + x₃^Z )exp(t) , x 2 =0

X

x з=0 д X

(1 + xZ + xZ) exp(t) ,

dx, ^x1=1

= Z exp(t) ,

X

X

= (2 + xZ + x 3) exp( t), x₂=1

= (2 + xZ + xZ ) exp(t). x 3 =1

Легко видеть, что задача (21), (22) имеет точное явное решение

X = (1+E xa) exp( t)• a=1

Использовались следующие значения входных данных: h_m = h, v_m = 1 (m = 1, 2, 3), Z = 6, g = 10, т = 0,002. При исследовании сходимости разностных уравнений ИИМ фиксировалась максимальная относительная погрешность

IX - JXI100%

• £ = max^J1-------,

• x, t^eQ        X

4.    Выводы

где X — явное точное решение (23), X~ — приближенное численное решение задачи (21), (22) при значении t_k = 5 . Из таблицы видно, что оно сходится к точному при измельчении шагов разностной сетки (t₀ — расчётное время в минутах на ПЭВМ).

Таблица. Максимальная относительная погрешность s при различных т и h

т	0,002	0,002	0,002	0,002	0,01	0,02
h	0,2	0,1	0,05	0,025	0,05	0,05
s, %	5,96	1,96	0,458	0,128	0,658	2,67
t0	0,2	0,4	1,15	8,1	0,35	0,12

• 1.    На основе ИИМ получена разностная схема для решения уравнения параболического типа с погрешностью аппроксимации на равномерной сетке O[Δξ + Δϕ +(Δζ)²].

• 2.    Показано, что при постоянных (или гладких) коэффициентах переноса разностные уравнения безусловно устойчивы [7].

• 3.    Проведен тестовый расчет, и анализируется практическая сходимость разностной схемы.

• 4.    Построенные разностные схемы использованы для решения уравнения типа пограничного слоя [10] с учетом различных режимов течения и в широком диапазоне чисел Рейнольдса и Маха.

Отметим, что положительные свойства разностных схем, получаемых с помощью ИИМ (сплайн-аппроксимация искомых решений и др., возможность получения схем повышенной точности) для одномерных краевых задач, сохраняются и при решении многомерных задач [7]. В частности, если в качестве начального приближения в логической схеме из [7] вместо линейной функции взять параболу, то для уравнения (21) при h_m = h, g = v_m =1, по алгоритму ИИМ [7] можно получить дифференциально-3

разностную задачу с погрешностью аппроксимации O( ∑ h_m⁴ ).

m =1

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект № 08-01-00496).