Обобщение итерационно-интерполяционного метода для решения трехмерного параболического уравнения общего вида Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 4, № 2, 1999

ОБОБЩЕНИЕ ИТЕРАЦИОННО-ИНТЕРПОЛЯЦИОННОГО МЕТОДА ДЛЯ РЕШЕНИЯ ТРЕХМЕРНОГО ПАРАБОЛИЧЕСКОГО УРАВНЕНИЯ ОБЩЕГО ВИДА *

А. М. Гришин, А. С. Якимов Томский государственный университет Томский университет систем управления и радиоэлектроники, Россия

Making use of iterational interpolation method absolutely stable difference schemes have been obtained for the solution of non-linear boundary problems of heat conduction. The approximation error has been found. The algorithm of the method is presented and the test calculation example is given for a three-dimensional parabolic equation in a cube under boundary conditions of the first and second kind.

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

В [4] на основе ИИМ и поочередного интегрирования по пространственным переменным получены однородные устойчивые явно-неявные разностные схемы с погрешностью аппроксимации в обычном смысле на равномерных сетках по пространству второго порядка точности. В последнем случае разностная схема для решения трехмерного уравнения теплопроводности условно устойчивая с ограничением на шаг по времени т < hj/(6x) (h3 — шаг по неявному направлению x3, X = Рср), X — коэффициент температуропроводности).

Часто при решении стационарных задач методом установления и особенно трехмерных уравнений математической физики необходимо, чтобы разностная схема была безусловно устойчивой, а также безытерационной, если коэффициенты переноса — постоянные величины. Цель данной работы показать, что последним требованиям отвечают разностные

* Работа выполнена в рамках Федеральной целевой программы "Государственная поддержка интеграции высшего образования и фундаментальной науки на 1997-2000 г.", а также при финансовой поддержке Российского фонда фундаментальных исследований, гранты №96-01-00964 и 96-01-00011.

© А. М. Гришин, А. С. Якимов, 1999.

схемы, полученные на основе алгоритма ИИМ [2, 4]. Отметим, что при этом неявно использовался подход Дугласа — Рэкфорда для трехмерного уравнения теплопроводности, приведенный в монографиях [6, 7].

1. Алгоритм метода

C целью простоты дальнейших выкладок логическую схему метода приведем для уравнения теплопроводности с постоянными теплофизическими коэффициентами внутри параллелепипеда R = {x = (xl,x2,x3), 0 < xm < dm, m = 1, 2, 3} при 0 <t < tk:

dT 3 02T

-Ш = + f(x-t) ■ (i-i)

m= 1 m

Для построения разностного аналога уравнения (1.1) введем сетку с координатами узлов xM = ih1 (i = 0,1, 2,...,N), x2J = jh2 (j = 0,1, 2,...,M), x3k = kh3 (k = 0,1, 2,...,L), tn = пт (n = 0,1, 2,...).

При написании операторов разностных схем вводятся п+1/3, п+2/3 — промежуточные (формальные) слои и n, п + 1 — начальный и конечный (фактические) слои по времени t, при этом шаг на каждом слое предполагается целым и равным т.

Используя результаты работы [4], введем следующие обозначения (точка сверху означает частную производную по времени):

• д2Т д2Т

a(t,x) = Т - С2dxj - С3Щ, (1-2)

д2Т

b(t,x) = Т — сзщ , (1-3)

hi,i-i = xM — xi)¿-i, h2,j-i = x2j — x2,j-i, hi,i = 0, 5(hi)j-i + hi,i),

h1,i = xi,i+i — xi,i , h33,k-l = x3,k — x33,k-l .

Разностные операторы внутри области определения R взяты из [4] и имеют вид

Aiai = [(hia)i-i + 4(hia)i + hi,iai+i]/6fti,i,

= Cl[(Ti+l — Т)/hi,i + (Т— — Ti)/hM-i]/frM,

i = 1, . . . , N — 1 , 0 < xm < dm , m = 2, 3,

A2bj = [(h2bi)j-i + 4(^2bi) j + h2,jbi,j+l]/6^2;j , Л2Ti,j = c2 [(Ti,j+l — Ti,j )/h2,j + (Ti,j-l — Ti,j )/h2,j-l]/^2,j ,

i = 1,...,N — 1, j = 1,...,M — 1, 0 < x3 <d3 ,

A3T ,j,k = [h3,k-lTi,j,k-l + 4^3,k'Ti,j,k + h3,kTi,j,k+l]/6^3;k , Л3Ti,j,k = C3[(Ti,j,k+l — Ti,j,k)/h3,k + (Ti,j,k-l — Ti,j,k)/h3,k-l]/^3,k ,

i =1,...,N — 1, j = 1,...,M — 1, k =1,...,L — 1. (1.4)

Разностные операторы внутри области определения Я для нового, предлагаемого ниже, подхода имеют вид

АТ = [ьм_1Т

к]/6ЯМ ,

Л1Т = С1 [(Тг+1,.,к — Тг,.,к)/^1,г + (Тг—1,.,к — Тг.,к)/^1,г_1]/^1,г ,

А2Т = [^

Л2Т = с2 [(Тг,.+1,к — Тг,.,к)/^2,. + (Тг,. — 1,к — Тг.,к)/^2,. —1]/^2,. , А3Т = А3Тг,.,к , Л3Т = ЛзТг).,к ,

г = 1,... - 1, ] = 1,... ,М - 1, к =!,...,£ - 1.

(1.5)

Следуя алгоритму из [4] и используя обозначения (1.2), (1.4) для уравнения (1.1) последовательно получим:

по направлению координаты х1

д2Т + /

а =С1 дХ2 + /

разностные схемы поочередного интегрирования по пространственным переменным на основе ИИМ имеют вид

^га+1/3

Ла = Л1Т^^ + А /Г, 1,..., N - 1, 0 < хт , т = 2, 3;

(1.6)

по направлению оси 0Х2

С2

д 2Т 5x2

+ «г

1,...,#- 1

разностные уравнения записываются как

АЬ,. = Л2ТГ+2/3 + ^2«,

