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

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

Том 10, № 6, 2005

УСЛОВИЯ ДИССИПАТИВНОСТИ И М-МАТРИЧНОСТИ РАЗНОСТНОГО ОПЕРАТОРА КОНВЕКЦИИ-ДИФФУЗИИ С ГРАНИЧНЫМИ УСЛОВИЯМИ ТРЕТЬЕГО РОДА*

Л. Г. ЧИКИНА, И. Н. ШАБАС ЮГИНФО Ростовского государственного университета, Ростов-на-Дону, Россия e-mail: [email protected], [email protected]

Positiveness of the real part and M-matrix property conditions for the convection — diffusion difference operator with boundary conditions of the third type are obtained. Positiveness of the real part and M-matrix properties depend on the central or upwind approximation of the convective terms and on the specific representation of these terms in the symmetric or nondivergent form.

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

В области П х T, П = П U Г, рассматривается система трехмерных уравнений, описывающая процессы переноса вещества в несжимаемой среде:

dS Л д ( dSN Л д , m dS £

St- £ dXi ("'sXj +1 £ aX (v‘s) +(1 - Y) £v'al, + es =1 (1)

t=l 4 7 '=! ,= 1

divv = 0, (2)

где П = {х = (х, у, г)} — трехмерная область расчета (акватория водоема); Б — концентрация вещества; ^,^2,^3 — коэффициенты турбулентной диффузии вещества; vl,v2,v3 — скорости движения среды по направлениям х, у и г соответственно; V = {и, V, /ш} — вектор скорости; в = в(х,у,г,Ь) > 0 описывает консервативность рассматриваемого вещества; 7 £ [0,1] — параметр записи уравнения (1). Условие несжимаемости среды (2) позволяет

записать систему (1) в недивергентном виде (7 = 1), в дивергентном виде (7 = 0) и в

эквивалентной симметричной форме (7 = 1/2).

Система (1) замыкается начальными

Б(х, 0) = Б0(х), х € П, (3)

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 03-01-00005 и № 05-01-00096-а), а также Российского фонда фундаментальных исследований и администрации Ростовской области (№ 04-01-96807) и программы Университеты России (УР.03.01.024). © Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

и краевыми условиями

Мх)Х(Х)*(Х) = Г(Х)’ Х е Г’ (4)

где ^(х), х(х), г(х) — кусочно-гладкие функции. В зависимости от равенства нулю функций ^(х), х(х), г(х) допускается возможность постановки условий 1, 2 и 3-го рода.

2. Разностная аппроксимация

В области Й вводится равномерная по всем направлениям разностная сетка Йь — и Гь с векторным параметром Н — (Нх, Ну, Н%), где Нх,Ну,Н% — соответствующие шаги сетки. Здесь — множество внутренних узлов сетки; Гь — множество граничных узлов.

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

При аппроксимации первых производных в граничных условиях используются правые либо левые разности в зависимости от того, какая граница рассматривается [1].

Рассмотрим, например, разностный аналог на семиточечном шаблоне на оси ШЕ (рис. 1) краевых условий 3-го рода. Он выглядит следующим образом:

—

------к----------+ 80]к — гш,

Нх

&Ызк — -!зк ,

№е Н + Хе SNjk ,

преобразованный вид этих условий следующий:

(^ш + Нх)

Нх Н

(№е + Хе Нх ) №е

S0jk--------Г~ s1jk —

Ж 5

°х Нх

SNjk------Т~ SN — 1 jk — г

Нх Н

е •

я

'Т(юМ\) У'у Х/У(и+и)

Шг-ЧЮ Е(}+Чк)

О (и, к) X

>Вш,к-г>

Рис. 1. Семиточечный шаблон.

Таким образом, в общем виде аппроксимация 3-го граничного условия будет иметь вид

дв 5 + рв ? = ге,

где —, — — значения концентрации вещества соответственно в некотором граничном и приграничном узле. Для коэффициентов дв ,ре ,те справедливо

