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

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

Том 3, № 6, 1998

НЕЯВНЫЙ МЕТОД ДРОБНЫХ ШАГОВ

С РАСЩЕПЛЕНИЕМ СМЕШАННЫХ

ДИФФЕРЕНЦИАЛЬНЫХ ОПЕРАТОРОВ

В.Ф. Формалев, О. А. Тюкин Московский государственным авиационный институт, Россия e-mail: [email protected]

A new resource-sparing, totally stable method of fractional steps for the numerical solution of parabolic problems with mixed differential operators which are split on an asymmetric pattern into coordinate directions is suggested. Theorems of approximation are proved as well as those of stability and convergence. The results of comparative numerical experiments are presented in comparison with the well-known classical methods.

Идеи академика Н. Н. Яненко по разработке полностью неявных и в то же время экономичных схем численного решения задач механики сплошных сред привели его к созданию метода дробных шагов [1], широко используемого в настоящее время. Однако он неоднократно высказывал идеи о более полном расщеплении задач механики сплошных сред не только по независимым переменным и физическим процессам, но и о расщеплении отдельных дифференциальных операторов [2].

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

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

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

1. Введение

3

а,в=1

© В.Ф. Формалев, О. А. Тюкин, 1998.

u(x, О) = u0(x), x Є Q, t = О, u(x, t) = О, x Є dQ, О < t < T,

3

kав — ^ ^ kr ^та^тв, а, в — 1 2 3, kав — k,

r=1

3

3 3

C1 п0 < kавПаПв < C2 п0•

(2)

(3)

(4)

(5)

la'lp —

a=1 а,в=1 a=1

Здесь n — произвольный вектор nT = (n1, П2, Пз)> c1, C2 = const > 0.

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

2. Метод полного расщепления

На пространственно-временной сетке = {х = (ііЛ-і, і2^2, і3Л<3), іа = 0,1, ... , Да =

/а, а =1, 2, 3} , й)т = {£ = ^ + ат/3, tj = Іт, І = 0,1,... , І — 1; Іт = Т, а =1, 2, 3} рассматривается следующая конечно-разностная схема, аппроксимирующая дифференциальную задачу (1) — (5):

ыа =ЛаУ(а), x Є Wh; y(x, О) = uo(x), x Є Uh,

У(а) = О, x Є dwh, t Є WT, а =1, 2, З,

(б)

где

У*а

Ла

У(а) - У(а—1)

т

а—1

У(а) = yj+°/3 = y(x,tj і ат/3),

в=1

і

і

в—а+1

Лаау кааУхаха 5 а 15 3 , ^ав кав 5

Л +вУ = кавухвХа 5 Л+вУ = кавухвХа 5

Л«вУ = кавуХ'вХа 5 Л«вУ = кавуХ'вХа 5

а, в = 1 5 3 5 а = в-

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

= Л„у^ + (^ЛГ2 + !±£»Л«Л г/+‘/>.

т

2

2

3

у,+1 7,+1/2 =Л22у>+' + (^Л+ + 1+^-Л+) (7)

где

*-аау 111ааухаха} ■‘■42 У к12ух2х1 ,

Лаау кааухаха) Л12 ^+1/2 =

Л12У'+1/2 = к12уХ+х\/2, Л+^'+1 = ^21уХ2х\ ,

Л+1 ^+1 = к21уХ+х\; Став = кав, а, в = 1, 2.

£ ,

т

Рис. 1. Шаблон схемы метода полного расщепления в двумерном случае.

При этом, если скалярная прогонка осуществляется по переменной х1, то перемещение шаблона осуществляется вдоль другой переменной х2 и наоборот (на рис. 1 в первой подсхеме перемещение шаблона осуществляется в направлении от нижней границы к верхней, а во второй подсхеме — от правой границы к левой, показано стрелкой).

3. Аппроксимация

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

Будем считать, что задача (1)-(5) имеет единственное решение — функцию и(ж,£) £ С2 (3т), где С!т(<3т) — класс функций, имеющих непрерывные производные т-го порядка по I и п-го по х1, х2, х3 в области 3т = 33 х [0,Т]. Тогда верна

Теорема 1. Схема (6) обладает свойством суммарной аппроксимации на решении и(ж,£) задачи (1) - (5) с погрешностью аппроксимации 0(т + + Л| + Л|).

Доказательство. Характеристикой точности схемы является разность

г'+а/3 = у'+а/3 - и'+і. (8)

Подставляя (8) в схему (6), получим конечно-разностную задачу для нахождения погрешности

2£а = Лаг(а) + ^(х, 0) = 0, X Є Ой, (9)

£(а) = 0, X Є до^, і Є от,

где

^а = Лам'^1 — ^а>1м*, $ад — символ Кронекера.

Используя разложения Тейлора, нетрудно показать, что:

3

