О построении численного метода для задачи Стокса с разрывным коэффициентом вязкости Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 14, № 2, 2009

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

А. В. Руклвишников Хабаровское отделение Института прикладной математики ДВО РАН, Россия

e-mail: [email protected]

Для задачи Стокса с разрывным коэффициентом вязкости применен метод декомпозиции области в сочетании с аппроксимацией задачи на подобластях при помощи неконформных конечных элементов. Для полученной системы линейных уравнений построен итерационный метод с блочным переобусловливанием матрицы системы.

Ключевые слова: метод декомпозиции области, неконформные конечные элементы, задача Стокса, разрывные коэффициенты, переобусловливатель.

Введение

Использование неконформного метода конечных элементов [1, 2] на декомпозиционной области в сочетании с мортарными склейками на интерфейсе между подобластями может быть весьма эффективным при вычислении течения вязкой несжимаемой жидкости с неоднородным (кусочно-постоянным) коэффициентом в эллиптической части уравнения [3].

Метод мортарных элементов впервые был представлен в [4] для решения уравнения Пуассона. Основные идеи метода заключены в следующем:

1) в качестве условий согласования на интерфейсе используются условия слабой непрерывности;

2) необходимо правильно выбирать конечно-элементные пространства на общей границе подобластей.

Об оценках скорости сходимости приближенного решения к точному для задач эллиптического типа можно узнать, например, из работы [5]. Теоретическому исследованию задачи Стокса с однородным (постоянным) коэффициентом кинематической вязкости на декомпозиционной области с нестыкующимися сетками на интерфейсе между подобластями посвящены работы [6, 7].

В отличие от вышеупомянутых работ в настоящей статье рассмотрен случай разрывного (кусочно-постоянного) коэффициента кинематической вязкости в эллиптической части уравнения. Исходная область, в нашем случае, разбивается на подобласти так,

*Работа выполнена при финансовой поддержке Президиума ДВО РАН (проект № 06-111-А-01-001), Российского фонда фундаментальных исследований (грант № 07-01-00210-а) и гранта Президента РФ МК-2092.2007.1.

© ИВТ СО РАН, 2009.

чтобы на каждой из них коэффициент вязкости был постоянный. Для такой декомпозиции области предложена новая вариационная постановка задачи, учитывающая следующие условия согласования на линии разрыва множителя:

1) условие слабой непрерывности компонент вектор-функции скорости;

2) равенство потоков скоростей с давлением на функционалах.

Подход настоящей статьи более общий по сравнению с [6], так как позволяет качественно учитывать особенности решения задачи Стокса на интерфейсе между подобластями в случае разрывного коэффициента кинематической вязкости. В работах [8, 9] получены априорные оценки скорости сходимости приближенного по методу конечных элементов решения к точному решению в различных нормах специальных пространств. В настоящей статье построен и реализован итерационный процесс, состоящий из четырех этапов, решения системы линейных алгебраических уравнений

А В " С "

Вт 0 ъ

с переобусловливанием матриц А и 5 = Вт А-1 В.

Предложенную в работе методику, а именно: 1) разбиение области на подобласти, на которых коэффициент вязкости постоянен; 2) использование сеток на подобластях, не стыкующихся на общем интерфейсе, в сочетании с мортарными склейками решения на общей границе с помощью условий слабой непрерывности; 3) построение блочного переобусловливания, можно перенести на задачу течения двухфазной вязкой несжимаемой жидкости (нелинейные уравнения Навье—Стокса), когда интерфейс между фазами представляет собой кривую линию [10-12].

Работа состоит из четырех разделов и заключения. В разд. 1 обсуждается постановка задачи Стокса с разрывным коэффициентом кинематической вязкости на декомпозиционной области. Раздел 2 посвящен изложению алгоритма численного решения задачи Стокса с использованием неконформных конечных элементов и функций склейки на границе между подобластями (мортарные элементы). Построение итерационного метода решения полученной системы линейных алгебраических уравнений описано в разд. 3. В разд. 4 приведены результаты тестовых расчетов модельной задачи. В заключении делаются выводы об эффективности использованного подхода.

1. Постановка задачи

Пусть П С Я2 — прямоугольник

П = П и П2 = {(х, у) : 0 < х < 2п, 0 < у < п}

4