— (^е + Хе На)

в Н ’

11а

№е

Ре = — ~ , те — те ,

На

где 0 заменяется на символы Ш,Б,Б,Т,Ы,Е, если граница области приходится соответственно на (г — 1), (] — 1), (к — 1), (к + 1), (] + 1), (г + 1)-й узел семиточечного шаблона,

а — х,у, г.

В результате задаче (1), (3), (4) поставлен в соответствие дифференциально-разностный аналог

д—н

+ Ьнвн — ¡н(х,£), х,г е Пн х П*,

де5 + ре 5 — те, х,г е Гн X П*, (5)

вН — ин(х, 0), х,Ь еПн х {0},

П* — П* и{0},

где Ьн — Ьон + Lcн + Ь^н, Ь^н — разностный оператор диффузионного переноса; Ьсн — разностный оператор конвективного переноса; Ь^н — разностный аналог функции взаимодействия веществ.

Исключив решение в граничных точках области Пн и учитывая разностные краевые условия, перейдем к неявной операторно-разностной схеме:

„п+1 — -п

-------н + ьн(-п)-п+1 — л;, -н — и0, п > о. (6)

т

Здесь Ьн — это оператор Ьн, в который уже включены граничные условия.

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

Теорема 1. (Таусски) (достаточное условие положительной определенности) [2, 3]. Пусть симметричная матрица А е Мп:

1) неразложима;

2) матрица с диагональным преобладанием, т. е.

п

\аи\ > ^ ] |\; j=i

3) матрица, у которой хотя бы для одного і

П

\аа\ > 'У ] \аіі\;

і=і

4) матрица, все диагональные элементы которой положительны

ац > 0, г,] — 1, 2,..., п.

Тогда все собственные значения матрицы А строго положительны.

Существует несколько эквивалентных определений М-матричности. Для наших исследований необходимо следующее

Определение 1. Невырожденная матрица А называется М-матрицей, если ее вне-диагональные элементы неположительны, а обратная матрица поэлементно неотрицательна [4].

Теорема 2. Пусть матрица А е Мп и

а^ < 0, г — ], ац > 0, г,] — 1, 2,...,п

имеет строгое диагональное преобладание или же А является неразложимой и имеет диагональное преобладание. Тогда матрица А является М-матрицей [3].

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

3. Диссипативность рассматриваемого разностного оператора

Определение 2. Невырожденная матрица А называется диссипативной, если ее реальная часть А0 — ^ (А + АТ) положительно определена [1].

Для разностных операторов Ьн, получаемых в результате пространственной аппроксимации уравнений системы (1) при центрально-разностной аппроксимации конвективных членов в уравнении конвекции-диффузии (1), (2) с 7 — 1/2 при наличии граничных условий 3-го рода и постоянного коэффициента консервативности рассматриваемого вещества, доказано достаточное условие диссипативности рассматриваемого разностного оператора.

Теорема 3. Пусть в уравнении конвекции-диффузии (1), (2), записанном в симметричной форме (ч — 1/2), с краевыми условиями 3-го рода (5) и в — в(х,у, г,Ь) > 0 конвективные члены аппроксимируются центральными разностями.

Тогда для того чтобы оператор Ьсн — разностный аналог стационарной задачи конвекции-диффузии — был диссипативной матрицей, достаточно выполнение неравенств

в^к — ^ ¿е (1 + Р^) 0С{О)^к > 0, (7)

в=Ш,Я,Б,Т,М,Е ' де '

г — 1, ...,ЫХ — 1, ] — 1,..., Ыу — 1, к — 1, ...,ЫХ — 1,

где хотя бы одно из них является строгим. Здесь 0^к — элементы матрицы А0, симметричной составляющей оператора Ьсн в (6) в приграничном узле для соответствующей границы; 5в — символ Кронекера для соответствующей границы; дв ,рв — коэффициенты из (5), дв = 0.

