Неконформные схемы МКЭ для гиперболических систем линейных уравнений Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

Том 152, кн. 1

Физико-математические пауки

2010

УДК 519.63

НЕКОНФОРМНЫЕ СХЕМЫ МКЭ ДЛЯ ГИПЕРБОЛИЧЕСКИХ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ

Е.М. Федотов

Аннотация

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

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

Ключевые слова: гиперболические системы уравнений, сеточные схемы, метод конечных элементов, разностные схемы с весами.

Методам построения сеточных схем для линейных и нелинейных гиперболических систем уравнений и их решению уделяется большое внимание. Это связано с большим количеством приложений для таких уравнений. При конструировании численных методов используются различные подходы (см.. например. [1] и имеющуюся там библиографию), которые позволяют строить сеточные схемы, воспроизводящие различные особенности решений этих уравнений.

В настоящей работе предлагается метод построения схем сквозного счёта для многомерных гиперболических систем линейных уравнений первого порядка в предположении их симметрии. Метод построения сеточных схем для таких систем уравнений основан на использовании предложенного в [2] метода аппроксимации скалярного уравнения конвекции-диффузии. Использование предельной аппроксимации Галёркина Петрова смешанной постановки исходной задачи позволило учесть направление конвективного переноса. При аппроксимации систем уравнений такой подход естественным образом позволяет учитывать локальные направления характеристик, а также сохранять основные свойства пространственного оператора исходной задачи. Предложенные здесь сеточные схемы являются вариантом неконформного метода конечных элементов, в котором приближённое решение ищется в пространстве кусочно-полиномиальных функций, допускающих разрывы на границах элементов триангуляции.

Построенные в работе сеточные схемы включают в себя как частный случай известные противоиотоковые схемы первого порядка точности и являются естественным обобщением на случай гиперболических симметричных систем известных разрывных локальных схем метода Галёркина (см. [3 5] и цитированную там литературу).

1. Постановка задачи и построение полудискретной схемы для задачи Коши

Рассмотрим задачу об определении вектор-функции и : Д" х [0, Т] ^ Дт, т > 1) удовлетворяющей уравнению

о ТЭ " Л