с границей Г = У Г^, где Г1 и Г3 — нижняя и верхняя стороны прямоугольника, а Г2 и

г=1 _

Г4 — правая и левая соответственно. Подобласти П1 и П2 = П\П1 являются квадратами, где П1 = {(х, у): 0 <х<п, 0 < у < п}. Через Г12 обозначим общую часть границы (интерфейс) П1 и П2.

Требуется найти вектор-функцию скорости w = (и, у) и скалярную функцию давления р, удовлетворяющие следующей системе уравнений и граничных условий (задача Стокса):

— ^у (а0(х, у) У^ + Ур = f, ^у w = 0 в П, w • п = 0 на Г, w • т = дг на Гг, г =1, 2, 3, 4,

/ рйхйу = 0,

(1)

где f = (/1, /2) и дг (г = 1, 2, 3, 4) — заданные функции, а п и т — единичные нормальный и касательный векторы к границе.

Предполагается, что положительный коэффициент а0(х,у) кусочно-постоянен:

ао(х, у)

ах, (х, у) € Пь а2, (х, у) € П2.

Через У(П) обозначим пространство вектор-функций, каждая из компонент которых принадлежит ¿2(П), их сужение на к-ю подобласть принадлежит Н1 (Пк), к = 1, 2, и удовлетворяет граничным условиям (1). Через Ь2/Л(П) обозначим пространство функций из Ь2(П), ортогональных единице, а через V0(П) определим пространство вектор-функций, каждая компонента которых принадлежит Ь2(П), их сужение на каждую подобласть Пк принадлежит Н1 (Пк), к = 1, 2, и удовлетворяет однородным условиям Дирихле на Г. Кроме этого, введем пространство М(Г12) как множество вектор-функций, каждая компонента которых принадлежит пространству, дуальному к Н1/2 (Г 12) (относительно ^(Г^)).

Через вк обозначим сужение функции в на к-ю подобласть, а через [в]|г12 — функцию разности двух следов (скачок) на общей границе двух подобластей, т.е. [в]|г12 := в1|г12Пй1 — в2|г12пп2. Аналогично вводятся обозначения для вектор-функции ! [й]|Г12 .

Б и

Пусть даны вектор-функции ъ = (г1, г2) и q = (д1, д2), тогда

( 5^1 5^1 + 5^1 5^1 \

Уъ : Уq

дх дх 5^2 5^2

+

ду ду д<?2

и ъ • q

дх дх ду ду

Определим обобщенное решение системы (1) как решение следующей вариационной задачи: найти тройку Л) € V(П) х Ь2/Д(П) х М(Г12) такую, что для любых

(ф, ф, V) € V0(П) х Ь2(П) х М(Г12) выполнены соотношения

a(w, ф) + 6(ф, р) + ^(ф, Л) =< f, ф >, 6^, ф) = 0, V) = 0,

где билинейные и линейная формы в (2) имеют вид

a(w, ф) = ^^ак(w, ф), ак(w, ф)= ак Уwk : Уфк ^х^у,

(2)

к=1

Пк

6(w, ф) = ^^6к(w, ф), 6к(w, ф) = J —фк diу wk ^х^у,

Пк

к=1

2

¿К V) = ^4К V), 4К V) = I (-1)к+1 wk|г12пйк ■ V ¿у, к= Г12

< {, ф > = ^ 1к (ф), 1к (ф) = I {к ■ фк ¿хйу. к=1

0,к

Для согласования решения в (2) на Г12 использованы следующие условия:

1) слабой непрерывности вектора скоростей

У Мг^п^ - ^г12пп2) ■ Д¿У = 0 длялюбой Д € М(Г12); (3)

Г12

2) равенства потоков скоростей с давлением на функционалах

/д С д^^2