Доказательство. Оператор Ьсн = Ьон+ЬСн+Ь^н, полученный в результате включения в него граничных условий 3-го рода, сохранил кососимметричность оператора ЬСн. Так как вклады граничных условий 3-го рода вошли в диагональ оператора Ьон, нам достаточно найти условия диссипативности оператора ЬСн + Ь^н ■ Оператору ЬН соответствует матрица Ас = А0 + А1, где кососимметрическая составляющая А\ = ЬСн, а симметричная составляющая А0 = Ьон + Ь^н. Матрица А0 является неразложимой [3, 5]. Так как матрица А0 является неразложимой, по теореме 1 достаточно найти условия, при которых матрица А0 имеет диагональное преобладание, а ее диагональные элементы будут положительными:

1 агг1 ^ 1 ац \ 1

агг > 0.

Сумма норм — это неотрицательная величина, поэтому достаточно найти условия выполнения неравенства

агг ^ \ агЦ \1

г=3 агг = 0.

Эти условия обеспечат положительность диагональных элементов. Симметричная матрица А0 имеет вид

\

Ас

Ао

где Шг,Ог,Ег — матрицы размером (Му — 1) х (Му — 1):

(8)

( °1 Е1 0 ■■■ 0

W2 О е2

0 Wз Оз ■■■ 0

■■■ ENx—2

0 0 WNX-1 °^-1

(9)

Wi

Ег

Wгl 0

00

0

0 0 ■■■ WiNv-1

Ег1 0 ■■■

0

Егм.-1

( Ог1 Nil 0 ■■■ 0 \

Бг2 Ог2 N2

0 Бг3 Огз ■■■ 0

■ ■■ -2

0 0 1- гГ г 1- Бг /

А

Здесь Wij, Бц, Оц, N,¿3 ,Ец — матрицы размером — 1) х — 1) Распишем поэлементно матрицы Ог, Wi, Ег:

О

( П^С 0 ■■■ 0 \

Щ2 Щ2 ГПС ■■■

0 Щз Щз ■■■ 0

ГПС 1ijNz -2

0 0 ^^N2 -1 пс DijNz -1 )

Б] 0

0 0

ОС

°г]Мг -1

N

г]

г]

Е

г]

Щ]1 0 ■ ■ 0

0 0 N С ■ ■ ^г]Мг-1

Щг 0 ■ ■ 0

0 0 ■ ■ Щм2-1

С] С Е 0 0

0 0 7 : ^ Е

Коэффициенты матрицы А0 для внутренних точек области имеют следующие значения:

— элементы первой, второй и третьей поддиагоналей соответственно

т0 ■ к =

(0) г] к

а

г]к

кХ

О С0 О(0)]к

г]к

Щ ’

в Со

(0)г]к

а

диагональные элементы (без вкладов граничных условий)

ВСо

(0)г]к

ах _|_ ах

иг+1]к ' игщк

+

а

г]+1к

+ а

у ах + а2

г]к . г]к+1 ' гщк . п

+ ——7^-----------— + вг]к;

кХ Щ к2

элементы первой, второй и третьей наддиагоналей соответственно

(10)

Т С0 (о)г]к

а

гщк+1

к2

NСо. к

(0)г]к

1г]+1к

ц

Е Со

(0)г]к

а

г+1]к

кх ■

В этих обозначениях в из (6) будет иметь вид

ЬН8 = тг°к 8г-1]к + Б] 8г]-1к + вЩк $г]к-\ + ОЩк 8гщк

+ТгЩк вг]к+1 + ^^гщ+гк + ЕгЩк вг+1]к,

где

тЩк

(1 - Ъ Кгк вЩ

(1 - бв )в;

(о)г]к;

Щ = (1 - ¿т)Т0щ; $°к = (1 - К )^:щ; ЕС0к = (1 - ¿е)ЕС0щ,

'(0)г]к’

ь0

