close

Вход

Забыли?

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

?

О построении численного метода для задачи Стокса с разрывным коэффициентом вязкости.

код для вставкиСкачать
Вычислительные технологии
Том 14, № 2, 2009
О построении численного метода для задачи Стокса
с разрывным коэффициентом вязкости∗
А. В. Рукавишников
Хабаровское отделение Института прикладной математики ДВО РАН, Россия
e-mail: alexeyruk@mail.ru
Для задачи Стокса с разрывным коэффициентом вязкости применен метод декомпозиции области в сочетании с аппроксимацией задачи на подобластях при
помощи неконформных конечных элементов. Для полученной системы линейных
уравнений построен итерационный метод с блочным переобусловливанием матрицы системы.
Ключевые слова: метод декомпозиции области, неконформные конечные элементы, задача Стокса, разрывные коэффициенты, переобусловливатель.
Введение
Использование неконформного метода конечных элементов [1, 2] на декомпозиционной области в сочетании с мортарными склейками на интерфейсе между подобластями
может быть весьма эффективным при вычислении течения вязкой несжимаемой жидкости с неоднородным (кусочно-постоянным) коэффициентом в эллиптической части
уравнения [3].
Метод мортарных элементов впервые был представлен в [4] для решения уравнения
Пуассона. Основные идеи метода заключены в следующем:
1) в качестве условий согласования на интерфейсе используются условия слабой
непрерывности;
2) необходимо правильно выбирать конечно-элементные пространства на общей границе подобластей.
Об оценках скорости сходимости приближенного решения к точному для задач эллиптического типа можно узнать, например, из работы [5]. Теоретическому исследованию задачи Стокса с однородным (постоянным) коэффициентом кинематической вязкости на декомпозиционной области с нестыкующимися сетками на интерфейсе между
подобластями посвящены работы [6, 7].
В отличие от вышеупомянутых работ в настоящей статье рассмотрен случай разрывного (кусочно-постоянного) коэффициента кинематической вязкости в эллиптической
части уравнения. Исходная область, в нашем случае, разбивается на подобласти так,
Работа выполнена при финансовой поддержке Президиума ДВО РАН (проект № 06-III-А-01-001),
Российского фонда фундаментальных исследований (грант № 07-01-00210-a) и гранта Президента РФ
МК-2092.2007.1.
c ИВТ СО РАН, 2009.
°
∗
110
111
А. В. Рукавишников
чтобы на каждой из них коэффициент вязкости был постоянный. Для такой декомпозиции области предложена новая вариационная постановка задачи, учитывающая
следующие условия согласования на линии разрыва множителя:
1) условие слабой непрерывности компонент вектор-функции скорости;
2) равенство потоков скоростей с давлением на функционалах.
Подход настоящей статьи более общий по сравнению с [6], так как позволяет качественно учитывать особенности решения задачи Стокса на интерфейсе между подобластями в случае разрывного коэффициента кинематической вязкости. В работах [8, 9]
получены априорные оценки скорости сходимости приближенного по методу конечных
элементов решения к точному решению в различных нормах специальных пространств.
В настоящей статье построен и реализован итерационный процесс, состоящий из четырех этапов, решения системы линейных алгебраических уравнений
¸
¸· ¸ ·
·
A B
ω
ζ
=
z
η
BT 0
с переобусловливанием матриц A и S = B T A−1 B.
Предложенную в работе методику, а именно: 1) разбиение области на подобласти,
на которых коэффициент вязкости постоянен; 2) использование сеток на подобластях,
не стыкующихся на общем интерфейсе, в сочетании с мортарными склейками решения
на общей границе с помощью условий слабой непрерывности; 3) построение блочного
переобусловливания, можно перенести на задачу течения двухфазной вязкой несжимаемой жидкости (нелинейные уравнения Навье—Стокса), когда интерфейс между фазами
представляет собой кривую линию [10 – 12].
Работа состоит из четырех разделов и заключения. В разд. 1 обсуждается постановка
задачи Стокса с разрывным коэффициентом кинематической вязкости на декомпозиционной области. Раздел 2 посвящен изложению алгоритма численного решения задачи
Стокса с использованием неконформных конечных элементов и функций склейки на
границе между подобластями (мортарные элементы). Построение итерационного метода решения полученной системы линейных алгебраических уравнений описано в разд. 3.
В разд. 4 приведены результаты тестовых расчетов модельной задачи. В заключении
делаются выводы об эффективности использованного подхода.
1. Постановка задачи
Пусть Ω ⊂ R2 — прямоугольник
Ω = Ω1 ∪ Ω2 = {(x, y) : 0 ≤ x ≤ 2π, 0 ≤ y ≤ π}
с границей Γ =
4
S
Γi , где Γ1 и Γ3 — нижняя и верхняя стороны прямоугольника, а Γ2 и
i=1
Γ4 — правая и левая соответственно. Подобласти Ω1 и Ω2 = Ω\Ω1 являются квадратами,
где Ω1 = {(x, y) : 0 < x < π, 0 < y < π}. Через Γ12 обозначим общую часть границы
(интерфейс) Ω1 и Ω2 .
Требуется найти вектор-функцию скорости w = (u, v) и скалярную функцию давления p, удовлетворяющие следующей системе уравнений и граничных условий (задача
Стокса):
О построении численного метода для задачи Стокса с разрывным...
− div (a0 (x, y) ∇w) + ∇p = f ,
w · n = 0 на Γ,
112
div w = 0 в Ω,
w · τ = gi на Γi ,
Z
p dx dy = 0,
i = 1, 2, 3, 4,
(1)
Ω
где f = (f1 , f2 ) и gi (i = 1, 2, 3, 4) — заданные функции, а n и τ — единичные нормальный
и касательный векторы к границе.
Предполагается, что положительный коэффициент a0 (x, y) кусочно-постоянен:
(
a1 , (x, y) ∈ Ω1 ,
a0 (x, y) =
a2 , (x, y) ∈ Ω2 .
Через V(Ω) обозначим пространство вектор-функций, каждая из компонент которых принадлежит L2 (Ω), их сужение на k-ю подобласть принадлежит H 1 (Ωk ), k =
1, 2, и удовлетворяет граничным условиям (1). Через L2 /R(Ω) обозначим пространство функций из L2 (Ω), ортогональных единице, а через V0 (Ω) определим пространство вектор-функций, каждая компонента которых принадлежит L2 (Ω), их сужение
на каждую подобласть Ωk принадлежит H 1 (Ωk ), k = 1, 2, и удовлетворяет однородным условиям Дирихле на Γ. Кроме этого, введем пространство M(Γ12 ) как множество
вектор-функций, каждая компонента которых принадлежит пространству, дуальному
к H 1/2 (Γ12 ) (относительно L2 (Γ12 )).
Через sk обозначим сужение функции s на k-ю подобласть, а через [s]|Γ12 — функцию разности двух следов (скачок) на общей границе двух подобластей, т. е. [s]|Γ12 :=
s1 |Γ12 ∩Ω1 − s2 |Γ12 ∩Ω2 . Аналогично вводятся обозначения для вектор-функции s : sk и
[s]|Γ12 .
Пусть даны вектор-функции z = (z1 , z2 ) и q = (q1 , q2 ), тогда


∂z1 ∂q1 ∂z1 ∂q1
¶
µ
 ∂x ∂x + ∂y ∂y 
z
q
1
1

∇z : ∇q = 
 ∂z2 ∂q2 ∂z2 ∂q2  и z · q = z2 q2 .
+
∂x ∂x
∂y ∂y
Определим обобщенное решение системы (1) как решение следующей вариационной
задачи: найти тройку (w, p, λ) ∈ V(Ω) × L2 /R(Ω) × M(Γ12 ) такую, что для любых
(φ, ψ, ν) ∈ V0 (Ω) × L2 (Ω) × M(Γ12 ) выполнены соотношения
a(w, φ) + b(φ, p) + d(φ, λ) =< f , φ >,
b(w, ψ) = 0, d(w, ν) = 0,
где билинейные и линейная формы в (2) имеют вид
a(w, φ) =
2
X
ak (w, φ),
k=1
b(w, ψ) =
2
X
k=1
ak (w, φ) =
Z
ak ∇wk : ∇φk dxdy,
Ωk
bk (w, ψ),
bk (w, ψ) =
Z
Ωk
−ψ k div wk dxdy,
(2)
113
А. В. Рукавишников
d(w, ν) =
2
X
dk (w, ν),
dk (w, ν) =
k=1
Z
(−1)k+1 wk |Γ12 ∩Ωk · ν dy,
Γ12
< f , φ >=
2
X
k=1
lk (φ),
lk (φ) =
Z
f k · φk dxdy.
Ωk
Для согласования решения в (2) на Γ12 использованы следующие условия:
1) слабой непрерывности вектора скоростей
Z
(w1 |Γ12 ∩Ω1 − w2 |Γ12 ∩Ω2 ) · µ dy = 0 для любой µ ∈ M(Γ12 );
(3)
Γ12
2) равенства потоков скоростей с давлением на функционалах
Z
Z
∂w2
∂w1
1
−(−a2
(−a1
+ p n1 ) · ϕ dy =
+ p2 n2 ) · ϕ dy для любой ϕ ∈ L2 (Γ12 ), (4)
∂n1
∂n2
Γ12
Γ12
nk — внешняя нормаль к Γ12 — куску границы области Ωk , k = 1, 2.
Для того чтобы замкнуть с помощью равенства (4) полученные на подобластях Ω1
и Ω2 уравнения, в (2) введена вспомогательная вектор-функция λ, определенная из
соотношения
Z
Z
∂w1
+ p1 n1 ) · ϕ dy для любой ϕ ∈ L2 (Γ12 ).
(5)
λ · ϕ dy =
(−a1
∂n1
Γ12
Γ12
2. Численное решение
Выполним триангуляцию Th области Ω. Сначала каждую из подобластей Ω1 и Ω2 вертикальными и горизонтальными линиями разбиваем на элементарные квадраты со стороной 2 h1 и 2 h2 соответственно. Затем каждый из квадратов диагоналями делим на
четыре треугольника, которые обозначим символом e и назовем конечными элементами.
S
(1)
(2)
e,
Далее, введем Th и Th — триангуляции подобластей Ω1 и Ω2 , а через Ωh =
e∈Th
S
S
Ω1h =
e и Ω2h =
e обозначим разбиения Ω, Ω1 и Ω2 соответственно. Отметим,
(1)
e∈Th
(2)
e∈Th
что построенная триангуляция квазиравномерна на подобластях (см. [2]). В качестве
(k)
узлов аппроксимации {aj } для каждой компоненты вектор-функции скорости на Ωk
выберем середины сторон конечных элементов e, совокупность которых обозначим чеS
(k)
рез S . Подмножество узлов, принадлежащих Ωk Γ12 , определим как S (k) , k = 1, 2.
Выбор узлов аппроксимации для скалярной функции давления не принципиален.
Не ограничивая общности, будем полагать h = 2 h1 = 4 h2 ; тогда подмножества
узлов
T
T
аппроксимации S (1) и S (2) для компонент вектора скоростей w на Γ12 Ω1 и Γ12 Ω2
не совпадают, т. е. сетки не стыкуются на интерфейсе Γ12 .
Введем конечно-элементные пространства на Ωkh .
(k)
1. Vh (Ωkh ) — подпространство кусочно-линейных непрерывных на e конечно-элементных функций из пространства L2 (Ωkh ), базисные функции в котором определены
114
О построении численного метода для задачи Стокса с разрывным...
(k)
(k)
естественным образом: базисная функция φi , соответствующая узлу ai ∈ S (k) , рав(k)
на единице и нулю в остальных узлах. При этом φi тождественно равна нулю вне
(k)
конечных элементов, содержащих узел ai .
(k)
2. Xh (Ωkh ) — подпространство кусочно-постоянных конечно-элементных функций
(k)
из пространства L2 (Ωkh ), базисные функции ψi которого равны единице на одном из
треугольников ei и нулю — на остальных.
Сеточное решение системы (2) на подобластях Ωk ищем в виде
X
X
(k)
(k)
ukh (x, y) =
uki φi (x, y), vhk (x, y) =
vik φi (x, y),
(k)
(k)
{i;ai ∈S (k) }
{i;ai ∈S (k) }
X
pkh (x, y) =
(k)
pki ψi (x, y),
k = 1, 2.
(k)
{i;ei ∈Th }
Введем мортарное конечно-элементное пространство на Γ12 . Отрезки
∆0 = [(π, 0), (π, h1 )],
∆j = [(π, h1 + 2 (j − 1) h1 ), (π, h1 + 2 j h1 )],
∆N = [(π, h1 + 2 (N − 1) h1 ), (π, π)],
образуют разбиение Γ12h =
N
S
j = 1, ..., N − 1,
π
N = ∈ N,
h
∆i и называются мортарными конечными элементами
i=0
(см. [4]). В качестве узлов аппроксимации выберем подмножество узлов S (1) на Γ12 :
N −1
S̃ = {b̃j }j=0
= {(π, h1 + 2 j h1 ), j = 0, ..., N − 1}.
Обозначим через Λh (Γ12h ) подпространство линейных на ∆j , j = 1, ..., N −1, и постоянных на ∆0 , ∆N непрерывных конечно-элементных функций из пространства L2 (Γ12h ),
базисные функции νj которого имеют следующий вид:
1) ν0 (b̃0 ) = 1, ν0 (b̃j ) = 0 ∀j = 1, ..., N − 1, ν0 (π, y) = 1 при y ∈ ∆0 , ν0 (π, y)
при y ∈ ∆1 — линейная, ν0 (π, y) = 0 при y ∈ ∆j ∀j = 2, ..., N ;
2) νj (b̃j ) = 1, νj (b̃i ) = 0 ∀i 6= j, νj (π, y) при y ∈ ∆j , νj (π, y) при y ∈ ∆j+1 — линейная,
νj (π, y) = 0 при y ∈ ∆i ∀i 6= j, j + 1, ∀j = 1, ..., N − 2;
3) νN −1 (b̃N −1 ) = 1, νN −1 (b̃j ) = 0 ∀j = 0, ..., N − 2, νN −1 (π, y) = 1 при y ∈ ∆N ,
νN −1 (π, y) при y ∈ ∆N −1 — линейная, νN −1 (π, y) = 0 при y ∈ ∆j ∀j = 0, ..., N − 2.
Функции склейки в пространстве Λh (Γ12h ) ищем в виде
λkh (y)
=
N
−1
X
i=0
λki
νi (π, y),
βhk (y)
=
N
−1
X
βik νi (π, y).
i=0
Как и в работе [13], полагаем λh = λ1h = −λ2h , βh = βh1 = −βh2 , которые являются
сеточными аналогами компонент λ (см. (5)) и определяют конечно-элементную векторфункцию склейки λh = (λh , βh ), а Mh (Γ12h ) = Λh (Γ12h ) × Λh (Γ12h ).
Далее, определим сеточные пространства на Ωh .
1. Для компонент вектор-функции скорости (u и v): пространство Vh (Ωh ) — подпространство функций из пространства L2 (Ωh ), удовлетворяющих краевым условиям задачи (1) в узлах на границе Γ, сужение каждого представителя определяемого простран(k)
ства на k-ю подобласть принадлежит Vh (Ωkh ), k = 1, 2, Vh (Ωh ) = Vh (Ωh ) × Vh (Ωh ).
115
А. В. Рукавишников
2. Для функции давления (p): пространство Xh (Ωh ) — подпространство функций
из пространства L2 (Ωh ), сужение каждого представителя которого на k-ю подобласть
(k)
принадлежит Xh (Ωkh ).
Кроме этого, введем пространство Vh0 (Ωh ) как подпространство функций из пространства L2 (Ωh ), удовлетворяющих однородному условию Дирихле в узлах на грани(k)
це Γ и таких, что их сужение на k-ю подобласть принадлежит Vh (Ωkh ), Vh0 (Ωh ) =
Vh0 (Ωh ) × Vh0 (Ωh ).
Поскольку пространство Vh (Ωh ) не вложено в исходное пространство V(Ω), то, следуя [2], метод нахождения решения задачи (2) будем называть неконформным методом
конечных элементов.
Приближенным решением задачи (1) по методу конечных элементов будем называть
тройку (wh , ph , λh ) ∈ Vh × Xh × Mh , удовлетворяющую системе уравнений
ah (wh , φh ) + bh (φh , ph ) + dh (φh , λh ) =< f , φh >h ,
bh (wh , ψh ) = 0, dh (wh , ν h ) = 0
Z
для любых (φh , ψh , ν h ) ∈ Vh × Xh × Mh и условию ph dxdy = 0. Здесь
(6)
Ω
ah (wh , φh ) =
2
X
akh (wh , φh ) =
k=1
bh (wh , ψh ) =
X Z
ak ∇whk : ∇φkh dxdy,
k=1 e∈T (k) e
h
2
X
bkh (wh , ψh ) =
k=1
dh (φh , λh ) =
2
X
2
X
2
X
X Z
k=1 e∈T (k) e
h
dkh (φh , λh ) =
k=1
< f , φh >h =
−ψhk div whk dxdy,
2 Z
X
k
λh · φkh |Γ12 ∩Ωk dy,
k=1 Γ
12
2
X
k=1
lkh (φh ) =
2 Z
X
k=1 Ω
f k · φkh dxdy.
k
Конечно-элементная задача (6) порождает систему линейных алгебраических уравнений, имеющую седловую матрицу:
¸
·
¸· ¸ ·
A B
ω
ζ
,
(7)
=
z
η
BT 0
где
 (u)


 (u)