(-«1 дП" + ^1) ■ ^ ¿у = ( «2 "дП" + Р2П2) ■ ^ ¿у длялюбой ^ € Ь2(Г12), (4)

Г12 Г12

Пк — внешняя нормаль к Г12 — куску границы области Пк, к = 1, 2.

Для того чтобы замкнуть с помощью равенства (4) полученные на подобластях и П2 уравнения, в (2) введена вспомогательная вектор-функция Л, определенная из соотношения

__г дwl

Л ■ ^ ¿у = (-Й1 "дп- + Р1П1) ■ ^ ¿у длялюбой ^ € Ь"(Г12). (5)

Г12 Г12

2. Численное решение

Выполним триангуляцию Т^ области П. Сначала каждую из подобластей и П2 вертикальными и горизонтальными линиями разбиваем на элементарные квадраты со стороной 2 к1 и 2 к2 соответственно. Затем каждый из квадратов диагоналями делим на четыре треугольника, которые обозначим символом е и назовем конечными элементами. Далее, введем Т^ и Т^2) — триангуляции подобластей П1 и П2, а через П = и е,

е&Тк

П1н = и е и П2ь = и е обозначим разбиения П, П1 и П2 соответственно. Отметим, вет« вет®

что построенная триангуляция квазиравномерна на подобластях (см. [2]). В качестве узлов аппроксимации {«к)} для каждой компоненты вектор-функции скорости на Пк выберем середины сторон конечных элементов е, совокупность которых обозначим через Б( к. Подмножество узлов, принадлежащих Пк и Г12, определим как Б(к), к = 1, 2. Выбор узлов аппроксимации для скалярной функции давления не принципиален.

Не ограничивая общности, будем полагать к = 2 к1 = 4 к2; тогда подмножества узлов аппроксимации Б(1) и Б(2) для компонент вектора скоростей w на Г12 Р| П1 и Г12 Р| П2 не совпадают, т.е. сетки не стыкуются на интерфейсе Г12. Введем конечно-элементные пространства на П^.

1. ^(к)(П^) — подпространство кусочно-линейных непрерывных на е конечно-элементных функций из пространства Ь2(Пкъ), базисные функции в котором определены

естественным образом: базисная функция ф(к), соответствующая узлу а(к) € равна единице и нулю в остальных узлах. При этом ф(к) тождественно равна нулю вне

(к)

конечных элементов, содержащих узел а .

2. Х^к)(П^) — подпространство кусочно-постоянных конечно-элементных функций из пространства базисные функции которого равны единице на одном из

треугольников е и нулю — на остальных.

Сеточное решение системы (2) на подобластях П ищем в виде

«£(*,у)= £ ик #(х,у), ^(х,у) = £ и* ф(к)(х,у),

{г;а(к)€5(к)} {г;а(к) ей1« }

к „/,(к)/

РЙ(х,у)= £ рк Л^у), к =1, 2.

{ад ет«}

Введем мортарное конечно-элементное пространство на Г12. Отрезки

Ас = [(п, 0), (п, М,

А,- = [(п, Д + 2 (э - 1) (п, Д + 2э &!)], э = 1,..., N - 1,

п

Ам = [(п, Д + 2 ^ - 1) (п, п)], N = - € N

N

образуют разбиение Г12^ = и А^ и называются мортарными конечными элементами

г=С

(см. [4]). В качестве узлов аппроксимации выберем подмножество узлов 5(1) на Г12: 5 = {6,Ы-1 = {(п, Д + 2эЛх), э = 0,..., N - 1}.

Обозначим через Л^(Г12^) подпространство линейных на А,, э = 1,..., N - 1, и постоянных на Ас, АN непрерывных конечно-элементных функций из пространства Ь2(Г12^), базисные функции V, которого имеют следующий вид:

1) ^(6о) = 1, ^(6,-) = 0 У? = 1,..., N - 1, ^(п, у) = 1 при у € Ас, ^(п, у) при у € А1 — линейная, ^(п, у) = 0 при у € А, У; = 2,..., N;

2) V, (6,) = 1, V, (6^) = 0 Уг = э, V, (п, у) при у € А,, V, (п, у) при у € А,+1 — линейная, V,(п, у) = 0 при у € Аг Уг = + 1, У; = 1,..., N - 2;

3) VN-l(bN-1) = 1, VN-l(6j) = 0 Уэ = 0,..., N - 2, VN-l(п, у) = 1 при у € АN, vN-1(п, у) при у € АN-1 — линейная, vN-1(п, у) = 0 при у € А, У? = 0,..., N - 2.

Функции склейки в пространстве Л^(Г12^) ищем в виде

N-1 N-1

г=С г=С

АЙ(у) = £ ^(п,у), $(у) = £ вк ^(п,у).

Как и в работе [13], полагаем А^ = А^ = -А|, = вЙ = -вЙ, которые являются сеточными аналогами компонент Л (см. (5)) и определяют конечно-элементную вектор-функцию склейки Лн = (ай, вд а М^(Гт) = Л^(Гт) х Л^(Гт).

Далее, определим сеточные пространства на П^.

1. Для компонент вектор-функции скорости (и и «): пространство ^(П^) — подпространство функций из пространства Ь2 (П^), удовлетворяющих краевым условиям задачи (1) в узлах на границе Г, сужение каждого представителя определяемого пространства на к-ю подобласть принадлежит (ПкЙ), к =1, 2, ) = ) х ^(ПЙ).