-дГ + Т,д^(л<-')и)+Си = ^ .г-еД'\^е(0 ,Т] (1)

¿=1 ¿

с начальным условием

и(х, 0) = и0(х), X е Д", (2)

Будем полагать, что В(х, £), С (х,Ь), А(г) (х,Ь), % = 1 ,...,п, являются вещественными непрерывно дифференцируемыми т х т матрицами, причём матрица В предполагается симметричной и положительно определённой. Будем также полагать, что при любом £ е Д" матрица

"

А(£)=£ А« & = А о £ (3)

¿=1

является симметричной, что обеспечивает гиперболичность системы уравнений (1).

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

J ■ п (1х ~ А(<)и • ¿т = - Си) • п с1х, (4)

ш ш ¿=1 ш

справедливого для любой ограниченной области ш С Д" и любой непрерывно дифференцируемой в ш вектор-функции п такой, что п|дш =0.

Для сокращения записи при выводе сеточных уравнений будем полагать, что матрицы А^ , % = 1,..., п постоянны, а матрица С и вектор f равны нулю.

Частным случаем изучаемой системы уравнений является риманов случай, когда система уравнений (1), благодаря введению новых переменных, распадается на систему независимых скалярных или слабосвязанных уравнений. Естественно требовать, чтобы сеточная схема, построенная для систем уравнений общего вида, при этом также приводила к устойчивым аппроксимациям соответствующих уравнений.

Указанный случай соответствует попарной перестановочности матриц А». % = 1,... ,п. Рассмотрим его. Известно, что в таком случае существует ортогональная матрица Q, что

Л« = QA(¿) QT, % =1 ,...,п, (5)

где Л^ - диагональные вещественные матрицы. Их диагональные элементы А^ , а = 1,... ,т, являются собственными числами матриц

Учитывая (5), запишем интегральное тождество (4) следующим образом:

= (б)

ш ш ¿=1

Введём обозначения г = Qu для переменных Римана, а также q = Qn для пробных функций. Тогда (6) запишется в виде

[г;У; л,,;,!,■ II. V, = (Л^, . . . , А<">). (7)

ш ^=1 ш

Аппроксимируем тождество (7). Заметим, что каждое слагаемое под знаком суммы в этом тождестве соответствует конвективному слагаемому в скалярном уравнении (1) [2] (см. там же (4)). Воспользуемся методикой раздельной аппроксимации слагаемых, предложенной в указанной работе.

Аналогично [8] построим в Д" конформную триангуляцию 7н. Через К условимся обозначать элементы триангуляции, через ек — границы элементов К. Приближённое решение ин задачи будем искать в иространстве Ут = Ун х Ун х ■ ■ ■ х Ун кусочно-полиномиальных, вообще говоря, разрывных функций, где

Ун = Уик = {пн е Ьто(Дп) : пн|к е Рк(К), К е Тн},

Рк (К) - пространство полиномов степени не выше к по каждой переменной на элементе К. Через ^н будем обозначать произвольную конечную часть триангуляции

о

7Н, через Ун - подпространство функций из Ун, равных нулю вне .

В соответствии с [2] тождество, аппроксимирующее (7). будет иметь вид

гдВи г т

1/!

т л

■ Пн ¿х + - / г^н^- ■ dx+

¿=1 к

т „

+ ЩЩ У ■ ^ + - «¿,Н|К')

ек ¿=1ек

= 0 V пн еУ т, V ^н с 7Н, (8)

где К' - элемент о>н, смежный с К через сторону ек, V - внешняя нормаль к ек со стороны элемента К, = (Н ± «)/2 - положительная и отрицательная части функции V.

Учитывая связь гн = ^ин и qн = QпН, первые слагаемые в квадратных скобках (8) запишем в виде

тт л л п

X У • / X • 777" (9)

¿=1 к к ¿=

Введём в рассмотрение матрицу

А+ = Qт Л+Q,

где Л+ = diag((v1 ■ v)+,..., ■ v) + ), тогда второе слагаемое в квадратных скобках (8) примет вид

т ,, ,,

/ ■ *)+ (д^н - |к') ¿х = ^ / А+ин ■ (пн - Пн|к') ¿х. (10)

ек ¿=1ек ек ек

Для симметричной т х т матрицы С = Ст построим спектральное разложение С = ТтЛТ, где Л = diag(Лl(C),..., Лт(С)), Л^(С) — собственные числа матрицы С Тт = Т-1. Введём в рассмотрение матрицы С± = ТтЛ±Т, где Л± = = diag(Л±(C),..., ЛТ(С )). Соответственно, определим матрицу |С | = С+ + С Ясно, что построенные таким образом матрицы неотрицательны, симметричны и С = С + - С- .

Нетрудно видеть, что

п

А+ = АМ+ = ¿(¿О + = (Ло^+. 4=1 '

С учётом (9). (10) тождество (8) запишется в виде

3Buh ^

• щ dx + ^

А® uh ■ dx + j Cuh ■ щ dx+

Th Ke7hL K i=1 K

^ (Aov)+Uh • (nh - Vh |k') dx

= J f • nh dx V nh GVm, V "h cTh. (11) Th

Под решением сеточной схемы метода прямых для задачи Коши для уравнения (1) с начальным условием (5) будем понимать вектор-функцию Uh G Vm при

о

любых nh GVm, С Th, удовлетворяющую тождествам (11) и

У (uh - uo) • nh dx = 0 V^h Gl°m, ^h CTh, t = 0. (12)

Th

Такую приближённую постановку будем использовать как в римановом случае, когда матрицы перестановочны, так и в общем случае, когда относительно матрицы A(£), определённой равенством (3), предполагается лишь симметричность при любом £ G Rm, включая случаи, когда матрицы B, , i = 1, 2,..., n, зависят от переменных x и t.

Отметим, что наиболее близкими к предлагаемым здесь схемам являются варианты схем разрывного метода Галёркина для систем Фридрихса [6], а также схемы с использованием численных потоков Виджайясундарама (Vijayasuridaram) (см., например, [7, § 2.4] н приведённую там библиографию), применяемые для численного решения нелинейных гиперболических систем уравнений.

2. Схема метода прямых для смешанной краевой задачи

Пусть О - многоугольная область в Rn с границей Г. Рассмотрим задачу о нахождении вектор-функции u : О х [0, T] ^ Rm, являющейся решением уравнения

О г) n о

+ = *е(о,т] (13)

i=1 г

с граничным условием

(Aov)-u = 0, x G r-(t) = {x G Г : (Aov)- = 0}, t G (0, T], (14)

где v - внешняя нормаль к г ранице Г,ис начальным условием

u(x, 0) = u0(x), x G О. (15)

Построение схемы метода прямых проведём указанным выше способом, путём аппроксимации интегрального тождества

/ ~§Г 'пdx ~ / ' Ъ7 dx + /(Аог/)+и' dx = У(f " С") ' Vdx, (16)

Q ш i=1 Г Q

справедливого для любой достаточно гладкой вектор-функции n G Rm •

Пусть, как и раньше, Th - конформная триангуляция области О с границей rh, согласованной с границей Г K - элементы триангуляции, eK - границы элементов,

Vh = Vh,fc = {nh G LTO(0) : nh|k G Pfc(K), K G Th}.

Под приближённым решением схемы метода прямых для задачи (13) (15) будем понимать вектор-функцию и^ € удовлетворяющую при любой вектор-функции п € Ут интегральным тождествам (ср. с (11), (12)):

Г дВиь ■ щ с1х + 53 [- / Е . ^ +

dt

Th

кеЪ

h к

¿=i

• (nh - Пь|к') dx = J(f - CUh) • nh dx, (17)

eK ек Th

j(Aov)-Uh • nh dx = 0, V nh G Vhm, J(uh - uo) • nh dx = 0, V nh G Vh", t = 0.

(18) (19)

В тождестве (17) для однородности записи принято, что пл.|К' = 0 для сторон ек, являющихся частью грапицы Г^, приграничных элементов К.

Введём в рассмотрение полунормы в Дт : = (|Ао^|

\i/2 , , v^v) и ,v,A±

1/2

= ((Ао^)±у^у) • Обозначим через 7 множество сторон элементов разбиения , лежащих в области П, символом [|и^|] обозначим скачок функции и^ вдоль нормали к стороне е € 7. Справедлива

Теорема 1. Сеточная схема .метода прямых (17) (19) устойчива и справедлива оценка

t t Uh(t)||| + / / ,[|uh|],A dxdt' + / / ,Uh,A+ dxdt' <

0 еет e ^ 0 r\

< ec

|uh(0)|B + / ||fd0

./o

, c > 0, t G (0,T], (20)

где ||v|B±i = J B±1v • v dx.

Th

nh = uh

1 д ^ f n du /*

Buft • uft ¿ж + ^ _ / X! ' "o- dx + / • иЛ

2 dt

Th

KeT,

dx,-

h к

к

+ XIУ (Аог/)+ил • (ил - иЛ |А-/) йж = J Г иЛ с1ж - \ I ' ил (21)

еК ек Тг Тг

Преобразуем сначала первое слагаемое в квадратных скобках в (21), пользуясь следующим представлением подынтегрального выражения, учитывающим симметрию матриц

1

2 дж,

dx,- 2 dx,- 0 -1-

Имеем:

/п ди. 1 С п С 1

ЕА^пн' с1х = 2 У 'ил ~ ^У 2(Лог/)ил'ил

к

Тогда, принимая во внимание равенство Ло^ = (Ло^)+ — (Ло^) , второе слагаемое в (21) запишем в виде

1 Р п дА.(*)

•/=2У Е Е

Т 4=1 КЕТ, ек (

^(Аог/)+иЛ -Щ +

+ ^(Аог/) иЛ • и/) - (А о г/) + иЛ •

¿ж.

Учитывая далее кратность вхождения сторон элементов триангуляции ек, лежащих в О, граничное условие (18) и отсутствие смежных элементов для сторон, лежащих на границе области, получим

1д_ 2~д1

1=\ [ Ё ■ил с1х+\ Е /1 Иил |] +^ / 1ил Iа+1..

Т- еЕ^е " Г-

Подставляя полученное для J выражение в (21), придём к равенству

J Виь ■ иЛ &с+У | [|иЛ |] 3,х + ^ J\uh 3,х =

Т- еЕ7 е Г-

/1 Г (дВ п дА^*) \ ¡■иксЬ-~ + + 2С )ил ' ил (23)

Т- Т- *=1 *

откуда, пользуясь непрерывной дифференцируемостью коэффициентов матриц В, С и А^, а также неравенством Гронуолла, придём к утверждению теоремы. □

3. Сеточная схема с весами

Аппроксимируем схему метода прямых (11), (12) схемой с весами. Для этого построим па отрезке [0, Т] сетку с!т = , ] = 0,..., N¿1, ст = с!т \ {Т}. Определим пространство сеточных функций Х т следующим образом:

Х,т = М*) е Улт, * е Ст}.

Под решением сеточной схемы с весами будем понимать сеточную функцию у е Х,т, удовлетворяющую тождествам

)« •

(Ву)г ■ щ ¿х + Е " / Е А('°У(<Т) ' "Г" ' Е / (Ао^)+у(<т) •

Т- КЕТь К

4=1

— Пь|к') ¿ж

= У^ — Су(<т)) • пь ¿ж, V пь е К™, а > 0, (24)

п

У (Ло^) у • пь ¿х = 0, V пь е УГ, * е , У (у - ио) • пь ¿х = 0, V пь е УТ, * = 0.

(25)

(26)

Здесь и далее для средних и разностных отношений будем пользоваться обозначениями. принятыми в [9]:

(*) = + т) + (1 - ст)«(*), «»(*) = («(* + т) - «(*))/т.

Исследуем устойчивость сеточной схемы (24) (26). Для этого определим скалярное произведение в V™ формой:

(V, п) = У V • п ¿х, V V, п е V™.

Тн

Сеточные операторы А и В ^ е X,т, у0 е V™ определим формами

(Лу,»?) = Е "У У У СУ-^йс

к

(Ву, п) = J Ву • пь^ж, (^,п)^У Ж • пь^х, У (уо - ио) • пь ¿х, V V, п е V™.

Тн Тн Тн

Схему (24) запишем в операторном виде

(Ву)» + Ау(<7) = * е Шт, у(0) = уо.

Очевидно, оператор В является самосопряжённым и положительно определён-А

матриц коэффициентов, может быть оценен следующим образом:

(Лу, у) > \ ^ УI [|у|] 1а,, + ^ У |у 1а+ Аг' - С1

з, С1 > 0.

Тогда, следуя [10, с. 62 69], нетрудно доказать, что при значениях весового параметра 1/2 < а < 1 сеточная схема с весами (24)-(26) абсолютно устойчива и верна оценка

Цу(*)11в + Е тЕ/М*)С ¿X + £ т /|у(0|А+ ¿X <

»'£"т ^ »'С" -

< м

1|у(0)ПВ + Е т||г(*')

II-

, м > 0, * е шт

г

е

"г

4. Однородная неконформная аппроксимация

Для решения исходной задачи (13) (15) можно также воспользоваться однородной по пространству и времени сеточной аппроксимацией, предложенной в [6]. Для этого по аналогии с [5, § 6.3.2] построим па отрезке [0, Т] сетку с!т = {¿¿, ] = 0, ... , . Определим пространство сеточных функций

Хь,т = {V е ЬТО[0,Т): V е р(г,,¿¿+1), «(*) е УЛт, 3 =0,1,...,^ — 1}, 1 > 0.

Под решением однородной неконформной сеточной схемы для задачи (13) (15) будем понимать сеточную функцию у е , удовлетворяющую при всех значениях п е тождествам

(Ву • пь)? ¿ж +

Т-

^ кеТ-

-

1+1 п

¿3 к

1+1

+ Е / /(Ло^)+у^(пь —пь|к') ¿ж ¿г

¿3 ек

1

Т,

(Г — Су) • пь ¿ж ¿г, ^ = 0,1,...^ — 1, (27)

¿3 Т-

(Ло^) у • п^ ¿ж & = 0 Vпh е

¿3 Г-

У (у — ио) • пь ¿ж = 0 Vпh е г = 0.

Т-

(28)

(29)

где (> • = ((v • ^)|4,+1-о — -о • +о^/т,, Т, = (г,+ 1 — )-

Отметим, что сеточная схема (27)—(29) является неявной. При 1 = 0 она совпадает с чисто неявной сеточной схемой с весами (а = 1), изученной в § 3..

Характерной особенностью схемы (27) (29) является её консервативность. Действительно, полагая в тождестве (27) п^ = Т,е(я), в = 0,1,...,ш, в(я) = (0,..., 0,1,0,..., 0) е Дт (единица — в позиции в) и суммируя полученные равенства по ], придём к сеточному аналогу закона сохранения:

г г

J Ву о ¿ж + уу(Ло^)+у ¿ж ¿г = у Ву о ¿ж + JJ(Í — Су) ¿ж ¿г, г е ст \{0}.

0Г-

о Т-

Обозначим через В единичную т х т матрицу, через Со = (С + Ст)/2 -симметричную часть матрицы С. Скобками (| • |) обозначим скачок функции по переменной

Воспользовавшись полученными в [6] для этой схемы результатами, придём к следующему утверждению.

Теорема 2. Пусть при всех (ж, г) е О х [0, г] выполнено неравенство

дВ

~дГ

+Е

4=1

дж*

+ 2Со > 2соВ, со > 0,

(30)

¿

3+1

тогда сеточная схема (27) (29) абсолютно устойчива и справедлива оценка

j—1 j _ iiy(tj -o)nb+EKiyD(ti)lB + / E/|[|y(t)iilAvdx+J|y(i)iA+ dx

dt+

£o 2

E dt < M

l|y(o)ll2

||f(t)||E dt

M> 0, j = 1,...,Nt.

Заметим, что условия устойчивости схемы (27) (29) являются более жёсткими по сравнению с условиями для схемы метода прямых и схемы с весами. При I > 0 условие (30) обеспечивает однозначную разрешимость уравнений на каждом временном слое.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проекты Х- 09-01-97015. 10-01-00728).

г

Summary

E.M. Fedotov. Nonconformal Finite Element Schemes for Hyperbolic Linear Systems of Equations.

In this paper we propose a variant of nonconformal finite element method of approximation of the multidimensional linear first order hyperbolic system. The approach is used that was suggested earlier for the scalar convection-diffusion equation, based on Galerkin Pet.rov approximation for the mixed formulation of the original problem, taking into account the direction of convection. Using this approach for the approximation of symmetric systems of equations allows naturally to take into account the local direction of the characteristics, as well as preserve the basic properties of the spatial operator of the original problem.

Unconditional stability of the semidiscret.e scheme, implicit two-layer difference schemes with weights is proved.

Key words: linear hyperbolic systems, mesh schemes, non-conformal finite element methods, difference schemes with weights.

Литература

1. Куликовский А.Г., Погорелой H.B., Семёнов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001. 608 с.

2. Ляшко А.Д., Федотов Е.М. Предельные схемы Галёркипа Петрова для уравнения конвекции-диффузии // Дифферепц. уравнения. 2009. Т. 45, Л' 7. С. 1042 1052.

3. Cockburn В., Shu G.-W. The local discontinuous Galerkin method for time-dependent, convection-diffusion systems // SIAM J. Numer. Anal. 1998. V. 35, No 6. P. 2440 2463.

4. Cockburn B. Discontinuous Galerkin methods for convection dominated problems // HighOrder Methods for Computational Physics. Lecture Notes in Computational Science and Engineering / Eds. T. Bart.li, H. Deconink. Springer-Verlag, 1999. V. 9. P. 69 224.

5. Em A,, Gucrmond J.-L. Theory and Practice of Finite Elements. V. 159 of Applied Mathematical Sciences. N. Y.: Springer-Verlag, 2004. 520 p.

6. Jensen M. Discontinuous Galerkin Methods for Friedriclis Systems with Irregular Solutions: PhD Thesis. Oxford: Oxford University, 2004. URL: http://web.comlab.ox.ac.uk/oucl/research/na/theses.html.

7. Hartmann R. Discontinuous Galerkin methods for compressible Hows: higher order accuracy, error estimation and adaptivity // VKI LS 2006-01: CFD-Higlier Order Discretization Methods , Nov. 14 18, 2005 / Eds. H. Deconinck. M. Ricchiut.o. Rhode Saint Genèse. Belgium: Von Karman Institute for Fluid Dynamics, 2005.

8. Сья'рле. Ф. Метод конечных элементов для эллиптических задач. М.: Мир. 1980. 512 с.

9. Самарский A.A. Теория разностных схем. М.: Мир. 1977. 656 с.

10. Ляшко А.Д., Федотов Е.М. Разностные схемы для нелинейных нестационарных задач. Казань: Изд-во Казан, ун-та, 2008. 199 с.

Поступила в редакцию 18.01.10

Федотов Евгений Михайлович доктор физико-математических паук, профессор кафедры вычислительной математики Казанского государственного университета. Е-шаП: eugeMy.fedotovQksu.ru