'г] к

о

(0)г]к;

г]к

(0)гЗк’

Щк = °Со г]к -—тСог]к - 6~—Осо г]к - ¿в—вС0г]к - ¿Т—тСоог]к-

(0)г]к ^ д (0)г]к Б д (0)г]к

'в д (0)г]к Т д (0)г]к

-а ЕкNС0 - а Ее ес0

°м дк 1У(0)г]к °Е дЕ Е(0)г]к ■

0

У

х

У

У

2

С

0

Для выполнения (8) необходимо, чтобы

- 5,^- 5 1Рэ£е°,,, - 5в1РвБс°;,, - Тс°;,,-

Г)Со _ 5 _ - 5 ..

( 0)4'^ ш д ( о)цк э д ( о)цк

'т дт" (0)іік

_ 5 Ек м с° _ 5 — Е с° >

5^ дк < о )іік °Е дЕ Е(°)іік >

Ре

(1 - )w;°iіk + (1 - 53)С

( °)іік

+

+

' ( о )іік

(1 - 5в )в(°)іік + (1 - 5т)Т(°)іік + (1 - 5к)м(°°)іік + (1 - 5е)Е

( о )іік

Так как Щ(°)г]к, , в] Т^к > N0]> ЕС(0)г]к являются отрицательными величинами по

построению, а (1 - ¿©) > 0 по определению символа Кронекера и значения (1 - А©)©^] являются неположительными,

1(1 - 5е)ОІіік| = -(1 - 5©)0

(о)іік

и неравенство (8) в нашем случае имеет вид

Г)со _ 5 РШ шсо _ 5 РЭ ясо

и(Р)іік 0ш дш ш(о)іік °э д С(о)Цк

- 501РвВV,, - 5^рттс°и-

"в дв (о)іік

'т д/ (о)іік

-5»-м1ик - 5ер-Ес:,іік > -(і - к- (і - 5эКі-

Е

дк

-(1 - 5в )В(о)іік _ (1 - 5т)т(°)іік _ (1 - 5кЖ°)іік _ (1 - 5Е)Е2і

■ ( о)іік

о

( о)іік

(о)іік'

(11)

Из (10) видно, что диагональные элементы без вкладов граничных условий можно записать в виде

вСо

( о)Цк

тд/со Ссо Всо Тсо Nсо Есо + в •

_ ( о)цк °(о)цк (о )іік (о)цк (о )іік (о)іік + віік‘

Тогда после приведения подобных предыдущее неравенство (11) примет вид

1 _ 5ш _ 1 _ 5ш 1дШ) Ш(о)Цк + ^ _ 5э _ 1 _ 5эіік+

+ ( 1 _ 5 в _ 1 _ 5 в — } В ,1іИк + ( 1 _ 5т _ 1 _ 5т ~ ) Т,°.іИк +

дв

( о)іік

'т дт ! "( о )іік

+ ( 1 _ 5к _ 1 _ 5к рк) Кііік + _ 5Е _ 1 _ 5е Р^) Е\0))іік + віік > 0

или после элементарных преобразований

Рш

дш

( о)іік

Рэ

_5ш 1 + — ш;оіік _ 5э 1 + ^ хо^к - 5в 1 + — Віік -

дэ

( о)іік

Рв

дв

(о)іік

-5т 0+р£) т^к _ 5~ I1+д,)_ 5е І1+іе) Е°о”к+віік >0

следовательно, в компактной форме записи

Е 5в( 1+Е) 0^к+вик >

« В Т N Е V У в /

©=Ш, Б, В, Т, М, Е

о

о

с

о

Таким образом, имеем (7). Теорема доказана. □

Для операторов нестационарного уравнения конвекции-диффузии с краевыми условиями 3-го рода в случае неявной схемы при центрально-разностной аппроксимации конвективных членов показана устойчивость используемой разностной схемы при выполнении условий теоремы 3.

Запишем схему (6) в каноническом виде