г = 1,... - 1, ] = 1,... ,М - 1, 0 < х3 < 4

(1.7)

Окончательно по направлению координаты х3

Т

52Тг

г,.

С3-

г.

5x3

+ , г = 1,...,Ж - 1, ] = 1,

,М- 1

разностные схемы ИИМ принимают вид

,.,к = Л3Тг,+к + А3Ьг,.,к ,

г = 1,...,Ж - 1, ^ = 1,...,М - 1, к = 1,...,Ь - 1. (1.8)

Отметим, что разностные уравнения вида, например (1.6), получаются непосредственно из логической схемы ИИМ [2]. Для получения явной зависимости разностного соотношения (1.6) от неизвестного Т подставим в него значение а из (1.2); тогда получим

А (Т - - С3= Л1Т"+1/3 + /

i = 1,..., N - 1, 0 <Xm <dm , m = 2, 3 . (1.9)

Из выражений (1.5), (1.9) видно, что матрица 1 дискретного оператора AiT действует в направлении оси 1. Производные д2Т^/дж^, д2Ti/dx2 в (1.9) берутся от неизвестного T по координатам x2, x3, которые изменяются непрерывно. Поэтому эти производные можно вынести за знак матрицы 1 и аппроксимировать независимым образом через операторы À2Tij , A3Ti,j,k из формул (1.4). В результате, воспользовавшись обозначениями (1.5), имеем вместо (1.9) по координатному направлению x1 дифференциально-разностную схему (A1T = AiTijk, Л1T = A1Ti)j)fc, так как x2, x3 в интервале 0 < xm < dm, m = 2, 3 изменяются непрерывным образом в (1.4)):

A1T = Л1 Tn+1/3 + Л2Тn + Л3Тn + 1 fn ,

i = 1,...,N - 1, j = 1,...,M - 1 , k = 1,...,L - 1. (1.10)

Далее найдем i из уравнений (1.6), умножая обе его части на A-1 слева. Это можно сделать, так как матрица, соответствующая оператору, трехдиагональна: 1/6 4/6 1/6 при hm = const, m = 1, 2, 3 и с диагональным преобладанием. Тогда, замечая, что A-11 = E, получим

ai = Ar%Tn+1/3 + fn . (1.11)

Подставим в левую часть операторного уравнения (1.7) значение b из (1.3), а в правую — значение i из (1.11). В результате (1.7) перепишется в виде

2 (Tij - С3 ij = Л2ТГ3 + A2(A-1 Л1Т71/3 + fn ) ,

i = 1,... ,N - 1, j = 1,... ,M - 1, 0 < X3 <4 . (1.12)

Используя обозначения (1.5) и соотношение 2A-1 = E, аналогично предыдущему (1.9), (1.10) получим из (1.12) (A2T = A2Tij;fc, Л1Т = Л^,^, Л2Т = Л2Ti)j, так как x3 из промежутка 0 < x3 < d3 изменяется непрерывно в (1.4)):

A2T = Л1Тn+1/3 + Л2Тп+2/3 + ЛЗТп + A2fn ,

i = 1,...,N - 1, j = 1,...,M - 1, k = 1,...,L - 1. (1.13)

Наконец, подставим в правую часть уравнения (1.8) значение bi;j, полученное из разностного соотношения (1.7) и предварительно умноженное слева на A-1. Тогда разностная схема (1.8) перепишется в виде

A3T, = ЛзТ,1 + A3(A-^^,2^ + A-^1Tn+fc1/3 + fj). (1.14)

Окончательно, замечая, что A3A-1 = E, A3A-1 = E, из операторного уравнения (1.14) получим при использовании выражений A3T = A3TiJ)k, Л3Т = Л^. д. из (1.5):

A3T = Л^^3 + Л2Tn+2/3 + Л3Tn+1 + 3fn ,

i = 1,...,N - 1, j = 1,...,M - 1, k = 1,...,L - 1. (1.15)

Таким образом, для решения трехмерного параболического уравнения (1.1) внутри области определения R найдены разностные уравнения (1.10), (1.13), (1.15), аппроксимация и устойчивость которых будет рассмотрена ниже в разделах 3 и 4.

2. Постановка задачи и метод построения разностного решения

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

дТ А д ( дТ\ А д 2Т А дТ

д Т

ж = ^

m=i

dxr

dxr

+

Ur

m=1 k = m

r

dxmdxk л dxm

m=1

' dxr

+ f

(2.1)

в параллелепипеде R = {x = (xi,x2,x3), 0 <xm < dm, m = 1, 2, 3} при 0 < t < tk, w = const с начальным условием

T U = po(x) (2.2)

и с граничными значениями вида

T |xi=di T |X2=d2

T|x3=d3 = P3(t,d3,xi,x2,), (2.3)

Pl(t,di,x2 ,x3) , P2(t,d2,xi ,x3) , P3(t,d3,xi,x2,)

дТ

vi( dxi)

дТ

v2( дд~)

dx2 —)

dx3

x1=0

X2=0

хз=0

-qi(t,x2 ,x3) -q2(t,xi ,x3) -q3(t,xi, x2)

(2.4)

Уравнения типа (2.1) применяются в механике реагирующих сред [8] и в математической теории лесных пожаров [9]. В общем случае при наличии смешанных производных имеем 27-точечный шаблон в пространстве, а при отсутствии — 7-точечный крест в параллелепипеде К.

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

Разностные операторы внутри области определения К записываются в виде

ЛТ = т-1[Нц-1(дТ)г-и,к + 4Нц(дТ+ кц(дТ)г+ц,к]/6Ям ,

АТ = т-1[Н2,]-\{дТ)г,]-1,к + (дТ)^к + Н2,](дТ)г,]+1,к]/6Н2,] ,

(дТ) (дТ)

Л{Т = [(VI ,г + (ъ1,г-1,],к + у1,г,],к) х

х (Т

Л2Т = [(Ь2 к )(Тг )/Л2] + (У2 ,г,]-1,к + Ъ2,г,],к) х

х (Тг,]-1,к — Тг,],к,