^а = 0(1), ^ = ^ ^а = 0(т + й? + Л| + Л|).

а=1

Следовательно, схема (6) обладает суммарной аппроксимацией и является аддитивной. Теорема доказана.

4. Устойчивость

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

3

(и, V) = м(х)^(х)Н, Н = ^ йа, ||м|| = а/(и, и),

хЄ^н а=1

получим гильбертово пространство сеточных функций. В этом случае операторы Лаа, Л,

Л-в, Л +в, Л+в можно рассматривать как линейные операторы, отображающие на (их

область определения и область значений совпадают со всем пространством С^), поскольку значения функции на границе равны нулю. Такая интерпретация разностных операторов позволяет вместо схемы (6) рассмотреть операторно-разностную схему

(Е — тЛа)у'+а/3 = у'+(а“1)/3, у(х, 0) = Мо(х), X Є оЛ,

исключая в которой промежуточные значения сеточной функции, получим

Ву'+? = у', у(х, 0) = мо(ж), X Є о^, (10)

где

В = (Е — тЛ?)(Е — тЛ2)(Е — тЛ3).

Используя тождество у'+? — у' = ту* , запишем схему (10) в каноническом виде [3]:

Ву* + Ау' = 0, у(ж, 0) = мо(х), х Є о^, (11)

где

Теорема 2. Пусть компоненты ^ав тензора задачи (1)-(5) удовлетворяют условию сильной эллиптичности (5). Тогда схема (11) устойчива по начальным данным, так что для решения задачи с начальным условием у(ж, 0) = м0(х) справедлива оценка

Цу'+І < ||у'|| < ... < ||у(х,0)| = ||мо(х)||. (12)

Доказательство. Поскольку выполнено условие сильной эллиптичности для коэффициентов ^ав, то оператор А > 0 и В = Е + тА > 0. Следовательно, существует оператор В-1. Следуя [3], умножим (11) слева на В-1, получим

у* + (Е + т А)-1Ау' = 0,

или, что то же самое,

у'+? = 5у', 5 = Е — т (Е + тА)-1А. (13)

Устойчивость схемы (11) эквивалентна выполнению операторного неравенства 5*5 < Е

(5* — оператор, сопряженный 5), означающего, что норма оператора 5 не превышает единицы. Неравенство 5*5 < Е имеет в данном случае вид

(Е + тА*)-1 А* + А(Е + тА)-1 > т(Е + тА*)-1А*А(Е + тА)-1

и эквивалентно операторному неравенству

А*(Е + тА) + (Е + тА*)А > тА*А,

которое, в свою очередь, эквивалентно неравенству

А + -А*А > 0. (14)

2 - V у

Поскольку А > 0 и А*А > 0, то (14) выполнено и ||$|| < 1. Поэтому

Цу'+І < ||£IIНу'II < ||у'! < ... < ||у(х,0)| = 1Ых)||.

Теорема доказана.

Перейдем теперь к исследованию устойчивости по правой части. Для этого рассмотрим схему

Ву* + Ау' = /', у(х, 0) = мо(ж), х Є оЛ. (15)

В работе [3] доказана теорема о том, что из устойчивости схемы (11) по начальным данным

в некоторой норме || • || 1 следует устойчивость схемы (15) по правой части, взятой в норме

II/' 112 = ||В-1/'|| 1. Отсюда следует

Теорема 3. Схема (15) устойчива по правой части, если выполнены условия теоремы 2, и для решения справедлива оценка

'

Иу'+ЧКЫ! + £ т ||В-1 /‘ и. (16)

к=о

Поскольку В-1 < 5 < Е, вместо (16) можно получить оценку

'

1'||<М + £т/(17)

к=о

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

3

задачи при т ^ 0, ^ ^ 0 сходится к решению дифференциальной задачи со скоростью 0(т + + Л| + Л3).

5. Результаты численного эксперимента

Экспериментальное исследование рассматриваемой в работе разностной схемы проводилось на примерах многочисленных задач анизотропной теплопроводности: линейных, квазилинейных, нелинейных, двумерных, трехмерных, а также путем сравнения с полученными авторами аналитическими решениями задач анизотропной теплопроводности и с результатами, полученными с помощью метода дробных шагов Н. Н. Яненко (МДШ), метода переменных направлений (МПН) Писмена — Рэчфорда, явной схемы. Эти исследования проводились в широком диапазоне чисел Куранта от 1 до 150 (Ки = кт/Л2, к = 4||к||/ср, Л = у/Щ + Л2 + Л2, ||к|| = шах(к^, , &£)), где — компоненты главного

тензора теплопроводности, ср — объемная теплоемкость). Теоретическое и экспериментальное обоснование метода приведено в [4].

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

В работах [5, 6] получено и исследовано аналитическое решение следующей задачи анизотропной теплопроводности в полосе —то < х < +то, 0 < у < /:

ди

ср

, д2и д2и д 2и

Кц 2 + 2 «12 Ті— -----------+ «22

дг дх2

— сю < х < +оо

дхду

0 < у < I,

и(х, 0,і) : и(х, /, і)

и1,

0,

И2,

0,

|х| < Ьь |х| > Ь1, |х| < І2, |х| > ^2,

У

У

У

У

дУ2 і > 0; о, г > о,

о, г > о,

/, г > о,

/, г > о,

м(±то, у, г) = о,

ди(±то,у, г)

дх

о,

и(х,у, 0) = 0, —то < х < +то, Решением этой задачи является функция

о < у < I, г = о.

(18)

(19)

(20)

(21)

(22)

(23)

(24)

и(х, у,г) = |

£\(у,т)(егі+ еі'Г)^т + — I ^2(у,г)(егГ+ еі'ГЬ2)^г,

(25)

где

±

Ь1 ± ту ^ х

2^? ;

Ь±

I

г

т =

^12

^22

5 =

кп срк22 ’

к2 =

ср

22

0.08

0.04

\2

:±__

0.8

0.4

У 13 4

2 __

/ 1

0.016

0.08

II

_ 4

2

/ "Г

0 0.1 0.2 0.3 0.4 0.5 г, с

0.8

0.4

4 /

3 / 1

0.08

0.04

з/ 4 2

/ 1^ / / У

/ А/ 1У2 //

О 0.1 0.2 0.3 0.4 0.5 с 0 0.1 0.2 0.3 0.4 0.5 с

Рис. 2. Зависимость погрешности е от времени для чисел Куранта

I - Ки = 1.0 (а), 10 (б), 150 (в), II - Ки = 0.25 (а), 1.0 (б).

2п

ГО

^1(ут) = <

к212 ^ П=1

. /ппул

П • 81П ^ ехр

22 п2п2 т

к212

0 < у <

5(т) = 5(т + 0), у = 0;

0, у = 1.

— 2п

В-ч

к2/2

n=1

о, у = о,

, ^(г) = £(т + 0), у = /.

• /ППУ\ n ■ sin ^ —— j exp

n2 п2т k2/2

0 < у < /,

Решение (25) автоматически удовлетворяет уравнению (18), а для удовлетворения краевым условиям необходимо использовать следующие свойства дельта-функции и функции ошибок:

t

J /(тЖт + 0)dT = /(0); erf(^) = 10

Поскольку задача (18)-(24) имеет аналитическое решение, то для оценки точности и эффективности предлагаемой конечно-разностной схемы используется критерий

e(t) = max |u(xi,yj,t) - ufc(xi,yj,t)|,

где u,uk — соответственно точное и конечно-разностное решения.

Для входных данных T = 0.6 с; / = 0.02 м; cp = 103 кДж/м3К; u1 = u2 = 1; L1 = L2 = 0.02 м; = 0.02 кВт/мК; An = 0.08 кВт/мК; ^ = 45 °; u(x, y, 0) = 0; h1 = h2 = 0.0025 м на рис. 2, I представлены сравнительные результаты зависимости критерия e(t) от времени и чисел Куранта: а) Ku =1; т = 0.01 с; к = 0.625 ■ 10-3 м2/с; б) Ku =10; т = 0.1 с; к = 0.625 ■ 10-3 м2/с; в) Ku = 150; т = 0.1 с; к = 10-3 ■ 9.375 м2/с.

Во всех трех вариантах кривая 1 соответствует предлагаемому методу полного расщепления, кривая 2 — методу дробных шагов Н. Н. Яненко, кривая 3 — методу переменных направлений Писмена — Рэчфорда. Из рисунка видно, что при малых числах Куранта (Ku < 1) погрешности перечисленных схем одного порядка. При увеличении чисел Ku на порядок (Ku ~ 10) явная схема становится неустойчивой (погрешность ее резко возрастает), остальные схемы остаются устойчивыми, хотя у схемы МПН погрешность примерно вдвое выше, чем у схем МДШ и полного расщепления. При увеличении чисел Ku на два порядка (Ku ~ 100), схема МПН становится неустойчивой (кривая 4 резко возрастает при t = 0.4 с) и только МДШ и предлагаемая схема остаются устойчивыми; при этом погрешность у МДШ выше, чем у схемы полного расщепления.

В табл. 1 и 2 приведены результаты решения двумерной начально-краевой задачи теплопроводности в квадрате двумя методами: методом полного расщепления (табл. 1) и методом переменных направлений (табл. 2) для т = 0.01 с; t = 0.02 с.

Таблица 1

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0000

0.2712

0.1477

0.1317

0.1299

0.1306

0.1401

0.2367

1.0000

1.0000 1.0000 1.0000

0.1458 0.1312 0.1298

0.0265 0.0143 0.0134