2. Для функции давления (р): пространство Хн(Пн) — подпространство функций из пространства ¿2(Пн), сужение каждого представителя которого на к-ю подобласть принадлежит X^ (Пкн)•

Кроме этого, введем пространство УЦ (Пн) как подпространство функций из пространства ¿2(Пн), удовлетворяющих однородному условию Дирихле в узлах на границе Г и таких, что их сужение на к-ю подобласть принадлежит Ун(к)(Пкн), V00(Пн) = У0(Пн) х у0(Пн).

Поскольку пространство Vн(Пн) не вложено в исходное пространство У'(П), то, следуя [2], метод нахождения решения задачи (2) будем называть неконформным методом конечных элементов.

Приближенным решением задачи (1) по методу конечных элементов будем называть тройку Ан) € Vн х Хн х Мн, удовлетворяющую системе уравнений

ан(^н, Фн) + Ьн(фн,Рн) + 4(фн, Ан) =< f, Фн >н, bн(wн,^н) = 0, Vн ) = 0

для любых (фн, фн, ин) € Vн х Хн х Мн и условие У рн ¿хйу = 0. Здесь

п

2 2 „

ан(^ фн) = X] akH(wH, фн) = X X] ак : Уфн dxdУ,

к=1 _1 (к) "

(6)

к=1 ееТ(к)

2 2 „ ЬнМн) = 52 Ькн^н,фн) = x x /

к=1 к=1 еет(к) е

—фк Wk ¿хйу,

Пк ^,

4(фн Ан) = X Ан) = X Ан • фкн\г12ппк

к= 1 к=1-р

2 2 л

< f, фн >н= X 1кн(фн) = X / fk • фн ^хйу.

к=1 к=1Пк

Конечно-элементная задача (6) порождает систему линейных алгебраических уравнений, имеющую седловую матрицу:

А В " С "

В т 0 2

где

А

А(и) А1н 0 0

0 А(и) А2н 0

0 0 А (V) А1н

0 0 0

2н

В

В

(и)

1н

В

(V)

1н

0 Ан 0

В2н ^2н 0

0 0 Ан

В (V) В2н 0 Ан

(7)

(8)

" ин ' " рн"

ин ^н , п = Рн Ан

. вн _

(и)

1н

р(и)

р2н

Р (V)

Р (V) р2н

0 0 0 0

0

0

С

2

3. Построение итерационного процесса

Для нахождения решения системы (7), (8) построим итерационный процесс с переобусловливанием матрицы системы, состоящий из четырех этапов. Отметим, что А — симметричная, положительно определенная матрица, а В — прямоугольная (неквадратная) матрица полного ранга.

На первом этапе находим решение системы А£с = ш по формуле

С+1 = СП + аП А-1(ш - АСП). (9)

На втором этапе вычисляется решение системы Бп = Вт£с - 2:

