CC BY
49
Вычислительные технологии
Том 15, № 5, 2010
Экономичный полностью неявный метод численного решения параболических уравнений, содержащих смешанные производные*
E.J1. Кузнецова, В.Ф. Формален Московский авиационный институт (государственный технический университет), Россия e-mail: [email protected]
Предложен экономичный полностью неявный абсолютно устойчивый метод численного решения уравнений параболического типа, содержащих смешанные дифференциальные операторы, основанный на более глубоком расщеплении конечно-разностной задачи по сравнению с классическими методами.
Ключевые слова: неявный абсолютно устойчивый метод, численные методы, уравнения параболического типа, двумерная нестационарная задача, анизотропная теплопроводность.
Введение
Большинство известных методов численного решения многомерных задач математической физики, основанных на процедуре расщепления конечно-разностных операторов по координатным направлениям, используют явную аппроксимацию смешанных дифференциальных операторов, что при определенных условиях, накладываемых на матрицу характеристик переноса, может приводить к неустойчивости численного решения [1, 2]. Кроме того, наличие смешанных дифференциальных операторов в уравнениях существенно затрудняет процедуру улучшения точности аппроксимации, особенно для краевых условий, содержащих производные. Так, в задачах анизотропной теории теплопроводности практически невозможно построить экономичную полностью неявную конечно-разностную схему расщепления по координатным направлениям без аппроксимирования смешанных дифференциальных операторов на нижнем временном слое и использования дополнительной информации о решении.
В данной работе предложен экономичный полностью неявный конечно-разностный метод расщепления по координатным направлениям численного решения нестационарных задач математической физики, содержащих смешанные дифференциальные операторы. Метод основан на расщеплении не только по координатным направлениям, но и по отдельным дифференциальным операторам, входящим в дифференциальную задачу, а, также на использовании апостериорной информации, полученной на верхнем временном, слое, что является новым, подходом, к построению конечно-разностных схем,.
* Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00655-а) и грантов Президента РФ МК-1184.2009.8, МК-8460.2010.8, НШ-8460.2010.8.
Метод рассматривается на примере двумерной нестационарной задачи анизотропной теплопроводности с граничными условиями первого рода. Аналогичной проблеме были посвящены работы [3 - 5], Предлагаемый метод использовался при численном моделировании тепловой защиты гиперзвуковых летательных аппаратов.
1. Постановка задачи и метод численного решения
Рассмотрим следующую первую начально-краевую задачу теории теплопроводности в анизотропном прямоугольнике П = (0, /1) х (0, /2) + Г:
ди д2и д2и д2и д2и ,
ср— = Ли— + Л12т-^Г + л21т-т- + А22тг^, {х,у)еЪ1, ¿>0, (1)
дЬ дх2 дхду дудх ду2
и (х,У,Ь) |г = ф (х,У,Ь), (х,у) € Г, Ь> ° (2)
и(х,у,0) = -ф(х,у), (х,у)еП, ¿ = 0, (3)
где Ли, Л12, Л21, Л22 — компоненты тензора теплопроводности, П = П + Г, Ли > 0, Л22 > 0 Ли ■ Л22 — Л12 > 0, причем Л12, Л21 могут быть как положительными, так и отрицательными.
На расчетную пространственно-временную область накладывается конечно-разностная сетка
= {хг = гЬъ г = 0, /, ^ = ]Н2, ] = 0, J, Г = кт, к = 0,1, 2...} , (4)
где к1 = 11/1, к2 = 12/3, на которой задача (1)-(3) аппроксимируется следующей полностью неявной конечно-разностной схемой [атп = \тп/ср\ т,п = 1,2):
к+1/2 — к
_
т
_ «11 / кЦ _ кЦ кЦ \ 1 - О сц2 ( кЦ _ кЦ _ к+1 кЦ \
1 + а а\2 / к+± к+± кЦ к+± \ , ,
«Г — «Г2
а22 ( к+1 о к+1 к+1 с 1 а а21 ( к+1 к+1 к+1 к+1с
\иг,]+1 ~ + иг,]-1) ^ 2 ~ ,3 ~ 1 + )
I 1 + (Т "21 С-/+1 _ _ I ук+1 Л /-ч
Шаблоны этой схемы представлены на рисунке.
Подсхема (5) содержит три неизвестных и1к+^/2, ик+1/2, и^^2, определяемых скалярными прогонками вдоль координатных линий, параллельных оси Ох, а значения ик+1/2_ 1, ик+_1(2 при а = —1 (Л12 > 0) и значения , и^1^ при а = 1 (Л12 < 0)
уже определены на верхнем временном полуслое Ьк+1/2 = (к + 1/2) т методом прогонки вдоль координатной линии yj_1 = (] — 1) к2. В подсхеме (6) также три неизвестных:
О-
У
ч
о-
7'+ 1,7
—О-
<7=1
У
7 _ 1,7 7 7 О—
/,7-1 / + 1,7-1 а
0~
а = -1
7 + 1,7 —о—
—□
/-1,7-1 и-1 б
/'-1,7 + 1
.V'
7'- 1,7
Г
□-о 7,7 + 1 7,7 + 1 о-□
/ + 1,7 +1
<7 = — 1
о-
77
/,7- 1о
7'+ 1,7
г
□ — Узлы с известными значениями О — Узлы с искомыми значениями
Шаблоны схемы глубокого расщепления: а, б— подсхема (5); в, г — подсхема (6)
мк?+15 определяемых в скалярных прогонках вдоль оси Оу а значения и^-^-,
+1 при а = — 1 (Л21 > 0) и значения ик—ик_^ при а = 1 (Л21 > 0) уже определены на верхнем временном слое Ьк+1 = (к + 1)т методом прогонки вдоль координатных линий = (г + 1)к1 при а = — 1 и = (г — 1)^ при а = 1.
Коэффициент а = —1 в случае, когда тепловой поток имеет направление "третий—
а=1
Если наложить друг на друга шаблоны б и г и а и в (см. рисунок), то получим известную центрально-симметричную схему Самарского [2]. В ней в отличие от предлагаемой схемы смешанные дифференциальные операторы аппроксимируются на нижнем временном слое (явно).
Схема (о), (6) является экономичной, поскольку использует только скалярные прогонки вдоль координатных направлений. Покажем, что данная схема имеет диагональное преобладание, для чего запишем ее в виде следующих систем линейных алгебраических уравнений с трехдиагональными матрицами, например, при а = — 1 :
к+1/2 (<Н1_ _ 0'У2 \ _ к+1/2 (2а 11 , ^ _ а12 \ к+1 (^п | _
^ Щ ~ ыы) ~ \Щ + Г " ыы) 1 "
к
и%3 а> 12 ( к+1/2 к+1/2
= + г = 1,1-1 ^ = 1,7-1, (7)
,,А;+1 („,к+1 (2а22 _ °21 \ к+1 (^22 °21 _
~г , , Т И;
\ щ ) \щ т нгн2) " \ щ нгн2) ~ к+1/2
В выражениях (7), (8) выполняются соответственно неравенства
2«ц 1 «12
Т /1/2
2а22 1 «21
Т /1/2
>
>
«11 «12
«22
«11
н\
/2
«22 «21 /2 н1н2
Схема (5), (6) обладает частичной аппроксимацией па каждом дробном временном шаге и полной аппроксимацией па целом временном шаге с порядком О (т + (/ + /2)2) . Действительно, запишем схему (5), (6) в следующей векторно-операторной форме при а = — 1 (при а = 1 рассуждения аналогичны):
и
к+1/2 — ик
Т
ик+1 — ик+1/2
Т
Лпик+1/2 + Л12ик+1/2, = Л22ик+1 +Л21ик+1,
(9) (10)
где
А2 2 Д2 ,
д _ „ _А — п 111112 А — п _-2. А — п 112111
Ли — а11~Го~) А12 — «12 ~—;-, ^22 — «22 ' А21 ~~ Й21
,2 > —12 -12 , , , —22 -22 , 2 ' —21 -21 , , /¿1 ^1^2 ^2 ^2^1
Исключая из (9), (10) векторно-матричные операторы на промежуточном временном
слое, получим эквивалентную схему
ик+1 — ик
Т
= (Лц + Л22 + Л12 + Л21) ик+1 + О (т)
(Н)
откуда следует первый порядок аппроксимации по времени. Далее имеем
д2и
Л11ик+1 = «п
Л22и
к+1
^2 |Г + О (Л?) ,
д2и
«22
дУ
к+1
2 I«.?
О(Ы
(12)
Рассмотрим операторы (Л12 + Л21) ик+1. Разлагая па точном решении в ряды Тейлора сеточную функцию на шаблоне в окрестности центрального узла получим (в случае «12 = «21)
(Л12 + Л21) ик+1 =
П1П2 а21 ( к+1
бик+1 — ик+1 — и к+1 + и к+1 \ + 1« ,?-1 и«-1 , з + и«-1 , з-и +
«12 /1/2
аг,3+1^ агЗ
1 _ —и — М _ _Ё!_М О (К
ду 2 <9у2 2 <9у3 6 ^
1 д ^ д2 /2 д3 ^ ^ 4 1 - — Ы + —- —т1 + О
ик+1-
ик+1—
дх дх2 2 дх3 6
ик+1+
д , д \ ^ ^ д . Т-Л-1 + Т-Л-2 + - Т-Л-1 + Т-Л-2 I -дх ду у 2 удж ду
<+' > +
«21
5 5 \ 1 ( д д
1 + ( —кг + —1г2 + - —кг + — Ь2 дх ду ) 2 \дх ду
+
д2 й2 д3 й3
д д д 4 ^ дх^1 дх2 2 ^ дх3 6 ^ ^
ик+1—
ик+1—
1 д , д2 ^2 д3 ^2 ^,7 1 + —Ьъ + —■+ —^ + О {¡г,
ду ду2 2 ду3 6
«12 «21
{"ху№ + О (^4) + О (й4) + О ((^1 + й2)4) }
4^ ч к+1
г.
+ О ((^1 + + О (й4) + О (й*) } = «12"ху |к+1 + «21иух|к+1 + О ((^1 + Л^) .
к+1 г.
(13)
Из (11)—(13) следует аппроксимация схемы (5), (6) с порядком О (г + + Л,2) ) .
Для исследования устойчивости схемы (5), (6) энергетическим методом умножим екалярно эквивалентную схему (11) на вектор и = ("¿+1 — ик) /т, используя тождество
и
к+1
ик+1 + ик т
УСЛОВИЯ ЭЛЛИПТИЧНОСТИ
Лц > 0, Л22 > 0, Л11Л22 — Л^2 > 0,
(14)
(15)
а также положительность и самосопряженность конечно-разностного оператора А = — Л = — (Лц + Л12 + Л21 + Л22) [2]
(Аи,г>) = (и, Аг>), А > 0.
(16)
В результате получим
(щ,щ)=(Аик+1,щ) = (л ^ откуда следуют равенства
ик+1 + ик т
+ —Щ ) , щ
Л
ик+1 + ик
т
) + - (Ащ,щ)
((е- 1л) «,.«,) = (А ( ((в-1л)И,,и,) = (л(
ик+1 + ик
,и(
ик+1 + ик\ ик+1 — ик
— (Л 2т
ик+1,ик+1
--(Аи\ик)--(А
2т
ик+1,ик
+ _1 (ли*,«*+1).
2т
(17)
2
2
2
2
2
2
В силу (16) и коммутативности скалярного произведения получаем энергетическое тождество
(Аик+1,ик+1) +2т + = (Аик,ик) , (18)
Т
в котором оператор Е + —А > 0, откуда следует принцип максимума
бАик+1,ик+1) < ^ик,ик) (19)
ИЛИ
1|ик+11| а < 1|ик || а < ... < ||и° || а = 11Ф (*,У)Н , (20)
где ф (х,у) — начальное распределение функции в условии (3), а ||и||А _ норма в энергетическом пространстве На сеточных функций и2 со скалярным произведением = = (Аи2,^2) ДЛЯ € На И нормой
|и2|А = (Аи2,и2)1/2 . (21)
Поскольку не принималось никаких ограничений на сеточные характеристики т, , /2 и коэффициенты переноса, кроме условий эллиптичности (15), которые для тензоров второго ранга выполняются всегда, можно заключить, что схема (5), (6) абсолютно устойчива.
Таким образом, схема (5), (6) глубокого расщепления обладает следующими достоинствами при решении нестационарных задач, содержащих смешанные производные:
1) экономичность;
2) неявная аппроксимация всех дифференциальных операторов, включая смешанные и, как следствие, абсолютная устойчивость;
3) применимость к задачам любой размерности по пространственным переменным;
4) отсутствие ограничений на величину коэффициентов «11, «12, «21, «22 за исключением естественных ограничений вида «11 > 0, «22 > 0, «11«22 — «22 > 0.
Как и большинство схем расщепления, схема (5), (6) имеет первый порядок аппроксимации по времени и второй по пространственным переменным.
Поскольку на возрастающем по времени решении неявная конечно-разностная схема аппроксимирует точное решение сверху, а явная снизу, а на убывающем — наоборот: явная схема — сверху, а неявная — снизу [4], то для увеличения порядка аппроксимации по времени можно использовать неявно-явную схему с весом в (0 ^ в ^ 1) при неявной части и 1 — в при явной. При этом следует ожидать условной устойчивости при 0 ^ в ^ 0.5 и безусловной при 0.5 ^ в ^ 1, Второй порядок аппроксимации по времени достигается при в = 0.5, однако при этом значительно снижается запас устойчивости, заложенный в неявную аппроксимацию.
Для предложенной конечно-разностной схемы (5), (6) схема с весом в = 1 — в = 0.5 имеет вид
к
—г-=КАи+^+Тль)"'+1/2+КЛи++1,1
иШ -"" 1 'л* + + 1±^л+У+Ч^л22 + + 1±£л*
2 V 2 2 21 2 21 / 2 V 2 2
2. Сравнительный анализ схемы глубокого расщепления
При анализе эффективности предлагаемой схемы, как и других конечно-разностных схем, целесообразно использовать максимальную абсолютную погрешность
е (£) = тах (£) — МуТ(£)|
(22)
где Му (£) — сеточная функция, МуТ(£) — точное решение тестовой задачи, В качестве такой задачи рассмотрим задачу для квазилинейного параболического уравнения
ср
дм д д£ дх
<9и Л12 / ди 11 дх 2 уду
+
д_
дУ
Л12 / <9гЛ ди
3 \ дх) 22 ду
+ 2ср£ — 10,
0 < х < /1, 0 < у < /2, £ > 0,
(23)
с начальным условием
, . 2х 3 2 м(ж,у,0) = — + —у , А11 Л22
и краевыми условиями первого рода
3
0 < х < /1, 0 < у < /2, £ = 0,
п (0, ¿) = — уг + 0 <у<12, t> 0, А22 23
и{1иу,г) = —121 + 1— у2 + Л о <у<к, ¿>о,
А11 А22 2
и(х, 0, = -—ж2 + ¿2, 0 < х < ¿>0, А11
м (х,/2,£)
2
х
А11 А22 Задача (23)-(25) имеет частное решение
+ —0 < ж < /ь ¿>0.
(24)
(25)
и (х, у, ¿) = -—ж2 А11
А
22
(26)
На этом решении протестируем следующие конечно-разностные схемы:
1) предлагаемую схему метода глубокого расщепления (МП'):
2) схему метода дробных шагов Яненко;
3) схему метода переменных направлений Писмена — Рэчфорда;
4) центрально-симметричную схему Самарского;
5) классическую явную схему,
приняв нижеприведенные входные данные. Компоненты тензора теплопроводности:
А11 = А^ осе2 ф + Ап ■ вт2 ф, А22 = А^ вт2 ф + Ап ■ осе2 ф, А12 = (А^ — Ап) осе ф ■ вт ф,
(27)
где А^, Ап — главные коэффициенты теплопроводности, ф — угол ориентации главной оси тензора теплопроводности относительно оси Ох декартовой системы координат,
2
1х = 12 = 0.08 м, конечное время ¿к = 0.6 с, теплоемкость с = 0.4 кДж/(кг-К), р = 2500 кг/м3, ф = 45°, кх = к2 = к.
Варьируя главные коэффициенты теплопроводности А^, Ап и шаг и тик пространственно-временной сетки с помощью перечисленных методов были проведены расчеты при различных значениях чисел г, где
ат 4 ЧЛЧ мхи м \ ^
г = —, ||Л|| = тах(Л?,Лч).
В табл. 1, 2 приведены результаты для фиксированных момента времени и значения переменной ж, по которым можно судить об абсолютных погрешностях различных методов, Из табл. 1 видно, что при г < 1, а именно в области устойчивости классической явной схемы, и т = 0.025 с погрешности перечисленных схем имеют практически один порядок, С увеличением г до значения ~ 2.0 (см, табл. 2) резкий скачок абсолютной
Таблица 1. Значения температур в сечении х = 0.05 м в момент времени Ь = 0.6 с при т = 0.025 с, Н = 0.01 м, а = 0.001 м2/с, г = 0.25
у, м
Схема 0.01 0.02 0.03 0.04 0.05 0.06 0.07
Точное решение 0.4583 0.5033 0.5783 0.6833 0.8183 0.9833 1.1783
Схема глубокого
расщепления 0.4725 0.5183 0.5934 0.6983 0.8333 0.9950 1.1975
Схема метода
дробных шагов 0.4713 0.5205 0.5913 0.6975 0.8314 0.9936 1.1955
Схема метода
переменных
направлений 0.4654 0.5109 0.5858 0.6918 0.8260 0.9913 1.1853
Центрально-сим-
метричная схема 0.4725 0.5183 0.5933 0.6993 0.8335 0.9984 1.1923
Явная схема 0.4825 0.5284 0.5932 0.6963 0.8427 0.9987 1.1923
Таблица 2. Значения температур в сечении х = 0.05 м в момент времени Ь = 0.6 с при т = 0.2 с, Н = 0.01 м, а = 0.001 м2/с, г = 2.0
у, м
Схема 0.01 0.02 0.03 0.04 0.05 0.06 0.07
Точное решение 0.4583 0.5033 0.5783 0.6833 0.8183 0.9833 1.1783
Схема глубокого
расщепления 0.4961 0.5557 0.6460 0.7616 0.8945 1.0426 1.2100
Схема метода
дробных шагов 0.4956 0.5617 0.5849 0.7448 1.3592 -0.300 2.9474
Схема метода
переменных
направлений 0.4999 0.8148 3.4225 7.8926 98.754 89.546 -72.34
Центрально-сим-
метричная схема 0.5062 0.5634 0.6134 0.7240 1.4283 -0.5634 3.7213
Явная схема 0.4877 0.7158 0.2567 0.2894 2.4564 -13.84 -11.34
погрешности e(i) явной схемы, схемы метода переменных направлений, центрально-симметричной и схемы метода дробных шагов иллюстрирует их неустойчивость, в то время как схема глубокого расщепления остается устойчивой, хотя его абсолютная погрешность возрастает по сравнению со случаем r = 0.25 (см, табл. 1) в соответствии с порядком точности O (т + (hi + h2)2).
Возникновение неустойчивости в рассматриваемых схемах уже при умеренных чнс-r
емешанных дифференциальных операторов на нижних временных слоях (явно).
Следует отметить, что нелинейности в коэффициентах задачи (23)-(25) во всех конечно-разностных схемах, кроме явной, учитывались в итерационном процессе Ньютона, при этом всюду достаточно было не более одной итерации, не считая нулевой.
Выводы
Предложен новый экономичный абсолютно устойчивый метод глубокого расщепления численного решения задач для уравнений параболического типа, содержащих смешанные дифференциальные операторы, основанный на более глубоком, чем в существующих методах, расщеплении смешанных производных и их неявной аппроксимации.
Исследования метода глубокого расщепления на примерах решения задач анизотропной теплопроводности показали зависимость вида конечно-разностной схемы, а следовательно и устойчивости, от знаков коэффициентов теплопроводности при смешанных производных. Показано, что абсолютная устойчивость схемы метода глубокого расщепления имеет место в случаях, когда смешанный конечно-разностный оператор ориентирован по вектору теплового потока, В связи с этим разработан специальный алгоритм, сохраняющий абсолютную устойчивость метода.
Представленный метод является достаточно общим и может быть применен не только к задачам анизотропной теплопроводности, но и к другим задачам математической физики.
Список литературы
[1] Яненко H.H. Метод дробных шагов решения задач математической физики. Новосибирск: Наука, 1967. 196 с.
[2] Самарский A.A. Теория разностных схем. М.: Наука, 1983. 616 с.
[3] Формалев В.Ф. Метод переменных направлений с экстраполяцией по времени для параболических задач со смешанными производными // Вычисл. технологии. 1996. Т. 1, № 2. С. 99-103.
[4] Формалев В.Ф., Ревизников Д.Л. Численные методы. М.: Физматлит, 2004. 400 с.
[5] Формалев В.Ф., Тюкин O.A. Неявный метод дробных шагов с расщеплением смешанных дифференциальных операторов // Вычисл. технологии. 1998. Т. 3, № 6. С. 82-91.
Поступила в редакцию 28 марта 2010 г., с доработки — 28 мая 2010 г.