ЛзТ = [(ь3 к )(Тг

х (Тг,],к-1 — Тг,],к)/Лз,к-1]/2^3,к , Л12Т = 2т[(Т+1 ,]+1,к — Тг+1],к — Тг-1,]+1,к + Тг- 1,],к )/Л2] +

+ (Тг+1 1,],к + Тг-1]-1,к) /Л2]—1}/4^1,г ,

g

v

предполагается, что краевая задача 2.1-2.4 корректна, т.е. ее решение существует, единственно и непрерывно зависит от данных задачи.

Л1,зТ = 2Ц(Тт + Т

+ (Тг+1 — Тг+1,.7>—1 — Тг—+ Тг—1 1)1] /4^-1,г ,

Л2,зТ = 2Ц(Тг

+ Т

+ (Т

Л1-3Т

А1Т А2Т А3Т

X

X

(Л1,2 +Л2,3 + Л1,3)Т, [(и1 №

[(и2,г,^+1,к + 2и2,г,^,к)(Тг,.+1,к — ) + (и2,г,. —+ 2и2,.) х

№

г = 1,..., N - 1, 3 = 1,...,М - 1, к = - 1.

(2.5)

Разностные операторы на плоскости х1 = 0 имеют вид

А1ТЖ1 =0

А2Тжх =0

А3ТХ1 =0

Л1ТХ1 =0

Л2ТЖ1 =0

Л3ТХ1 =0

Л1,2ТХ1 =0

+Л1,3ТЖ1 =0 Л2,3ТЖ1 =0

А1ТХ1 =0

А2ТЖ1 =0

А3ТХ1 =0

= т-1[2(дТ+ (¿О^/З ,

= («1

= [(«2

= [(«3

= 2^[(Т1;^+1;к — Г. — Т0)^+1)к + +

= 2^[(Т1.^+1 — — Т0,.,к+1 + +

+ — — + Т0,.,к-1)/^3,к-1 ]/2^1,0 ,

= 2^[(Т0.+1;&+1 — Т0,.+1,к — Т0,.-1,к+1 + Т0,.-1,к)/

= (и1,+ 2и1,0,.,к)(Т1.,к — )/3^1,0 ,

= [(и2,0,.+1,к + 2и2,0,.,к )(Т0,.+1,к —

+ 2и2,0,.,к)(Т0,.,к — Т0.—)]/6^2. ,

= [(и3,0,.,к+1 + 2и3,0,.,к )(Т0,.,к+1 —

,0,.'—

+

— 1 +

М- 1, к = 1,

- 1.

(2.6)

1

3

На оси ОХ3 (0 < х3 < ^3) пересечения плоскостей х1 = 0 и х2 = 0 при использовании

обозначений Тх3 =0 = Т(х = 0, х2 = 0) операторы записываются в виде

= т-1[2($Т)п,п,к + (^т)1,п,к]/3 ,

= т-1[2($Т)п,п,к + (^Т)п,1,к]/3 , -1

хз =0

А1Тж3=0 А2Тжз=П А3Тж3=П

Л1Тх3=П Л2Тж3=П Л3Тжз=П

= Т (^Т)п,п,к-1 + 4Яз,к (¿Т )п,п,к + (^Т)п,п,к+1]/6Яз,к

= («1

1,0

= («2,п,1,к + «2,0,0,к)(Т0,1,к — ТП,П,к)/Л1,0 ,

= [(«3

+ «3,0,0,к)(Т0,0,к-1 — Т0,0,к)/Л3,к-1]/2^3,к ,

Л1,2Тж3=0 = 2М(Т1,2,к — Т1,1,к — Т0,2,к + Т0,1,к )/Л2,1 +

Л1,3Тж3=0 Л2,3ТЖ3=0 ^1Тх3=0

Д2ТЖ3=0

^3Тж3=0

+ (Т1,1,к — Т1,0,к — Т0,1,к + Т0,0,к)/Л2,0]/2Л1,0 , = 2Ч(Т1,0,к+1 — Т1,0,к — Т0,0,к+1 + Т0,0,к )/Л3,к +

+ (Т1,0,к — Т1,0,к-1 — Т0,0,к + Т0,0,к-1)/Л3,к-1 ]/2Л1,0 ,

Т0

+ (Т0,1,к — Т0,1,к-1 — Т0,0,к + Т0,0,к-1)/Л3,к-1 ]/2Л2,0 , = (и1,1,0,к + 2и1,0,0,к)(Т1,0,к — Т0,0,к)/3Л1,0 , = (и2,0,1,к + 2и2,0,0,к)(Т0,1,к — Т0,0,к)/3Л2,0 , = [(и3,0,0,к+1 + 2и3,0,0,к )(Т0,0,к+1 — Т0,0,к) + (и3,0,0,к-1 +

+ 2м3,п,п,к)(Тп,п,к - Тп,п,к-1)]/6^3,к , к = 1,

- 1

(2.7)

На оси 0Х1 (0 < х1 < пересечения плоскостей х3 = 0 и х2 = 0 подобно операторным обозначениям (2.7) получим

А1ТЖ1=0 = Т [Л,М-1(дТ )г—1,0,0 + 4ЯМ(#Т)гД0 + Л-М (^Т)г+1,п,п]/6^1,

А2ТЖ1=0

А3ТХ1=0

Л1ТХ1=0

= т [2(#Т)г,0,0 + (^т)г,1,0]/3 , = т - 1[2(#Т)г,0,0 + (^т)г,0,1]/3 ,

= [(«1 ,г+1,0,0 + «1,1,0,0 )(Тг+1,0,0 — Тг,0,0 )/Л1,г + («1,г-1,0,0 + -1,0,0 — ТгД0)/Л1,г-

Л2ТХ1=0

Л3ТЖ1=0

Л1,2ТХ1=0

Л1,3ТХ1=0 Л2,3ТЖ1=0

Д1Тж1=0

^2Тж1=0 ^3Тж1=0

,1,0 ТгД0

2,0

= («3, г,0,1 + «3,гД0 ,0,1 — Тг,0,0)/Л2,0 ) = 2Ц(Тт ,2,0 — Тг+1 ,1,0 Тг- 1,2,0 + Тг-1,1,0)/^2,1 +

+ (Тг+1,1,0 — Тг+1,0,0 — Тг-1,1,0 + Тг-1,0,0)/Л2,0]/4Я1,г ,

= 2Ц(Тт

,0,2 — Тг+1,0 ,1 Тг- 1,0,2 + Тг-1Д1)/Л3,1 + + (Тг+1,0,1 — Тг+1,0,0 — Тг-1,0,1 + Тг-1,0,0)/Л3,0]/4Я1,г ,

= 2Ц(Тц,2 - Тг,1,1 - Тг,0,2 + Тг,0,1)/^3,1 +

+ (Тг

,1,1 Тг,1,0 Тг, 0,1 + ТгД0)/Л3,0]/2Л2,0 , = [(и1,г+1,0,0 + 2и1,г,0,0) (Тг+1,0,0 - Тг,0,0) + (и1,г-1,0,0 + + 2и1,г,0,0) (Тг,0,0 - Тг-1,0,0)]/6^1,г ) = (и2,г,1,0 + 2и2,гД0)(Тг

,1,0 Тг,0,0

= («3,гД1 + 2М3,г,0,0)(Тг,0,1 - Тг,0,0)/3Л-3,0 , г = 1, . . . , N - 1

(2.8)

2

В вершине параллелепипеда Я, совпадающей с началом координат: хт = 0, т = 1, 2, 3,

имеем

А1Тжт=0 А2Тжт=0 А3Тжт=0 Л1Тжт=0 Л2Тжт=0 Л3Тжт=0 Л1,2Тжт=0

т_1[2(^Т)0,0,0 + (^Т) 1,0,0]/3 , т_1[2(^Т>0,0,0 + (£Т)0,1,0]/3 ,

т

_1г

[2(^Т)0,0,0 + (^Т)0,0,1]/3 ,

(«1 ,1,0,0 + «1,0,0,0)(Т1,0,0 - Т0,0,0)/^10 («2,0,1,0 + «2,0,0,0)(Т0,1,0 - Т0,0,0)/^2

2,0

= («3,0,0,1 + «3,0,0,0)(Т0,0,1 - Т0,0,0

3,0

2^[(Т1,2,0 - Т110 - Т0,2,0 + Т0,1,0)/^2,1 +

+ (Т1,1,0 - Т1,0,0 - Т0,1,0 + Т0,0,0)/^2,0]/2^1,0

Л1,3ТЯт=0 = 2Ц(Т1,0,2 - Т1,0,1 - Т0,0,2 + Т0,0,1)/Л-3,1 +

+ (Т1,0,1 - Т1,0,0 - Т0,0,1 + Т0,0,0)/^3,0]/2^1,0

Л2,3ТЖт=0 = 2^[(Т0,1,2 - Т0,1,1 - Т0,0,2 + Т),0,1)/Л.3,1 +

+ (Т0,1,1 - Т0,1,0 - Т0,0,1 + Т0,0,0)/^3,0]/2^2,0

А1Тжт=0 = (и1,1,0,0 + 2и1,0,0,0)(Т1,0,0 - Т0,0,0)/3^1,0 ,

^2Тжт=0 = (и2,0,1,0 + 2и2,0,0,0)(Т0,1,0 - Т0,0,0)/3^2,0 ,

Д3Тжт=0 = (и3,0,0,1 + 2и3,0,0,0)(Т0,0,1 - Т0,0,0)/3^3,0 .

Наконец, на оси ОХ2 (0 < х2 < ^2) пересечения плоскостей х3 (2.7), (2.8) получим

(2.9)

0 и х1 = 0 аналогично

А1ТЖ2=0

А2Тж2=0 А3Тж2=0

Л1ТЖ2=0

Л2Тж2=0

Л3Тж2=0 Л1,2Тж2=0

-_1[

т ~[2(#Т)0,.,0 + (^Т)1,.-,0]/3 ,

= т_1[^2,._1 (^Т)0,._1,0 + 4^2,. (¿Т)0,.,0 + ^2,. (^Т)0,.+1,0]/6^2,.

= т_1[2(^Т)0,.,0 + (^Т)0,.,1]/3 , = («1

= [(«2,0,.+1,0 + «2,0,.,0)(Т0,.+1,0 - Т0,.,0)/^2,. + («2,0,._1,0 + + «2,0,.,0)(Т0,._1,0 - Т0,.,0)/^2,._ 1 = («3,0,7,1 + «3,0,.,0)(Т0,.,1 - Т0,.,0)/^2,0 , = 2^[(Т1,.+1,0 - Т1,.,0 - Т0 ,.+1,0 + Т0,.,0)/^2,. + + (Т1,.,0 - Т1,._1,0 - Т0,.,0 + Т0,._1,0)/^2,._1 ]/2^1,0 ,

Л1,3Тх2=0 = 2^[(Т1,.,2 - Т. - Тз,.,2 + Т0,.,1)/^3,1 +

+ (Т1,.,1 - Т1,.,0 - Т0,.,1 + Т0,.,0)/^3,0]/2^1,0 , Л2,3ТЖ2=0 = 2М(Т0,.+1,2 - Т0,.+1,1 - Т0,._1,2 + Т0,._1,1 )/^3,1 +

^1Тж2=0

^2ТЖ2=0

^3Тж2=0

+ Т0,.+1,1 - Т0,.+1,0 - Т0,._1,1 + Т0,._1,0)/^3,0]/4^2,. ,

= (и1,1,.,0 + 2и1,0,.,0)(Т1,.,0 - Т0,.,0)/3^1,0 ,

= [(и2,0,.+1,0 + 2и2,0,.,0)(Т0,.+1,0 - Т0,.,0) + (и2,0,._1,0 +

+ 2и2,0,.,0)(Т0,.,0 - Т0,._1,0)]/6^2,. ,

= («3,0,.,1 + 2М3,0,.,0)(Т0,.,1 - Т0,.,0)/3^3,0 , ] = 1,... М - 1.

(2.10)

Не представляет трудности выписать по аналогии с (2.6) разностные операторы для плоскостей х3 = 0 и х2 = 0.

2

В случае переменных коэффициентов vm = vm(x,t), um = um(x,t), и hm = const, m = 1, 2, 3, используя обозначения (2.5), получим следующую систему разностных уравнений внутри параллелепипеда R:

Ai(Tп+1/3 - Tn) = (А! + Д1)Тга+1/3 + (Л2 + A2)Tn +

+ (Лз + Дз)Тп + (Л1,2 + Л2,3 + Л1,з)Тп + TAifn , (2.11)

А2(тп+2/3 - Тп) = (Л1 + Д1)Тп+1/3 + (Л2 + Д2)Тп+2/3 +

+ (Л3 + Д3)Тп + Л1-3Тп + тА2fn , (2.12)

A3 (Тп+1 - Тп) = (Л1 + Д1)Тп+1/3 + (Л2 + Д2)Тп+2/3 +

+ (Л3 + Д3)Тп+1 +Л1-3Тп + ТА3Р ,

i = 1,...,N - 1, j = 1,...,M - 1, k =1,...,L - 1. (2.13)

Практически при исследовании устойчивости разностных схем (2.11) - (2.13) и для написания программы расчета удобно пользоваться вместо двух последних уравнений (2.12), (2.13) более простыми разностными уравнениями [6, 7]. Вычтем из второго уравнения первое, а из третьего второе, тогда вместо двух последних уравнений системы (2.11) - (2.13) получим

Л2Тп+2/з - Л{Тп+1/з = (Л2 + А2)(Тп+2/3 - Тп) + (Л2 - А)(Т + т!)п , (2.14) АзТп+1 - л2Тп+2/3 = (Лз + Дз) (Тп+1 - Тп) + (Лз - Л2) (Т + т!)п ,

I = 1)...)Ы - 1, 3 = 1,...,М - 1, к = 1,...,1 - 1. (2.15)

На плоскости х1 = 0 разностные уравнения при использовании обозначений (2.6) имеют вид (все разностные уравнения на границах получаются из логической схемы ИИМ [2, 4] и алгоритма из раздела 1)

А!^3 - Тп1=о) = (Л1 + ДОТ^3 + (Л2 + Д2)Тжп1=0+ + (Лз + Дз)Т,п о + (Л1,2 + Л2,з + Л1,з)Тжп1=о + тЛ1!1=о + 2д1М,к/к^ ,

x 1=0 + (Л1,2 + Л2,3 + Л1,3) Т xi=0 + Т A1Jxi=0 + 241,0,j,k/ h1,0

А2Тп+03 - AX^T = (Л2 + Д2)(Тп+03 - Тп=0) + (A2 - A1) X (Т™1=0 + ТРХ1=0), Aзтnlt10 - A2Tnx=* = (Л3 + ДзХТП+о - Тжп1=0) + (A3 - A2) х (^=0 + tJU) ,

j = 1,...,M - 1, k =1,...,L - 1. (2.16)

На оси OX1 (0 < x1 < d1) при помощи операторов (2.8) получим

A^f - ТХ1=0) = (Л1 + Д1)2Х503 + (Л2 + Д2)Тх™=0+

+ (Л3 + Д3)ТХ1=0 + (Л1,2 + Л2,3 + Л1,3)Тхп1=0 + TA1fX1=0 + 2^,0,0/h2fl + q3,i,0,0^3,0), A2T:+X3 - AX+f = (Л2 + Д2)(Тхп1+=20/3 - Тхп1=0) + (A2 - A1) х (Тхп1=0 + Tfx™=0) , A3Txn1+=0 - A2TX1+=20/3 = (Л3 + Д3)(Тхп+10 - TX1=0) + (A3 - A2) х (Т^ + тД^) ,

1,..., N - 1. (2.17)

В начале координат = 0, т = 1, 2, 3 при использовании операторов (2.9) имеем

л^3 - т>т=о) = (Л1+дот^о3 + (л2+Д2)т?т=0 + (л+Дз)т?т=0+

+ (Л1,2 + Л2,з + Л1,з)ТПт=о + ТА1 /Пт + 2 5т,о,о,о/Лто,о ,

т=1

А2ТП+:?/О3 - ЛТ:^3 = (Л2 + Д2)(ТГт+==/03 - Т") + (А - ^(Т™ + / )

АзТП+Ь, - ^Т^03 = (Л3 + Дз)СС=> - тхт=о) + (А - Л2)СС=о + /). (2.18)

Не представляет труда записать разностные уравнения для осей ОХ2, ОХ3 и для плоскостей х2 = 0, х3 = 0, используя операторные обозначения (2.6) - (2.10) и разностные схемы (2.16), (2.17).

На граничных плоскостях = т = 1, 2, 3 имеем

= РТ+М^ , 3 ^...М к = 0,._Ь

ТЙ = рПДЦ^, * = 0,...,^, к = 0,...,ь,

Т+1 = 1хз=^з, * = 0,...,^,3 = 0,...,М. (2.19)

В заключение приведем типичную схему ИИМ (2.20) - (2.23) на плоскости х1 = (0<хт< ^т, т = 2, 3) при заданном тепловом потоке и ит = т = 0, так как она используется при решении тестовой задачи (6.1), (6.2):

«/3 = Р^и^ , т =1, 2, 3 , к = 0,..., Ь , т;йт/3 = Р+и^ , т =1, 2, 3 ,3 = 0,..., М, (2.20)

№ —1 + /т ] =

ж N—1 - —+/т] + —1+

+ {3[(^2 + + —+ )х

- + 1 + 1 - )/^3,к—1 ]/2^3,к +

+2[(дТ+ ] + (дТ— +/т + /м — }п ,

3 = 1,..., М - 1, к = 1,..., Ь - 1; (2.21)

Х[3(^2,М,+ —+ )/^2,+—1 + + + 4^2+ /т] +

=

+ ^2,3 /т + )]П + 2^2,3 - )/т - +

— + (T.n-1l^73)fc - ^— + )/т - /м—+] - 3[(v2,^Vj+1,fc + )(TW,3+1,fc - )/^2,+ +

3

+ (V2 ,N,j-i,k + V2,N,j,k) (TN,j-i,k TN,j,k)/h2,

j = 1,...,M - 1, к = 1,... ,L - 1;

(2.22)

+ [3(V3 ,N,j,k+i + V3,N,j,k) x (TN,j,k+i TN ,j,k )/h33,k + 3(v3,N,j,k-i + V3,N,j,k )(Tn ,j,k-i - TN ,j,k)/h3,k-i}n ,

j = 1,...,M - 1, к = 1,... ,L - 1.

(2.23)

На остальных пяти гранях параллелепипеда R значения искомых величин T считаются известными: например, они заданы граничными условиями первого рода в виде функций от t и xm, m = 1, 2, 3.

3. Об экономичности метода

Представляет интерес число арифметических действий для нахождения разностного решения и аппроксимация. Внутри параллелепипеда R получается Si = (N- 1)(M- 1)(L-1) разностных уравнений (2.11) - (2.13) или (2.11), (2.14), (2.15) для определения Si неизвестных T.

На граничных плоскостях Xm, m =1, 2, 3 выписываются порядка S2 = ^-^(M-^^^-1)(L-1) + (M- 1)(L- 1)+N+M-1+L-1 разностных схем для определения S2 неизвестных. На граничных поверхностях xm = dm, m = 1, 2, 3 имеем также S3 = NM + NL + ML конечных алгебраических выражений (2.19) для определения S3 известных функций (2.3).

Легко видеть, что расчет по однородным разностным схемам (2.11) - (2.13) экономичен (с учетом условий абсолютной устойчивости (4.9), (4.11)), так как для получения решения в узлах области определения понадобится O(Si) арифметических действий, пропорциональное числу узлов R.

Если существуют ограниченные вплоть до четвертого порядка производные по пространству и второго порядка по времени, то для уравнения (2.1) локальная погрешность аппроксимации разностной схемы (2.11) - (2.13) выражается формулой O[^2^a=i(hav -ha,v+ г] (v = i при а = 1, v = j при а = 2, v = к при а = 3). Для hm = cm (cm = const, m = 1, 2, 3) локальная погрешность аппроксимации тех же схем — (т + hi + h2 + Щ) ; для w = Um = 0 и Vm = Cm на решении задачи — 0[т + Y;3m=i(h'm,q - hm,qhm,q-i + h2mq-i)} при hi,i - hi,i-i = h2,j - h2,j-i = h3,k - h3,k-i.

На граничной плоскости xi = 0 (0 < xm < dm, m = 2, 3) локальная погрешность аппроксимации разностных уравнений (2.16) при hm = cm — 0(т + hi + h2 + h3) (имеется в виду первое уравнение системы (2.16), так как остальные два уравнения необходимы для устойчивости схемы [6]). При um = w = 0 и Vm = cm погрешность аппроксимации этой схемы на решении задачи 0[т + hi,0(hi,0 + h2,j - h2,j-i + h3,k - h3,k-i)].

На оси ОХ1 (0 < х1 < локальная погрешность аппроксимации разностных схем (2.17) при Нт = ст — 0(т + Н + Н2 + Н3), а при ит = и> = 0 и = ст погрешность аппроксимации этой схемы на решении задачи — 0[т + ^¿(Н^ — 1 + Н2,0 + Н3,0)].

В начале координат = 0, т = 1, 2, 3 локальная погрешность аппроксимации разностных уравнений (2.18) — 0(т + Л1)0 + Н2,0 + Н3,0).

Не представляет труда найти, что локальная погрешность аппроксимации разностных схем (2.21) - (2.23) при = ст на решении задачи имеет второй порядок точности по пространству и первый по времени.

Начальное условие (2.2) удовлетворяется точно. Граничные уравнения (2.19) на соответствующей граничной плоскости из (2.3) аппроксимируются точно, так как они заданы от времени явно на целом слое.

4. Устойчивость по начальным данным

Для исследования устойчивости разностных схем (2.11) - (2.13) выпишем типичную схему ИИМ [1, 2] для внутренних узлов параллелепипеда Л, используя соотношения (2.5), (2.11), (2.14), (2.15) и результаты пункта 1:

Т—1 [Лм— 1#г— 1.,к/т — 3(«1,г— 1.,к + )/^1,г—1 + 2« к ] +

+Т.1/3 [4йм£г.,к/т + 3(«М— 1,.,к + «м.,к)/НМ—1 + 3(«1)г+1,.)к + «М.,к)/^1,г + «1,.+1 ,.,к - «1,.— 1,.,к] +

+Т+11?(3[НМ Х ^г+1,.,к/т - 3(«1,г+1,.,к + «1,г,.,к)/Н1,. - 2и1,г,.,к - и1,г+1,.,к] =

= Лм_1(5Т/т + / )П— 1.,к + 4Й1,.(5Т/т + /). + Нм(уТ/т + / )П+1,.,к+

+6Я1,г(Л2 + Л3 + А2 + А3 + Л1—3)ТП ,

г = 1,... - 1, 3 = 1,... ,М - 1, к = 1,... ,Ь - 1; (4.1)

ТГ,-— 1,к [Н2

—1^.,.—1,к/т 3(«2,г,, —1,к + «2,.,., к )/Л2 1,к ] +

+Т.к2/3[4^2,.£г.,к/т + 3(«2,.,. —1,к + «2.)/Н2,.— 1 + 3(«2,.,.+1,к + «2,г.,к)/Н2. + «2,г.-+1,к -«2,.,. — 1,к] +

+ТГ+21/к[Л2,. Х ^»,.+1,к/т - 3(«2,.,.+1,к + «2,г.,к)/Н2. - 2и2,г.,к - и2,г.+1,к] =

= {[Л1,.— 1(£ТГ+1/3)г— 1,.,к/т + 4^1,г(^Тга+1/3)г.,к/т + НМ X х(£ТГ+1/3)г+1.,к/т]/6^1,. - (Л2 + А2)ТГ + (А2 - Л)(Т + т/)Г},

г = 1,..., N - 1, 3 = 1,..., М - 1, к = 1,..., Ь - 1; (4.2)

1 [3(«3,г.,к—1 + «3,.,.,к)/Н3,к—1 - Н3,к—1^г,.,к—1/т] +

Т. X [3(«3,.,.,к—1 + «3,.,.,к)/Н3,к—1 + 3(«3,.,.,к+1 + «3,.,.,к)/Н3,к + 4^3,к х /т] +

+TÍJ+fc1+1[3(v3,г,j,k+1 + «3,.,.,к)/Н3,к - Н3,к^г,.,к+1/т] =

= 6^3,к {[Н^—^Т^3).. —1,к/т + 4^2,. ($Т га+2/3)г,.,к/т + Н2. X (^ТГ+2/3)г,.+1,к/т]^. -

-(Л3 + А3)ТГ + (А3 - А2)(Т + т/)Г}, г = 1,.. .N - 1, 3 = 1,... ,М - 1, к = 1,... ,Ь - 1. (4.3)

Методом Фурье в приближении замороженных коэффициентов [6, 10] можно найти условие устойчивости явно-неявных разностных уравнений (2.11) - (2.13) или (4.1) - (4.3) по начальным данным при vm = um = hm = const, m = 1, 2, 3, g = 1 и f = 0. Представляет интерес оценки величины п — множителя роста разностных схем ИИМ (4.1) - (4.3). Сделаем стандартную подстановку [10]:

Tn+1 = nTn , Tj = exp[/ (gi xM + g2X2,j + g3X3,fc )] , (4.4)

где / = \J — 1, gm = 0, ±1, ±2,..., m = 1, 2, 3 — гармоники при переходе со слоя на слой. Подставим (4.4) в (4.1) — (4.3) и для сокращения дальнейших выкладок введем обозначения:

= 4rvmhm2 sin2 gmhm/2 , 6m = rumh^1 sin gmhm,

rm =(1 + 2 cos2 gmhm/2)/3 , m = 1, 2, 3 ,

Ci = T(hih2)-1 sin gihi sin g2^2 , C2 = T(^3^2)-1 sin g3^3 sin g2^2 ,

C3 = T(hiЛ.3)-1 singihi sing3^,3 , c = 2w(ci + C2 + C3). (4.5)

Тогда по каждому из направлений получим

Tn+1/3(ri + «i — /61) = Tn(ri — «2 — «3 + /62 + /63 — c), (4.6)

Tn+2/3(r2 + «2 — /62) = riTn+1/3 + Tn(r2 — ri + «2 — /62) , (4.7)

Tn+1(r3 + «3 — /63) = r2Tn+2/3 + Tn(r3 — r2 + «3 — /63). (4.8)

Выразим Tn+1/3 из (4.6) и подставим в (4.7), затем найдем Tn+2/3 из полученного соотношения и подставим в (4.8). В результате окончательно имеем уравнение для определения п

П = (A — /B)/(C — /D), |п| = [(A2 + B2)/(C2 + D2)]0'5, (4.9)

где

A = F + s — c, C = F + a, D = E + 6, B = E + q, /2 = —1,

F = r3r2ri + г1(й2й3 — 6263)+ г2(й1й3 — 6163)+ ^(«2 «1 — — 62 6i) + «^«3 — «16263 — «26163 — «36261 ,

s = Й1Г2(г3 — ri) + Г1Й2(г3 — Ы , ü = Г3Г2Й1 + r3Й2Г1 + Й3Г2Г1 ,

E = ri(62«3 + «263)+ r2 (6i «3 + «163)+ ^(62^ + «261) +

+ «2«36i + ^«362 + ^«263 — 616263 ,

q = 61 r2(r3 — ri) + ri62(r3 — r2) ,

6 = r3r26i + r362ri + 63r2ri. (4.10)

Для двухслойных разностных схем необходимое условие устойчивости по начальным данным записывается [6, 10] в виде |п| < 1. В силу того, что | singaha| < 1, | cosgaha| < 1, найдем оценки сверху для п. В частности, знаменатель будет минимален в (4.9) при gmhm/2 = kn (m = 1, 2,3, k = 1, 2,...). Тогда из выражений (4.5), (4.9), (4.10) имеем rm = 1, singmhm = singmhm/2 = 0,m = 1, 2, 3, s = q = c = 6 = « = 0, E = 0, F = 1, A = C = F и

|п| = F/F = 1. (4.11)

Из (4.9), (4.11) видно, что разностная схема (2.11) - (2.13) или (4.1) - (4.3) безусловно устойчивая т > 0. Аналогично доказывается безусловная устойчивость остальных разностных схем (2.16), (2.18), (2.21) - (2.23). При отсутствии конвективных и смешанных производных (um = w = 0, m = 1, 2, 3) множитель роста (4.9) совпадает с множителем роста разностной схемы полученным в [6] для трехмерного параболического уравнения при vm = const, m =1, 2, 3, rm =1, g = 1, f = 0

£ = (1 + ai«2 + «3^2 + «1^3 + «1 а2«э)/[(1 + «i)(1 + «2)(1 + «з)] •

Как отмечено в монографии [10], способ "замораживания" коэффициентов схемы обоснован для некоторых классов параболических уравнений с гладкими коэффициентами (в ряде случаев достаточно непрерывности коэффициентов). Отметим, что критерий устойчивости, полученный этим способом, хорошо согласуется с результатами численных расчетов, приведенных ниже при решении неодномерного параболического дифференциального уравнения с постоянными коэффициентами переноса.

5. Устойчивость прогонки

Так как матрица коэффициентов уравнений (4.1) - (4.3) трехдиагональна, то для нахождения неизвестных воспользуемся методом скалярной прогонки [5]. Рассмотрим условия устойчивости прогонки для типичной схемы ИИМ (4.1) - (4.3), для которой имеет место диагональное преобладание. Воспользуемся сокращенной записью уравнения(4.1) при g = vm = um = hm = const и f = 0

—i + BTi + CTi+1

= F, (5.1)

где A = —6v1 /h2 + g/т + 3u1/h1, B = 4g/T + 12v1/h? , C = -6v1/h2 + g/т - 3u1/h1, F = A1Tn + б(Л2 + Лз + Д2 + Аз + Л1—з)Тп; тогда оценки [5]: A > 0, C > 0, B > A + C -для устойчивости прогонки в уравнении (5.1) дают условия

а) C = —6v1/h2 + g/т — 3u1/h1 > 0 и T<g/(6v1/h2 + 3u1/h1), (5.2)

поэтому тем более т < g/(6v1/h2 — 3u1/h1) при A > 0;

б) если решение уравнения (2.1) обладает свойством монотонности, то важно сохранение этого свойства в разностных схемах [5]. Поэтому при A< 0, C< 0 и B > 0 имеем

т > g/(6v1/h2 — 3VM • (5.3)

Если величина конвективного переноса значительно превышает кондуктивный перенос, то для выполнения условия (5.3) целесообразно применить принцип монотонизации [5, 11] разностной схемы (5.1).

Условие (5.2) может быть практически выполнено, когда граничные условия второго рода из (2.4) ненулевые. Например, задан тепловой поток большой интенсивности: > 106 Вт/м2. Последнее может приводить к частому дроблению шага по времени т при решении задачи с заданной точностью.

Кроме того, это условие для уравнения теплопроводности не является обременительным, так как g = pcp ~ (1 — 3, 8)106 Дж/(м3- К), v1 = 386 Вт/(м- К) для тела, например, из меди, а шаг по пространству в трехмерном случае не удается взять менее 0, 01 м при

dm = 1 м, m = 1, 2, 3 и двойной точности (сетка 100 х 100 х 100) даже в современных ПЭВМ (Pentium 16 Мбайт) из-за наличия в общем случае 10 трехмерных массивов (vm, um, g, f, Tn+1, Tn, m = 1, 2, 3) в системе разностных уравнений (2.11) - (2.13). Аналогичные оценки (5.2), (5.3) имеют место для разностных уравнений (4.2), (4.3).

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

Сходимость разностных схем ИИМ (2.11) - (2.13) установим практически на последовательности сгущающихся сеток по т и Кт, т = 1, 2, 3 при решении частного случая трехмерного уравнения (6.1) со смешанными условиями (6.2)

дТ

дж

Т U

Т |Х2=0 Т 1жз=0 Т 1Х2 = 1

д 2Т

3

Va дх2а

а=1 a

+ [1 + Y;[XZa - Z(Z - 1)VaxZa Z]} exp(t)

(6.1)

a=1

1 + Y, xa, т|

xi=0

(1 + x2 + xZ) exp(t)

= (1 + x\ + xZ)exp(t),

a=1

Z)

= (1 + x\ + xZ)exp(t),

дТ

dx1

= z exp(t)

xi = 1

(2 + x\ + xZZ) exp(t), T|хз=1 = (2 + x\ + xz2) exp(t)

(6.2)

Для решения тестового примера воспользуемся разностными уравнениями (4.1) - (4.3), (2.20) - (2.23). Точное решение задачи (6.1), (6.2) Т = rexp(t) в кубе Q: [0 < xm < 1,m = 1, 2, 3] известно при 0 < t < tk. Были взяты следующие значения входных данных: hm = h, vm = 1, m = 1, 2, 3, Z = 6, g =10, т = 0, 002. Программа составлена на языке Фортран-77, расчет проводился на ПЭВМ Pentium (Транслятор Power Station 1) с двойной точностью. В таблице дается максимальная относительная погрешность е = |Т—Т|100%/Т (Т - точное, Т - приближенное численное решение) в момент времени tk = 1. Видно, что численное решение задачи сходится к точному при измельчении шагов разностной сетки (t0 — расчетное время в минутах на ПЭВМ).

Таблица

з

з

т 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

е, % 5, 23 1, 75 0, 458 0,128 0, 658 2, 67

t0 0, 9 1, 9 4, 75 30, 0 0, 75 0, 37

7. Заключение

На основании вышеизложенного можно сделать следующие выводы.

1. Для численного решения трехмерного нелинейного параболического уравнения общего вида на основе ИИМ получены разностные схемы с погрешностью аппроксимации в обычном смысле на равномерных сетках внутри параллелепипеда 0(т + Н\ + + К2).

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

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

Отметим, что положительные свойства разностных схем, получаемых с помощью ИИМ (сплайн-аппроксимация искомых решений и др., возможность получения схем повышенной точности) для одномерных краевых задач, сохраняются и при решении многомерных задач. В частности, если в качестве начального приближения в логической схеме из [2] вместо линейной функции взять параболу, то для уравнения (1.1) по алгоритму ИИМ [2, 4] можно получить внутри области определения дифференциально-разностную схему с погрешностью аппроксимации O(h^ + h2, + h33).

Список литературы

[1] Гришин А. М. Об одном видоизменении метода М. Е. Швеца. Инж.-физ. журн., 19, №1, 1970, 84-93.

[2] Гришин А. М., Берцун В.Н., ЗинчЕнко В. И. Итерационно-интерполяционный метод и его приложения. Изд-во Томского ун-та, 1981.

[3] Алберг Д., Нильсон Э., Уолш Д. Теория сплайнов и ее приложения. Мир, М., 1974.

[4] Якимов А. С. Об одном методе расщепления. Численные методы механики сплошной ореды, 16, №2, 1985, 144-161.

[5] Самарский А. А. Введение в теорию разностных схем. Наука, М., 1971.

[6] ЯнЕнко Н. Н. Метод дробных шагов решения многомерных задач математической физики. Наука, Новосибирск, 1967.

[7] РихтмАйЕр Р., Мортон К. Разностные методы решения краевых задач. Мир, М., 1972.

[8] Гришин А. М., Фомин В. М. Сопряженные и нестационарные задачи механики реагирующих сред. Наука, Новосибирск, 1984.

[9] Гришин А. М. Математическое моделирование лесных пожаров и новые способы борьбы с ними. Наука, Новосибирск, 1992.

[10] КАлиткин Н. Н. Численные методы. Наука, М., 1978.

[11] Самарский А. А., Гулин А. В. Численные методы. Наука, М., 1989.

Поступила в редакцию 18 декабря 1997 г., в переработанном виде 14 декабря 1998 г.