ПП+1 = Пп + аП 5-1(ВтСс - 2 - (10)

Здесь Б = ВтА-1В — дополнение по Шуру матрицы системы (7).

Далее, на третьем этапе, ищем поправку т к £с с помощью найденного вектора п как решение системы Ат = Вп:

тП+1 = тП + аП А-1(Вп - АтП). (11)

На последнем этапе вычисляем вектор £:

С = С с - т.

Здесь пара векторов (£, п) — решение системы (7); аП, аП, аП — параметры процессов (9)-(11); СоП+1, п"+1, тга+1 и СП, п", тП — значения векторов на (п + 1)-й и п-й итерации (9)—(11) соответственно; £пс, тс — задаваемые начальные приближения процессов (9)—(11). А и Б — переобусловливающие матрицы для А и Б соответственно.

На каждой итерации для всех из перечисленных этапов, за исключением последнего, требуется вычислить вектор невязки гП и вектор К-1 гП. Невязку гП для абстрактного итерационного процесса

ЯП+1 = яП + а„ К-1 гП (12)

решения системы Кq = х определим равенством гП = Х-Кqra. Вектор К-1 гП — результат умножения обратной переобусловливающей матрицы системы на вектор невязки гП.

Для экономичного нахождения qra+1 в (12) используем разложение вектора К-1 гП в т-мерном подпространстве Крылова (см., например, [14, 15]). Отметим зависимость: чем больше т, тем быстрее будет сходиться итерационный процесс. Приближенное решение ищем с помощью ОМИЕВ-метода (см. [14]). Алгоритм проведения одной итерации (12) с использованием переобусловливающей матрицы К имеет следующий вид:

1) пусть qс — начальное приближение;

2) вычисляем невязку гс = х - Кqс и вектор г^ = К-1 гс;

3) находим его евклидову норму д = ||гЦ|;

4) в результате нормирования получаем новый вектор 01 = г^/д;

5) для э = 1,..., т;

6) находим вспомогательный вектор р: р' = К0,, р = К-1 р';

7) для г = 1,...,э;

8) вычисляем скалярное произведение векторов = (р, 0^);

9) ищем поправку вектора р = р - 0^;

10) окончание цикла по г, шаг 7;

11) вычисляем норму а7+1 = ||р||, если а7+1 j = 0, то выход из цикла по шаг 5;

12) иначе, вычисляем новый вектор = р/а7+1;

13) окончание цикла по шаг 5;

14) в результате получаем матрицу вт, составленную из векторов г = 1,...,т (столбцов), и матрицу = (ст^-}, которая имеет "почти" верхнетреугольный вид, диагональ под главной диагональю также ненулевая, более того, она не содержит нулевых элементов, кроме, может быть, последнего;

15) находим q1: q1 = q0 + (с101 + ... + где вектор ^ = [с ]т: ^ = ащшт ||д е1 — ^0||, е1 = [1,0,..., 0]Т, д — норма из шага 3;

16) вычисляем новое значение вектора q0 = q1 и переходим к следующей итерации (возвращаемся к шагу 1).

Для построения переобусловливающей матрицы А будем использовать идею (см. [15]) неполного (в некотором смысле приближенного) разложения Холецкого матрицы на треугольные множители ЮТ(0), т.е. А = Ь ■ Ьт, где Ь — верхняя треугольная матрица.

В силу того, что матрица А блочно-диагональная, А имеет вид

r(«) L1 h 0 0 0 " (L(1ih))T 0 0 0

0 ?(«) L2 h 0 0 0 f f (u)\ T (L2 h ) 0 0

0 0 L1 h 0 0 0 (Lfh)T 0

0 0 0 L(v) L2 h . 0 0 0 f f (v)\ T (L2 h )

Следует отметить, что процесс получения матрицы A в зарубежной литературе (см., например, [16]) называется Row-oriented level-0 процедурой.

В алгоритме 1-16, в его шагах 2 и 6, имеем r0 = Ar°, р' = Ар. Для того чтобы вычислить векторы r^ и р, необходимо сначала найти векторы г^т, р'ьт: г^т = LT r°, p'LT = LTр как решения систем r° = L г^т , р' = L р'^т. Эти системы линейных алгебраических уравнений легко разрешимы благодаря тому, что L и LT соответственно нижняя и верхняя треугольные матрицы. Таким образом, сначала находим векторы rLT, р^т: rLT = L-1r°, р'^ = L-V, а затем — векторы r°, р: r° = (LT)-1rLT, р = (LT)-1р'1^т, используя построчные исключения сверху вниз и снизу вверх, в первом и во втором случаях соответственно.

Процесс построения переобусловливающей матрицы S и эффективное решение системы с ней, в шагах 2 и 6 алгоритма 1-16, несколько сложнее. Рассмотрим вспомогательную матрицу S:

S = вт А-1 в = (вT (LT )-1)(L-1 в) = (L-1 b)t (L-1 в) = XT X.

