CC BY
51
Вычислительные технологии
Том 5, № 5, 2000
ОБ ОДНОМ МЕТОДЕ РЕШЕНИЯ НЕКОТОРЫХ ТРЕХМЕРНЫХ УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ *
А. М. Гришин, А. С. Якимов Томский государственный университет, Россия e-mail: [email protected]
The finite-difference schemes for a three-dimensional partial differential equation of the first order have been obtained. This schemes based on the iteration-interpolational method and Lax-Wendroff approach are absolutely stable in the ordinary sense and characterized by dissipative properties. The algorithm of the method is presented and the convergence is demonstrated by test examples for some partial differential equations.
Известно [1, 2], что явления диссипации, подобные вязкости и теплопроводности, оказывают на скачок (разрыв решения) сглаживающее влияние. Разностные схемы, использующие искусственную вязкость [3, 4], успешно применяются при газодинамических расчетах. Их недостатком является “размазывание” разрывов и “подгонка” скачков.
Разностные схемы, предложенные Лаксом — Вендроффом [5, 6] для решения газодинамических задач, также обладали диссипативными свойствами. Однако аппроксимацион-ная вязкость в разностных схемах появлялась не искусственно, а из условия повышения аппроксимации некоторой исходной разностной схемы. Этот принцип гибридизации двух линейных разностных схем [5, 6] в итоге приводил к затуханию высокочастотных компонент решения в процессе вычислений. В отечественной [1] и зарубежной [4] литературе есть сведения, что для успешного решения этого класса задач перспективны разностные схемы с более высоким порядком (выше второго) аппроксимации по пространству.
Вопрос о свойствах многоточечных аналогов производных повышенного порядка аппроксимации изучался, в частности, в теоретических работах [7, 8]. В [7] исследованы неявные схемы максимального порядка точности для гиперболических систем. В статье
[8] рассмотрены разностные задачи четвертого порядка точности на пятиточечном шаблоне для уравнений Пуассона и Гельмгольца. В работе [9] для многомерных уравнений с первыми и вторыми производными получены неявные высокоточные схемы и приводятся некоторые результаты расчета течения вязкой несжимаемой жидкости.
В дальнейшем в статье [10] на основе второго приближения к искомому решению итерационно-интерполяционного метода (ИИМ) [11], поочередного интегрирования по пространственным переменным и подхода Лакса — Вендроффа [2] для трехмерного уравнения
* Работа выполнена в рамках Федеральной целевой программы: Государственная поддержка интеграции высшего образования и фундаментальной науки на 1997 - 2000 г., программы “Университеты России", а также при финансовой поддержке Российского фонда фундаментальных исследований, гранты №98-0103005, 99-01-00363.
© А. М. Гришин, А. С. Якимов, 2000.
неразрывности получены явно-неявные условно устойчивые разностные схемы, обладающие диссипативными свойствами и внутри области определения погрешностью аппроксимации в обычном смысле — 0(т2 + ^2т=1 к^)- При решении трехмерных уравнений математической физики для экономии времени расчета необходимо, чтобы разностная схема была абсолютно устойчивой по начальным данным-
Цель данной работы показать, что последнему требованию отвечают разностные схемы, полученные на основе алгоритма ИИМ [11, 12], в котором неявно присутствует подход Дугласа — Рэкфорда (схема стабилизирующей поправки [13]).
1. Постановка задачи и выбор метода решения
Рассмотрим уравнение в частных производных из [1]
ди д (и и) д(иу) д(и,ш)
dt + дх\
+
дх2
+
дхз
0, U > 0
с начальным условием
U |t=o = p(x).
:i.2)
Пусть для определенности величины и(Ь,х), у(Ь,х), т(Ь,х) положительны и заранее известны при 0 < Ь < Ьк - Тогда граничные условия в общем случае задаются в виде [1]
U |xi=0 = Pl(t,X2,X3), U |x2=0 = P2(t,xi,x3), U |хз=0 = P3(t,Xi,X2).
:i.3)
Выпишем разностную схему ИИМ внутри параллелепипеда Q : х = (х1,х2,хз), (0 < хт < dm, m =1, 2, 3) при 0 <t < tk.
Для построения разностного аналога уравнения (1.1) введем сетку с координатами узлов х1,і = ihl (i = 0,1, 2, ..., N), х2j = jh2 (j = 0,1, 2, ..., M), х3,к = kh3 (k = 0,1, 2, ..., L),
tn = nr(n = 0,1, 2, ..., в) и обозначим T = Uu, R = Uv, S = Uw.
Воспользуемся операторами из статьи [10] при hm = const, m = 1, 2, 3
AiU =
Ui-l,j,k + 4 Ui,j,k + Ui+l,j,k 6 !
AiT
Ti-l,j,k - Ti
i+l,j,k
2hl
1, ..., N - 1, j = 1,...,M - 1, k =1,...,L - 1.
:i.4)
Операторы Л2и, Л3и, Д2Я, Д3$ по другим координатным осям 0Х2, 0Х3 получаются циклической заменой индексов. В результате разностные схемы для уравнения (1-1) согласно [12] имеют вид (точка “сверху” приписывается производной по времени)
Л^ = Д1Тп+1/3 + Д2Кп + Д3Бп, Л2(У = Д1Тп+1/3 + Д2Кп+2/3 + ДзБп,
A3U = Al Tn+l/3 + A2Rn+2/3 + A3Sn+l,
i = 1,...,N - 1, j = 1,...,M - 1, k =1,...,L - 1.
(i.5)
2. Метод построения разностной схемы на основе подхода Лакса — Вендроффа
Подход Лакса — Вендроффа [2, 5, 6], как было отмечено в [2], обобщается на многомерный случай. Он может иногда применяться и при наличии скачков (ударных волн), если не требуется очень тонкое пространственное разбиение. Тогда, вводя обозначения р = Цп, Р2 = Цу, Рз = Ц'ш, для уравнения (1.1) имеем
Цп+1 _ ЦТ
дхт
дх* + дх
т=1
т=1
дхт
(2.1)
где г,] = 1, 2, 3, т = г = ], т. е. при т =1 : г т =3: г =1, ] = 2,
2, ] = 3, при т =2: г = 1, ] = 3, и при
Дт
-Рт
ди
(2.2)
Чтобы воспользоваться формулами высокого порядка точности по пространству, сумму правой части уравнения (2.1) необходимо записать в развернутом виде относительно производных от функции Дт и Рт, т = 1, 2, 3. В результате в (2.1) дополнительно будет 18 слагаемых. Используя (2.2), получим, например, по направлению координаты ОХ1 (Д1 = п, Д2 = у, Дз = ы):
д
дх1
дх2 дхз
дх1
дп дР1 пд2Р1 дп дР2 пд2Р2 дп дРз пд2Рз
2.3)
дх1 дх1 дх1 ' дх1 дх2 ' дх1дх2 ' дх1 дхз ' дх1дхз
Расписав аналогично (2.3) производные по осям ОХ2, ОХз, имеем для уравнения (2.1
ит+1 — ит
дТ дК дБ
дхл дх2 дхз
+ 0.5т
дрг
¿г^
дп ду ды + тт- +
дхт дх2 дхз
Г
пд
Г
дРт
<т=1
дхг
уд
дРт
дхл
Г
\т=1
дхт
ыд
дх2
Г
дРт
<т=1
дхт
дхз
(2.4)
Разностные операторы для первых, вторых и смешанных производных в квадратной скобке правой части уравнения (2.4) внутри области определения Q записываются на четырех-, пяти- и восьмиточечном шаблоне соответственно. В плоскостях, отстоящих на шаг по пространству (кт, т = 1, 2, 3) от граничных плоскостей хт = 0, хт = dm, т = 1, 2, 3, смешанные производные записываются на десятиточечном шаблоне. Для разрешения пятиточечной скалярной прогонки на граничных плоскостях хт = dm, т = 1, 2, 3 записываются дополнительные разностные соотношения, которые представляют собой аппроксимацию исходного уравнения (1.1) [14].
Формулы высокого порядка аппроксимации по пространству для первых, вторых и смешанных производных можно легко получить методом неопределенных коэффициентов на четырех-, пяти- и т.д. точечном шаблоне, записывая эти производные через ряд Тэйлора.
У искомых или заданных функций и, V, /ш, не входящих под знак частных производных, в разностях отсутствуют операторные значки Дт, Лт, ... и поэтому они аппроксимируются внутри области определения независимым образом, например, в центре любого (четырех-, пяти-, и т. д.) шаблона как , vij,k, ... При написании операторов смешанных
д2” д2
производных предполагалось, что —------- — = ——-— (Лтр = Лр т), т = р, т, р = 1, 2, 3.
дхтдхр дхрдхт
Разностные операторы внутри области определения для первых производных на четырехточечном шаблоне записываются внутри параллелепипеда Q следующим образом:
■г-; гр Ti—2,3,к 8Ti—1,3,к + 8Ti+1 ,3,к Т^2,3,к
У’Т =-------------------12^1----------------■
г = 2, ..., N - 2, з = 2,...,М - 2, к = 2,...,Ь - 2. (2.5)
Операторы У2Я, У3Б в (2.5); vЛ2R, 1иЛ3Б и ниже в (2.7); vЛ2Ri, 1 ,к, 'шЛ3Б^> 1 в (2.8); ,шЛ1,3Т, гшЛ2, 3R в (2.9); гшЛ1 , 3Т1 ^,к в (2.10); У3Бм,з,к в (2.12) получаются один из другого круговой заменой индексов.
Разностные операторы для первых производных в приграничных плоскостях хт = Нт, хт = Ат - кт, т =1, 2, 3 на четерехточечном шаблоне (рис. 1) имеют вид
-2Т0,3,к — 3Т1,3,к + 6Т2,3,к — Т3,3,к
УіТі
VI-і =
6Н1 ’
2ТМ,з,к + 3ТМ-1,з,к - 6ТМ-2,з,к + ТМ-3,з,к
6Н1
з = 2, ..,М - 2, к = 2,...,Ь - 2. (2.6)
Разностный оператор для вторых производных внутри параллелепипеда Q записывается на пятиточечном шаблоне
иЛіТ
иі,3,к (-Ті-2,з,к + 16Ті-1,з,к - 30Ті,з,к + 16Ті+1,з,к — Ті+2,з,к)
12Н\
і = 2, ..., N - 2, і = 2,...,М - 2, к = 2,...,Ь - 2. (2.7)
Разностный оператор в приграничной плоскости х1 = Н1 имеет вид
д ГТ1 _ и1,3,к (11Т0,з,к - 20Т1,з,к + 6Т2,з,к + 4Т3,з,к - Т4,3,к)
мЛ1Т1 =------------------------------------------------------Щ-’
і = 2, ..., N - 2, і = 2,...,М - 2, к = 2,...,Ь - 2. (2.8)
Разностный оператор для смешанной производной на восьмиточечном шаблоне внутри области определения Q записывается
УЛ1,2Т = Vі,з,к [(Ті+1,з+1,к - Ті-1,з+1,к - Ті+1,з-1,к + Ті-1,з-1,к)/3-
- (Ті+2,з+2,к - Ті-2,з+2,к - Ті+2,з-2,к + Ті-2,з-2,к)/48]/(^1^2), і = 2, ..., N - 2, і = 2,...,М - 2, к = 2,..., І - 2. (2.9)
Рис. 1. Система координат и обозначения.
Разностный оператор для смешанной производной
д 2Т дх1дх2
в приграничной плоскости
х1 = Н1 на двенадцатиточечном шаблоне по направлению оси координат 0Х1 имеет вид
^Л1,2Т1 = г°1,з,к [5(Т2,з+1,к - Т0,з+1,к - Т2,з-1,к + Т0,з-1,к)/24 - (Т2,з+2,к-
-Т0,з+2,к - Т2,з-2,к + Т0,з-2,к)/24 + (Т3,з + 1,к - Т0,з + 1,к - Т3,з-1,к + Т0,з-1,к)/6 -(Т4,з + 1,к - Т0,з + 1,к - Т4,з-1,к + Т0,з-1,к)/16] /(^1^2),
і = 2,...,М - 2, к = 2,...,Ь - 2. (2.10)
Вблизи угла хт = 0, т = 1, 2, 3 в узле пересечения плоскостей хт = Нт, т =1, 2, 3 (на
і /і\ д 2^1,1,1 д 2Т^-1,1,1
рис. 1 узел А) разностные операторы для ——-—, —-— ---------------- записываются следующим
образом:
дх2дх3 ’ дх1дх2
Л2,3^1,1,1 = ^1,1,1[(^1,2,2 - ^1,2,0 - ^1,0,2 + ^1,0,0)/8 +
+ (^1,3,2 - ^1,3,0 - ^1,0,2 + К-1,0,0)/6 - (^1,4,2 - ^1,4,0 - ^1,0,2 + ^1,0,0)/8+
+ (^1,4,3 - ^1,4,0 - ^1,0,3 + ^1,0,0)/12 - (^1,4,4 - ^1,4,0 - ^1,0,4 + ^1,0,0)/32]/(^2^3), Л1,2ТЫ-1,1,1 = -1,1,1 [(Т^,2,1 - ТЫ-2,2,1 - Т^,0,1 + ТЫ-2,0,1 )/8+
+ (Т^,2,1 - ТЫ-3,2,1 - Т^,0,1 + ТЫ-3,0,1)/6 - (ТЫ,2,1 - ТЫ-4,2,1 - Т^,0,1 + ТЫ-4,0,1)/8+
+ (Т^,3,1 - ТЫ-4,3,1 - Т^,0,1 + ТЫ-4,0,1)/12 - (Т^,4,1 - ТЫ-4,4,1 - ТЫ,0,1 +
+Т^-4,0д)/32]/(^1^2). (2.11)
Аналогично выписываются разностные операторы для смешанных производных в остальных приграничных плоскостях х2 = к2, х3 = к3, хт = ¿т - кт, т = 1, 2, 3, линиях и узлах вблизи вершин параллелепипеда Q.
Предполагается, что на поверхностях, ограничивающих область Q : хт = йт, т =
1 , 2, 3, справедливо исходное уравнение (1.1). Аппроксимируем это уравнение независимым образом.
На плоскости х1 = ¿1 (0 < хт < йт, т = 2, 3) разностные операторы записываются
для первых производных в виде
У1Т
11Тм,1,к - 18^-1,п,к + 9Ты-2,1,к - 2^-
N
- N-3,3,к
V2R
6к1
RN.i—2.k 8RN.i-1.k + 8R
*'N¿+1$
- R
■N^+2^
12^2
3 = 2, ..., М - 2, к = 2, ..., Ь - 2.
(2.12)
Так как в логической схеме ИИМ [12] для решения трехмерного уравнения (1.1) неявно присутствует подход Дугласа Рэкфорда [13], то, используя обозначения (1.4), (2.5), (2.7), (2.9), для уравнений Лакса — Вендроффа (2.4) можно выписать разностные соотношения, аналогичные (1.5). В целях краткости дальнейших выкладок введем обозначения
f = 0.5т [и(Л1^ + Л^3 Б) + v(Лl)2T + Л2,3Б) + и^Л^Т + Л2)зR)],
С = У1и + У^ + У3'ш, тогда получим разностные схемы внутри области определения Q
А1(цп+1/3 - цп)/т = (Д1 + 0.5тиЛ1 )Тп+1/3 + (Д2 + 0.5™Л2^п+
+ (Д3 + 0.5™Л3)Бп + 0.5т (У1Тп+1/3 + V2Rn + У3 Бп)Сп + Р, А2(ип+2/3 - ип)/т = (Д1 + 0.5тиЛ1)Тп+1/3 + (Д2 + 0.5тvЛ2)Rn+2/3+
2)
+(Д3 + 0.5™Л3)Бп + 0.5т (У1Тп+1/3 + У2Rn+2/3 + У3Бп)Сп + [п, А3(ип+1 - ип)/т = (Д1 + 0.5тиЛ1)Тп+1/3 + (Д2 + 0.5^Л2^п+2/3+ + (Д3 + 0.5т^Л3)Бп+1 + 0.5т (У1Тп+1/3 + У2Rn+2/3 + У3Бп+1)Сп + ¡п г = 2, ..., N - 2, з = 2,...,М - 2, к = 2,...,Ь - 2.
(2.13)
(2.14)
(2.15)
Практически при исследовании устойчивости разностных схем (2.13) - (2.15) и для написания программы расчета удобно пользоваться вместо двух последних уравнений (2.14), (2.15) более простыми разностными уравнениями [2]. Вычтем из второго уравнения первое, а из третьего второе, тогда вместо двух последних уравнений системы (2.13) - (2.15)
получим
(А2ип+2/3 - Ахип+1/3)/т = (Д2 + 0.5^Л2)^п+2/3 - Rn)+ +0.5тУ2^п+2/3 - Rn)Gn + (А2ип - Ахип)/т, (А3ип+1 - А2ип+2/3)/т = (Д3 + 0.5т^Л3)(Бп+1 - Бп)+ +0.5тУ3(Бп+1 - Бп)Сп + (А3ип - А2Ип)/т, г = 2, ..., N - 2, з = 2,...,М - 2, к = 2,...,Ь - 2.
(2.16)
(2.17)
Разностные схемы в приграничной плоскости х 1 = Н 1, при использовании операторов из (1.4), (2.6), (2.8), (2.10), перепишутся
X = У 1и 1 + У 2V 1 + У3^ 1,
f1 = u(Л1,2R1 + Л1,3Б1) + v (Л1,2 Т1 + Л2,3Б1) + ^(Л1,3Т1 + ^,3^^
0 < хт < йт, т = 2, 3,
А1(ип+1/3 - Щ)/т = (Д1 + 0.5тиЛ1 )Тп+1/3 + (Д2 + 0.5^Л2^п+
+ (Д3 + 0.5т^Л3)БП + 0.5т (У1Тп+1/3 + У2Rn + У3Бп)Хп + 0.5тf1n,
(А2ип+2/3 - Ахип+1Г3)/т = (Д2 + 0.5^Л2)^п+2/3 - Rn)+
+0.5тУ2^п+2/3 - Rn)Xn + (А2ип - Ахип)/т,
(А3иГ+1 - А2ип1+2/3)/т = (Д3 + 0.5т^Л3)(БП+1 - Бп)+
+0.5тУ3(Бп+1 - БГ)Хп + (А3ип - А2!ЛГ)/т,
з = 2,...,М - 2, к = 2,...,Ь - 2. (2.18)
Аналогично выписываются разностные схемы для приграничных плоскостей хт = Нт, т = 2, 3 и хт = йт - Нт, т = 1, 2, 3. Разностные уравнения на линии пересечения плоскостей х1 = Н1, х3 = Н3, при использовании операторов (1.4), (2.5)-(2.8), (2.10), имеют
вид
Х1,-,1 = У1и1,-,1 + У 2V1,j,1 + У 3 ^1,^,1, f1,j,1 = u(Л1,2R1,j,1 + Л1,3Б1,j,1) + v(Л1,2T1,j,1 + Л2,3Б1,j, 1) + ^(Л1,3Т1-,1 + Л2,3R1,j,1) ,
А1(ип+1/3 - — /т = (Д1 + 0.5тиЛ1)Т1п+1/3 + (Д2 + 0.5тvЛ2)Rnj 1 +
+ (Д3 + 0.5тшЛ3)Б1^. 1 + 0.5т (У^+У3 + У 2#— + У3 Б-д)Х1— + 0.5тЯ-д, (А2ип+2/3 - А1ип+1/3)/т = (Д2 + 0.5™Л2)(:йГ+У - ДГ-д) +
+0.5тУ2^п,+,1/3 - + №¿1 - А1 и—)/т,
^3^+1 - А2иШ+1/3)/т = (Д3 + 0.5тшЛ3)(Бг+1 - БГл)+
+0.5тУ3(Б^1 - Б^-д)Х1^,1 + (А3ЦТ,,. 1 - А^^/т,
з = 2,...,М - 2. (2.19)
Аналогично можно записать разностные схемы для остальных одиннадцати линий пересечения плоскостей (х1 = Н1, х2 = Н2), (х3 = Н3, х2 = Н2) и (хт = йт - Нт, хр = йр - Нр), т, р =1, 2, 3, т = р. Разностные уравнения вблизи начала координат: в узле пересечения плоскостей хт = Нт, т = 1, 2, 3, при использовании операторов (1.4), (2.6), (2.8), (2.10), (2.11), записываются
Х1,1,1 = У1и1,1,1 + У2V1,1,1 + У3W1,1,1, f1,1,1 = и (^^1,2 ^^1,1,1 + Л1,3Б1,1,1) + +v(Л1,2T1,1,1 + Л2,3Б1,1,1) + ^(Л1,3Т1,1,1 + ^^2,3^^1,1,1) ,
А^ип+У3 - и^д)/т = (Д1 + 0.5тиЛ1)ТЩ^+1/3 + (Д2 + 0.5тvЛ2)Rn,1,1+
+ (Д3 + 0.5т^Л3)Бп 1 1 + 0.5т (У1Т1п+1/3 + У 2^1 1 1 + У3БГ1 ^1 1 + 0.5тЛ\ 1,
(А2ип+1/3 - А^Г+^/т = (Д2 + 0.5тvЛ2)(Rn,^2/3 - RnlЛЛ)+
+0.5тУ2^?,+2/3 - Д?д,1 )Х£М + (А2^^1,1 - А^д )/т,
(А3ипЦ - А2ип+1/3)/т = (Д3 + 0.5™Л3)(Бп,+1 - Бп,1,1)+
+0.5тУ3(БГ+,1 - БП,1,1)ХПм + (А3и^1,1 - А2и£м )/т. (2.20)
Аналогично разностным уравнениям (2.20) можно записать разностные схемы для остальных семи узлов пересечения плоскостей хт = Нт, хт = йт - Нт, т = 1, 2, 3.
Разностные уравнения на ребре (0 < х1 < й1) пересечения граничных плоскостей х2 = й2, х3 = й3, а также на линиях пересечения плоскости х1 = й1, плоскостей хт = Нт, хт = йт - Нт, т = 2, 3 ив вершине хт = йт, т = 1, 2, 3 при помощи операторов из (2.5), (2.6), (2.12) имеют вид
п+1/3 п п+1/3 п п
(ицг - г-Мх)/т = -У- у2апмх - У3Бих, (ипм2/3 - ии+£/3)/т = -У2(яМ+!/3 - Л^ь).
(ии+1 - иМ+2/3)/т = -У3(вМ^1 - = 2, N - 2,
(опих - и^_мх)/т = -У1ТГМ1/3 - У2ДГ,м,£ - У3Б,мх
(цГМ2Х3 - и;:м,1Х3)/т = -у2(лгм,Х3 - я;.м,ь). (ипм1^ - цг+М!ь3)/г = -у^ми - б'Гм.х),
(и£1/'м.ь - их —1,м,х)/т = -У1ТГ-1/- У2ДТ_1м,х -
(ип+2/3 ип+1/3 ) У7 ( пп+2/3 рп )
(UN—1,U,L и N —1,м,Ь)I т = У 2(RN—1,м,Ь RN—1,U,L),
(их—\мл - и£1мЛ)/т = -У3(ЯГ—1,м,х - —1,м,х),
п
,м,х - UN — 1,м,х)/ 1 = У 3(SN—1,м,х - ^ — 1,м,х)
(иК5 - ип,м,х)/т = -У^+м/Х - У2Лймг - У3б’У
^,м,Х UN,U,L)/т = y11N,U,L У 2RN,U,L У3SN,U,L,
п+2/3 п+1/3 п+2/3 п
(UN,U,L UN,U,L)|т = У2(RN,U,L RN,U,L),
(Цйи - )/т = - У3(5Т+Д1 - ,м,ь)■ (2.21)
Разностные схемы на плоскости х1 = й1 при помощи операторов (2.12) записываются
(и^1/3 - Щ )/т = -УlTN+1/3 - У2RnN - У3БГ, (UN+2/3 - и^1/3)/т = -У2№+2/3 - RnN),
(UN+1 - UN+2/3)/т = ^(Б^1 - БnN),
з = 2,...,М - 2, к = 2,...,Ь - 2. (2.22)
Не представляет труда записать разностные уравнения для плоскостей хт = йт, т =
2, 3, для ребра (0 < х2 < й2) пересечения плоскостей х1 = й1, х3 = й3 и ребра (0 < х3 < й3) пересечения плоскостей х1 = й1, х2 = й2, используя операторные выражения, подобные (2.5), (2.6), (2.12).
Наконец, на граничных плоскостях хт = 0, т = 1, 2, 3 имеем
и-1 = —„-о, з = 0,...,М, к = 0,...,Ь,
4,0+1 = Р231и=о, г = 0, к =0,...,Ь,
и—О = Р^о. г = 0, ..., N. з = 0,...,М. (2.23)
3. О погрешности аппроксимации и устойчивости разностной схемы
Внутри параллелепипеда Q: x £ (0 <xm<dm, m=1, 2, 3) получается B1=(N — 1)(M — 1)(L — 1) разностных уравнений (2.13) - (2.15), (2.18), (2.19) плюс пять систем уравнений вида (2.18), плюс одиннадцать разностных схем системы уравнений типа (2.19) и плюс семь разностных уравнений, подобных уравнениям (2.20), для определения Bi неизвестных U.
На граничных плоскостях xm = dm, m = 1, 2, 3 выписываются B2 = (N — 1)(M — 1) + (N — 1)(L — 1) + (M — 1)(L — 1) + N — 1 + M — 1+L — 1 разностных уравнений вида (2.21), (2.22) для определения B2 неизвестных. На граничных поверхностях xm = 0, m = 1, 2, 3 имеем также B3 = NM + NL + ML конечных алгебраических выражений (2.23) для определения B3 заданных функций из (1.3).
Легко заметить, что расчет по однородным разностным схемам (2.13)-(2.15), (2.18) — (2.23) экономичен (с учетом условий абсолютной устойчивости (3.6), (3.8)), так как для получения решения в узлах области определения понадобится O(Bi) арифметических действий, пропорциональное числу узлов области Q.
Если существуют ограниченные вплоть до шестого порядка производные по пространству и третьего порядка по времени, то для уравнения (1.1) локальная погрешность аппроксимации разностной схемы (2.13) — (2.15) при hm = h, m = 1,2, 3 внутри области определения Q на решении задачи — 0(т2 + h4). Условие hm = h, m = 1, 2, 3 для ПЭВМ Pentium-2 (ОП — 128 Мб) не является обременительным.
Локальная погрешность аппроксимации разностных схем (2.21), (2.22) на уравнении (1.1) — 0(т + h3), а схем (2.18) — (2.20) на решении задачи — 0(т2 + h3). Имеются в виду первые уравнения этих систем по направлению координаты 0Xi, так как вторые уравнения типа (2.16), (2.17) по направлению координат 0X2, 0X3 необходимы только для устойчивости разностных схем.
Начальное условие (1.2) удовлетворяется точно. Граничные уравнения (2.23) на соответствующей граничной плоскости из (1.3) аппроксимируются точно, так как они заданы от времени явно на целом слое.
Методом Фурье в приближении “замороженных” коэффициентов [1, 13, 14] можно найти условие устойчивости явно-неявных разностных уравнений (2.13) — (2.15) по начальным данным при u = v = w = c (c = const). Принцип “замороженных” коэффициентов [14] позволяет ориентироваться на эвристическом уровне строгости и при исследовании устойчивости нелинейных задач.
Если Z — множитель роста разностных схем (2.13), (2.16), (2.17), то спектральное условие Неймана имеет вид (3.6). Сделаем стандартную подстановку [1]
Un+i = ZUn, Ujk = exp [I (gixii + g2X2,j + 93X3,1,)], (3.1)
где I = \J — 1, gm = 0, ±1, ±2, ..., m = 1, 2, 3 — гармоники при переходе со слоя на слой. Подставим (3.1) в (2.13), (2.16), (2.17) и для сокращения дальнейших выкладок введем обозначения:
f2 = v|t=o = c, /3 = w|t=o = c, fi = u|t=o = c, bm = fmhm1 sin 9mhm, Tm = (1 + 2 COS2 gmhm/2)/3т,
am T fmhm (16 sin gmhm/2 sin gmhm)/6, m 1, 2, 3
di,j = т(hihj)-ififj sin gihi singjhj, qid = т(hihj)-1 fifj sin 2g¿hi sin 2gjhj,
i,j = 1, 2, 3, i = j,
Y = 2 (d1,2 + d1,3 — qi,2 — ?1,э) + 2 (d2,3 — ?2,з)- (3-2)
Тогда по каждому из направлений получим, замечая, что G|t=0 = 0
Un+1/3(ri + ai + Ibi) = Un(ri - a2 - аз - Ib2 - Ib:i - Y), (3.3)
Un+2/3(r-2 + a2 + Ib2) = riUn+1/3 + Un(r2 - ri + a2 + Ib2), (3.4)
Un+1(r3 + a3 + Ib3) = r2Un+2/3 + Un(r3 - r2 + a3 + Ib3). (3.5)
Выразим Un+1/3 из (3.3) и подставим в (3.4), затем найдем Un+2/3 из полученного соотношения и подставим в (3.5). В результате окончательно имеем уравнение для определения Z:
z = (H + IB)/(C + ID), |Z | = [(H2 + B 2)/(C2 + D2)]0-5, (3.6)
где
H = F + s - r1 r2Y, C = F + a, D = E + b, B = E + q, I2 = -1,
F = r3r2r 1 + ri(a2a3 - b2b3) + r2(aia3 - bib3) + r3^ai - b2bi) + a^a3 - aib2b3 - a2bib3 - a3b2bi,
s = air2(r3 - ri) + ria2(r3 - r2), a = r3r2ai + r3a2ri + a3r2ri, E = ri(b2a3 + a2b3) + r2(bia3 + aib3) + r3(b2ai + a2bi) + a2a3bi + aia3b2 + aia2b3 - bib2b3, q = bir2(r3 - ri) + rib2(r3 - r2),
b = r3r2bi + r 3 b2 r 1 + b3r2ri. (3.7)
Для двухслойных разностных схем необходимое условие устойчивости по начальным данным записывается [1, 13]: |Z| < 1. В силу того, что | singmhml < 1, | cosgmhml ^ 1, найдем оценку сверху для Z в (3.6). В частности, знаменатель будет минимален в (3.6) при gmhm/2 = kn ( m = 1, 2, 3, k = 1, 2,...). Тогда из выражений (3.2), (3.6), (3.7) имеем rm = 1/т, sin gm hm = sin gm hm/2 = 0, m =1, 2, 3, a = b = s = q = E = 0, C = F = 1/т3. При этом числитель в (3.6) равен H = F = 1/т3, так как s = Y = B = 0. В результате окончательно получим
|Z | = H/C = 1. (3.8)
Из (3.6), (3.8) видно, что разностные схемы (2.13), (2.16), (2.17) безусловно устойчивые
при т > 0. Аналогично доказывается абсолютная устойчивость остальных разностных
схем (2.18) - (2.22).
Рассмотрим разностные схемы в направлении координаты OX1. Это будет разностная схема (2.13) с пятидиагональной матрицей, алгебраические выражения из первого соотношения типа (2.18) при x1 = h1, и Xn-1 = dn-1 - h^-i, а также (2.22) при x1 = d1.
Воспользуемся сокращенной записью этих уравнений:
Hi Ui-2 + BiUi-i + CiUi + DiUi+i + EiUi+2 = Fi, (3.9)
i = 2, ..., N - 2,
U0 = P1|Xl=0, (3.10)
H1U0 + B1U1 + C1U2 + D1U3 + E1U4 = Fi, (3.11)
Нм-1им-4 + Вм-\Un-3 + См-1им-2 + -1и«-1 + Ем-1им — рм-1, (3-12)
Вм им-3 + См им -2 + ^м им-1 + Ем им — рм • (3-13)
Так как матрица коэффициентов уравнений (3-9)-(3-13) пятидиагональна, то для на-
хождения неизвестных воспользуемся методом скалярной прогонки [15]- Используя методику работы [15], нетрудно доказать устойчивость пятиточечной прогонки для разностных схем (2.13), (2-18)-(2-23)-
Для разностных схем вида (2-22), например, на плоскости х1 — ¿]^, ограничивающей область определения задачи Q, условие устойчивости прогонки имеет вид
т < К/(1.5и), (3-14)
где и — тах им^^, 3 — 1, •••, М — 1, к — 1, Ь — 1-
Нетрудно показать [12, 16], что разностные схемы (2-13) - (2-15), (2-18) - (2-22), получаемые при помощи поочередного интегрирования по пространственным переменным ИИМ, обладают свойством слабой консервативности (консервативность удовлетворяется с порядком аппроксимации 0(т + К2)) в смысле результатов работы [17]-
4. Сходимость и примеры применения метода
Сходимость разностных схем (1.5), (2.21) - (2.23) установим практически на последовательности сгущающихся сеток по пространству при решении трехмерной дифференциальной задачи (4.1) с граничными условиями (4.2)
f + ± £ = F. UU = f, (4.1)
m=1
f = 1 + x? + x3 + (X2 - 1)z, F = 4ft3 + z(1 + t4)[x1-1 + x3-1 + (X2 - 1)z-1],
J-4\fl I (™ 1 \Z I TT\ _ (Л 14- \Tl I ( 1 \Z I I
U|xi=0 = (1 + t )[1 + (x2 — 1)z + x3]. U 1x2=0 = (1 + ¿4)[1 + (—1)z + X1 + x3]
=0
U U3=0 = (1 + t4)[1 + X1 + (x2 — 1)z ]- (4.2)
Точное решение задачи (4.1), (4.2) : U = f (1 + t4) в кубе Q : [0 < xm < 1, m =
1, 2, 3] известно при 0 < t < tk. Были взяты следующие значения входных данных: hm = h, m = 1, 2, 3, z = 6, tk = 7, т = 0.005. Программа составлена на языке Фортран-77, расчет проводился на ПЭВМ Pentium-2 (200 МГ, Транслятор Power Station 4) с двойной точностью.
На рис. 2 изображена максимальная относительная погрешность е = I—щ—I 100 %
(U — точное, U — приближенное численное решение) от времени. Кривые 1, 2, 3 отвечают пространственной сетке h = 0.1 (N = 11, M = 11, L = 11), h = 0.05 и h = 0.025
соответственно. Штриховые кривые получены по предлагаемому методу и время, затраченное для их расчета (t0): t0>1 =2.1 мин (h = 0.1), t0,2 = 15 мин (h = 0.05), t0,3=75 мин (h = 0.025). Сплошные кривые на рис. 2 отвечают расчету по неявной разностной схеме бегущего счета (29) на стр. 345 из [1]:
Ui,+k — Ui>',j,k + (Ui,j,k — Ui-1,j,k)n+1 + (Ui,j,k — Ui,j-1,k)П+1 + (Ui,j,k — Ui,j,k-1)n+1 = Fn
T + h1 + h2 + h3 = j ’
1-75 3.5 5.25 г
Рис. 2. Максимальная относительная погрешность решения задачи (4.1), (4.2) от времени. Кривые 1, 2 и 3 соответствуют пространственным сеткам к = 0.1, к = 0.05 и к = 0.025. Сплошные кривые — разностная схема (4.3), штриховые кривые — двухслойная разностная схема по времени из алгоритма ИИМ (1.5), (2.21) - (2.23).
г = 1, ..., N ] = 1, ..., М, к =1, ..., Ь, (4.3)
которая в трехмерном случае имеет погрешность аппроксимации 0(т + Л + Л2 + Л3), а время расчета их составило ¿0>1 = 1 мин, ¿0,2 = 8 мин, ¿0,3 = 46 мин соответственно. Из сравнения кривых видно, что на конкретном примере точность расчета по предлагаемому методу лучше в 5 раз, а время ¿0 3 у них различается только в 1.6 раза.
Уравнение Лакса — Вендроффа (2.4) записывается только для внутренних узлов области определения задачи и имеет в общем случае параболический вид. Поэтому сходимость разностных схем (2.13), (2.16)-(2.20) проверим на решении дифференциального уравнения (4.4) с граничными условиями первого рода (4.5)
ди Л ди Л д2и Л д2и п тт.
С4_дГ + С12-, дХ- =+Сз ^ дх дх + ^ и|<=0 = / (4.4)
—=1 —=1 — —=1,р=— ^
^ = 4С4/^3 + (1 + ¿4)[х1-1 + ХТ1 + (Х2 - 1)2-1] - С2^(^ - 1)(1+ ¿4)[х1-2 + Х3-2 + (Х2 - 1)^-2],
и|Х1=0 = (1 + ¿4)[1 + (х2 - 1)2 + Х3]5 и|х2=0 = (1 + ¿4)[1 + (-1)2 + Х1 + Х3]5
и|х3=0 = (1 + ¿4)[1 + Х1 + (х2 - 1)2], и|х1 = 1 = (1 + ¿4)[2 + (х2 - 1)2 + Х3],
и|х2 = 1 = (1 + ¿4)[1 + Х1 + Х3^ и|х3 = 1 = (1 + ¿4)[2 + Х1 + (х2 - 1)2]. (4.5)
Использовались следующие значения входных данных: = Л, т = 1, 2, 3, = 7,
£ = 6, с4 = 100, с1 = 800, с2 = 1, с3 = 2. Из таблицы видно, что численное решение задачи
сходится к точному при измельчении шагов разностной сетки, так как є уменьшается (¿о — расчетное время в минутах на ПЭВМ). Отметим, что погрешность аппроксимации разностных схем, при помощи которых решена задача (4.4), (4.5), — 0(т + Л2).
Т а б л и ц а
Т 0.004 0.004 0.004 0.007
Н 0.1 0.05 0.025 0.05
є 18.2 3.7 1.17 4.35
Іо 5 19 62 11
Правильность предложенного алгоритма ИИМ (2.13), (2.16) - (2.23) целесообразно также проверить, например, при решении квазилинейного уравнения переноса [1], которое записывается в дивергентной форме и для трехмерного случая имеет вид
ди Л ди2
иг + °-5£
т=1
дх
0, и > 0
(4.6)
с начальным и граничным условиями
и|4=о = а : хз < Жз,о, 0 < жт < 1, т = 1, 2,
и|/=о = Ь : хз > Хз,о, ° < хт < 1, т = 1, 2,
и |Хт=о = а, т = 1, 2, 3, а>Ь> 0.
Тогда уравнение Лакса — Вендроффа (2.4) перепишется
(4.7)
(ип+1 - ип) 3 ди2
^+0*5 £
т=1
дх
0.5т
А 1 д 2и3 2 +
з дхт
т= 1 т
д2 и3 д2и3 д2 и3
+ т—+
3 \ дх1дх2 дх1дх3 дх2дх3
(4.8)
Находим решение, имеющее разрыв. Пусть наклон линии разрыва соответствует скорости О = (а + Ь)/2, тогда искомое решение имеет вид
и = а : х3 — х3,о < Ді, 0 < < 1, т =1, 2,
и = Ь : х3 — х3 о > Ді, 0 < < 1, т = 1, 2.
(4.9)
На рис. 3 изображен результат решения задачи (4.6) - (4.9) в сечении х = Ж1>0, х2 = ж2>0 при а = 2, Ь =1, ж—>0 = 0.5, = Л, т = 1, 2, 3 по разностным уравнениям ИИМ, по схеме
Лакса — Вендроффа [2] с погрешностью аппроксимации 0(т2 + ^—= Л—) и по уравнениям (4.3) (сплошные, штриховые и штрихпунктирные кривые соответственно).
В дальнейшем для сокращения назовем их подходами 1, 2 и 3. Шаг по времени т выбирался из условия устойчивости пятиточечной прогонки (3.14) для Л = 0.0125.
Кривые 1, 2, 3 отвечают моментам времени ^ = 0.05, ¿2 = 0.1, ¿3 = 0.2, жирные линии соответствуют точному решению задачи (4.9) (х3 = + х3,0, т = 1, 2, 3). Цифры с
двумя штрихами приписываются результату решения при Л = 0.05, с одним штрихом — Л = 0.025, а без штриха — Л = 0.0125. Расчет проводился до момента времени ¿3 = 0.2, когда граница Х3 = 1 практически не влияет на результат решения задачи.
При Л = 0. 05 имеются выбросы в численных решениях по предлагаемому подходу 1 (в одномерном случае [10] колебаний нет) и подходу 2, а подход 3 размазывает точное
Рис. 3. Результат решения задачи (4.6)-(4.9) в сечении х = Х1,0, Х2 = £2,0 вдоль координаты Х3 для моментов времени 1 — 0.05, 2 — 0.1, 3 — 0.2. Цифры с двумя штрихами отвечают Ь = 0.05, с одним штрихом — Ь = 0.025, а без штриха — Ь = 0.0125. Сплошные кривые — разностные схемы ИИМ, штриховые кривые — схема Лакса — Вендроффа, штрихпунктирные —
разностные уравнения (4.3).
решение (угол вдоль координаты х3) в силу большой схемной вязкости. Для Л = 0.05 при подходе 2 штриховые кривые 1" - 3" на рис. 3 не приводятся, так как амплитуда колебаний их еще выше, чем в одномерном расчете из [10]. При Л = 0.025 колебания в решении в подходе 2 остаются и исчезают только при Л = 0.0125, в то время как при подходе 1 точное решение (ступенька) с уменьшением шага по пространству отслеживается лучше (см. на рис. 3 сплошные кривые 1 - 3' и 1 - 3). Отметим, что подход 3 для Л = 0.025, Л = 0, 0125 при I = ¿2 и особенно при I = ¿3 вообще не воспроизводит точное решение (на рис. 3 см. штрихпунктирные кривые 2 и 3).
Заключение
1. Для численного решения трехмерного уравнения в частных производных первого порядка на основе ИИМ получены разностные схемы с погрешностью аппроксимации в обычном смысле внутри области определения на решении задачи 0(т2 + Л4) при = Л, т = 1, 2, 3.
2. При гладких коэффициентах переноса и подходе Лакса — Вендроффа разностные уравнения ИИМ — абсолютно устойчивы по начальным данным и диссипативны.
3. Лучший результат, полученный по разностным уравнениям ИИМ настоящей работы для тестового примера, подтверждает их эффективность при решении задачи с разрывом решения.
Список литературы
[1] КАЛИТКИН Н. Н. Численные методы. Наука, М., 1978.
[2] РИХТМАЙЕР Р., МОРТОН К. Разностные методы решения краевых задач. Мир, М., 1972.
[3] Войнович П. А., ЖмАкин А. И., Попов Ф. Д., Фурсенко А. А. О расчете разрывных течений газа. Препр. ФТИ им. А. Ф. Иоффе АН СССР. Л., №561, 1977.
[4] ОРАН Е., Борис Д. Численное моделирование реагирующих потоков. Мир, М., 1990.
[5] Lax P. D., WendrüFF B. Systems of Conservation Laws. Comm. Pure Appl. Math., 13, 1960, 217-237.
[6] Lax P. D., WendrüFF B. Difference Schemes for hyperbolic equations with high-order of accuracy. Ibid., 17, 1964, 381-398.
[7] ПААСоНЕН В. И. Абсолютно устойчивые разностные схемы повышенной точности для систем гиперболического типа. В “Численные методы механики сплошной среды". Новосибирск, 3, №3, 1972, 82-91.
[8] ВАЛИУЛЛИН А. Н. Многопараметрические итерационные процессы для разностных задач Дирихле повышенного порядка точности. Там же, 4, №2, 1973, 32-42.
[9] ПААСоНЕН В. И. Об одном методе построения высокоточных разностных схем и их применении в механике жидкости. Там же, 7, №6, 1976, 111-126.
10] Якимов А. С. Об одном методе решения линейного уравнения переноса. В “Вычислительные технологии". ИВТ СО РАН, Новосибирск, 4, №10, 1995, 322-332.
11] ГРИШИН А. М. Об одном видоизменении метода М. Е. Швеца. Инженернофизический журнал, 19, №1, 1970, 84-93.
12] ГРИШИН А. М., ЯКИМов А. С. Обобщение итерационно-интерполяционнго метода для решения трехмерного параболического уравнения общего вида. Вычисл. технологии. 4, №2, 1999, 26-41.
13] ЯНЕНКо Н. Н. Метод дробных шагов решения многомерных задач математической физики. Наука, Новосибирск, 1967.
14] Годунов С. К., РЯБЕНЬКИЙ В. С. Разностные схемы. Наука, М., 1973.
15] ШАМАНСКИЙ В. Е. Методы численного решения краевых задач. Изд-во АН УССР, Киев, 1963.
16] ЯКИМов А. С. Об одном методе расщепления. Численные методы механики сплошной среды, 16, №2, 1985, 144-161.
17] Войнович П. А., Жмакин А. И., Попов Ф. Д., Фурсенко А. А. Численное исследование течений газа с разрывами сложных конфигураций. Журн. вычисл. математ. и мат. физики, 19, №2, 1979, 1608-1614.
Поступила в редакцию 25 февраля 1999 г., в переработанном виде 10 мая 2000 г.
