close

Вход

Забыли?

вход по аккаунту

?

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

код для вставкиСкачать
Вычислительные технологии
Том 15, ќ 5, 2010
Экономичный полностью неявный метод
численного решения параболических уравнений,
?
содержащих смешанные производные
Е. Л. Кузнецова, В. Ф. Формалев
Московский авиационный институт
(государственный технический университет), оссия
e-mail:
lareynamail.ru
Предложен экономичный полностью неявный абсолютно устойчивый метод численного решения уравнений параболического типа, содержащих смешанные диеренциальные операторы, основанный на более глубоком расщеплении конечноразностной задачи по сравнению с классическими методами.
Ключевые слова : неявный абсолютно устойчивый метод, численные методы,
уравнения параболического типа, двумерная нестационарная задача, анизотропная теплопроводность.
Введение
Большинство известных методов численного решения многомерных задач математической изики, основанных на процедуре расщепления конечно-разностных операторов
по координатным направлениям, используют явную аппроксимацию смешанных диеренциальных операторов, что при определенных условиях, накладываемых на матрицу характеристик переноса, может приводить к неустойчивости численного решения
[1, 2?. Кроме того, наличие смешанных диеренциальных операторов в уравнениях
существенно затрудняет процедуру улучшения точности аппроксимации, особенно для
краевых условий, содержащих производные. Так, в задачах анизотропной теории теплопроводности практически невозможно построить экономичную полностью неявную
конечно-разностную схему расщепления по координатным направлениям без аппроксимирования смешанных диеренциальных операторов на нижнем временном слое и
использования дополнительной инормации о решении.
В данной работе предложен экономичный полностью неявный конечно-разностный
метод расщепления по координатным направлениям численного решения нестационарных задач математической изики, содержащих смешанные диеренциальные операторы. Метод основан на расщеплении не только по координатным направлениям,
но и по отдельным диеренциальным операторам, входящим в диеренциальную
задачу,
а также на использовании апостериорной инормации, полученной на верх-
нем временном слое, что является новым подходом к построению конечно-разностных
схем.
?
абота выполнена при инансовой поддержке ФФИ (грант ќ 11-01-00655-а) и грантов Прези-
дента Ф МК-1184.2009.8, МК-8460.2010.8, НШ-8460.2010.8.
72
73
Экономичный полностью неявный метод численного решения...
Метод рассматривается на примере двумерной нестационарной задачи анизотропной теплопроводности с граничными условиями первого рода. Аналогичной проблеме
были посвящены работы [3 5?. Предлагаемый метод использовался при численном моделировании тепловой защиты гиперзвуковых летательных аппаратов.
1. Постановка задачи и метод численного решения
ассмотрим следующую первую начально-краевую задачу теории теплопроводности в
анизотропном прямоугольнике
c?
? = (0, l1 ) Ч (0, l2 ) + ?:
?2u
?2u
?2u
?2u
?u
= ?11 2 + ?12
+ ?21
+ ?22 2 ,
?t
?x
?x?y
?y?x
?y
u (x, y, t) |? = ? (x, y, t) ,
u (x, y, 0) = ? (x, y) ,
где
?22
(x, y) ? ?,
(x, y) ? ?,
?11 , ?12 , ?21 , ?22 компоненты тензора
> 0, ?11 · ?22 ? ?212 > 0, причем ?12 , ?21
(x, y) ? ?,
t > 0,
t > 0,
(2)
t = 0,
теплопроводности,
(1)
(3)
? = ? + ?, ?11 > 0,
могут быть как положительными, так и
отрицательными.
На расчетную пространственно-временную область накладывается конечно-разностная сетка
?h? = {xi = ih1 , i = 0, I, yj = jh2 , j = 0, J, tk = k?, k = 0, 1, 2...} ,
где
h1 = l1 /I, h2 = l2 /J,
(4)
на которой задача (1)(3) аппроксимируется следующей пол-
ностью неявной конечно-разностной схемой
k+1/2
uij
amn = ?mn /c?; m, n = 1, 2
? ukij
:
=
?
1?? a a11 k+ 21
k+ 1
k+ 21
k+ 1
k+ 21
k+ 21
k+12
12
? 2uij 2 + ui?1,j
+
? ui?1,j
+ ui?1,j?1
+
= 2 ui+1,j
uij 2 ? ui,j?1
h1
2 h1 h2
1 + ? a12 k+ 21
k+ 21
k+ 1
k+ 21
+
? uij 2 + ui,j?1
,
ui+1,j ? ui+1,j?1
2 h1 h2
(5)
k+1/2
=
a22
h22
k+1
uij
? uij
=
?
1 ? ? a21
k+1
k+1
k+1
k+1
k+1
uk+1
+ ui,j?1
+
uk+1
+
i,j+1 ? 2uij
i+1,j+1 ? ui+1,j ? ui,j+1 + uij
2 h1 h2
+
1 + ? a21
k+1
k+1
uk+1
? uk+1
i,j+1 ? uij
i?1,j+1 + ui?1,j .
2 h1 h2
(6)
Шаблоны этой схемы представлены на рисунке.
k+1/2
k+1/2
k+1/2
Подсхема (5) содержит три неизвестных ui?1,j , uij
, ui+1,j , определяемых скалярными прогонками вдоль координатных линий, параллельных оси Ox, а значения
k+1/2
k+1/2
k+1/2
k+1/2
ui?1,j?1, ui,j?1 при ? = ?1 (?12 > 0) и значения ui,j?1 , ui+1,j?1 при ? = 1 (?12 < 0)
k+1/2
уже определены на верхнем временном полуслое t
= (k + 1/2) ? методом прогонки вдоль координатной линии
yj?1 = (j ? 1) h2 .
В подсхеме (6) также три неизвестных:
74
Е. Л. Кузнецова, В. Ф. Формалев
Шаблоны схемы глубокого расщепления: а,
б подсхема (5); в, г подсхема (6)
k+1
k+1
k+1
uk+1
i,j?1 , uij , ui,j+1, определяемых в скалярных прогонках вдоль оси Oy , а значения ui+1,j ,
k+1
k+1
uk+1
i+1,j+1 при ? = ?1 (?21 > 0) и значения ui?1,j+1 , ui?1,j при ? = 1 (?21 > 0) уже опредеk+1
лены на верхнем временном слое t
= (k + 1)? методом прогонки вдоль координатных
линий xi+1 = (i + 1)h1 при ? = ?1 и xi?1 = (i ? 1)h1 при ? = 1.
Коэициент ? = ?1 в случае, когда тепловой поток имеет направление третий
первый квадранты, и наоборот ? = 1 при направлении четвертыйвторой квадранты.
Если наложить друг на друга шаблоны
б и г и а и в (см. рисунок), то получим из-
вестную центрально-симметричную схему Самарского [2?. В ней в отличие от предлагаемой схемы смешанные диеренциальные операторы аппроксимируются на нижнем
временном слое (явно).
Схема (5), (6) является экономичной, поскольку использует только скалярные прогонки вдоль координатных направлений. Покажем, что данная схема имеет диагональное преобладание, для чего запишем ее в виде следующих систем линейных алгебраических уравнений с трехдиагональными матрицами, например, при
k+1/2
ui?1,j
a11
a12
?
2
h1
h1 h2
?
k+1/2
uij
2a11 1
a12
+ ?
2
h1
?
h1 h2
+
uk+1
i+1,j
? = ?1 :
a11
h21
=
ukij
a12 k+1/2
k+1/2
+
ui,j?1 + ui?1,j?1 , i = 1, I ? 1 j = 1, J ? 1,
?
h1 h2
a22
2a22 1
a21
a22
a21
k+1
k+1
k+1
ui,j?1
? uij
+ ?
+ ui,j+1
?
=
h22
h22
?
h1 h2
h22
h1 h2
=?
k+1/2
uij
=?
?
+
a21
k+1
uk+1
i+1,j ? ui+1,j+1 ,
h1 h2
i = 1, I ? 1, j = 1, J ? 1.
(7)
(8)
75
Экономичный полностью неявный метод численного решения...
В выражениях (7), (8) выполняются соответственно неравенства
2a11 1
a11
a11 a
a
12
12
h2 + ? ? h1 h2 > h2 ? h1 h2 + h2 ,
1
1
1
2a22 1
a21 a22 a22
a21 h2 + ? ? h1 h2 > h2 + h2 ? h1 h2 .
2
2
2
Схема (5), (6) обладает частичной аппроксимацией на каждом дробном временном
O ? + (h1 + h2 )2 .
шаге и полной аппроксимацией на целом временном шаге с порядком
Действительно, запишем схему (5), (6) в следующей векторно-операторной орме при
? = ?1
(при
?=1
рассуждения аналогичны):
uk+1/2 ? uk
= ?11 uk+1/2 + ?12 uk+1/2 ,
?
uk+1 ? uk+1/2
= ?22 uk+1 + ?21 uk+1,
?
где
?11 = a11
?2h2
1
h21
,
?12
?2
= a12 h1 h2 ,
h1 h2
?22 = a22
?2h2
2
h22
,
?21 = a21
(9)
(10)
?2h2 h1
.
h2 h1
Исключая из (9), (10) векторно-матричные операторы на промежуточном временном
слое, получим эквивалентную схему
uk+1 ? uk
= (?11 + ?22 + ?12 + ?21 ) uk+1 + O (? ) ,
?
(11)
откуда следует первый порядок аппроксимации по времени. Далее имеем
? 2 u k+1
+ O h21 ,
ij
2
?x
?2u = a22 2 k+1
+ O h22 .
ij
?y
?11 uk+1 = a11
?22 uk+1
(?12 + ?21 ) uk+1 . азлагая на точном решении в ряды Тейлоункцию на шаблоне в окрестности центрального узла ij, получим (в слу-
ассмотрим операторы
ра сеточную
чае
(12)
a12 = a21 )
a12
k+1
k+1
uk+1
? uk+1
ij
i,j?1 ? ui?1,j + ui?1,j?1 +
h1 h2
a21
k+1
k+1
k+1
+
uk+1
=
i+1,j+1 ? ui+1,j ? ui,j+1 + uij
h1 h2
k+1
a12
?
? 2 h22
? 3 h32
k+1
4
=
uij ? 1 ?
h2 + 2
? 3
+ O h2 uij ?
h1 h2
?y
?y 2
?y 6
k+1
?
? 2 h21
? 3 h31
4
? 1?
h1 + 2
? 3
+ O h1 uij +
?x
?x 2
?x 6
"
2
?
?
1 ?
?
+ 1?
h1 +
h2 +
h1 +
h2 ?
?x
?y
2 ?x
?y
(?12 + ?21 ) uk+1 =
76
Е. Л. Кузнецова, В. Ф. Формалев
)
#
3
1 ?
?
4
?
h1 +
h2 + O (h1 + h2 )
uk+1
+
ij
6 ?x
?y
("
2
?
?
1 ?
?
a21
+
1+
h1 +
h2 +
h1 +
h2 +
h1 h2
?x
?y
2 ?x
?y
#
3
?
1 ?
4
+
h1 +
h2 + O (h1 + h2 )
uk+1
ij ?
6 ?x
?y
k+1
? 2 h21
? 3 h31
?
4
? 1+
h1 + 2
+ 3
+ O h1 uij ?
?x
?x 2
?x 6
k+1
?
? 2 h22
? 3 h32
4
k+1
? 1+
h2 + 2
+ 3
+ O h2 uij + uij
=
?y
?y 2
?y 6
k+1
a12 =
uxy h1 h2 + O h42 + O h41 + O (h1 + h2 )4 ij +
h1 h2
k+1
a21 +
uyx h1 h2 + O (h1 + h2 )4 + O h41 + O h42 ij =
h1 h2
k+1
= a12 uxy k+1 + a21 uyx + O (h1 + h2 )2 .
ij
(13)
ij
Из (11)(13) следует аппроксимация схемы (5), (6) с порядком
O ? + (h1 + h2 )2 .
Для исследования устойчивости схемы (5), (6) энергетическим методом умножим
k+1
скалярно эквивалентную схему (11) на вектор ut = u
? uk /?, используя тождество
uk+1 =
uk+1 + uk ?
+ ut ,
2
2
(14)
условия эллиптичности
?11 > 0,
а
также
положительность
и
?22 > 0,
?11 ?22 ? ?212 > 0,
самосопряженность
A = ?? = ? (?11 + ?12 + ?21 + ?22 )
конечно-разностного
(15)
оператора
[2?
(Au, v) = (u, Av) ,
A > 0.
(16)
В результате получим
(ut , ut ) = ?u
k+1
k+1
k+1
u
+ uk ?
u
+ uk
?
, ut = ?
+ ut , ut = ?
, ut + (?ut , ut ) ,
2
2
2
2
откуда следуют равенства
uk+1 + uk ? E ? ? ut , ut = ?
, ut ,
2
2
uk+1 + uk uk+1 ? uk ? E ? ? ut , ut = ?
,
,
2
2
?
? E ? ? ut , ut =
2
1
1
1
1
=
?uk+1, uk+1 ?
?uk , uk ?
?uk+1, uk +
?uk , uk+1 .
2?
2?
2?
2?
(17)
77
Экономичный полностью неявный метод численного решения...
В силу (16) и коммутативности скалярного произведения получаем энергетическое
тождество
Au
в котором оператор
k+1
,u
k+1
? + 2? E + A ut , ut = Auk , uk ,
2
?
E + A > 0,
2
(18)
откуда следует принцип максимума
Auk+1 , uk+1 ? Auk , uk
или
где
? (x, y)
(19)
k+1
u ? uk ? ... ? u0 = k? (x, y)k ,
A
A
A
начальное распределение ункции в условии (3), а
HA сеточных ункций uh
uh , vh ? HA и нормой
энергетическом пространстве
(uh , vh )A = = (Auh , vh )
для
(20)
kukA
норма в
со скалярным произведением
kuh kA = (Auh , uh )1/2 .
(21)
Поскольку не принималось никаких ограничений на сеточные характеристики
?, h1 ,
h2 и коэициенты переноса, кроме условий эллиптичности (15), которые для тензоров
второго ранга выполняются всегда, можно заключить, что схема (5), (6) абсолютно
устойчива.
Таким образом, схема (5), (6) глубокого расщепления обладает следующими достоинствами при решении нестационарных задач, содержащих смешанные производные:
1) экономичность;
2) неявная аппроксимация всех диеренциальных операторов, включая смешанные и, как следствие, абсолютная устойчивость;
3) применимость к задачам любой размерности по пространственным переменным;
a11 , a12 , a21 , a22
a11 > 0, a22 > 0, a11 a22 ? a212 > 0.
4) отсутствие ограничений на величину коэициентов
чением естественных ограничений вида
за исклю-
Как и большинство схем расщепления, схема (5), (6) имеет первый порядок аппроксимации по времени и второй по пространственным переменным.
Поскольку на возрастающем по времени решении неявная конечно-разностная схема
аппроксимирует точное решение сверху, а явная снизу, а на убывающем наоборот: явная схема сверху, а неявная снизу [4?, то для увеличения порядка аппроксимации по
? (0 6 ? 6 1) при неявной части и 1?? при явной. При этом следует ожидать условной устойчивости при 0 6 ? 6 0.5
и безусловной при 0.5 6 ? 6 1. Второй порядок аппроксимации по времени достигается
при ? = 0.5, однако при этом значительно снижается запас устойчивости, заложенный
времени можно использовать неявно-явную схему с весом
в неявную аппроксимацию.
Для предложенной конечно-разностной схемы (5), (6) схема с весом
? = 1 ? ? = 0.5
имеет вид
uk+1/2 ? uk
1
1 ? ? ? 1 + ? + k+1/2 1
1?? ?
1+? + k
=
?11 +
?12 +
?12 u
+
?11 +
?12 +
?12 u ,
?
2
2
2
2
2
2
uk+1 ? uk
1
1?? ?
1 + ? + k+1 1
1 ? ? ? 1 + ? + k+1/2
=
?22 +
?21 +
?21 u +
?22 +
?21 +
?21 u
.
?
2
2
2
2
2
2
78
Е. Л. Кузнецова, В. Ф. Формалев
2. Сравнительный анализ схемы глубокого расщепления
При анализе эективности предлагаемой схемы, как и других конечно-разностных
схем, целесообразно использовать максимальную абсолютную погрешность
? (t) = max |uij (t) ? uijT (t)| ,
(22)
ij
где
uij (t)
сеточная ункция,
uijT (t)
точное решение тестовой задачи. В качестве
такой задачи рассмотрим задачу для квазилинейного параболического уравнения
"
"
#
2 #
3
?
?u ?12 ?u
? ?12 ?u
?u
?u
c?
=
?11
+
+
+ ?22
+ 2c?t ? 10,
?t
?x
?x
2
?y
?y 3
?x
?y
0 < x < l1 ,
0 < y < l2 ,
t > 0,
2x2
3 2
+
y ,
?11
?22
0 ? x ? l1 ,
0 ? y ? l2 ,
(23)
с начальным условием
u (x, y, 0) =
t = 0,
(24)
и краевыми условиями первого рода
3 2
y + t2 , 0 ? y ? l2 , t > 0,
?22
2 2
3 2
l1 +
y + t2 , 0 ? y ? l2 , t > 0,
u (l1 , y, t) =
?11
?22
2 2
x + t2 , 0 ? x ? l1 , t > 0,
u (x, 0, t) =
?11
2 2
3 2
u (x, l2 , t) =
x +
l + t2 , 0 < x < l1 , t > 0.
?11
?22 2
u (0, y, t) =
(25)
Задача (23)(25) имеет частное решение
u (x, y, t) =
2 2
3 2
x +
y + t2 .
?11
?22
(26)
На этом решении протестируем следующие конечно-разностные схемы:
1) предлагаемую схему метода глубокого расщепления (М);
2) схему метода дробных шагов Яненко;
3) схему метода переменных направлений Писмена эчорда;
4) центрально-симметричную схему Самарского;
5) классическую явную схему,
приняв нижеприведенные входные данные. Компоненты тензора теплопроводности:
?11 = ?? cos2 ? + ?? · sin2 ?,
?22 = ?? sin2 ? + ?? · cos2 ?,
?12 = (?? ? ?? ) cos ? · sin ?,
где
?? , ??
главные коэициенты теплопроводности,
оси тензора теплопроводности
O?
относительно оси
(27)
?
угол ориентации главной
Ox декартовой
системы координат,
79
Экономичный полностью неявный метод численного решения...
l1 = l2 = 0.08 м, конечное время tk = 0.6
2500 кг/м3 , ? = 45? , h1 = h2 = h.
с, теплоемкость
Варьируя главные коэициенты теплопроводности
c = 0.4
?? , ??
кДж/(кг·К),
и шаги
?
и
h
? =
простран-
ственно-временной сетки с помощью перечисленных методов были проведены расчеты
при различных значениях чисел
r=
a?
,
h2
r,
где
a=
4 k?k
,
c?
k?k = max (?? , ?? ) .
В табл. 1, 2 приведены результаты для иксированных момента времени и значения
переменной
x,
по которым можно судить об абсолютных погрешностях различных ме-
тодов. Из табл. 1 видно, что при
явной схемы, и
? = 0.025
порядок. С увеличением
r < 1,
а именно в области устойчивости классической
погрешности перечисленных схем имеют практически один
r до значения
? 2.0
(см. табл. 2) резкий скачок абсолютной
Т а б л и ц а 1. Значения температур в сечении x = 0.05 м в момент времени t = 0.6 с при
? = 0.025 с, h = 0.01 м, a = 0.001 м2 /с, r = 0.25
Схема
Точное решение
Схема глубокого
расщепления
Схема метода
дробных шагов
Схема метода
переменных
направлений
Центрально-симметричная схема
Явная схема
0.01
0.4583
0.02
0.5033
0.03
0.5783
y, м
0.04
0.6833
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.4825
0.5183
0.5284
0.5933
0.5932
0.6993
0.6963
0.8335
0.8427
0.9984
0.9987
1.1923
1.1923
0.05
0.8183
0.06
0.9833
0.07
1.1783
Т а б л и ц а 2. Значения температур в сечении x = 0.05 м в момент времени t = 0.6 с при
? = 0.2 с, h = 0.01 м, a = 0.001 м2 /с, r = 2.0
Схема
Точное решение
Схема глубокого
расщепления
Схема метода
дробных шагов
Схема метода
переменных
направлений
Центрально-симметричная схема
Явная схема
0.01
0.4583
0.02
0.5033
0.03
0.5783
y, м
0.04
0.6833
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.4877
0.5634
0.7158
0.6134
0.2567
0.7240
0.2894
1.4283
2.4564
?0.5634
?13.84
3.7213
?11.34
0.05
0.8183
0.06
0.9833
0.07
1.1783
80
Е. Л. Кузнецова, В. Ф. Формалев
погрешности
?(t)
явной схемы, схемы метода переменных направлений, центрально-
симметричной и схемы метода дробных шагов иллюстрирует их неустойчивость, в то
время как схема глубокого расщепления остается устойчивой, хотя его абсолютная погрешность возрастает по сравнению со случаем
2
порядком точности O ? + (h1 + h2 ) .
r = 0.25
(см. табл. 1) в соответствии с
Возникновение неустойчивости в рассматриваемых схемах уже при умеренных чис-
лах
r
является следствием использования в них конечно-разностной аппроксимации
смешанных диеренциальных операторов на нижних временных слоях (явно).
Следует отметить, что нелинейности в коэициентах задачи (23)(25) во всех
конечно-разностных схемах, кроме явной, учитывались в итерационном процессе Ньютона, при этом всюду достаточно было не более одной итерации, не считая нулевой.
Выводы
Предложен новый экономичный абсолютно устойчивый метод глубокого расщепления
численного решения задач для уравнений параболического типа, содержащих смешанные диеренциальные операторы, основанный на более глубоком, чем в существующих методах, расщеплении смешанных производных и их неявной аппроксимации.
Исследования метода глубокого расщепления на примерах решения задач анизотропной теплопроводности показали зависимость вида конечно-разностной схемы, а
следовательно и устойчивости, от знаков коэициентов теплопроводности при смешанных производных. Показано, что абсолютная устойчивость схемы метода глубокого
расщепления имеет место в случаях, когда смешанный конечно-разностный оператор
ориентирован по вектору теплового потока. В связи с этим разработан специальный
алгоритм, сохраняющий абсолютную устойчивость метода.
Представленный метод является достаточно общим и может быть применен не только к задачам анизотропной теплопроводности, но и к другим задачам математической
изики.
Список литературы
[1?
Яненко Н.Н.
[2?
Самарский А.А.
[3?
[4?
[5?
Метод дробных шагов решения задач математической изики. Новосибирск: Наука, 1967. 196 с.
Теория разностных схем. М.: Наука, 1983. 616 с.
Метод переменных направлений с экстраполяцией по времени для параболических задач со смешанными производными // Вычисл. технологии. 1996. Т. 1, ќ 2.
С. 99103.
Формалев В.Ф.
Формалев В.Ф., евизников Д.Л.
Численные методы. М.: Физматлит, 2004. 400 с.
Неявный метод дробных шагов с расщеплением смешанных диеренциальных операторов // Вычисл. технологии. 1998. Т. 3, ќ 6. С. 8291.
Формалев В.Ф., Тюкин О.А.
Поступила в редакцию 23 марта 2010 г.,
с доработки 28 мая 2010 г.
1/--страниц
Пожаловаться на содержимое документа