(Е + тАс)вг + Ас5га = ¡. (13)

В работе [6] для случая несамосопряженного оператора Ас сформулирован критерий устойчивости. Приведем формулировку этой теоремы для неявной схемы (13).

Теорема 4. [6] Если оператор А-1 в схеме

(Е + тА)в1 + Ав = 0

существует, то для устойчивости схемы в гильбертовом пространстве Ь2н, т. е. выполнения оценки

||вга+1|| < М, п > 0, необходимо и достаточно выполнения операторного неравенства

А + 0,5тА*А > 0. (14)

Неравенство (14) означает выполнение условия

(Ах, х) + 0,5т|| Ах||2 > 0.

Поэтому имеет место

Теорема 5. Для устойчивости неявной схемы (6) при центрально-разностной аппроксимации конвективных членов нестационарного уравнения конвекции-диффузии (1),

(2), записанного в симметричной форме ('у = 1/2), с краевыми условиями 3-го рода (5) достаточно выполнения неравенств (7), в котором хотя бы одно из неравенств является строгим.

Доказательство. Если выполнены условия теоремы 4, то для устойчивости схемы (6) достаточно диссипативности оператора Ас.

Теорема доказана. □

4. М-матричность рассматриваемого разностного оператора

Для разностных операторов Ьрн, получаемых в результате пространственной аппроксимации уравнений системы (1) при противопотоковой аппроксимации конвективных членов в уравнении конвекции-диффузии (1), (2) с 7 _ 1 при наличии граничных условий 3-го рода и постоянного коэффициента консервативности рассматриваемого вещества, доказано достаточное условие М-матричности получаемого разностного оператора.

Теорема 6. Пусть в уравнении конвекции-диффузии (1), (2), записанном в недивергентной форме (ч =1), с краевыми условиями 3-го рода (5) и в = в(х,У,г,^) > 0 конвективные члены аппроксимируются разностями “против потока”.

Тогда для того чтобы оператор Ьрн — разностный аналог стационарной задачи конвекции-диффузии — был невырожденной М-матрицей, достаточно выполнение неравенств

в=Ш,Я,Б,Т,М,Е

1 + £.) «и * °

(15)

3

1, ■■■уіїу — 1,

к = 1, ■■■,Ы7 — 1,

в котором хотя бы одно из неравенств (15) является строгим. Здесь ®р0)іцк — элементы матрицы Ар оператора Ьсн в (6) в приграничном узле для соответствующей границы; 8в — символ Кронекера для соответствующей границы; дв ,рв — коэффициенты из (5), 9@ = °.

Доказательство. В случае противопотоковой аппроксимации матрица Ар имеет квадратно-блочно-трехдиагональную структуру вида (9). Такая матрица является неразложимой [3, 5]. Кроме того, внедиагональные элементы матрицы Ар отрицательны. Так как вклады граничных условий 3-го рода вошли в диагональ оператора Ь^н, нам в силу теоремы 2 для доказательства свойства М-матричности достаточно показать выполнение диагонального преобладания матрицы Ар и положительность элементов матрицы, лежащих на главной диагонали. Таким образом, необходимо показать выполнение

|а«| >5] \ац|,

і=і ац > 0-

Сумма норм — это неотрицательная величина, поэтому достаточно найти условия выполнения неравенства

ац > ^ \ац |, і=і аіі = °

(16)

Эти условия обеспечат нам положительность диагональных элементов.

В результате противопотоковой аппроксимации пространственных производных получена разностная схема на семиточечном шаблоне

8і-1Цк + ‘С(0)іЦк вЧ-1к + ВР(0)іік вчк-! + ^(0)іЦк 8цк + Т(0)Цк вцк+1+

+ Що)г]квЧ+1к + Е\о)1]к8г+1]к 1^к, (17)

коэффициенты которой для внутренних точек сеточной области имеют следующий вид:

шр

ш(о)цк

іЦк