Построим матрицу X. Поскольку L — невырожденная квадратная матрица, а B — прямоугольная матрица полного ранга, то X — прямоугольная матрица, размерность которой совпадает с размерностью матрицы B. Обозначим через Xj, Bj и Lj — i-ю строку матриц X, в и L соответственно. Общий алгоритм вычисления матрицы X = L-1B заключается в следующем:

1) для i = 1,..., M;

2) присваиваем i-й строке матрицы X значение i-й строки матрицы B : Xj = Bj;

3) для k = 1,..., i — 1, и I= 0;

4) выполняем преобразование над i-й строкой матрицы X: Xj = Xj — • X^;

5) окончание цикла по к, шаг 3;

6) в заключение поделим Х^ на диагональный элемент г-й строки матрицы Ь : Х^ =

7) окончание цикла по г, шаг 1.

Оказывается, использование данного алгоритма неприемлемо. Дело в том, что если первые строки матрицы X не содержат большого числа нулей, то она в результате преобразований алгоритма 1-7 будет заполненной (значения элементов г-й строки передаются элементам ]-й строки, если ] > г). В таком случае, во-первых, необходимо хранить в памяти компьютера огромное число элементов матрицы X, а во-вторых, потребуются большие временные затраты умножения вектора на эту матрицу, и как следствие — нецелесообразность использования такого подхода для построения переобусловливающей матрицы системы.

Рассмотрим модернизированный [14] алгоритм построения матрицы X, который преодолевает отмеченные выше трудности. Через 1еу(тг;^) определим уровень заполнения (1еуе1-оГ-Ш1) произвольного элемента матрицы М = {т^-}:

Тогда в процедуре исключения Гаусса, определенной по формуле т^- = ^—тъ,к, уровень заполнения элемента т^- вычисляется следующим образом:

Здесь представляет собой элементы нижней треугольной матрицы, если i > j, и верхней треугольной матрицы, если i < j.

Таким образом, ICT(p)-преобразование исходной матрицы состоит из процедуры исключения Гаусса с одной поправкой: все элементы, чей уровень заполнения превосходит p, исключаются, т. е. принимают значение нуль. Алгоритм построения матрицы X = L-1B, если уровень заполнения равен p (Row-oriented Level-p процедура вычисления матрицы X), будет следующей:

1) для i = 1,..., M;

2) присваиваем i-й строке матрицы X значение i-й строки матрицы B : Xj = B^;

3) вычисляем уровень заполнения элементов i-й строки матрицы X = {Xi;j}:

где Ь^- — элемент г-й строки ^-го столбца матрицы В;

4) для к = 1,...,г — 1, и /ъ,к = 0;

5) выполняем преобразование над г-й строкой матрицы X: XI = XI — /¿,к • Xk;

6) находим уровень заполнения каждого элемента г-й строки матрицы X (^ = 1,..., г): 1еу(Х^-) = шт{1еу(Х^), 1еу(Х^к) + 1еу(Хк,^) + 1};

7) окончание цикла по к, шаг 4;

8) если 1еу(Х^-) > р, то Хъ, ^ = 0;

9) в заключение поделим XI на диагональный элемент г-й строки матрицы Ь: XI =

= 0, = 0.

lev(mi;j) = min{lev(mi;j), lev(mi)k) + lev(mk)j-) + 1}.

10) окончание цикла по i, шаг 1.

Следует отметить (см. [16]), что если при построении переобусловливающей матрицы к А был выбран алгоритм с уровнем заполнения равным г, то при построении переобусловливающей матрицы к Б необходимо выбирать алгоритм с уровнем заполнения, равным 2 г + 1.

В качестве параметра г, для того чтобы матрица X, а значит, и матрица Б были разреженными и легко умножались на вектор, необходимо выбрать нуль. Отметим, что приведенный алгоритм сохраняет блочную структуру матрицы В для матрицы X:

X

(Ь1Н) Н

О

(Ь1 Н) В1Н

О

О

( г 2л)-1

( г й-1

(¿й)-1 аН (¿2Ц))-1 ^Н

О

О

(¿й)-1 Ан

Ан

Прямоугольная матрица X полного ранга не является квадратной, размерность которой совпадает с размерностью матрицы В. В связи с этим следует отметить, что эффективное решение системы с матрицей Б невозможно в отличие от эффективного умножения на эту матрицу. Тем не менее матрица $ служит вспомогательной матрицей при построении переобусловливателя 5\

