CC BY
67
УДК 519.62
ЧИСЛЕННОЕ РЕШЕНИЕ НЕЛИНЕЙНОГО УРАВНЕНИЯ В ЧАСТНЫХ ПРОИЗВОДНЫХ ПАРАБОЛИЧЕСКОГО ТИПА
© 2007 г А.Б. Усов
This paper is devoted to the construction of the difference scheme of the finally difference method for the numerical solving of the nonlinear parabolic equation. The coordinate transformation was built for the receipt of the even differences net in the calculation plane. The physical plane was dynamic and curved.
Уравнения в частных производных параболического типа используются для описания разнообразных физических процессов: распространения тепла в твердых средах, загрязняющих веществ в жидкости и многих других. Поэтому весьма актуальными в настоящее время являются вопросы разработки численных методов решения этих уравнений. Несмотря на большое количество работ по этому вопросу, например [1-5], до сих пор открытым остается вопрос о разработке новых численных методов решения нелинейных уравнений параболического типа в переменной по времени области криволинейной формы. Нестационарные пространственные области создают значительные сложно -сти при численном решении уравнений параболического типа. В статье предлагаются преобразования координат, позволившие получить равномерную разностную сетку в вычислительной плоскости, несмотря на изменчивость во времени, криволинейность физической области, в которой ведется рассмотрение. Задача решается методом конечным разностей.
Математическая постановка задачи
Введем декартову прямоугольную систему координат Oxz (рис. 1). В двумерной области с криволинейными границами рассмотрим нелинейное уравнение параболического типа вида
dB n3B SB 1
-+ B-+ B-= —
д t д x д z A
+ F (x, z, t), (1)
где t - временная координата (0 < t < Д); Д =const; x, z -пространственные координаты (0 < x < L; h < z < H); L = const; функции H=H (x,t); h=h(x,t); E=E(x,z,t), A=A(x,z,t), F=F(x,z,t) считаются заданными.
Уравнение (1) рассматривается со следующими граничными и начальными условиями:
k -dBcosn,z) + k1 -dBcosn,x) + B=f1, z = h (x,t); dz dx
д B д B
k2—cosn,z) + k2—cos(n,x) + B=f2, z = H(x,t); (2) д дс
B(0, z, t) = B1(z, t); B(L, z, t) = B2(z, t), x=0 и x=L; B(x, z,0) = B3 (x, z), t=0. (3)
Здесь функции B1, B2,B3 , kj(x,t), k2(x,t), f\(x,t), f2(x,t) заданы; cos(n,x), cos(n,z) - проекции вектора внешней нормали к границам физической области в точке границы с координатой x в момент времени t.
" 3 ' B' д Г ,д B1
EA- EA-
д x д x д z
Построение вычислительной области
Задача (1)-(3) рассматривается в сложной области (рис. 1). Решение задачи значительно упростится, если использовать специально построенную расчетную сетку.
Z=h (x,t)
X
Рис. 1. Физическая область
Расчетная область в физической плоскости зависит от времени, поэтому при переходе с одного временного слоя на другой необходимо пересчитывать значения искомой функции на пространственных слоях. Избежать такого пересчета позволяет переход к переменным (хь х1), которые переводят переменную по времени криволинейную физическую область (рис. 1) в постоянную вычислительную область р (рис. 2). 21 = (г-о/(#-С); Х1 = х/Ь. (4)
Физическая область £ < г < Н; 0 < х< Ь переходит в вычислительную 0 < х1 < 1; 0 < <1.
Метрика преобразования (4) определяется как
д zi д x
z-Z
(H -Z)2
д H dZ
д x д x
i z .
H-Zдx '
д xi _ 1 .
д x L
д xi д xi 0 д zi
(5)
д t
д zi
~дГ
д z z-Z
(H -Z)2
д z H-g' (дH ддЛ
д t
дд
д t H-Zдt
J
Вблизи границ происходит резкое изменение концентраций загрязняющих веществ. Чтобы уловить эти изменения, необходима мелкая пространственная сетка. Но считать во всей вычислительной области на такой сетке нельзя из-за ограниченности памяти ЭВМ, поэтому счет проводится на неравномерной разностной сетке. Введем в вычислительной области
i
i
Р неравномерную сетку (хь г1), которая сгущается вблизи границ (рис. 2).
О
1 Х1
* Zi
Рис. 2. Вычислительная область Р с разбиением по пространственным переменным
Для измельчения разностной сетки вблизи границ сделаем замену переменных, которая область Р (рис.
2) переведет в область Я (рис. 3).
О 1Х
Z0
Рис. 3. Вычислительная область Я с равномерной сеткой
В области Я разностная сетка будет равномерной. Замена переменных имеет вид
z0 = ln
X0 = ln
ß1 + 2 z1 -1 ßl - 2 zi -1
ß2 + 2 Xi -1
/ln
ß +1
/ ln
ß1 -1
ß2 +1" ß2 -1_
ß2-2 X1 -1 Обратно
(ß + 1)((ß1 + 1)/(ß1 -1))2 z 0-1 -ß
(6)
z1 = 0,5
1 + ((ß1 + 1)/(ß1 -1))
2 z0-1
x, = 0,5
(ß2 + 1)((ß2 + 1)/(ß2 - 1))2X0-1 -ß2
1 + ((ß2 + 1)/(ß2 -1))
2 X 0 -1
Заметим, что 0 < x0 < 1; 0 < x1 < 1; 0 < z0 < 1; 0 < z1 < 1 .
Метрика преобразований (6)
d x 0
d X,
= 2 ß2 / ( - (2Xi -l)2 ]ln [(ß2 + 1)/(ß2 -1)])
д z0 д z0 д z1 д x
д X1 д z1 д X1 д X1
= 2 ß/ ( - (2z1 -1)2
д z 1
(7)
ln [[ + 1)/(ß1 -1))
я 0 о 0 я
д z _ д z д z1 д t д z1 д t
Параметры Д, /в задаются вычислителем и влияют на густоту разностной сетки вблизи границ рассматриваемой области. Чем они ближе к единице, тем сильнее сгущается сетка вблизи границ.
Таким образом, благодаря указанным преобразованиям координат удалось получить равномерную, не зависящую от времени сетку в вычислительной плоскости (рис. 3), хотя узлы сетки в физической плоскости подвижны и расположены неравномерно (рис. 1, 2).
Переходя в уравнении (1) к координатам x0z0, получим
+
дt дг0
dB дz0 дB
д t = J_
A
+B
д x 0 д B д x дг0
+B
д z 0 д B дx дг0
+B
д z 0 д B д z Cz 0
0 d х d EA (я 0 d x d B 0 d z d B 1
дх 0 cX0 +- dz 0 J
d х d x V d x
д z 0 д
д x dzс
EA
д x 0 д B д z0 д B Л
дx дх0 дX dzC
d£_
д z
V J
Л2
dz0
EA
д B dz 0
+ F.
k
Граничные условия примут вид
дz0 д B
cos(n,z)+k
( dz0 d B дх0 d B Л
dz dz0 +B = f1 при z0 = 1;
dz0 d B
к
dz dz0
C0s(n, z) + к2
дх dz0 dx CX0
( dz0 d B дх0 d B Л
дх dz0 дх дх0
+b = f2 при z0 = 0;
cos(n, х) + (8) cos(n,х)+
(9)
B(0, z0, t) = B1; B(L, z0, t) = B2 при х0 = 0 и х0 = 1; B(x0, z0,0) = B3(x0, z0) при t=0. (10)
Численная схема метода конечных разностей
Используется разностная сетка
{хг°,z0,i = 0,1,...,Nx;у = 0,1,...,N.;к = 0,1,...,Nt; 4 =
= к Д1; х10 = i / Ых; г/ = (I - 0,5) / Ы2; Дt = Д / Ы Дх =1 / Ых; д. =1/N
где Ых, Ыю N - количество узлов сетки; Дх, Дг, Дt -шаги разностной сетки по х0, г0, t.
Для удовлетворения граничных условий по г граничные узлы разностной сетки должны быть расположены на расстояниях 1/(2 Ы) от горизонтальных
границ и введены фиктивные узлы (х0, г N +1),
(х0, гЦ); i, I = 0,1,..., Ых +1.
Учитываются соотношения 0 < х0 < 1; 0 < х1 < 1;
0 < г0 < 1; 0 < г! < 1 0,5С/£/+1/2 -/к}-1/2), имеющие второй порядок аппроксимации по пространственной переменной. Такой подход позволяет удовлетворить условия на границах, не вводя дополнительных граничных условий для метода конечных разностей. Наличие адвективных слагаемых (1)
+
+
+
д
+
/
/
( B д B B д B
(B-; B-) заставляет использовать разностные
д x д z
схемы, которые аппроксимируют эти слагаемые «против потока». В результате получим конечно-разностный аналог задачи (8)-(10), который для простоты выпишем в случае A, E = const.
Уравнение (8) в разностной форме примет вид (i = 1, 2, ..., j = 1, 2, ..., Ni; k = 0, 1, ..., Nt ).
nk+i Bk
Bi,J - B',J + Bk -+ Bi, J
At
д x0 д x
k -nk+8 -nk+8
Bi, j - Bi-i,J
+B
и j
д z0 д x
J
nk+8 nk+8
Bi, J - Bi, J-i
Az
i, J
Ax
д z0 д t
_ fi- (ii)
nk+8 nk+8
Bi, J - Bi, J-i + Bk -+ Bi, J
_ E
Az
д 2 z 0 д x 2
д z 0 д z
1, J
k nk+8 nk+8
Bi, J - Bi, J-i
1, J
k k+8 k+8 Bi, J - Bi, J-i
Az
i, J
Az
д x 0 д x
и J
д z 0
д x
д x 0
k Bk+Y Bk+Y Bk+Y + Bk+Y
Bi+i,J+i - Bi+i,J-i - Bi-i,J+i + Bi-i,J-i
д x
д z 0 д x
3 x
B
i+i, J
4 AzAx
О n/t + M . n/t + M - 2Bi, j + Bi-i, J
Л Гд z0 ]
+ J д x
д x
. . (Ax)2 _ .. . .
i,j v 7 J L Ji,j
k Bk+ß 2Bk+ß + Rk+ß Bi, j+i - 2Bi, j + Bi, j-i +
. . (Az)2
Uj
k r>k+Y T>k+Y T>k+Y . Rk+/ Л Bi+i,J+i - Bi+i,J-i - Bi-i,J+i + Bi-i,J-i
Ь J
( 3 0 Л 2
д z
д z
v J
4 AzAx
Rk+ß - Rk+ß Rk+ß I Bi,M - 2B.,j + Bi,j-1 1+
(Az)
2
h j
Здесь величины а, в, у, 3 постоянны, равны нулю или единице. Условие (8) примет вид (г = 1, 2, ..., Гх )
д z0
д z
cos(n, z) + ki
i,i/2
д z0
dx
cos(n, x)
.,i/ 2
(i2)
Bk+ß Bk+ß Bk+ß Bk+ß + Bk+ß Bk+ß + —;-;-+ —;-;-;-— ki x
д x 0
д x
Az
k
<
i,i/ 2
2Ax
s(n,x) + 0.5^ + Buß)_ fi-
Nx)
Первое условие (9) запишется в виде (i = i, 2, ...,
¿2 cos(«, x)
3 z 0 3 x
Bk+ß _ Bk+ß Bi,Nz+i Bi,Nz
k 2 cos(«, z)
0
о z
"sT
i,Nz+i/2 k
Az
Bk+ß - Bk+ß i,Nz+i П1,М2
Az
(i3)
B
k+ß
+ Bk+ß Bk+ß
+ Bi,Nz - Bi,Nz +
-B
k+ß
k2 cos(w, x)
д x0 д x
2 Ax
k
+0,5(ß + Biki+ß)_ f2.
i, Nz +i/2
При использовании начальных условий (10) и всей разностной схемы следует помнить, что в них входят функции, зависящие от переменных х0, а в постановке задачи функции зависели от переменных х, г.
Заметим, что при а = в = У = 3 = 0 из (11)—(13) получим явную схему метода конечных разностей, при а = в = У = 3 = 1 - неявную. Другие комбинации параметров а, в, У, 3 дают набор полунеявных разностных схем, часть слагаемых в которых берется на текущем временном слое, часть - на следующем. Каждая разностная схема аппроксимирует исходную задачу с порядком О (А/ + Дх + Дг), что доказывается аналогично [6, 7].
Опишем, для примера, алгоритм счета по предложенной разностной схеме в случае а = в = У = 3 = 1.
Алгоритм счета по неявной схеме метода конечных разностей.
1) Задается число узлов сетки по пространственным переменным х0, г0 и шаг по времени А/.
2) Переход от узлов {г)у- ^0; {(х) }}=х0 к узлам
{г/ Г-^; {х° }=0 во всех входных функциях, и в начальных условиях (3) по формулам (4)-(7), расчет сеточных функций
( . „\к
(gz0^k
д z
(g z 0
Ji, j
д x
Ji, J
д 2 z 0
д x'
2
(gx0^k
д x
Ji, J
Ji, J
Nx
3) Разностные уравнения (11) (г, ] = 1, 2, (Г2)), условия (10), (12), (13) (г = 1, 2, Гх) служат для определения функции В. Все слагаемые в этих уравнениях берутся на (к+1)-м слое по времени (а = =в = У = 3 = 1). В результате для определения функции В на каждом слое по времени получается система Гг+2) (Гх+2) уравнений с (Г2+2) (Гх+2) неизвестными, которые решаются одним из численных методов, изложенных, например, в [6].
5) Если не достигнут последний момент времени, то осуществляется переход на следующий шаг по времени (к = к + 1) и возврат на п.3 алгоритма. В противном случае следует перейти на п. 6. При расчете на первом временном слое используются начальные условия задачи.
6) Переход от
у- {fc. {}
1 x
к узлам
{(г)у 0; {(х)г }Гх0 в рассчитанных сеточных функциях.
В случае неявной разностной схемы возникают системы алгебраических уравнений большого порядка (1000 и более уравнений), решение которых сопряжено со значительными вычислительными трудно -стями. Для полунеявных схем получаются разреженные матрицы меньшего размера, например, в случае а = 0; в = 1; У = 3 = 0 задача распадается на подзадачи.
+
x
x
+
k
k
+
k
+
x
+
x
k
+
x
x
+
+
k
+
k
k
ki
+
x
k
+
+
+
На каждом слое по времени и по пространственной переменной х решается система (Ы2+2) уравнений с (Ы2+2) неизвестными. Основная матрица системы является трехдиагональной. Для ее решения можно использовать метод прогонки [7]. Накопления погрешности при вычислении прогоночных коэффициентов в силу диагонального преобладания основной матрицы не происходит. Условия устойчивости выписываются на основе спектрального принципа устойчивости фон Неймана и принципа Куранта-Фридрихса-Леви.
Выводы
Предложен набор разностных схем метода конечных разностей для расчета нелинейного уравнения параболического типа в переменной области криволинейной формы. Выбор используемой для вычислений разностной схемы определяется задачами, стоящими перед исследователем, входными данными, промежутком времени, на котором проводится исследование, математической моделью, в состав которой входят решаемые уравнения, и многими другими факторами. Использование явной схемы позволяет проводить вычисления по простым явным формулам, но
Ростовский государственный университет_
накладывает жесткие ограничения на шаг по времени, неявная разностная схема не накладывает ограничений на шаг по времени, но вычисления по ней связаны со значительными вычислительными сложностями. Использование полунеявной схемы позволяет сгладить недостатки явной и неявной схем, ослабляет ограничения на шаг по времени по сравнению с неявной схемой и не сильно осложняет счет по сравнению с явной.
Литература
1. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислитель-
ная гидромеханика и теплообмен. М., 1990. Т. 1, 2.
2. Бим Р.М., Уорминг Р.Ф. // Ракетная техника и космос.
1978. Т. 16. № 4. С. 145-156.
3. Маккормак Р.В. // Аэрокосмическая техника. 1983. Т. 1.
№ 4. С. 114-123.
4. Роуч П. Вычислительная гидродинамика. М., 1980.
5. Белоцерковский О.М. Численное исследование совре-
менных задач газовой динамики. М., 1974.
6. Марчук Г.И. Методы вычислительной математики. М.,
1989.
7. Самарский А.А. Теория разностных схем. М., 1977.
20 ноября 2006 г.