ьХ

п.

+

іЦк

К

Ср

(0)іЦк

іік

ьУ

V,

+

Ц к Ьу, :

в р

В(0)іЦк

іі к ч

+ іі к

Ь,

Пр

0)іік

а

і+1Цк

+ аХц

ч

\пЦ к \

Ьх

+

у

1іЦ+1к

ьу

+ аЬк+і + аЬк+

Ы

ьу

ь,

(18)

у

X

,

+

тР _ _ аі]к+ї . 'Ші]к Nр

(0)і]к _ ^2 + Ь ’ (0)іІ&

Для ОРр0)^к имеет место

Що^к = --^{0)Цк - ^(О^к - Б(0)г7к - Т{0)Цк — ^(0)^ - Що^к + А?к• (19)

Аналогично доказательству теоремы 3 для выполнения (16) необходимо, чтобы

имеет место

Используя (19), а также то, что по построению , 5^к, Б^к , Щк ,Е^к являются

отрицательными величинами, а (1 — 5@) > 0 по определению символа Кронекера, имеем по аналогии с выкладками в доказательстве теоремы 3 и, используя определение модуля,

Для операторов нестационарного уравнения конвекции-диффузии с краевыми условиями 3-го рода в случае неявной схемы при противопотоковой аппроксимации конвективных членов с помощью принципа максимума получены оценки решения задачи через ее правую часть, начальные и граничные условия при выполнении условий теоремы 6.

Теорема 7. Для решения нестационарного уравнения конвекции-диффузии (1), (2), записанного в недивергентой форме ('у = 1), с краевыми условиями 3-го рода (5) в случае использования разностей “против потока” при выполнении условий (15) для неявной схемы (6) справедливы оценки:

или после элементарных преобразований

значит, в компактной форме записи

а следовательно, получим (15). Теорема доказана.

□

для консервативного вещества (ß(x,y,z,t) = 0)

◦ + max

C n>0

f n

e 1 1 + I eijk e=W,S,B,T,N,E \ 9&

Pe

E Sei 1 + - e?

+

с *

+

E Ser-eePjk

e=W,S,B,T,N,E ge

e=W,S,B,T,N ,E

£ Se( 1 + pel ej

re

C*

для неконсервативного вещества (ß(x,y,z,t) > 0)

fn

ß

+ max

° n>0

C

fn

ßijk - E 1 + — ) eijk

e=W,S,B,T,N,E

re

+

C*

+

E Sereejk

e=W,S,B,T,N,E ge

ßijk E

e=W,S,B,T,N ,E

Se I 1 + — I e

re

ijk

C*

Доказательство. Приведенные оценки получены с использованием идеи методики [7],

а также при условии выполнения условия M-матричности (15) из теоремы (6).

Теорема доказана. □

Случай гв = 0 рассмотрен в работе [6].

Список литературы

[1] Самарский А.А., Гулин А.В. Устойчивость разностных схем. М.: Наука, 1973.

[2] Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. М.: Мир, 1991. 367 с.

[3] Ортега Дж., РЕйнволдт В. Итерационные методы решения нелинейных систем уравнений со многими неизвестными. М.: Мир, 1975. 558 с.

[4] РихтМАйЕР Р., МОРТОН К. Разностные методы решения краевых задач. М.: Мир, 1972. 420 с.

[5] ХОРН Р., ДЖОНСОН Ч. Матричный анализ. М.: Мир, 1989. 655 с.

[6] САМАРСКИЙ А.А., Вавищевич П.Н. Численные методы решения задач конвекции-диффузии. М.: Изд-во УРСС, 1998. 272 с.

[7] КРУКИЕР Л.А. Неявные разностные схемы и итерационный метод их решения для одного класса систем квазилинейных уравнений // Изв. вузов. Математика. 1979. № 7. С. 41-52.

Поступила в редакцию 30 июня 2005 г., в переработанном виде — 8 августа 2005 г.