0.0143 0.0024 0.0014

0.0133 0.0014 0.0004

0.0143 0.0024 0.0014

0.0260 0.0143 0.0133

0.1400 0.1306 0.1299

1.0000 1.0000 1.0000

1.0000

0.1306

0.0143

0.0024

0.0014

0.0024

0.0143

0.1316

1.0000

1.0000

0.1405

0.0260

0.0143

0.0134

0.0143

0.0264

0.1472

1.0000

1.0000

0.2384

0.1402

0.1305

0.1298

0.1312

0.1454

0.2692

1.0000

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

П

В. Ф. Формалев, О. А. Тюкин Таблица 2

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0000

0.2435

0.1448

0.1381

0.1378

0.1379

0.1419

0.2147

1.0000

1.0000 1.0000 1.0000

0.1169 0.1066 0.1060

0.0167 0.0071 0.0065

0.0106 0.0009 0.0004

0.0102 0.0006 0.0001

0.0105 0.0009 0.0004

0.0164 0.0070 0.0065

0.1133 0.1063 0.1060

1.0000 1.0000 1.0000

1.0000

0.1063

0.0070

0.0009

0.0006

0.0009

0.0071

0.1066

1.0000

1.0000

0.1133

0.0164

0.0105

0.0102

0.0106

0.0167

0.1169

1.0000

1.0000

0.2147

0.1419

0.1379

0.1378

0.1381

0.1448

0.2435

1.0000

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

Входные данные принимали следующие значения: и|г = 1; и(ж,у, 0) = 0; ср = = 1000 кДж/м3К; кц = к22 = 0.05 кВт/мК; к12 = -0.03 кВт/мК. Угол ^ ориентации главных осей относительно осей декартовой системы координат составляет 45 °, так что распределение температуры должно быть симметричным относительно побочной диагонали квадрата, что можно использовать для контроля точности.

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

В качестве следующего теста метода полного расщепления рассматривается сравнительный анализ (по сравнению с классическими численными методами) поведения во времени абсолютной погрешностиотносительно аналитического решения следующей существенно нелинейной задачи для уравнения параболического типа:

ди

Срд£

д I к ди к12 / ди дж 1 11 дх1 + 2 1дж2

+

_д_

дж2

1 ( ди

3 \дх1

+

ди \

+к22 ^— I + 2ср£ — 10, 0 < х1 < /1, 0 < ж2 < /2, £ > 0,

дх2 )

2 2 3 2

и(х1,х2, 0) = -— х1 + -— х2, 0 < х1 < /1, 0 < ж2 < /2, £ = 0,

к11 к22

3

и(0,х2,£) = -— х2 + £2, х1 = 0, 0 < ж2 < /2, £ > 0,

к22 2

23

и(11 ,ж2,£) = -— ^ + -— х2 + £2, х1 = /1, 0 < ж2 < /2, £ > 0,

к11 1 к22 2

2

и(х1,0,£) = -— ж? + £2, 0 < х1 < /1, ж2 = 0, £ > 0,

к11 1

23

и(ж1,/2,£) = -— ж2 + -— + £2, 0 < ж1 < /1, ж2 = /2, £ > 0.

к11 1 к22 2

Для этой задачи подобрано частное аналитическое решение

и(х1, 2,£) = т~+ -т3-ж2 + £2, к11 1 к22 2

(26)

(27)

(28)

(29)

(30)

(31)

(32)

3

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

e(t) = max |yibj2 (t) - (t) |,

il ,i2

где u — точное решение (32), a y — сеточная функция, полученная в результате численного решения различными методами.

На рис. 2, II представлены результаты тестирования по критерию (33) предлагаемого метода полного расщепления (кривые 1), метода дробных шагов (кривые 2), явного метода (кривые 3) и метода переменных направлений (кривые 4) для чисел Ku = 0.25; т = 0.0025; h = 0.01; к = 0.001 (а) и Ku = 1.0; т = 0.1; h = 0.01; к = 0.001 (б).

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

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

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

[2] ЯНЕНКО Н. Н. Очерки. Статьи. Воспоминания. Сост. Н. Н. Бородина. Наука, Новосибирск, 1988.

[3] САМАРСКИЙ А. А. Теория разностных схем. Наука, M., 1983.

[4] Тюкин О. А. Метод полного расщепления численного решения параболических уравнений со смешанными производными: Дис. ... канд. физ.-мат. наук. МАИ, М., 1996.

[5] Формалев В. Ф., Москаленко А. А. Аналитическое решение трехмерной нестационарной задачи теплопроводности с тензором теплопроводности. Дифференциальные уравнения, 26, №7, 1990, 1277-1279.

[6] Формалев В. Ф., Тюкин О. А. Исследование температурных полей на основе аналитического решения двумерной задачи анизотропной теплопроводности, 32, №4, 1994, 518-523.

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