В общем случае, пусть для решения системы уравнений х' = Ну' построена переобусловливающая матрица Н, но она обладает неприятным свойством: эффективное решение системы с матрицей Н (нахождение обратной — Н-1) невозможно. Тогда следует построить такую переобусловливающую матрицу к Н — Н, для которой Н будет вспомогательной с той особенностью, что ее можно легко умножить на вектор. Процесс решения системы с матрицей Н будет выглядить как внутренняя итерационная процедура.

Пусть требуется найти вектор г* — решение системы Н г* = г, тогда внутренняя итерационная процедура будет следующей:

1) = 0;

2) и1 = и1-1 + в1 (Г — Ни1-1), I = 1,..., Ь;

3) г* =

где в1 — параметр процесса (в [13] — параметр чебышевского ускорения).

В настоящей работе в качестве такой процедуры при получении результата умножения вектора шагов 2 и 6 алгоритма 1-16 на матрицу Б-1 предлагается использовать ОМИЕВ-метод с вспомогательной матрицей Внутренняя итерационная процедура содержит в качестве своих итераций также алгоритм 1-16 с той особенностью, что

К = К

-1

Е.

Таким образом, построен итерационный метод решения системы (7) с переобусловливанием и использованием разложения вектора невязки в подпространстве Крылова (ОМИЕБ-метод). В численных экспериментах размерность подпространства Крылова варьировалась в пределах от 5 до 10. Начальное приближение (£°, п°, т°) = (0, 0, 0).

4. Результаты численного эксперимента

Рассмотрим тестовый пример задачи (1) с коэффициентом

(х, у) е пь

«°(х, у)

1,

Й2, (х, у) е П2-

О

О

О

При этом положительный параметр а2 изменяется таким образом, чтобы разрыв коэффициентов на подобластях увеличивался. Соотношение шагов сеток по каждому

направлению, как и ранее, Л = 2 Л2, число отрезков разбиения стороны П равно п

N

(к)

2 Лк

Разность между точным и приближенным решениями для компонент вектора скоростей будем определять в нормах ) и Н 1(е), а для функции давления — в

норме ).

Одна из рассмотренных характеристик — поведение погрешности в норме Ь2 на

полосах вблизи интерфейса Г12. Ширина полосы, примыкающей к интерфейсу, равна п

— с каждой стороны. Полосу в к

-й подобласти обозначим через ^^.

8

Пример. В качестве решения системы (1) выберем

эт х еоэ у + х еоэ у, (х, у) € П1,

—этх еоэ у +(2 п — х) еоэ у, (х, у) € П2;

тогда

у) =

—еоэ х эт у — эт у, (х, у) € П1, еоэх эту + эту, (х, у) € П2;

/l(x, у) =

р(х, у) = еоэ х эт у, (х, у) € П,

2 этх еоэ у + х еоэ у — эт х эт у, (х, у) € П1,

—а2 (2 эт х еоэ у — (2п — х) еоэ у) — этх эт у, (х, у) € П2;

ЛОх у) =

— 2 еоэ х эту — эту + еоэ х еоэ у, (х, у) € П1,

а2 (2 еоэ х эту + эту) + еоэх еоэ у, (х, у) € П2;

#1(х 0)

х + этх, (х, 0) € Г^ П1,

2п — х — эт х, (х, 0) € Г^ П2;

Таблица 1.

а2 ^(1) ни1 — инН ни2 — и£\\ н^1 — <н н^2 — ^11 нр1 — Нр2 — р1\\

10 10 10 16 32 64 0.0023431 0.0005929 0.0001489 0.0002180 0.0000550 0.0000137 0.0010125 0.0001985 0.0000469 0.0003786 0.0000920 0.0000230 0.0223022 0.0115002 0.0059003 0.0118944 0.0062433 0.0033015

100 100 100 16 32 64 0.0024603 0.0006296 0.0001581 0.0002250 0.0000561 0.0000140 0.0009873 0.0001885 0.0000445 0.0005175 0.0001269 0.0000321 0.0226011 0.0117994 0.0059376 0.0123778 0.0063928 0.0033667

