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

УДК 536.2

В. В. Г о р с к и й, Н. А. Г о р с к а я, В. Г. Р е ш

МЕТОДИКА РЕШЕНИЯ УРАВНЕНИЙ МАТЕМАТИЧЕСКОЙ ФИЗИКИ НА ОБЛАСТЯХ СЛОЖНОЙ ФОРМЫ МЕТОДОМ МАТРИЧНОЙ ПРОГОНКИ

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

E-mail: resch883@rambler.ru

Ключевые слова: уравнение Фурье, подвижная внешняя граница, нестационарный прогрев, метод матричной прогонки, неравномерная сетка.

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

Физико-математическая постановка задачи. Рассматривается задача нестационарного прогрева затупленного тела вращения, внешние контуры которого складываются из поверхности переменной кривизны, цилиндрического хвостовика и плоского заднего торца (рис. 1).

Для обозначения этих фрагментов поверхности будем применять подстрочные индексы w, cyl и flat, а для обозначения стыка первых

'cyl "

Zcy1 zflat Z

Рис. 1. Затупленное тело вращения

двух из них — надстрочный индекс *. При этом допускается перемещение во времени первого из указанных фрагментов поверхности параллельно оси симметрии тела. Взамен термина «поверхность переменной кривизны» будем использовать термин «стенка».

Температурное состояние рассматриваемого тела в цилиндрической системе координат описывается двумерным нестационарным уравнением Фурье вида [1]

2 дТ д ( . дТ ) д ( . дТ pcrs — = —I rAz— н--1 таг-

дт dz 1 dz ) dr l dr ,

(1)

Z G

(Zw(rZflat ), r

(o, rw (z,r)), z < / (r),

(o, rcyi), z > z* (r^ ге(a rend].

Граничные условия для уравнения (1) задаются следующим обра-

зом:

• на стенке

4 IT- Sin (3) + ЛГ |Т- COS (¿) = sQw (r,r), r g[0, rcyi ), rG(0,Tend ]; (2)

• на стыке стенки и поверхности цилиндра

дТ

= sQw (r, r)sin (.9),

5z

r = rcyb rG

(0,rend ];

(3)

, дТ Är — = s dr

Qw (r, r)cos (3) + Qcyi (z\r)],

r _ rcyl, r6( o

, ^"end I; (4)

• на плоском заднем торце дТ

4 = Фш (r, г е [О, гсу1 ], ге( 0, гепй ];

• на оси симметрии тела

дТ

— = 0, г е^0^ Гда I г£(0,Гепй ];

• на поверхности цилиндра

сгТ

Л -^Т = ^су1 (z,т), г е (^ (гсу1, ], 0, Гепй ].

В свою очередь, начальное условие задается в виде

(5)

(6)

(7)

Т{z, г,0) = ^ Z е[Zw(r,0), f], r

t0, rw М)], z < z (0);

Г "I * / ч (8)

[0, ^] , z > z (0).

В формулах (1)—(8): р — плотность, кг/м3; с — удельная теплоемкость, Дж/(кг • K); г, г — оси цилиндрической системы координат, направленные соответственно вдоль оси симметрии тела и по нормали к ней, в долях от характерного линейного размера тела; ^ — характерный линейный размер тела, м; Т(г, г, г) — температура, К;

г — временная координата, с; гсу1 — радиус цилиндра; гепй — верхняя граница рассматриваемого временного интервала, с; Л — теплопроводность, Вт/(м • К); 3 — угол между касательной к образующей поверхности тела и осью 2, рад; Q —удельный тепловой поток, подводимый по нормали к поверхности тела, Вт/м2; ТЬщ —

температура тела в начальный момент времени, К.

Ось г направлена внутрь тела, а ее нулевое значение совмещено с точкой пересечения оси симметрии тела с поверхностью его затупленной части в нулевой момент времени. Индекс теплопроводности указывает ось, вдоль которой она характеризует интенсивность кон-дуктивного теплопереноса.

Переход к подвижной системе координат. В подавляющем большинстве случаев численное решение задач рассматриваемого класса находят в исходной системе координат (назовем ее неподвижной) с использованием явных или неявных расчетных схем; в последнем случае обычно применяют метод переменных направлений [2]. К недостаткам использования таких подходов к решению задач, характеризующихся достаточно сложной формой расчетной области, относят:

• сложность формирования конечно-разностных операторов в окрестности границ этой области, обладающих вторым порядком точности аппроксимации (особенно при подвижности части из указанных границ);

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

Одним из возможных путей решения указанных проблем является переход от неподвижной системы координат | г, г, г] к подвижной