B1h
0 D1h 0
A1h
0
0
0




(u)
(u)
0
0 
 0 B2h D2h 0 
 0 A2h
,
B
=
(8)
A=


,
(v)
(v)
 B1h
 0
0
0 D1h 
0 A1h
0 
(v)
(v)
0
B2h
0 D2h
0
0
0 A2h
 (u) 
 1 
 1 
 
F1h
uh
ph
0
 (u) 
 u2 
2 




F


ph  , ω =  2h  , z =  0 
.
ζ =  1h  , η = 
 (v) 



0
 vh 
λh 
 F1h 
0
βh
(v)
vh2
F
2h
О построении численного метода для задачи Стокса с разрывным...
116
3. Построение итерационного процесса
Для нахождения решения системы (7), (8) построим итерационный процесс с переобусловливанием матрицы системы, состоящий из четырех этапов. Отметим, что A — симметричная, положительно определенная матрица, а B — прямоугольная (неквадратная)
матрица полного ранга.
На первом этапе находим решение системы Aζ 0 = ω по формуле
= ζ n0 + αn1 Ã−1 (ω − Aζ n0 ).
ζ n+1
0
(9)
На втором этапе вычисляется решение системы Sη = B T ζ 0 − z:
η n+1 = η n + αn2 S̃ −1 (B T ζ 0 − z − Sη n ).
(10)
Здесь S = B T A−1 B — дополнение по Шуру матрицы системы (7).
Далее, на третьем этапе, ищем поправку τ к ζ 0 с помощью найденного вектора η
как решение системы Aτ = Bη:
τ n+1 = τ n + αn3 Ã−1 (Bη − Aτ n ).
(11)
На последнем этапе вычисляем вектор ζ:
ζ = ζ0 − τ .
Здесь пара векторов (ζ, η) — решение системы (7); αn1 , αn2 , αn3 — параметры процессов
(9)–(11); ζ0 n+1 , η n+1 , τ n+1 и ζ n0 , η n , τ n — значения векторов на (n + 1)-й и n-й итерации
(9)–(11) соответственно; ζ 00 , η 0 , τ 0 — задаваемые начальные приближения процессов
(9)–(11). Ã и S̃ — переобусловливающие матрицы для A и S соответственно.
На каждой итерации для всех из перечисленных этапов, за исключением последнего,
требуется вычислить вектор невязки rn и вектор K̃ −1 rn . Невязку rn для абстрактного
итерационного процесса
qn+1 = qn + αn K̃ −1 rn
(12)
решения системы Kq = χ определим равенством rn = χ−Kqn . Вектор K̃ −1 rn — результат умножения обратной переобусловливающей матрицы системы на вектор невязки rn .
Для экономичного нахождения qn+1 в (12) используем разложение вектора K̃ −1 rn
в m-мерном подпространстве Крылова (см., например, [14, 15]). Отметим зависимость:
чем больше m, тем быстрее будет сходиться итерационный процесс. Приближенное решение ищем с помощью GMRES-метода (см. [14]). Алгоритм проведения одной итерации (12) с использованием переобусловливающей матрицы K̃ имеет следующий вид:
1) пусть q0 — начальное приближение;
2) вычисляем невязку r0 = χ − Kq0 и вектор r01 = K̃ −1 r0 ;
3) находим его евклидову норму ̺ = ||r01 ||;
4) в результате нормирования получаем новый вектор θ 1 = r01 /̺;
5) для j = 1, ..., m;
6) находим вспомогательный вектор ρ: ρ′ = Kθ j , ρ = K̃ −1 ρ′ ;
7) для i = 1, ..., j;
8) вычисляем скалярное произведение векторов σi j = (ρ, θ i );
9) ищем поправку вектора ρ = ρ − σi j θ i ;
10) окончание цикла по i, шаг 7;
117
А. В. Рукавишников
11) вычисляем норму σj+1 j = ||ρ||, если σj+1 j = 0, то выход из цикла по j, шаг 5;
12) иначе, вычисляем новый вектор θ j+1 = ρ/σj+1 j ;
13) окончание цикла по j, шаг 5;
14) в результате получаем матрицу Θm , составленную из векторов θ i , i = 1, ..., m
(столбцов), и матрицу Σm = {σi j }, которая имеет “почти” верхнетреугольный вид, диагональ под главной диагональю также ненулевая, более того, она не содержит нулевых
элементов, кроме, может быть, последнего;
15) находим q1 : q1 = q0 + (ς 1 θ 1 + ... + ς m θ m ), где вектор ς = [ς 1 , ..., ς m ]T :
ς = argmin ||̺ e1 − Σm ς 0 ||, e1 = [1, 0, ..., 0]T , ̺ — норма из шага 3;
16) вычисляем новое значение вектора q0 = q1 и переходим к следующей итерации
(возвращаемся к шагу 1).
Для построения переобусловливающей матрицы Ã будем использовать идею (см.
[15]) неполного (в некотором смысле приближенного) разложения Холецкого матрицы
на треугольные множители ICT (0), т. е. Ã = L̃ · L̃T , где L — верхняя треугольная
матрица.
В силу того, что матрица A блочно-диагональная, Ã имеет вид
  (u)
 (u)

