CC BY
12
Вычислительные технологии
Том 21, № 4, 2016
Применение разрывного метода Галёркина для решения сингулярно-возмущенных задач
Н.Б. ИткинА, С. И. Марков*
Новосибирский государственный технический университет, Россия *Контактный e-mail:[email protected]
Настоящая работа посвящена применению разрывного метода Галёркина для решения стационарного уравнения диффузии — конвекции — реакции в двумерном случае на прямоугольных конечных элементах. Приводятся результаты решения сингулярно-возмущенных задач при использовании базисов разных порядков, неструктурированных сеток и стабилизирующих членов.
Ключевые слова: разрывный метод Галёркина, сингулярно-возмущенная задача.
Введение
Большинство физических процессов описываются уравнениями в частных производных. Если уравнение содержит малый параметр, то говорят, что в задачу введено возмущение. Если нулевое значение данного параметра приводит к вырождению задачи, то возмущение называют сингулярным, в противном случае — регулярным. В качестве примеров сингулярно-возмущенных задач можно привести уравнение конвекции — диффузии с превалирующим конвективным членом и задачи с пограничным слоем.
Одним из эффективных методов решения данных задач является метод конечных элементов (МКЭ) [1]. Этот матод базируется на вариационных постановках, в которых решение исходной краевой задачи заменяется на минимизацию некоторого функционала. Областью определения этого функционала является гильбертово пространство функций, содержащее в качестве одного из своих элементов решение интересующей нас краевой задачи. Конечно-элементным пространством называют тройку К = {К, Р, Е}, где К С Rm — ограниченное замкнутое подпространство, непустое множество с кусочно-гладкой границей (геометрическое множество элементов), Р — конечномерное подпространство функционального пространства, определенного на K, dim Р = п (пространство функций формы); Е = {Li, L2...Ln} — базис подпространства в пространстве P* (сопряженное пространство к P) функциональных форм (пространство степеней свободы) [2, 3].
Можно выделить две тенденции развития конечно-элементных методов для решения задач конвекции — диффузии: первая — разработка специальных схем с дополнительными стабилизирующими операторами (SUPG, RFB); вторая — использование специальных алгоритмов построения адаптивных сеток. Основные идеи построения неравномерных сеток для конечно-разностных аппроксимаций можно найти в работах Н.С. Ба-хвалова, Г.И. Шишкина, В.Д. Лисейкина. Также существует направление адаптивных методов для построения конечно-элементных сеток (Adaptive FEM).
© ИВТ СО РАН, 2015
В зависимости от аппроксимирующих свойств пространств P и Е выделяют конформные и неконформные конечно-элементные методы. В конформных МКЭ предполагается, что конечно-элементные подпространства принадлежат тем пространствам, в которых определена вариационная задача, а билинейную форму можно вычислить точно на конечно-элементных подпространствах. Конформность аппроксимации в гильбертовом пространстве Нп, п > 1, нарушается, например, в случае неоднозначности следа функции на межэлементной границе.
Разрывный метод Галёркина (далее DG-метод, от английского Discontinuous Galerkin Method) является одним из наиболее перспективных численных методов для решения задач с пограничным слоем и внутренними слоями (см. работы D. Arnold, B. Cocburn, M.G. Larson, F. Brezzi). Оригинальная вычислительная схема на базе разрывного метода Галёркина предложена в 1973 г. Ридом и Ниллом [4] для решения устойчивого стационарного уравнения переноса нейтронов. В отличие от ранее разработанных вычислительных схем для решения физических задач, где непрерывность естественна, разрывный метод Галёркина допускает разрывное представление решения. Такая дополнительная свобода полезна при рассмотрении дифференциальных уравнений, решение которых включает в себя быстро меняющийся градиент или даже разрывы на межэлементной границе.
Отметим, что в параллельной компьютерной архитектуре применение традиционных методов, основанных на межэлементной интерполяции решения, крайне нежелательно, поскольку они обременяют канальную связь между компьютерными узлами и сложно распараллеливаются. Поэтому компактность DG-метода намного лучше подходит к сегодняшним компьютерным архитектурам при решении сложных задач с применением hp-технологий уточнения решения [3-8].
1. Постановка задачи
Рассмотрим трехмерную замкнутую область П С Д3 с границей дП, разделенной на границы с краевыми условиями Дирихле и Неймана дП = Гд иГ^, при этом Гд ПГ^ = 0. Рассмотрим в области П С Д3 краевую задачу
(ХУи) + а ■У и + 7и = ¡, (1)
Л (Уи ■ п)|г^ = , (3)
где и — решение задачи; /, д^, д^ € Ь2 (П); п — внешняя нормаль к дП; а — вектор конвективного переноса; Х,^ € К+ — коэффициенты диффузии и реакции.
2. Вариационная постановка
Сформулируем задачу (1) в виде задачи о седловой точке, введя дополнительную переменную а:
а = ХУи,
{
-V(Xa) + a ■ а + 7и = f. (4)
Введем разбиение области П : 5н = {К} с границей дП = и~ дК, где К — конечный элемент. Определим решение задачи в виде линейной комбинации полиномов на каждом конечном элементе. Для каждого К Е 5н обозначим через т наивысшую локальную степень полиномов. Таким образом, конечно-элементные пространства определим в виде
Ун = { V Е Ь2 (П) : ь\к = Рт(К), deg Рт(К) < т V К Е Ен} , (5)
£н = {т\ т Е [Ь2 (П)]3 : т\к = [Рт(К)]3, deg Рт(К) < т V К Е ~н} . (6)
Решение и внутри каждого конечного элемента представляется непрерывной функцией
п
и(Х,у,г) = ^2 Qг^г(X,У,Z), (7)
г=1
где Фг(ж, у, г) — базисная функция, ассоциированная с г-й степенью свободы конечного элемента. Чтобы связать след и значение функции на границе конечного элемента, введены два оператора — среднего и скачка. Операторы среднего:
М = + Ьу) оп еш,, {г>} = оп еый, (8)
{т} = 2(т 1 + Ту) оп {т} = Т1 оп еый. (9)
Операторы скачка:
[у] = ьп + У]п оп еш, [у] = ьп оп еЬп4, (10)
[т] = т 1 ■ п + т^ ■ п оп еш, [т] = т 1 ■ п оп еый, (11)
где е^ — внутреннее ребро, еьпа — внешнее ребро конечного элемента.
Запишем общую вариационную постановку разрывного метода Галёркина
У \VhUh ■ dП + У (а ■ а) V dП + J ^иь dП+ п п п
А
г
и -uh
■ [Vhv] - {■ М) dS+
+ j ^^ -Uhl [VhW] - ^ ["]) dS = / fV ^ (12) г0 п
В зависимости от вида численных потоков можно получить вычислительную схему с необходимыми свойствами. В работе рассматриваются две постановки DG-метода: NIPG-постановка и постановка в форме Басси.
NIPG-постановка со стабилизирующим параметром Арнольда
Uint = iuh] + n ■ [uh] , иD = n ■ [uh - 9d] , иN = u, n °int = n ■{^VhUh}- a, (K]), n oD = n ■ \VhUh - а, (KD , n °int = 9n,
у \VhUh ■У hvdП + J (а ■Униь) vdП + J с1П+ п п п
+ I (Ы ■ {Жhv} - {ХУьП} ■ М) dS+
Г0 иГд
+ Е / ^ь] ■ М = / (п ■ \Уhv) до ¿Б + (дм у) ¿Б.
^г°иг° { П г{ Г{
Постановка Басси
и™ = , = дв, им = щ °¿п( = {АУь^Ь} - аг ([^]) , сд = АУЬ^Ь - «г (Ы) , п а^ = дм,
J \yhUh ■Уь^ dП + J (а ■УьПь) V dП + J ^и^о dП-п п п
- I (Ы {АУь^} + {ХУьП} [V]) dS-Г°иГд
- ^ [{г (Ы)} ■ М = [ fvdП - I (п ■ АУь^) до dS + ( (дм у) (¡Б.
3. Вычисление стабилизирующих параметров
Стабилизирующие члены в вариационной постановке позволяют добиться устойчивости вычислительной схемы при использовании базисных функций первого порядка. Для базисных функций более высокого порядка применение стабилизирующих членов при решении задач с гладким решением не только бесполезно, но и оказывает отрицательный эффект на точность [9, 10].
Для несимметричной 1Р-постановки стабилизирующий параметр Арнольда а^ ([иь]) определяется в виде
ч(Ф,Ф)= Е М ■ ^ (13)
гег°иап I
VI
где и = —р , р — степень базиса, Щ — длина ребра, п > 0. Щ
Стабилизирующий параметр постановки Басси определяется с помощью лифтинг-оператора [2]
"г ([0]) = -VI {П ([0])} , / г1 ([0]) ■ тdП={
^ [0] ■{т} dS, Г} г (14)
- фп ■ тdS + дп п ■ тdS,
Го Гд
с учетом п. 1;
алгоритм вычисления которого можно записать в виде
к к
1) представить п ([ин]) = ^ Цф ^, ин = ^ д^ф ^;
3=1 3=1
- J К] ■ [т} ¿Б,
2) V I Е Г0 расписать / Г1 ([м^]) ■ тdQ = < г°л п
^ — и^п ■ тйБ + дп ■ т¿Б,
гп гп
3) решить СЛАУ относительно .
Приведем общую рекомендацию по выбору стабилизирующих членов:
1) для Басси- и ШРО-постановок параметр щ выбирается равным числу ребер (граней) или больше;
2) для симметричной 1Р-постановки параметр щ выбирается достаточно большим, индивидуально для каждой задачи.
4. Базисные функции
На одномерном мастер-элементе [0,1], который содержит к степеней свободы, определим функции
X, = ¡г (С), г = !Л (15)
где £ =-на конечном элементе [хк ,Хк+\]. Аналогично для оси Оу
%к+1 — %к
Гз = 9з (гп) , 3 = 1Л
:1б)
Таблица 1. Разрывный базис
Порядок базиса Степень свободы на ребре Степень свободы на грани Вид одномерных функций
1 2 4 1 Х1 = 1—е, Х2 = 2е,
2 3 9 = ^ — ^ (е — 1), = —— 1), *=4 - 9
3 4 16 *=—«—з) («-з)« -1) • 27 / 2\ ^2 = тек — 3ке — 1), 27 / 1\ ^э = тек — 3ке — 1), = — 3) — 3)
Таблица 2. Лагранжев базис
Порядок базиса Степень свободы на ребре Степень свободы на грани Вид одномерных функций
1 2 4 xi = i—е, Х2 = е,
2 3 9 = 2 ^ — ^ (£ — 1), Х2 = —4£(£ — 1), *=4—1)
3 4 16 *=— 2(«—3) («-2)(í—1) • 27 f 2\ Х2 = 27^U — 3) (^ — 1), 27 / Д xa = yde — 3 (е — 1), —3) —3)
У - У к Г 1 гр „
где ц = -на конечном элементе [ук ,Ук+1\. Тогда двумерный элемент с т сте-
Ук+1 - Ук
пенями свободы содержит базисные функции вида
где ^ (i) = (г — 1) modk +1, rq (i)
А = Хф)Yv{i), i =1,т, 'i — 1
+ 1.
:i7)
к
В табл. 1 и 2 приведены виды разрывных и лагранжевых базисных функций для разных порядков.
5. Задача с сингулярным возмущением
Рассмотрим задачу диффузии — реакции
—e2V (Vu(x, у)) + -/u(x, у) = f (х, у) in Q = [0,1] х [0,1],
, , Г 2, у < 0.5, | 7(Х,У)=\ 1, у > 0.5, U\Ш =0,
с аналитическим решением
и(х, у) = 1U(х, у) (sin 4пх + 2) — е—^ — е^ — е— еУ~-~
где U(х,у)
{
2у-1
3 — е 2е , у < 0.5, 1 + е^, у > 0.5.
:is)
Параметр е задает степень сингулярного возмущения. При е < 0.01 вычислительные схемы классического метода конечных элементов становятся полностью неустойчивыми, в то время как разрывный метод Галёркина позволяет получить приемлемое
Рис. 1. Аналитическое решение (а) и правая часть задачи (б) при е = 0.01
решение задачи даже на грубых сетках. На рис. 1 приведены аналитическое решение и правая часть задачи (18) для параметра возмущения е = 0.01. Видно, что решение и правая часть претерпевают разрыв. В табл. 3-7 приведены погрешности решений в норме пространства Ь2, полученных разными методами для разных параметров сингулярного возмущения и на вложенных сетках. Для решения задачи применялась регулярная сетка с шагом дискретизации кх = 0.1, ку = 0.1. Для разрывного метода использовались разрывные базисные функции первого, второго и третьего порядков, для непрерывного метода — лагранжевы базисы тех же порядков. Решатель СЛАУ — Б1СС81аЬ с ЬИ-факторизацией. Замер времени производился от начала сборки глобальной матрицы (с использованием численного интегрирования) до конца решения СЛАУ.
По результатам табл. 3 можно сделать вывод, что вычислительная схема непрерывного метода Галёркина перестает быть устойчивой при коэффициенте диффузии £2 ^ 10-8.
Таблица 3. Решение задачи непрерывным методом Галёркина (лагранжевы базисы)
£ = 10- э £ = 10- 4
Сетка Порядок Погрешность Время, Погрешность Время, Число Размер
базиса решения мс решения мс элементов СЛАУ
к 1 3.95е-2 31 3.78е+3 31 100 121
2 1.27е-2 141 1.13е+3 140 100 441
3 6.39е-3 406 6.77е+2 422 100 961
к/2 1 1.28е-2 78 1.37е+3 93 400 441
2 5.11е-3 547 4.07е+2 563 400 1681
3 2.97е-3 1656 2.39е+2 1672 400 3721
к/4 1 4.75е-3 375 4.99е+2 391 1600 1681
2 1.34е-3 2204 1.48е+2 2219 1600 6561
3 5.72е-4 6578 8.60е+1 6603 1600 14 641
Таблица 4. Решение задачи разрывным методом (стабилизатор Арнольда а^), е = 10 4
^ = 10 - 4 ^ = 0.1 ^ = 10
Сетка Порядок Число Погрешность Время, Погрешность Время, Погрешность Время,
базиса узлов решения мс решения мс решения мс
к 1 400 1.14е-3 265 1.14е-3 260 1.14е-3 266
2 900 8.08е-5 875 8.08е-5 860 8.08е-5 860
3 1600 4.82е-6 2360 4.82е-6 2344 4.82е-6 2344
к/2 1 1600 1.46е-4 1085 1.46е-4 1063 1.45е-4 1063
2 3600 1.61е-5 3570 5.13е-6 3570 5.23е-6 3579
3 6400 1.34е-5 9716 5.85е-7 9687 1.30е-6 9641
к/4 1 6400 1.81е-5 4250 1.81е-5 4250 1.81е-5 4266
2 14 400 4.97е-7 14 406 4.97е-7 14219 1.46е-6 14 438
3 25 600 1.42е-7 38 859 1.42е-7 38 859 1.76е-6 39 407
Таблица 5. Решение задачи разрывным методом (стабилизатор Басси аг), е = 10 4
^ = 10- 4 ^ = 0.1 ^ = 10
Сетка Порядок Число Погрешность Время, Погрешность Время, Погрешность Время,
базиса узлов решения мс решения мс решения мс
к 1 400 Стагнация невязки 3.29е-1 132 1.05е-1 91
2 900 5.02е-1 923 3.93е-1 895 6.19е-2 523
к/2 1 1600 Стагнация невязки 1.62е-1 436 5.77е-2 398
2 3600 1.94е-1 2186 1.21е-1 2128 1.82е-2 1642
к/4 1 6400 Стагнация невязки 9.81е-2 2469 2.74е-2 1984
2 14 400 8.90е-2 4594 8.85е-2 4198 4.73е-2 3749
Таблица 6. Решение задачи разрывным методом (стабилизатор Арнольда а^), е = 10 5
^ = 10- 3 ^ = 1 ^ = 100
Сетка Порядок Число Погрешность Время, Погрешность Время, Погрешность Время,
базиса узлов решения мс решения мс решения мс
к 1 400 Стагнация невязки 3.61е-1 134 1.21е-1 103
2 900 4.58е-1 853 3.14е-1 924 6.01е-2 574
к/2 1 1600 Стагнация невязки 1.81е-1 462 6.03е-2 356
2 3600 1.29е-1 2158 1.03е-1 2205 1.75е-2 1131
к/4 1 6400 Стагнация невязки 9.83е-2 2589 1.98е-2 1736
2 14 400 3.59е-2 4368 3.59е-2 4098 7.78е-3 3562
Таблица 7. Решение задачи разрывным методом (стабилизатор Басси аг), е = 10 5
^ = 10- 3 ^ = 1 ^ = 100
Сетка Порядок Число Погрешность Время, Погрешность Время, Погрешность Время,
базиса узлов решения мс решения мс решения мс
к 1 400 Стагнация невязки 3.15е-1 135 1.19е-1 101
2 900 5.15е-1 923 3.03е-1 903 5.89е-2 506
к/2 1 1600 Стагнация невязки 1.51е-1 431 6.14е-2 402
2 3600 1.89е-1 2201 8.58е-2 2209 1.59е-2 1625
к/4 1 6400 Стагнация невязки 8.18е-2 2426 3.52е-2 1864
2 14 400 8.81е-2 4559 2.49е-2 4264 4.08е-2 3647
По приведенным в табл. 4-7 данным можно сделать вывод, что при решении задачи постановки ШРО и Басси дают схожий результат. Параметр стабилизации влияет на устойчивость вычислительной схемы при использовании базисных функций первого порядка. Неправильный выбор параметра стабилизации, как правило, приводит к стагнации невязки при решении СЛАУ итерационными методами.
На рис. 2 представлены поверхности численных решений с помощью непрерывного (СОРЕМ) и разрывного (БОРЕМ) метода Галёркина. Использованы сетка Нх = 0.1, Ну = 0.1 и полиномы первого порядка. Несмотря на то, что в норме простран-
Рис. 2. Решение задачи с параметром £ = 0.001: а — СОРЕМ; б — БСРЕМ
Рис. 3. Сечение плоскостью у = 0.5: аналитическое решение, численное решение непрерывным и разрывным (стабилизатор Басси) методами, базисы второго порядка
ства Ь2 непрерывный метод показал хорошие аппроксимирующие свойства на грубой сетке для параметра г = 0.001 (см. табл. 3), реальную картину решения с помощью него без специальных технологий уточнения получить невозможно. Так, используя грубое разбиение в непрерывном методе, можно потерять существенную часть решения.
На рис. 3 показано сечение аналитического и численных решений плоскостью у = 0.5. Как можно видеть, разрывный метод позволяет получить не уступающее по точности непрерывному методу решение при использовании достаточно грубого разбиения. Так, сетка с 2500 конечными элементами соответствует СЛАУ с 22 500 неизвестными для разрывного метода (БОРЕМ), а сетка с 40 000 конечными элементами соответствует СЛАУ с 160 801 неизвестным для непрерывного метода (СОРЕМ). Это позволяет уменьшить размерность СЛАУ почти в 7 раз, что ускорит время решения задачи.
6. Задача с пограничным слоем
Задачи, связанные с течением "прилипающих" жидкостей в каналах, задачи аэродинамики, как правило, требуют детальной дискретизации расчетной области вблизи зон с быстроменяющимся решением или его градиента. Такие зоны будем называть пограничным слоем.
Рассмотрим задачу конвекции — диффузии
- V (Чи) + а •Чи = / 2 П = [-1,1] х [-1,1],
и\дп = СОй + у^1 - е^) - е^) , а =(2,1), е< 0.1,
где г — параметр крутизны пограничного слоя. Пограничные слои порядка О(е) находятся на правой и верхней границах области (рис. 4).
В табл. 8 представлены погрешности решений, полученных разными методами для разных параметров пограничного слоя. Для решения задачи использовалась сетка
Рис. 4. Градиент решения задачи с параметром £ = 0.1 (а) и £ = 0.01 (б)
Таблица 8. Погрешности решений задачи, полученные методами ССРЕМ и БСРЕМ в №РС-постановке (оператор а^), 1600 КЭ, базисы первого порядка
Погрешность £ = 0.1 £ = 0.01 £ = 0.001
П — П h и аналит UCGFEM ¿2 6.43e-2 3.54e+0 9.83e+1
Ни — П h II || ^аналит и DGFEMWlß 8.22e-2 9.13e-1 1.19e+0
Таблица 9. Повторные исследования на вложенных сетках (оператор а^), 6400 КЭ
Погрешность £ = 0.1 e = 0.01 e = 0.001
и - иh/2 V аналит u CGFEM L2 3.16e-2 3.98e-1 7.18e+1
и - иh/2 u аналит u DGFEM L2 4.10e-2 4.23e-1 6.22e-1
со следующими параметрами: Нх = 0.025, Ну = 0.025, всего 1600 конечных элементов, размерность СЛАУ для БОРЕМ - 6400, для СОРЕМ - 1681.
Как можно видеть, при параметре пограничного слоя е < 0.01 СЛАУ для непрерывного метода Галёркина не решена, фактически останов алгоритма происходит по максимальному числу итераций. Для е > 0.1 разрывный подход уступает по своей точности непрерывному. Это обусловлено главным образом накоплением машинной ошибки в связи с большим размером конечно-элементной матрицы СЛАУ1.
Покажем устойчивость методов на вложенных сетках. В табл. 9 представлены погрешности решений, полученных разными методами для разных параметров пограничного слоя на вложенных сетках. Из таблицы следует, что порядки аппроксимации (как и порядки сходимости) для методов с линейными базисными функциями для параметра £ = 0.1 совпадают,
log2
\U —Пh II
I Vаналит и DGFEM\\i,2
Ua
-U'.
h/2
DGFEM
1 log2
L2
\U —TJh II
I ^аналит uCGFEM\\i,2
Ua
-u;
h/2
CGFEM
L2
что соответствует теоретическим сведениям о порядке аппроксимации решения разрывным и непрерывным методом Галёркина.
Варьируя параметры стабилизации, можно получить достаточно устойчивую вычислительную схему БО-метода. В общем случае выбор параметра постановки Басси совпадает с выбором параметра для ШРО-постановки. Оптимальным выбором для параметра стабилизатора является выражение [6]
Vi
Уе'рУ? ||а|| h
:i9)
где — число ребер конечного элемента в двумерном случае или число граней в трехмерном случае; р — порядок базиса; А, 7, ||а|| — параметры уравнения; Не — длина ребра в двумерном случае или площадь грани в трехмерном случае. Влияние параметра стабилизации в постановке Басси показано в табл. 10.
1
1 Размерность матрицы СЛАУ для БО-метода в 4 раза больше, чем для ОС-метода при использовании линейных базисных функций на носителе.
Таблица 10. Влияние параметра стабилизации щ (оператор ar), 6400 КЭ, £ = 0.001
Vi TT — TTh и аналит и DG FEM L2
0.01 СЛАУ не решена
0.05 СЛАУ не решена
0.10 1.64e+0
0.72 6.22e-1
1.00 7.26e+0
1.50 СЛАУ не решена
2.00 СЛАУ не решена
Необходимо обратить внимание на высокую чувствительность вычислительной схемы к выбору параметра стабилизации. На линейных базисных функциях и неудачном выборе параметра стабилизации итерационные методы решения СЛАУ стагнируют.
На рис. 5 показано сечение плоскостью х = 0.1 полученного решения для параметра £ = 0.001 на грубых и вложенных сетках. Решения, полученные разрывным и непрерывным методом Галёркина, при дроблении сетки в пределе совпадают. Таким образом, необходимо сделать главный вывод: разрывный метод Галёркина позволяет находить решение задачи на более грубых сетках в отличие от непрерывного метода. На практике это выражается в уменьшении времени, которое тратится на решение СЛАУ.
Кроме того, возможность применять несогласованные грубые сетки в разрывном методе Галёркина позволяет сокращать размер матрицы СЛАУ в разы. При решении задачи с пограничным слоем установлено, что пограничные слои порядка О(е) находятся на правой и верхней границах области, наглядно пограничные слои можно рассмотреть на рис. 4. Целесообразно произвести адаптацию сетки (с применением Д-технологии уточнения решения) так, чтобы максимально подробно описать зоны сингулярного возмущения, исключая область, где решение ведет себя без резких осцилляций. Пример такой сетки приведен на рис. 6. В результате получена сетка с 1700 конечными элементами.
Рис. 5. Сечение аналитического и численных решений плоскостью х = 0.1
Рис. 6. Пример использования ^-технологии
Таблица 11. Вычислительные затраты при решении задачи с параметром слоя е = 0.001 классическим (СОРЕМ) и разрывным (БОРЕМ) методами
Метод решения Число конечных элементов Число неизвестных при решении СЛАУ Время решения СЛАУ (ВСО с ЬИ-факторизацией), с Время решения СЛАУ (вСОБ1аЬ с ЬИ-факторизацией), с Погрешность
СОРЕМ 12100 12 321 1.51 1.53 6.44е-1
БОРЕМ 1700 6800 0.83 0.83 6.03е-1
Таким образом, удалось сократить размер матрицы СЛАУ в 4 раза, при этом погрешность полученного решения не изменилась. Результаты исследования вычислительных затрат представлены в табл. 11.
7. Влияние сеточного числа Пекле
-, где a — скорость
Введем в рассмотрение специальное сеточное число Пекле Ре = конвекции; Но, — равномерный шаг по пространственной сетке; И — параметр диффе-
|a| hp 2D
Рис. 7. Влияние числа Пекле на погрешность а решений задачи с параметром е = 0.001
ренциального уравнения (теплопроводность, коэффициент диффузии и т.д.). Установлено, что для числа Пекле Pe > 105 прогнозирование поведения решения непрерывным методом Галёркина невозможно в силу нарушения устойчивости вычислительной схемы классического метода. Поэтому для ситуации D << 1 требуется адаптация сетки — измельчение шага по пространству, hn ^ 0. На рис. 7 приведены графики влияния числа Пекле на сходимость решений задачи с пограничным слоем (для разрывного метода использован лифтинг-оператор в NIPG-постановке).
Исходя из результатов решения задачи с пограничным слоем можно утверждать, что с помощью разрывного метода Галёркина можно получить более точное решение задачи при использовании достаточно грубых сеток, что является определяющим критерием эффективности метода при моделировании реальных физических процессов.
Заключение
По результатам проведенных исследований можно сделать вывод, что высокая гибкость разрывного метода Галёркина позволяет применять данную вычислительную схему для решения задач с сингулярными возмущениями и пограничными слоями. Установлено, что при оптимальном подборе параметра стабилизации можно эффективно применить h-технологию уточнения решения, тем самым снизив размерность матрицы СЛАУ в несколько раз. Благодаря локальной консервативности все вычислительные схемы на базе разрывного метода Галёркина дают адекватное представление решения задачи на достаточно грубых тесселяциях расчетной области.
Основные недостатки вычислительной схемы заключаются в резком увеличении размерности матрицы СЛАУ при использовании базисных функций высокого порядка и высокой чувствительности вычислительной схемы к выбору стабилизирующих параметров.
Список литературы / References
[1] Cocburn, B. Superconvergence of the numerical traces of discontinuous Galerkin and hybridized mixed methods for convection-diffusion problems in one space dimension // Math. Comput. 2007. No. 76. P. 67-96.
[2] Brenner, S.C., Scott, L.R. The mathematical theory of finite element methods. Texts in Applied Mathematics. Vol. 15: 3rd ed. Springer, 2007. 420 p.
[3] Solin, P., Segeth, K., Dolezel, I. High-order finite element methods. Hardcover: Chapman & Hall/CRC Press. 2003. 408 p.
[4] Arnold, D.N., Brezzi, F., Cocburn, B., Marini D. Unified analysis of discontinuous Galerkin methods for elliptic problems // SIAM J. Numer. Analusis. 2002. Vol. 39, No. 5. P. 1749-1779.
[5] Baumann, C.E., Oden, J.T. A discontinuous hp finite element method for convection-diffusion problems // Comput. Methods Appl. Mech. Eng. 2000. Vol. 175. P. 311-341.
[6] Oden, J.T., Babuska, I., Baumann, C.E. A discontinuous hp finite element method for diffusion problems // J. of Comput. Physics. 1998. No. 146. P. 491-519.
[7] Brezzi, F., Cockburn, B., Marini, L.D., Suli, E. Stabilization mechanisms in Discontinuous Galerkin finite element methods.
Available at: http://eprints.maths.ox.ac.uk/1162/1/NA-04-24.pdf.
[8] Sudirham, J.J., van der Vegt, J.J.W., van Damme, R.M.J. A study on discontinuous Galerkin finite element methods for elliptic problems: Memorandum No. 1690. Mathematics of Computational Science (MaCS). 2003. 20 p.
Available at: http://purl.utwente.nl/publications/65875
[9] Gilbarg, D., Trudinger, N.S. Elliptic Partial Differential Equations of Second Order. Berlin; New York: Springer, 2001. 500 p.
[10] Sangalii, G. Capturing small scales in elliptic problems using a residual free bubbles finite element method // Multiscale Model Simul. 2003. Vol. 1. P. 485-503.
Поступила в 'редакцию 2 декабря 2015 г., с доработки — 23 декабря 2015 г.
Discontinuous Galerkin Method for solution of singularly perturbed problems
Itkina, Natalya B., Markov, Sergey I.*
Novosibirsk State Technical University, Novosibirsk, 630073, Russia * Corresponding author: Markov, Sergey I., e-mail: [email protected]
The objective of our paper is to develop and verify a stable computational finite-element scheme for solution of singularly perturbed problems.
In this paper we present a nonconformal finite element method (the Discontinuous Galerkin Method) using special lifting-operators that enhance the stability of the computational scheme.
The computational scheme with special stabilizing terms and lifting-operators is proposed. The verification of the computational scheme is carried out in the model problem class. The optimal change range of the Bassi and Arnold stabilizers was defined. The efficacy of the ^-refinement technology application for solving the singularly perturbed problems was substantiated.
A high flexibility Discontinuous Galerkin Method allows you to apply this computational scheme for solving the problems with singular perturbations and boundary layers. It was found that the optimal choice of stabilization parameter can effectively apply the ^-refinement technology to find solutions and significantly reduce the dimension of the matrix linear systems. Because of the local conservatism, all the computing circuits based on the Discontinuous Galerkin Method give an adequate idea for solving the problem on a fairly coarse tessellation computational domain.
A sharp increase in the dimension of the SLAE matrix using the high-order basis functions and a high computational scheme sensitivity to the choice the stabilizing parameters is the main disadvantages of the computational scheme.
Keywords: the Discontinuous Galerkin Method, singularly perturbed problem purpose.
Received 2 December 2015 Received in revised form 23 December 2015
© ICT SB RAS, 2016