системе координат {х, у, г'), которая для рассматриваемой задачи имеет вид

г ~ ^ (г, ^ , х =----, у = г, г=г.

Тогда

аг _ т УУ -Т _ 1^x z

дт дт' дт дх дт ду дт дт'

дТ

а w,r ах

дТ _ дТ дт^+дТ_ дх + д£ су _ ±_ д£; dz дт' dz дх dz dy dz А dx'

dT _ dT dTf ^dT dx ^dT dy dr dT' dr dx dr dy dr

1 - x

dx dy'

(9) (10) (11)

rX Ю-iA + \ ik.'T^ + rk dTTy

dz ^ z dz J dT\ A dx J dz dx [ A öi J dz dy [ Д Si J dz

1 Af Ä T1 •

A dx l A dx )

(12)

4 —

dr l dr

d ( „ dT Л dT' df„ dT Л dx df

dT'

rl

dy

■ + dr dx

rXr--+

dy J dr dy

rl

dT_ ^

dy J dr

5r'lr r A zw,r ax J dr dxlr r A zw,r ax )dr

d ( „ 1 -x dT^dy --1 rAr-zwr- I —

dy v А ' dx J dr

d ( . dT> y\ —

dy

dy

1 - x . 1-x dT_

+ A zw,r & ly r A zw,r 5x

1 - x

d ( „ dTЛ о

dx

1 - x

dT

у r dy , dy ^у r A w,r dx

(13)

где zw r = ■

dzw

abl

dT sin (.9);

(V — линейная скорость абляции ма-

териала, из которого изготовлено тело, по направлению внутренней нормали к поверхности последнего, м/с); г, г А (у, г) =

dr

= z

flat

" zw (у , 7).

С учетом выражений (9)—(13) и очевидного равенства

dTЛ1-х

d_

dx

yÄr а zw,r ax

w,r

1 -X _д_( 1-х дГ\_ 1-х 2 дТ_

AZw,r I л Zw,r I . 2 Zw,r ^

cX I А сх J А2 дх

уравнение (1), граничные условия (2)—(7) и начальное условие (8) можно записать в виде

дТ д ( дТ Л д ( дТ ^ дТ " ' " + У^з

У%0 ~ У ~ I n I ^ Л.

от ох ^ ох) ду ^ ду

дх

-(1 - х

(

дх

%4 '

дТ

\

д ( дГ

Хл

дх

ду) ду

х е(0,1), у е(0, Усу!), r'e(0,rend] 2 1

Z0 = Pcs , Z1 = 77 А2

4 +Är (1 - х)

22 z

"О х)[ Z0Zw,t Л +^rZw,r ,2 I, Хл -

Z2 = у К,

уЛ

(14)

При х = 0 дТ 1

дТ

[Л; sin(3) + Лг cos(S)zw,r] +—Ar cos(5) = sQw, (15) дх A L ду

у G[0, усу!), 0,re'nd ];

1 дТ

~Лz А~дх = SQw у = ^У^ r'G(0,rend ]; (16)

Ar f = s к (уУ)

cos ($)- sin

+ 0су(х,г')к (17)

у = Уcyl, 0, ^end ]■

При х = 1

/L дТ

-Т-¡Г = öflat {у,т)s, у е[0, Уcyl ] , 0,<nd ].

А дх При у = 0

дТ

= 0, х е[0,1], r'E(0,rend ].

(18)

(19)

При у = Уcyl

дТ 1 - х дТ

Л- — = Qcyl (х,т)* + Аг—гЩг — > х е(0,1], г'е^,^]. (20)

При т' = 0

Т = ТЬе§, х е [0,1], у е [0, Усу1]. (21)

При преобразовании выражения (6) в вид (19) использовалось свойство затупленного тела вращения, согласно которому на его оси

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

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

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

Далее подробно рассматриваются вопросы, связанные с решением уравнения (14) именно методом матричной прогонки. Однако предварительно остановимся на вопросе использовании неравномерных сеток при интегрировании многомерных уравнений теплопроводности независимо от того, в какой форме эти уравнения записаны.

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

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

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

Суть этого подхода проиллюстрируем на примере уравнения (14), при решении которого целесообразность использования неравномерных сеток обусловлена следующими требованиями:

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

• шаг интегрирования по координате у желательно изменять таким образом, чтобы он был обратно пропорционален производной

по этой координате, которая имеет тенденцию к увеличению с ростом координаты у.

Указанным требованиям, в частности, соответствует преобразование системы координат вида

{х, у}^, , (22)

где

dx = p(x) = (1 + xf, nx < 0, x e[0,l], ^(0) = 0; (23)

d£ dy

= ¥{y) = {l + У Г, ny > 0, y e[0, ycyl ], £(0) = 0. (24)

Тогда

7](x ) =

ln(l + x), nx =-1, l

nx + 1

(1+x)"x+1 -1 , nx ^e[0,^fiat];

(25)

= ТГГ^+y)ny+1 -1], ^cyi ];

ny + 1

(26)

%at =^(1) =

in(2), nx =-1

1 1 [

2nx+1 _ 1

nx +

, nx * -1;

(27)

4yl =фсу! ) = n^^Y i1 + ycyl)

ny+1 _ 1

(28)

x

w=

exp (77)-1, nx =-1,

[^K +1) +1 ]1/("x+1)-1, nx *-1, T7e[0,^ ];

(29)

y Й = [Фу +1

+1

7(„y +1)

■1, 0, 4yl ].

(30)

Изменив теперь в уравнении (14), в граничных условиях (15)— (20) и в начальном условии (21) систему координат в соответствии с выражениями (22)—(30), а также заменив т' на г, получим окончательный вид физико-математической постановки задачи, решаемой ниже методом матричной прогонки:

дТ

yZ0^~ = yP-

от

о

(

дТ

\

о

(

ХгУ

дТ_

Л

+ УХзФ

дТ

дц

-(1 - х )(ру/

с

(

%4

дТ

Л

с

V дц)

%4 '

дТ

(31)

^G(0,%at ), £G(0, 4yl), Ге(0, Tend ].

При ri~ 0

\ tÄz Sin ^ + Är c0S ^ zw,r ] + c0S ^ = SQw ,

£e[0, 4yl), re(0,Tend]•

Ф dT

~Äz A^ = sQw r)sin (<9), £ = 4yi, r'e(0, rend];

XrWly = S |Qw ^ ^

СС8 (5)-8Ш

+ Qcyll *■)!

£ = 4yb 0, Tend ].

При 77 = T7fiat

Äzp dT A d^ При£ = 0

: Qflat (£ S £ G 4yl ], T G (0, ^"end I.

dT_

= a ^e[0,%at], a Tend].

При Z = l

cyl

Л ( dT 1 - x

v-

-zw.rV

dT_

d^

= Qcyl T), I0, ^flat ], G (0, ^end ].

А

При г = 0

Т(ъ £ г) = ТЬе§, £е[0,

Решение сформулированной задачи будем искать на сетке

Ч = М*"^ 1 = 1, hn =

f •

N„-1'

^=-1), j=1, n^, ^ =

1

cyl

NP -1

с использованием чисто неявной схемы аппроксимации пространственных дифференциальных операторов.

Метод матричной прогонки. Конечно-разностный аналог уравнения (31) на текущем шаге интегрирования по времени в окрестности расчетного узла области определения функции Т(г), £,,т) с координатами , ^ записывается в каноническом виде следующим образом:

¡2 _ _

Е Е »/, 1,1А-2,1+т-2 + Рг,} = 0, i = Щ, 1 = Щ, (32)

I =¡1 Ш=Ш1

где щ 11т, Д 1 — коэффициенты, рассчитанные по начальному приближению к искомому температурному профилю на текущей итерации.

В зависимости от месторасположения расчетного узла граничные значения параметров суммирования I и т принимают различные значения:

при / = 1 ¡1 = 2, ¡2 = 3; при 1 = 1 т1 = 2, т2 = 3;

при I = 2, 1 ¡1 = 1, ¡2 = 3; при 1 = 2, Ы^-1 т1 = 1, т2 = 3;

при I = Ы^ ¡1 = 1, ¡2 = 2; при 1 = Ы| т1 = 1, т2 = 2.

Согласно методу матричной прогонки, искомое решение рассматриваемой задачи представляется в виде [3]

Т,1 = Е 4,мТ+и + Ви, / = 1, ыч- 1, 1 = 1ГЫ7, (33)

к=1

где 4 1 к, Б{ 1 — прогоночные коэффициенты.

Если / = 2, Ыл -1, то конечно-разностный аналог (32) уравнения Фурье имеет вид

т^ _

Е {а', 1,1,тТ-1,1+т-2 + а1,1,2,тТ1,1+т-2 + Ц', 1,3,тТ+1,1+т-2 ) + Д, 1 = 0, 1 = 1, .

Исключим из этого уравнения Т-11+т_2 с применением (33). Тогда

(ы* \

m

Z

Ц, m

v k=1

Е 4-1, j+m-2,kTi, к + Bi-1, j+m-2

+ д. j = о, j = 1, N

+ а T T

^j,2,mJi,j+m-2 ^j,3,mJi+1,j+m-2

или

N i m 1

Z T,k \ Z lai,j,1,m4-1,j+m-2,k + «i,j,2,mÖU + m " 2 k)J f =

k=1 | m=m Nf m2

^ m2 ^

= "Z Ti+1, k Z ai, ] ,3, m^ (j + m - 2, k) - Д,] + Z aijXmBl-lj+m-2

k=1 m=m1 V m=mi

7 = 1,

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

Если г = 1, то конечно-разностный аналог (32) уравнения Фурье имеет вид

'"2

Z К j,2,mTi, j+m-2 + j,3,mTi+1,j+m-2

) + Д. j = 0, j = 1, Nü,

m=mj

что равносильно выражению

Nf

Z T,k Z ai, j,2,mÖU + m - 2, k) =

k=1 m=mj

N.

m2

= -Е Т+1,к Е 7,3,™^ (■/ + " 2 к) - Д,7 , 7 = 1, N4-

к=1 т=т1

Решение этого набора уравнений позволяет определить значения прогоночных коэффициентов А1]к и Б1]-.

Если г = Л^, то конечно-разностный аналог (32) уравнения Фурье

имеет вид

Z [ai, j,1,mTi-1, j+m-2 +ai, j,2,mTi, j+m-2 y

,) + Д. 7 = 0, 7 = 1, (34)

т=т1

Тогда, исключив Т+т_2 из (34) с применением (33), получим,

что

X

( N.

а

i, j ,1,m

Z A-1, j+m-2,kTi,k + Bi-1, j+m-2

k=1

+ a T

^ "i,j,2,m^i,j+m-2

+ ßu j =0,

J = 1, N^.

m

m

При этом искомая система линейных алгебраических уравнений, предназначенная для расчета температур на плоском торце,

Nf

Ё Ti ,k \ S , j XrnAi-1, j+m-2,k + 4, j XmÖ(j + m ~ 2 k )] \ = к=1 I m=mi

Pi, j + S j,1,mBi—1, j+m-2 , j = 1, Ni .

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

Исходные данные: S = 0,1 м; zcyl = 3; zflat = 5; Xz = 50 Вт/(м-К); Xr = = 50 Bt/(m-K); с = 2 000 Дж/(кг-К); p = 1 500 кг/м3; Qw(0) = 107 Вт/м2; Qw(0) = Qw(0) • sintf, Bt/m2; Qcyi = Qfiat = 0; Theg = 300 К; Tend = 50 c; N„ = 99; Nc = 99; nx = -5; ny = 0; h = 0,01 c.

Расчет выполнен для двух значений zw,r, не зависящих от координаты r.

На рис. 2—4 приведены результаты расчетов при zw,x = 2 • 10-3м/с. На рис. 5 и 6 показано изменение тех же показателей, что и на рис. 3 и 4, но при zw,r = 0.

^^ * ^

. • * г> -*

/ • / / / / < / / / / / - • -2

/ ' ■ / / f 1 | -4

1 ! i • 1 1 1

О 0,5 1 1,5 2 2,5 3 3,5 X

Рис. 2. Изменение геометрии тела в процессе обгара:

1 — т = 0; 2 — 10; 3 — 30; 4 — т= 50 с

1 \

* 1 'Л ч _ _ 1 2 3

1 » 1 \ \

1 1 \ 1 < * \ \

600 \ 1 * ч \

\ Ч ч ч \ s

о о, га о.ое о,ш 012 х

Рис. 3. Распределение температуры Т по оси тела при = 210-3 м/с:

1 — т = 1 с; 2 — 10; 3 — г= 50 с

----^

V V V V _\

----- " ■—, . — \

\ \ \\

0 0,3 0,6 0,9 1,2 у

Рис. 4. Распределение температуры Т на стенке при = 210-3 м/с:

1 — т = 1 с; 2 — 10; 3 — г= 50 с

\

\

\ V \ \ \ V ---2

\ \ \ \ i \

\ \ i \ 1

ч ч 1 1

\ \ \ \ N

\ \ \ \ \ ч ^

N. --------

О жв о® аов ас х

Рис. 5. Распределение температуры Т по оси тела при = 0:

1 — т = 1 с; 2 — 10; 3 — г= 50 с

О 0,3 0,6 0,9 1,2

Рис. 6. Распределение температуры Т на стенке при zw,r = 0: 1 — т = 1 с; 2 — 10; 3 — г= 50 с

Выводы. Разработанный численный метод решения широкого круга задач математической физики предназначен для повышения точности решения задач на областях со сложной формой граничных поверхностей. Подвижность части этих поверхностей не является принципиальным моментом данного подхода.

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

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

Работа выполнена при финансовой поддержке РФФИ (проект № 2010-01-00841а).

СПИСОК ЛИТЕРАТУРЫ

1. Беляев Н. М., Рядно А. А. Методы теории теплопроводности. - М.: Высш. шк., 1982. - Ч. 2. - 304 с.

2. Самарский А. А. Введение в теорию разностных схем. - М.: Наука, 1971. -550 с.

3. Пасконов В. М. Стандартная программа для решения задач пограничного слоя // Численные методы в газовой динамике / Под ред. Г.С. Рослякова и Л.А. Чудова. - М.: Изд-во Моск. ун-та, 1963. - С. 110-116.

4. Демидович Б. П., Марон И. А. Основы вычислительной математики. -М.: Физматгиз, 1963. - 659 с.

Статья поступила в редакцию 03.07.2012.

Авторы статей

Акинкин Дмитрий — студент Аэрокосмического факультета МГТУ им. Н.Э. Баумана.

Апельцин Виктор Филиппович — канд. физ.-мат. наук, доцент кафедры «Вычислительная математика и математическая физика» МГТУ им. Н. Э. Баумана.

Бабаев Ислам Акмурадович — инженер-испытатель ОАО «Ракетно-космическая корпорация «Энергия» им. С.П. Королева».

Боровин Геннадий Константинович — д-р физ.-мат. наук, профессор, зам. директора ИПМ РАН.

Бродский Юрий Игоревич — канд. физ.-мат. наук, доцент МГТУ им. Н.Э. Баумана, вед. науч. сотр. ВЦ РАН.

Бутина Татьяна Александровна — канд. физ.-мат. наук, доцент кафедры «Вычислительная математики и математическая физика» МГТУ им. Н.Э. Баумана.