(L̃1 h )T
0
0
0
L̃1 h
0
0
0
 


(u)
(u)
0
(L̃2 h )T
0
0
0
0  
 0 L̃2 h

à = 
 
.
(v) T
(v)
 0

0
0
(L̃1 h )
0
0 L̃1 h 0  
(v) T
(v)
0
0
0
(L̃2 h )
0
0
0 L̃2 h
Следует отметить, что процесс получения матрицы Ã в зарубежной литературе (см.,
например, [16]) называется Row-oriented level-0 процедурой.
В алгоритме 1–16, в его шагах 2 и 6, имеем r0 = Ãr01 , ρ′ = Ãρ. Для того чтобы
вычислить векторы r01 и ρ, необходимо сначала найти векторы r0LT , ρ′LT : r0LT = L̃T r01 ,
ρ′LT = L̃T ρ как решения систем r0 = L̃ r0LT , ρ′ = L̃ ρ′LT . Эти системы линейных алгебраических уравнений легко разрешимы благодаря тому, что L̃ и L̃T соответственно нижняя
и верхняя треугольные матрицы. Таким образом, сначала находим векторы r0LT , ρ′LT :
r0LT = L̃−1 r0 , ρ′LT = L̃−1 ρ′ , а затем — векторы r01 , ρ: r01 = (L̃T )−1 r0LT , ρ = (L̃T )−1 ρ′LT ,
используя построчные исключения сверху вниз и снизу вверх, в первом и во втором
случаях соответственно.
Процесс построения переобусловливающей матрицы S̃ и эффективное решение системы с ней, в шагах 2 и 6 алгоритма 1–16, несколько сложнее. Рассмотрим вспомогательную матрицу Ŝ:
Ŝ = B T Ã−1 B = (B T (L̃T )−1 )(L̃−1 B) = (L̃−1 B)T (L̃−1 B) = X̃ T X̃.
Построим матрицу X̃. Поскольку L̃ — невырожденная квадратная матрица, а B —
прямоугольная матрица полного ранга, то X̃ — прямоугольная матрица, размерность
которой совпадает с размерностью матрицы B. Обозначим через X̃i , Bi и L̃i — i-ю
строку матриц X̃, B и L̃ соответственно. Общий алгоритм вычисления матрицы X̃ =
L̃−1 B заключается в следующем:
1) для i = 1, ..., M;
2) присваиваем i-й строке матрицы X̃ значение i-й строки матрицы B : X̃i = Bi ;
3) для k = 1, ..., i − 1, и ˜li, k 6= 0;
4) выполняем преобразование над i-й строкой матрицы X̃: X̃i = X̃i − ˜li, k · X̃k ;
О построении численного метода для задачи Стокса с разрывным...
118
5) окончание цикла по k, шаг 3;
6) в заключение поделим X̃i на диагональный элемент i-й строки матрицы L̃ : X̃i =
X̃i /˜li, i ;
7) окончание цикла по i, шаг 1.
Оказывается, использование данного алгоритма неприемлемо. Дело в том, что если
первые строки матрицы X̃ не содержат большого числа нулей, то она в результате преобразований алгоритма 1–7 будет заполненной (значения элементов i-й строки передаются элементам j-й строки, если j > i). В таком случае, во-первых, необходимо хранить
в памяти компьютера огромное число элементов матрицы X̃, а во-вторых, потребуются
большие временны́е затраты умножения вектора на эту матрицу, и как следствие —
нецелесообразность использования такого подхода для построения переобусловливающей матрицы системы.
Рассмотрим модернизированный [14] алгоритм построения матрицы X̃, который
преодолевает отмеченные выше трудности. Через lev(mi, j ) определим уровень заполнения (level-of-fill) произвольного элемента матрицы M = {mi, j }:
(
∞, mi, j = 0,
lev(mi, j ) =
0, mi, j 6= 0.
Тогда в процедуре исключения Гаусса, определенной по формуле mi, j = mi, j −mi, k ·mk, j ,
уровень заполнения элемента mi, j вычисляется следующим образом:
lev(mi, j ) = min{lev(mi, j ), lev(mi, k ) + lev(mk, j ) + 1}.
Здесь mi, j представляет собой элементы нижней треугольной матрицы, если i ≥ j, и
верхней треугольной матрицы, если i < j.
Таким образом, ICT (p̄)-преобразование исходной матрицы состоит из процедуры
исключения Гаусса с одной поправкой: все элементы, чей уровень заполнения превосходит p̄, исключаются, т. е. принимают значение нуль. Алгоритм построения матрицы
X̃ = L̃−1 B, если уровень заполнения равен p̄ (Row-oriented Level-p̄ процедура вычисления матрицы X̃), будет следующей:
1) для i = 1, ..., M;
2) присваиваем i-й строке матрицы X̃ значение i-й строки матрицы B : X̃i = Bi ;
3) вычисляем уровень заполнения элементов i-й строки матрицы X̃ = {x̃i, j }:
(
0, bi, j 6= 0,
lev(x̃i, j ) =
∞, bi, j = 0,
где bi, j — элемент i-й строки j-го столбца матрицы B;
4) для k = 1, ..., i − 1, и ˜li, k 6= 0;
5) выполняем преобразование над i-й строкой матрицы X̃: X̃i = X̃i − ˜li, k · X̃k ;
6) находим уровень заполнения каждого элемента i-й строки матрицы X̃ (j = 1, ..., i):
lev(x̃i, j ) = min{lev(x̃i, j ), lev(x̃i, k ) + lev(x̃k, j ) + 1};
7) окончание цикла по k, шаг 4;
8) если lev(x̃i,j ) > p̄, то x̃i, j = 0;
9) в заключение поделим X̃i на диагональный элемент i-й строки матрицы L̃: X̃i =
X̃i /˜li, i ;
10) окончание цикла по i, шаг 1.
119
А. В. Рукавишников
Следует отметить (см. [16]), что если при построении переобусловливающей матрицы к A был выбран алгоритм с уровнем заполнения равным q̄, то при построении переобусловливающей матрицы к S необходимо выбирать алгоритм с уровнем заполнения,
равным 2 q̄ + 1.
В качестве параметра q̄, для того чтобы матрица X̃, а значит, и матрица Ŝ были
разреженными и легко умножались на вектор, необходимо выбрать нуль. Отметим, что
приведенный алгоритм сохраняет блочную структуру матрицы B для матрицы X̃:


(u)
(u)
(L̃1 h )−1 B1 h
(u)
(L̃1 h )−1 D1 h
0
0

(u)
(u)
(u)

0
(L̃2 h )−1 B2 h (L̃2 h )−1 D2 h
0

X̃ =  (v) −1 (v)
(v)
 (L̃1 h ) B1 h
0
0
(L̃1 h )−1 D1 h

0
(v)
(v)
(L̃2 h )−1 B2 h
0
(v)
(L̃2 h )−1 D2 h



.


Прямоугольная матрица X̃ полного ранга не является квадратной, размерность которой совпадает с размерностью матрицы B. В связи с этим следует отметить, что
эффективное решение системы с матрицей Ŝ невозможно в отличие от эффективного
умножения на эту матрицу. Тем не менее матрица Ŝ служит вспомогательной матрицей
при построении переобусловливателя S̃.
В общем случае, пусть для решения системы уравнений x′ = H y′ построена переобусловливающая матрица Ĥ, но она обладает неприятным свойством: эффективное
решение системы с матрицей Ĥ (нахождение обратной — Ĥ −1 ) невозможно. Тогда следует построить такую переобусловливающую матрицу к H — H̃, для которой Ĥ будет
вспомогательной с той особенностью, что ее можно легко умножить на вектор. Процесс решения системы с матрицей H̃ будет выглядить как внутренняя итерационная
процедура.
Пусть требуется найти вектор r∗ — решение системы H̃ r∗ = r̄, тогда внутренняя
итерационная процедура будет следующей:
1) ω 0 = 0;
2) ω l = ω l−1 + β l (r̄ − Ĥ ω l−1 ), l = 1, ..., L;
3) r∗ = ω L ,
где β l — параметр процесса (в [13] — параметр чебышевского ускорения).
В настоящей работе в качестве такой процедуры при получении результата умножения вектора шагов 2 и 6 алгоритма 1–16 на матрицу S̃ −1 предлагается использовать
GMRES-метод с вспомогательной матрицей Ŝ. Внутренняя итерационная процедура
содержит в качестве своих итераций также алгоритм 1–16 с той особенностью, что
K̃ = K̃ −1 = E.
Таким образом, построен итерационный метод решения системы (7) с переобусловливанием и использованием разложения вектора невязки в подпространстве Крылова
(GMRES-метод). В численных экспериментах размерность подпространства Крылова
варьировалась в пределах от 5 до 10. Начальное приближение (ζ 00 , η 0 , τ 0 ) = (0, 0, 0).
4. Результаты численного эксперимента
Рассмотрим тестовый пример задачи (1) с коэффициентом
(
1, (x, y) ∈ Ω1 ,
a0 (x, y) =
a2 , (x, y) ∈ Ω2 .
120
О построении численного метода для задачи Стокса с разрывным...
При этом положительный параметр a2 изменяется таким образом, чтобы разрыв коэффициентов на подобластях увеличивался. Соотношение шагов сеток по каждому
направлению, как и ранее, h1 = 2 h2 , число отрезков разбиения стороны Ωk равно
π
N(k) =
.
2 hk
Разность между точным и приближенным Q
решениями для компонент вектора скоH 1 (e), а для функции давления — в
ростей будем определять в нормах L2 (Ωk ) и
(k)
e∈Th
норме L2 (Ωk ).
Одна из рассмотренных характеристик — поведение погрешности в норме L2 на
полосах вблизи интерфейса Γ12 . Ширина полосы, примыкающей к интерфейсу, равна
π
с каждой стороны. Полосу в k-й подобласти обозначим через Qk .
8
Пример. В качестве решения системы (1) выберем
(
sin x cos y + x cos y,
(x, y) ∈ Ω1 ,
u(x, y) =
−sin x cos y + (2 π − x) cos y, (x, y) ∈ Ω2 ;
(
−cos x sin y − sin y, (x, y) ∈ Ω1 ,
v(x, y) =
cos x sin y + sin y,
(x, y) ∈ Ω2 ;
p(x, y) = cos x sin y,
(x, y) ∈ Ω,
тогда
(
2 sin x cos y + x cos y − sin x sin y,
(x, y) ∈ Ω1 ,
f1 (x, y) =
−a2 (2 sin x cos y − (2π − x) cos y) − sin x sin y, (x, y) ∈ Ω2 ;
(
−2 cos x sin y − sin y + cos x cos y,
(x, y) ∈ Ω1 ,
f2 (x, y) =
a2 (2 cos x sin y + sin y) + cos x cos y, (x, y) ∈ Ω2 ;
(
T
x + sin x,
(x, 0) ∈ Γ1 Ω̄1 ,
g1 (x, 0) =
T
2π − x − sin x, (x, 0) ∈ Γ1 Ω̄2 ;
Т а б л и ц а 1.
a2
10
10
10
100
100
100
0.1
0.1
0.1
0.01
0.01
0.01
N(1)
ku1 − u1h k
ku2 − u2h k
kv 1 − vh1 k
kv 2 − vh2 k
kp1 − p1h k
kp2 − p2h k
16
32
64
16
32
64
16
32
64
16
32
64
0.0023431
0.0005929
0.0001489
0.0024603
0.0006296
0.0001581
0.0010369
0.0002581
0.0000643
0.0009415
0.0002353
0.0000588
0.0002180
0.0000550
0.0000137
0.0002250
0.0000561
0.0000140
0.0018707
0.0004511
0.0001107
0.0019537
0.0004955
0.0001215
0.0010125
0.0001985
0.0000469
0.0009873
0.0001885
0.0000445
0.0021211
0.0004940
0.0001229
0.0022563
0.0005295
0.0001321
0.0003786
0.0000920
0.0000230
0.0005175
0.0001269
0.0000321
0.0011953
0.0003068
0.0000792
0.0013521
0.0003461
0.0000878
0.0223022
0.0115002
0.0059003
0.0226011
0.0117994
0.0059376
0.0205748
0.0105118
0.0054222
0.0202156
0.0103172
0.0053926
0.0118944
0.0062433
0.0033015
0.0123778
0.0063928
0.0033667
0.0129343
0.0067412
0.0034766
0.0135580
0.0067624
0.0034815
121
А. В. Рукавишников
Т а б л и ц а 2.
a2
10
10
10
100
100
100
0.1
0.1
0.1
0.01
0.01
0.01
N(1)
ku1 − u1h k
ku2 − u2h k
kv 1 − vh1 k
kv 2 − vh2 k
16
32
64
16
32
64
16
32
64
16
32
64
0.0506119
0.0269168
0.0140726
0.0521878
0.0277154
0.0144291
0.0420919
0.0210421
0.0105279
0.0419619
0.0209714
0.0104847
0.0219283
0.0110970
0.0056036
0.0210023
0.0104988
0.0052493
0.0767802
0.0402733
0.0216351
0.0794782
0.0435169
0.0236703
0.0314110
0.0142037
0.0066196
0.0318380
0.0143119
0.0066740
0.0321798
0.0146107
0.0068550
0.0324776
0.0147743
0.0069455
0.0142226
0.0066240
0.0031739
0.0147205
0.0068892
0.0033224
0.0246155
0.0121144
0.0060792
0.0266743
0.0131115
0.0065212
Т а б л и ц а 3.
a2
10
10
10
100
100
100
0.1
0.1
0.1
0.01
0.01
0.01
N(1)
ku1 − u1h k
ku2 − u2h k
kv 1 − vh1 k
kv 2 − vh2 k
16
32
64
16
32
64
16
32
64
16
32
64
0.0018073
0.0004592
0.0001155
0.0019234
0.0004946
0.0001242
0.0004027
0.0000988
0.0000244
0.0002491
0.0000624
0.0000155
0.0001034
0.0000265
0.0000066
0.0000947
0.0000190
0.0000067
0.0016287
0.0003758
0.0000897
0.0016974
0.0004116
0.0000983
0.0006844
0.0001374
0.0000306
0.0006879
0.0001491
0.0000345
0.0012020
0.0002999
0.0000762
0.0013127
0.0003281
0.0000834
0.0001755
0.0000397
0.0000094
0.0002778
0.0000637
0.0000154
0.0008459
0.0002023
0.0000500
0.0009524
0.0002272
0.0000553
(
T
−x − sin x,
(x, π) ∈ Γ3 Ω̄1 ,
g3 (x, π) =
T
x − 2π + sin x, (x, π) ∈ Γ3 Ω̄2 ;
g2 (2 π, y) = 2 sin y,
u(x, y) = 0,
(2 π, y) ∈ Γ2 ;
(x, y) ∈ Γ2 ∪ Γ4 ;
g4 (0, y) = −2 sin y,
v(x, y) = 0,
(0, y) ∈ Γ4 ;
(x, y) ∈ Γ1 ∪ Γ3 .
Для данной задачи определим значения основных параметров. Имеем (N(1) , N(2) ) =
(16, 32); (32, 64); (64, 128); a2 = 10; 100; 0.1; 0.01. Результаты расчетов приведены в
табл. 1–3: в табл. 1 — погрешность решения — компонент вектора скорости w = (u, v)
и функции давления p, на подобластях Ωk , в норме L2 (Ωk ); в табл. 2 — погрешность
Q
H 1 (e); в
решения — компонент вектора скоростей w, на подобластях Ωk , в норме
(k)
e∈Th
табл. 3 — на полосах Qk , в норме L2 (Qk ).
О построении численного метода для задачи Стокса с разрывным...
122
Заключение
Расчеты численного эксперимента показали, что погрешность не слишком сильно растет при увеличении
разрыва коэффициентов на подобластях в нормах пространств
Q
1
L2 (Ωk ) и
H (e), k = 1, 2. Большая часть погрешности (в норме пространства L2 )
(k)
e∈Th
сосредоточена вдоль полос Qk в окрестности интерфейса Γ12 между подобластями, т. е.
погрешность не распространяется во внутреннюю часть подобластей.
Легко видеть, что погрешностьQрешения для компонент вектор-функции скорости
H 1 (e) и скалярной функции давления p в норме
w = (u, v) в норме пространства
(k)
e∈Th
пространства L2 убывает пропорционально степени h, т. е. является, как и ожидалось,
величиной порядка O(h), а погрешность решения вектор-функции скорости w = (u, v)
в норме пространства L2 — пропорционально h2 , т. е. имеет порядок O(h2 ), что согласуется с априорными оценками скорости сходимости (см. [8, 9]).
Список литературы
[1] Crouzeix M., Raviart P.A. Comforming and non-conforming finite element methods for
solving the stationary Stokes equations // RAIRO Anal. Numer. 1973. Vol. 7. P. 33–76.
[2] Brezzi F., Fortin M. Mixed and Hybrid Finite Element Methods. N.Y.: Springer-Verlag,
1991.
[3] Рукавишников А.В. Численный метод решения задачи Стокса с разрывным коэффициентом // Вычисл. методы и программирование. 2005. Т. 6, № 1. C. 21–30.
[4] Bernardi C., Maday Y., Patera A. A New nonconforming approach to domain
decomposition: the mortar element method // Nonlinear Partial Differential Equations and
their Applications. P.: Pitman, 1989. P. 13–51.
[5] Braess D., Dahmen W., Wieners C. A multigrid algorithm for the mortar finite element
methods // SIAM J. Numer. Anal. 1999. Vol. 37, N 1. P. 48–69.
[6] Ben Belgacem F. The mixed mortar finite element method for the incompressible Stokes
problem: convergence analysis // SIAM J. Numer. Anal. 2000. Vol. 37, N 4. P. 1085–1100.
[7] Ben Belgacem F. A stabilized domain decomposition method with nonmatching grids for the
Stokes problem in three dimensions // SIAM J. Numer. Anal. 2004. Vol. 42, N 2. P. 667–685.
[8] Рукавишников А.В., Рукавишников В.А. Неконформный метод конечных элементов для задачи Стокса с разрывным коэффициентом // Сиб. журн. индустр. математики.
2007. Т. 10, № 4(32). C. 104–117.
[9] Рукавишников А.В. Об оценке точности в норме пространства L2 (Ω) // Тр. 5-й Международной научной конференции творческой молодежи, Хабаровск: Изд-во ДВГУПС,
2007. Т. 4(6). С. 105–109.
[10] Рукавишников А.В. Обобщенная постановка задачи течения двухфазной жидкости
с непрерывно изменяющимся интерфейсом // Матем. моделирование. 2008. Т. 20, № 3.
C. 3–8.
[11] Flemisch B., Melenk J., Wohlmuth B. Mortar methods with curved interfaces //
APNUM. 2005. Vol. 54. P. 339–361.
123
А. В. Рукавишников
[12] Рукавишников А.В. Построение численного метода решения задачи течения двухфазной жидкости // Инновационные технологии — транспорту и промышленности: Тр. 45-й
Международной научно-практической конференции ученых транспортных вузов, инженерных работников и представителей академической науки. Хабаровск: Изд-во ДВГУПС,
2007. Т 3(6). С. 53–58.
[13] Kuznetsov Yu.A. Efficient iterative solvers for elliptic finite elements problems on
nonmatching grids // Russ. J. Numer. Anal. Math. Modelling. 1995. Vol. 10, N 3.
P. 187–211.
[14] Saad Y. Iterative Methods for Sparse Linear Systems. New Jersey: PWS, 1996.
[15] Ильин В.П. Методы конечных разностей и конечных объемов для эллиптических уравнений. Новосибирск: Ин-т математики СО РАН, 2000.
[16] Little L., Saad Y. Block LU Preconditioners for Symmetric and Nonsymmetric Saddle
Point Problems. Minnesota: Minnesota Supercomputing Inst., 1999.
Поступила в редакцию 6 марта 2008 г.
Документ
Категория
Без категории
Просмотров
7
Размер файла
232 Кб
Теги
построение, метод, вязкости, стокса, коэффициента, разрывных, задачи, численного
1/--страниц
Пожаловаться на содержимое документа