0.1 0.1 0.1 16 32 64 0.0010369 0.0002581 0.0000643 0.0018707 0.0004511 0.0001107 0.0021211 0.0004940 0.0001229 0.0011953 0.0003068 0.0000792 0.0205748 0.0105118 0.0054222 0.0129343 0.0067412 0.0034766

0.01 0.01 0.01 16 32 64 0.0009415 0.0002353 0.0000588 0.0019537 0.0004955 0.0001215 0.0022563 0.0005295 0.0001321 0.0013521 0.0003461 0.0000878 0.0202156 0.0103172 0.0053926 0.0135580 0.0067624 0.0034815

Таблица 2.

a2 N(1) ||u1 — 411 ||u2 — ||v1 — у ||v2 — y

10 10 10 16 32 64 0.0506119 0.0269168 0.0140726 0.0219283 0.0110970 0.0056036 0.0314110 0.0142037 0.0066196 0.0142226 0.0066240 0.0031739

100 100 100 16 32 64 0.0521878 0.0277154 0.0144291 0.0210023 0.0104988 0.0052493 0.0318380 0.0143119 0.0066740 0.0147205 0.0068892 0.0033224

0.1 0.1 0.1 16 32 64 0.0420919 0.0210421 0.0105279 0.0767802 0.0402733 0.0216351 0.0321798 0.0146107 0.0068550 0.0246155 0.0121144 0.0060792

0.01 0.01 0.01 16 32 64 0.0419619 0.0209714 0.0104847 0.0794782 0.0435169 0.0236703 0.0324776 0.0147743 0.0069455 0.0266743 0.0131115 0.0065212

Таблица 3.

a2 N(1) ||u1 — ||u2 — ^ ||v1 — у ||v2 — у

10 10 10 16 32 64 0.0018073 0.0004592 0.0001155 0.0001034 0.0000265 0.0000066 0.0006844 0.0001374 0.0000306 0.0001755 0.0000397 0.0000094

100 100 100 16 32 64 0.0019234 0.0004946 0.0001242 0.0000947 0.0000190 0.0000067 0.0006879 0.0001491 0.0000345 0.0002778 0.0000637 0.0000154

0.1 0.1 0.1 16 32 64 0.0004027 0.0000988 0.0000244 0.0016287 0.0003758 0.0000897 0.0012020 0.0002999 0.0000762 0.0008459 0.0002023 0.0000500

0.01 0.01 0.01 16 32 64 0.0002491 0.0000624 0.0000155 0.0016974 0.0004116 0.0000983 0.0013127 0.0003281 0.0000834 0.0009524 0.0002272 0.0000553

. —x — sin x, (x, п) £ Г3 П П1, (x, п) = < -

x — 2п + sin x, (x, п) £ Г3 || П 2;

Ы2у) = 2 зту, (2 у) £ Г2; #4(0, у) = -2 зту, (0, у) £ А;

и (ж, у) = 0, (ж, у) £ Г2 и Г4; г>(ж, у) = 0, (ж, у) £ Г1 и Г3.

Для данной задачи определим значения основных параметров. Имеем (N(1), N(2)) = (16, 32); (32, 64); (64, 128); а2 = 10; 100; 0.1; 0.01. Результаты расчетов приведены в табл. 1-3: в табл. 1 — погрешность решения — компонент вектора скорости w = (и, V) и функции давления р, на подобластях , в норме ); в табл. 2 — погрешность

решения — компонент вектора скоростей w, на подобластях , в норме П Н 1(е); в

табл. 3 — на полосах ^, в норме ).

Заключение

Расчеты численного эксперимента показали, что погрешность не слишком сильно растет при увеличении разрыва коэффициентов на подобластях в нормах пространств ) и П H 1(e), k =1, 2. Большая часть погрешности (в норме пространства L2)

сосредоточена вдоль полос Qk в окрестности интерфейса Г12 между подобластями, т. е. погрешность не распространяется во внутреннюю часть подобластей.

Легко видеть, что погрешность решения для компонент вектор-функции скорости w = (u, v) в норме пространства H 1(e) и скалярной функции давления p в норме

eezf)

пространства 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 Belgaoem 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 Belgaoem 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] Flemisoh В., Melenk J., Wohlmuth B. Mortar methods with curved interfaces // APNUM. 2005. Vol. 54. P. 339-361.

[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 г.