CC BY
95
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
Том XLIV
2013
№ 3
УДК 519.63
ВЫБОР МЕТОДА АППРОКСИМАЦИИ ВЯЗКИХ ЧЛЕНОВ В МЕТОДЕ ГАЛЕРКИНА С РАЗРЫВНЫМИ БАЗИСНЫМИ ФУНКЦИЯМИ
Приведено описание неявного численного метода Галеркина с разрывными базисными функциями для решения двумерного модельного скалярного уравнения конвекции-диффузии. Рассмотрены три метода аппроксимации диффузионных членов — Bassi & Rebay 1, Cockbum & Shu и Bassi & Rebay 2. Описаны результаты тестовых расчетов, произведена оценка порядков точности. Проанализировано влияние разрывов граничных условий на скорость сходимости к точному решению. Даны рекомендации по выбору метода аппроксимации диффузионных членов.
Ключевые слова: метод Галеркина с разрывными базисными функциями, уравнение конвекции-диффузии, порядок точности.
Метод Галеркина с разрывными базисными функциями (РМГ) [1 — 5] является одним из наиболее перспективных подходов к конструированию схем для численного решения задач аэродинамики. Он сочетает в себе лучшие свойства методов конечного объема (возможность учета направления распространения информации при аппроксимации потоков на гранях ячеек и возможность сквозного счета разрывов) с таким важным свойством метода конечных элементов, как слабая зависимость от вида используемой расчетной сетки. Особенно привлекательным свойством метода РМГ является возможность построения схем произвольного порядка точности на компактном шаблоне, содержащем лишь текущую ячейку и ее ближайших соседей.
При решении уравнений Навье — Стокса и Рейнольдса возникает необходимость построения РМГ-аппроксимации диффузионных потоков, представляющих собой произведение коэффициента диффузии на производные от параметров газа (скорости, температуры и др.). РМГ-аппроксимация
В. В. ВЛАСЕНКО, А. В. ВОЛКОВ, А. И. ТРОШИН
1. ВВЕДЕНИЕ
ВЛАСЕНКО
ВОЛКОВ Андрей Викторович
доктор физико-математических наук, начальник отделения ЦАГИ
ТРОШИН
Владимир Викторович
кандидат физико-математических наук, начальник сектора ЦАГИ
Алексей Игоревич
инженер ЦАГИ
градиентов должна обеспечить заданный порядок аппроксимации уравнений и, по возможности, сохранить компактность шаблона схемы. Известно, что в методе Галеркина с разрывными базисными функциями использование градиентов, вычисленных непосредственным дифференцированием решения в ячейках, в диффузионных потоках приводит к потере аппроксимации [1]. Поэтому распространение получили методы, в которых градиенты (или связанные с ними величины) считаются дополнительными переменными и определяются путем решения дополнительных дифференциальных уравнений [4]. Первоначально наибольшую известность получил метод аппроксимации, предложенный F. Bassi и S. Rebay [5], в котором диффузионные потоки на гранях вычисляются по полусуммам значений переменных и их градиентов с двух сторон грани. Несмотря на простоту и физичность, этот подход имеет существенные недостатки — некомпактность и возможность возникновения четно-нечетного расщепления в решении [6]. Позднее F. Bassi и S. Rebay предложили другой вариант аппроксимации вязких членов [7]. Эта аппроксимация симметрична и, несмотря на более сложную формулировку, считается одной из перспективных [4]. B. Cockburn и C.-W. Shu предложили альтернативную аппроксимацию [8]. В ней при вычислении потоков значения переменных берутся с одной стороны грани, а их градиентов — с другой. Этот метод не имеет недостатков аппроксимации [5], но его конечные формулы несимметричны. В настоящей работе проведено сравнение трех перечисленных методов применительно к двумерному скалярному модельному уравнению и предложены рекомендации по их использованию при решении стационарных задач неявным методом РМГ.
Структура статьи следующая: в разделе 2 описан общий вид численного метода, аппроксимация явного и неявного операторов и метод решения системы линейных уравнений. В разделе 3 приведены результаты расчетов тестовых задач — одномерного скачка уплотнения, размытого вязкостью, задачи Вигтона [9] об искривленном слое смешения (в двух вариантах — с гладкими и разрывными граничными условиями) и модельного пограничного слоя на пластине (также в двух вариантах — с гладким решением и с решением, имеющим особенность). Сделанные выводы перечислены в заключительном разделе.
2. ЧИСЛЕННЫЙ МЕТОД
Рассматривается двумерное модельное скалярное уравнение для поля u (x, t):
du
—+ V-F(u,G) = 0. (1)
dt
Здесь F(u, G) = Fconv(u) + Fdiff(G) — сумма векторов конвективных и диффузионных потоков; G = Vu — градиент u. В тестовых задачах, которые будут описаны ниже, используются различные варианты векторов конвективных потоков. При построении численного метода будет рассматриваться вектор диффузионных потоков вида
Fdiff(G) = -v G.
В каждой ячейке расчетной сетки определяются локальные полиномиальные базисные функции фу (x), линейная комбинация которых аппроксимирует решение в ячейке:
Kf
u = z Vj ф j(x).
у=1
Коэффициенты разложения qj являются основными неизвестными величинами в методе РМГ. В качестве базисных функций используются полиномы вида фу-(х,у) = Xа3Ув3 , где аи
—х х ) —( )
Р — неотрицательные целые числа; X(х) = —Х——; У(у) = —У—Ус— ; (хс, ус) — координаты
Нх Ну
центра ячейки; кх — ширина ячейки; Ну — высота ячейки. Степень полинома (а 3 +р 3) меняется от 0 до К. Базисные функции должны образовывать полный набор, т. е. любой полином степе-
ни не выше К должен однозначно выражаться в виде линейной комбинации базисных функций. Число базисных функций равно К^ . В двумерном случае величина К^ связана с максимальной
(К + 1)(К + 2)
степенью базисного полинома К соотношением К^ =
2
Умножим уравнение (1) на базисную функцию ф. (х) и проинтегрируем по объему ячейки О:
+ У-Е(и, С) ]ф. (х)й О = 0.
О
Применив преобразование (У-Е )ф = У-( Е ф)-Е Уф и формулу Гаусса — Остроградского, получим уравнение
{—ф. (х)йО +| (Е - й" ф. (х)йЕ = { (Е - Уф. (х)) йО.
(2)
Здесь Е — поверхность ячейки. Интеграл ^ (Е - й^ф. (х)йЕ — это поток через поверхность
Е
ячейки; {(Е - Уф. (х)) й О — источниковый член, возникающий в методе РМГ из-за неоднородно-
О
сти базисных функций по пространству.
В настоящей работе строится неявный метод РМГ для получения стационарных решений модельного уравнения вида (1), основанный на тех же идеях, что и неявный метод конечного объема из работы [10]. Поскольку детали нестационарного процесса не представляют интереса, достаточно ограничиться первым порядком точности по времени. Для получения абсолютно устойчивой схемы потоки через поверхность аппроксимируются на неизвестном временном слое (п + 1). При этом они линеаризуются и представляются в виде суммы значения на известном временном слое п и приращения:
ЕС0ПУ(МП+1) Ес0пу(ип) +
дЕс
ди
Аи,
ЕЙЖ(Сп+1) = -УСп+1 = -УСп ^АС = ЕЙЖ(Сп) + ЕЙШ'(АС).
Здесь Аи = ип+1 - ип , АС = Сп+1 - Сп . Во всех вариантах аппроксимации градиентов, которые будут рассмотрены ниже, АС линейно выражается через значения Аи в текущей ячейке и некоторых ячейках вблизи нее.
Значения потоков на известном временном слое п аппроксимируются с той точностью, которая требуется от схемы. Неявная часть, связанная с приращениями, исчезает при приближении к стационарному решению и может быть аппроксимирована с пониженной точностью. Она выполняет роль «неявного сглаживателя» [10], который делает схему абсолютно устойчивой.
В источниковом члене {(Е-Уф.(х))йО используем то же представление Е(и,С), что и
О
в потоках, что дает полностью неявную аппроксимацию 1-го порядка по времени. Тогда уравнение (2) переписывается в виде
{—ф. (х)й О+|
О Т Е
= {
ЕсоШ(ип )-
дЕс
ди
Аи + Еда(С п) + ЕЙЖ(АС)
й ф. (х)й Е =
(ип )-
дЕс
ди
Аи + Еда(Сп) + Е да(АС)
-Уф. (х)й О,
где т — величина шага по времени.
О
Введем обозначения:
Kf
Au =ХЛ<?у Фj (x),
j=1
' dFconv \n y du j
Fnconv = Fconv • n, F°onv = Fconv ^ф-, An = A • n , A = A • Vфг, Fdiff = Fdiff • n , F;difr = Fdiff • Vфг, AGn =AG n , ЛОг = AG • Vcp;, Fn = Fnconv + Fndff , Ft = Fconv + Fdiff, MQ = J ф. (x)py (x)dQ,
Q
Су = J A (un)фу (x)dQ, Rur (un, Gn) = -J ф. (x)Fn (un, Gn) d2 +J F (un, Gn)dQ.
Q 2 Q
Здесь n — вектор внешней нормали к грани ячейки. Матрица MQ, составленная из элементов MQQ, далее будет называться «массовой матрицей».
В новых обозначениях схема (3) запишется в следующем виде:
1 Kf . Kf -ZMQAqj + £ J ф(x)AnAud2-£СуAqj -
T J=1 f 2f J=1 (4)
-vZ J ф- (x)AGnd2 + v J AG.dQ = Ru. (un,Gn).
f 2f Q
Потоки и их приращения суммируются по всем граням f текущей ячейки, 2f — поверхность грани f.
В правую часть уравнения (4) входят интегралы по граням ячеек J ф. (x)Fnconv(un) d2, где
2f
ф. (x) — это базисные функции для текущей ячейки. На каждой грани имеется два набора параметров: (uL, G L) — решения в ячейке «слева» от грани и (ur , G r ) — решения в ячейке «справа» от грани. Одна из этих двух ячеек — текущая, другая — соседняя. По этим данным для конвективных потоков решается точная задача о распаде произвольного разрыва:
((onv(un ) )f = Fconv (u^jecay ) , udnecay = u^y^, uR, n).
Интегралы по объему ячейки J F°onv(un )d Q берутся непосредственно по значениям un во
Q
внутренних точках ячейки.
Рассмотрим три метода аппроксимации градиентов G , используемых в диффузионных потоках Fdiff(Gn). Первые два — Bassi & Rebay 1 и Cockburn & Shu — заключаются в решении дополнительного дифференциального уравнения
G -Vu = 0.
Умножим его на ф. (x), проинтегрируем по объему ячейки Q и применим формулу Гаусса — Остроградского. В результате получим
cFxconv du
' dFconv \r>
U1 y
du
JGифг(x)dQ-^ J un фг(x)ndE + Jun Уфг(x)dQ = 0. (5)
Q f Ef Q
В силу линейности уравнения (5) по G и u, оно остается верным при замене Gn и un на AG и Au соответственно.
В обоих методах интегралы по объему JunVфг■ (x)dQ и J^dlff(Gn)dQ вычисляются по
значениям un и Gn во внутренних точках ячейки. Различия возникают в интегралах по поверхности J (x)ndE и J фг- (x)^ndlff(Gn)dE: в подходе Bassi & Rebay 1 значения un и Gn на гра-
Ef Ef
uL + uR GG L + GG —
ни берутся равными L — и —^—— соответственно, а в подходе Cockburn & Shu — ur[ и
Gn
R.
В третьем методе аппроксимации — Bassi & Rebay 2 — градиенты Gn в ячейке представляются в следующем виде:
Gn = (G£ )n R f, (6)
f
где СЕ — градиенты, полученные непосредственным дифференцированием поля ип в ячейке, а К/ — вклад грани /в поправку к СЕ , связанный с наличием разрыва ип на грани. Для R/ записываются интегральные уравнения
J R f ф (x)dQ = J (u4 - un, )ф (x) ndE, (7)
Q Ef
_ 1Л Т I 1Л п
где и у- на внутренних гранях равно -, а на границе расчетной области — иоаХ; и^ — ре-
и" + ип 2
шение в текущей ячейке; иоаХ — значение и, заданное на границе расчетной области. Как и в
уравнении (5), линейное уравнение (7) остается верным при замене Сп и ип на АС и Аи соответственно.
В подходе Ба881 & Rebay 2 в интеграл { ^^ (Сп О входят градиенты (6), а в интегралы
О
{ ф. (х)^^ (Сп) йЕ — полусумма «неполных» градиентов, содержащих поправку только для
Е/
текущей грани /
С п ( (С Е )п + К / )т +( (С Е )п + С / =-2-. (8)
Для аппроксимации приращений конвективных потоков { ф. (х)АпАийЕ используется реЕ/
шение задачи Римана по Роу, согласно которому на грани /
\Аиг, Ап > 0,
(Аи) / Т п (9)
}/ [АиЙ, Ап < 0.
Величина An вычисляется по осредненному значению u =
uL + u—
Рассмотрим аппроксимацию приращений диффузионных потоков в методе Bassi & Rebay 1. Разложим AG по базисным функциям в ячейке:
Kf
AG = Z Ag j ф j (x).
j=1
Подстановка разложений Au и AG в уравнение (5) дает
Kf
Kf f
ZMQAgj =Z Z-
j=1 k=1^ f
f I ^ Mf Aq{,L + mZr Aq{
Mk Mk
Aqk
Здесь M-L = J ф. (x)фfj,L (x)d2 , ф. (x) — i-я базисная функция текущей ячейки; ф^ (x)
2f
f L
j-я базисная функция ячейки, находящейся слева от граниf(индекс «L»), Aqj' — j-й коэффициент разложения в ячейке слева от грани f; аналогичные обозначения используются для величин с индексами «R» (M—R и Aqf,R). Одна из ячеек « f,L » и « f,R » является текущей. Кроме этого,
д (х) д х
Мх = [—--ф. (х)ёО и МУ = [—--ф. (х)ёО . Таким образом, Ag-■ выражается через коэф-
^дх 1 ду
О О ^
фициенты Aq текущей ячейки и ее соседей:
к - (
K
Ag. =Z j=1
(M Q) z
f - M%,l Aqf L + M f R Aq[,R
k=1
Z-
f
MXk Mjk
Aqk
(10)
Запишем теперь интеграл | ф- (x)AGndЕ , входящий в схему (4), через компоненты АО :
2f
J ф. (x)AGnd2 = J ф. (x)AG • n fd2 =
2f
2f
f <Z ф^,L (x)Agf,L +фfj 'R (x)Agf,R 1 Z l^jf A fL \/f2f A fR
J ф-j—^-— • nfd2 = -Е|_М&Agf + M2RAgj
2f j=1 и далее, с помощью (10):
■j=1
• nf,
1 Kf
J ф.(x)AGnd2 = -Z[M=jfLAgf,L + mIrAgf,R
2f 2 j =1
Kf
=Z
j=1
mil z
k=1
Kf f
((MQl ) Z Z
j l=1 ^ f'(L)
MkLLAq(',L + M2RAqfi ',R
Mx му
A
f ,L
+-
мf z
'j,R
2
k=1
Kf f
((qr ) Z Z
f ^ MrfL V'L + MjfR V
', R
Mx
Mykl
Aq,
f, R
•n f.
(11)
Здесь (mql ) и (mqr ) — обратные массовые матрицы в ячейках « f,L » и « f,R ».
>}к V 1 ,Я );к
Знаками ^ и ^ обозначено суммирование по граням /' ячеек « /, X » и « /, Л » соответст-
/'(¿) /'(Я)
венно.
В схему (4) входит сумма выражений (11) по всем граням / текущей ячейки. Как видно, аппроксимация Ба881 & ЯеЬау 1 имеет некомпактный шаблон: в (4) коэффициенты Ад текущей ячейки оказываются связанными не только с коэффициентами Ад ее соседей, но и с соседями ее соседей. Это значительно усложняет реализацию численного метода. Запишем теперь через Ад интеграл :
Kf
J AGidQ = J Уфг (x) • AGdQ = J Уф (x) • Z Ag y фу (x)dQ
Kf
=z
j=1
MX
My
ij
Kf
■z
k=1
Kf (
(M z z
Q
f - MkL Aq{,L + MkL Aqf'R
f
j=1
kl,R Aq
f
Mkl МУ
Aq
(12)
Итак, мы выразили приращения диффузионных членов J фг- (x)AGndЕ и J AGidQ через
If Q
векторы Aq в текущей ячейке (пусть ее индекс (i, j)) и ячейках (i ± 1, j), (i, j ± 1), (i ± 2, j), (i, j ± 2) . Вектор Aq в «заграничных» ячейках считается нулевым.
Повторим эту процедуру для случая аппроксимации Cockburn & Shu. Подстановка разложений Au и AG в уравнение (5) дает
Kf
Kf (
Z MQ Ag j =z Z Mf Aqk ,Ln f -j=1
к=1^ f
MX МУ
Aqk
откуда
Kf
Agi =Z j=1
(M Qt Z
Kf (
k=1
ZMXl Aqk,Lnf -
MXk M>k
Aqk
Интеграл J фг- (x)AGndЕ запишется в виде:
If
If
If
Kf
=z
j=1
K
j •z
f k=1
If _1 Kf (
Kf
)Z<
j=1
((Q,R )) Z Z. m£Aqf',Lnf
i=1V f'(R)
Kf
•n fd i=z
j=1
Mkl' Aqf ,R
My _
•n f
(13)
а интеграл J AGidQ — в виде:
Q
Kf
J AGidQ = J Уф- (x) • AGdQ = J Уф. (x) • Z Agx Фx (x)dQ =
j=1
Kf
=Z
j=1
MX му
Kf
■Z
k=1
Kf (
(M Q)-k Z ZMf Aqf,Ln f
Mki муш
4
По сравнению с методом Bassi & Rebay 1 метод Cockburn & Shu требует значительно меньшего количества операций и сохраняет компактность схемы. Однако у него есть и недостаток — несимметричность формулировки.
Рассмотрим теперь метод Bassi & Rebay 2. Разложим по базисным функциям величину
AG * =(G * )1 -(G * ) (приращение за шаг по времени градиентов, получаемых непосредст-
венным дифференцированием поля un в ячейке):
Kf
AG * =ZAge j ф j (x).
j=1
Заметим, что коэффициенты Age j можно выразить через коэффициенты разложения Au с помощью постоянных матриц:
Kf
Age j = Z
k=1
A4k,
так как дифференцирование полинома дает полином с коэффициентами, линейно связанными с коэффициентами исходного полинома.
Аналогично, разложим по базисным функциям величину АЯ/ = (к/ ) — (/ ) (приращение за шаг по времени вклада грани/ в поправку к градиентам ОЕ ):
к/
Ак/ =ЕАг/,3 ф3(х).
3=1
Подставим это разложение в уравнение (7):
K
f,OA KfM^out Aqj,out -a fMZ A%
Z Arf, j = Z-2-:
j=i k=i 2
(15)
где МЕои1 = | ф- (х)ф/^, ои1(х)а?Е — интеграл, в который входит --я базисная функция текущей Е/
ячейки и 3-я базисная функция соседней ячейки со стороны грани / М-/с = | ф- (х)фу (x)dЕ —
Е
интеграл, в который входят базисные функции только текущей ячейки. На внутренних гранях а/ = 1, на границах расчетной области а/ = — и Aqj Из(15)следует, что
7f,°ut = о.
Kf
Arf ,г =Z j=1
(M Q)-1 ZM Ж out Aqf,0ut-a fM f Aqk
k=1
Интеграл | ф- (x)AGndЕ (приращение диффузионного потока через грань /) запишется
в виде:
Ef
J ф. (x)AGndЕ = J Фг (x)(AG* + ARf )• nfdE =
Ef Ef
% ф?,ь (x) (Agef,L +Arffj ) + f (x) (Agef,R + Ar/f)
: J ф (x)Z—-- 2---^^
Ef j=1 2
• nfdE =
Kf
1 Kf
=1 Z
2 j=1
Kf
Mf Z
2Z m|,'l (Agef,L + Ar^j) + mIr (gef r + Дгfj ) (f ) jk
+Mf Z
fR l Q \-1 KfME(LR)Aqf,L -afME(RR)
Aq{,R + ((M Q,r ) ^ Z -2 f kl ,R
•n f.
1 к/ Ме/(X) Ад/,Я _а МЕ/(X) Ад/,х
л , Л/о Г^ Мк1,Я Ад/ а/Мк/Ад/
Г« [[(/), Гк ^12 М
[С*) 1Ад/-+() ^
к-1 _[(/) ] ]к1-1
Аг/, Аг/, ^^ , , ^/я , — это коэффициенты Аги матрицы Вх, в ячейках « /,£ » и « /,Я » соответственно. Также введены обозначения:
М/Ь) - | ф/,Х (х)ф/,Х (х^Е, М/- | ф/,х (х)Ф/,Я (х^Е ,
Е/ Е/
м/Я) - | ф/,Я (х)ф/,х (х)^Е, м/Я) - | ф/,Я (х)ф/,Я (х^Е .
Е/ Е/
Интеграл | АО^О (приращение источникового члена, связанного с диффузией) примет вид:
J AGidQ = J Уф. (x) • AGE + Z AR
Л
Kf (
f
f / Л
dQ =
JVфг (x) •Z Age j +Z Ar
Kf
=2:
j=1
ij
My
j
Kf
Z
k=1
DXk
j=1\ (
ф j (x) d Q =
(17)
Aqk+Z (м Q))1 Z
f V 1=1
1 K^M?out V,out -a fMlc Aqi
В отличие от схем Bassi & Rebay 1 и Cockburn & Shu полученная аппроксимация (16), (17) симметрична и компактна.
Последовательная подстановка в (4) всех базисных функций ф. (x) дает замкнутую систему из Kf уравнений относительно приращений коэффициентов Aqy-. Неизвестные величины в текущей ячейке связаны с неизвестными величинами в соседних ячейках (и, возможно, в более далеких). Эта связь осуществляется через приращения конвективных и диффузионных потоков.
Приведенные выше формулы позволяют записать (4) в общем виде:
Kf Г ( Л
Z-j [MQ -т Сij J Aq j + т Mpin Aq j +Z Mj™'out Aqf j=1 I V
out
-т
Л
out AqOUt
out / V
Mf inAq} + ZMdlff'out Aqout + т jin Aq} + Z jout Aq'
f
w
(18)
-idiff out Kqout
out /
= т Ru. (un, G n),
где
(
( \
Mco„v,m +Ц Mco„v,out д^ои1 out
Ч
J
это аппроксимация I J Фi (х)AnkudЕ с подстановкой (9),
f Ef
\
Mf Дц} +Z Mdlil'out Ддо
V
out
аппроксимация vI J Фi (х)ДС^Е (формулы (11), (13) и (16),
f Ef
просуммированные по граням f текущей ячейки),
С d
1 д^з +i cdifi;out ддо
out
аппрокси-
мация vJдGidQ, (формулы (12), (14), (17)). В каждое из этих слагаемых входит некоторая мат-
Q
рица (с индексом «in»), умноженная на вектор дд в текущей ячейке, и сумма матриц (с индексами «out»), умноженных на векторы дд в других ячейках, участвующих в аппроксимации. В методах Cockburn & Shu и Bassi & Rebay 2 индекс «out» относится только к соседним ячейкам, а в методе Bassi & Rebay 1 — еще и к ячейкам, расположенным вокруг соседних ячеек (за исключением текущей).
Введем для каждой ячейки с векторы размерностью Kf :
Aq(c) = (д^;...; дд^) , Ru(c) = (R^;...;RuKf) , и матрицы D и Hout размером Kf х Kf :
г\ / \ hjtO г* , a ^conv, in ъ ^diff, in , ^-idiff, in
Dij (c) = M°0 -тСз +хМз , -мз +тСз ,
conv, out — mdiff, out + сdiff, out
ij -Mj +Cij
Hout(c) = -(M, Тогда система (4) приводится к виду:
Ru(c) + £ H out(c)Aq (cout(c))
D(c)Aq(c) = т
out
где сои1(с) — любая из ячеек, участвующих в схеме (18), кроме текущей ячейки с. Система уравнений для всей расчетной области запишется в виде
D(c1)Aq(c1) = т
Ru(c1) + £ H out(c1)Aq (cout(c1))
out
D(C^Cells)Aq(C^ce„s) = T
cells >
^(CNcells ) + IHout (CNcells )^ (out ^ )) out
где — число ячеек в расчетной области.
Как и в неявном методе из работы [10], применим для решения этой системы итерационный поблочный метод Гаусса — Зейделя [11] с ограниченным числом итераций. Обозначим нижним индексом к номер текущей итерации метода Гаусса — Зейделя. Тогда
D(ci)(Aq(Cl)) = т Ru(Cl) + £Hout(Cl)(((c0Ut(Cl)))
out
k,k+1
D(CNcells ) ((^ ))+i = T F^U(CNce[ls ) + XH0Ut (CNce[ls ) (( (Cout (%се1к )))
out
к ,к+1
В этой системе уравнение с номером т используется для определения нового итерационного значения коэффициентов в ячейке ст — (Ад(ст ^. ^ Н ааХ(ст) (Ад (соиХ(ет))) ^ — сумма
out
по другим ячейкам, участвующим в схеме (18). Если в соседней ячейке Cout(cm) уже получено новое итерационное значение коэффициентов, то в этой сумме используется значение
(( out
(ст )))к+1. В пРотивном слУчае, используется значение (Aq(C out (ст )))к .
Для уменьшения зависимости решения от выбранной нумерации ячеек и для более равномерного распространения информации по всем пространственным направлениям, после каждой итерации производится перенумерация ячеек.
3. ТЕСТОВЫЕ РАСЧЕТЫ
Была написана программа, реализующая представленный численный метод для задач, описываемых уравнением вида (1), на одноблочных прямоугольных неравномерных сетках с возможностью задания неоднородных граничных условий.
Следует отметить, что в тестовых расчетах вместо формулы Fdiff (G) = -vG использовалось следующее выражение для диффузионных потоков:
F diff (G) = -v
0
Gy
Благодаря пренебрежению диффузионными потоками вдоль оси х в каждой из рассмотренных задач имелось аналитическое решение.
Первая рассмотренная задача — одномерная нестационарная непрерывная волна сжатия, переходящая в стационарный скачок уплотнения, размытый вязкостью. Решалось следующее уравнение:
= 0, v = 0.1,
ди + д сЛ ду
Требовалось получить стационарное решение, зависящее только от у, в области
и г
--vGy
2 y
-L < y < L, L = 4,
с граничными условиями
и (±L, t) = u0tanh |+ |, и0 = 1,
Gy (±L, t) =
—и,
0
2v cosh2I U°L
2v
и начальными условиями
u(y, 0) =
Mo, -L < y <-L2,
- 2uo y/L, - L/ 2 < y < L/ 2, L/2 < y < L.
-Un
Точное стационарное решение задачи:
и(у, t ^да) = и01апЬ
В расчетах этой задачи использовался набор равномерных сеток с одной ячейкой в х-направ-лении и 16, 32, 64, 128, 256 и 512 ячейками в ^-направлении. Расчеты доводились до состояния, максимально близкого (в рамках арифметики с двойной точностью) к стационарному.
Для нахождения порядков точности вычислялись С- и ¿2 -нормы разности между численным (ипит ) и точным (иех ) решениями:
С = таХ Ипит - иех |,
L2 =
уЦ(unum -Uex)2d^.
Пусть С(1) и ^ — значения С- и ¿2 -норм на наиболее подробной сетке, а С(2) и ¿(22) — значения С- и ¿2 -норм на сетке с вдвое большим шагом. Порядки точности Ыс и N1 вычислялись по этим значениям следующим образом:
с(2) l22)
NC = log2 ^ NL = log2 l2l)-
2
2
Таблица 1
Порядки точности по С- и Ь2 -нормам, полученные в задаче о скачке уплотнения, размытом вязкостью
Метод K = 0 K = 1 K = 2 K = 3
NC NL NC NL NC NL NC Nl
Bassi & Rebay 1 Cockburn & Shu Bassi & Rebay 2 1.07 0.99 0.03 0.98 0.99 -0.11 1.35 1.97 1.96 1.52 1.99 1.98 3.06 2.98 3.02 3.02 2.99 3.01 3.16 4 3.88 3.2 3.99 4
В табл. l приведены порядки точности, вычисленные описанным способом. На рис. l изображены графики зависимости С- и L2 -норм от числа ячеек сетки.
Из приведенных данных можно сделать следующие выводы:
1) в рассмотренной задаче теоретические порядки точности достигаются лишь на весьма подробных сетках, которые в практических задачах обычно не применяются;
2) порядки точности, вычисленные по норме С, при измельчении сетки приближаются к теоретическим медленнее, чем порядки, вычисленные по норме L2 ;
3) при K = 0 схемы Bassi & Rebay 1 и Cockburn & Shu имеют первый порядок точности, причем на одинаковых сетках Cockburn & Shu дает меньшие ошибки; решение по схеме Bassi & Rebay 2 не сходится к точному (порядок точности равен нулю);
4) при K = 1 схемы Cockburn & Shu и Bassi & Rebay 2 имеют второй порядок точности, причем на одинаковых сетках Bassi & Rebay 2 имеет меньшие ошибки по норме C; схема Bassi & Rebay 1 дает пониженный порядок точности — около 1.5;
5) при K = 2 все рассмотренные схемы имеют третий порядок точности, на подробных сетках наименьшие ошибки дает Bassi & Rebay 1, наибольшие — Cockburn & Shu, на грубых сетках — наоборот;
6) при K = 3 схемы Cockburn & Shu и Bassi & Rebay 2 имеют четвертый порядок точности, причем на одинаковых сетках Bassi & Rebay 2 имеет меньшие ошибки по норме C; схема Bassi & Rebay 1 дает пониженный порядок точности — около 3.2;
7) несмотря на пониженный порядок точности при нечетных K, в данной задаче на сетках средней густоты ошибки по схеме Bassi & Rebay 1 не превышают ошибок по другим рассмотренным схемам (имеющим порядок точности K + 1).
С практической точки зрения, интерес может представлять качество воспроизведения деталей решения на очень грубых сетках. На рис. 2 представлены профили скачка уплотнения, содержащего лишь около 2 ячеек в области размытия вязкостью. При K = 1 хорошо видна несимметричность решения по схеме Cockburn & Shu; более того, расчеты по этой схеме показывают, что при сближении границ расчетной области скачок уплотнения перестает быть стационарным и приобретает малую, но ненулевую скорость, перемещаясь в итоге на границу расчетной области. Однако этот недостаток проявляется лишь в тех задачах, в которых стационарное решение зависит от начальных условий. В типичных стационарных задачах аэродинамики после процесса установления решение зависит только от параметров, заданных на границе, а начальные условия «забываются». В этих случаях схема Cockburn & Shu не будет порождать нефизичных нестационарных процессов.
200 300 400600 N,
200 300 400500
Рис. 1. Зависимость С-нормы (а) и ¿2-нормы (б) от числа ячеек сетки в расчетах скачка уплотнения,
размытого вязкостью: ---метод Bassi & Rebay 1;----— Cockburn & Shu;------Bassi & Rebay 2
Рис. 2. Профили скачка уплотнения на грубой сетке:
столбцы: а — метод Bassi & Rebay 1; б — Cockburn & Shu; в — Bassi & Rebay 2; строки: I — K = 1, II — K = 2, III — K = 3
Сравнение решений при K = 2 и 3 показывает худшую работу схемы Bassi & Rebay 1 по сравнению с другими схемами — разрывы решения на гранях у нее более сильные. Интересно отметить, что при K = 3 (особенно по схеме Cockburn & Shu) полученные профили вполне удовлетворительно воспроизводят форму скачка, несмотря на то, что на него приходится всего 2 ячейки.
Обнаруженное в тестах отсутствие аппроксимации схемы Bassi & Rebay 2 при K = 0 легко можно вывести аналитически. Рассмотрим схему для одномерного уравнения теплопроводности на равномерной сетке с шагом h:
ди дО „ ди
--V-= 0, О =—.
дt дх дх
Проинтегрируем это уравнение по 7-й ячейке (считаем, что у нее есть соседи, а единственная базисная функция ф( х) тождественно равна 1):
д" Тudx-vG\x'^ = 0 ^ U-vGi+12 -Gi-12 = 0, (19)
dt J lxi-V2 dt h ' v '
xi-hl2
Xi +h/2
где U1 = J udx . Непосредственное дифференцирование поля u в ячейках дает тождественный
X- -Щ 2
нуль (в каждой ячейке u = U1 = const), так что при аппроксимации Bassi & Rebay 2 величины Gi±Xj2 состоят лишь из вкладов поправок R (см. (8)):
_ Rj—ght + RT _ Rfght + R-+t (20) Gi—1/2 =-2-' G'+V2 =-2-.
Величины R определяются по формуле (7), из которой
„right = „left = Ui - U1—1 „right = „left = Ui+1 - Ui 1-1 ^ 2h ' ^ ^+1 2h '
Подставляя эти выражения в (20), получаем
, = Ui — Ui—1 G = Ui+1 — U 1—12 ' Gi+12
Окончательно схема (19) запишется в виде:
уЦ+1 - 2Ц + и,-1 =0 дt 2й2 '
ди д^и ди 1 д2и т. е. вместо уравнения — = у- решается уравнение —= —V-. Это и означает потерю
дt дх2 дt 2 дх2
аппроксимации.
К сожалению, эта особенность схемы Bassi & Rebay 2 не позволяет использовать ее напрямую совместно с методом p-multigrid [12], где требуется получение решений при последовательно уменьшающихся К вплоть до нуля. Однако, решение при К = 0 можно найти по какой-либо другой схеме.
Второй тест — задача Вигтона [9] об искривленном слое смешения двух сред с параметрами и = —1 и и = 1. Среды движутся с одинаковой скоростью. х -компонента скорости постоянна и равна единице, а ^-компонента меняется по гармоническому закону. Смешение сред обусловлено диффузией в направлении оси у. Взаимодействие потоков начинается в точке (х,,; 0). Эта задача описывается следующим уравнением:
— + — +—|~1.25^ (5( х — х,) )и — V Оу ] = 0, у = 0.05.
дt дх дУ
Требовалось получить стационарное решение в области
0 < X < 1, —2 < y < 2,
Рис. 3. Общий вид решения задачи Вигтона: а — гладкие граничные условия; б — разрывные граничные условия
с граничными условиями
u(0, y, t) = uex (0, y), u(x, - 2, t) = Uex(x, - 2), u(x, 2, t) = Uex(x, 2),
Gy (x, - 2, t) =
duex(x, y)
dy
Gv (x, 2, t) =
duex(x, y)
y=-2
dy
y=2
Здесь uex(x, y) — точное стационарное решение:
uex(x, y) = u(x, y, t ^ да) = erf
f y - 0.25 sin(5(x - xs)) ^
С математической точки зрения, данная задача не требует постановки граничного условия на границе х = 1. В расчете параметры на этой границе брались изнутри расчетной области. Задавались следующие начальные условия:
u(x, y, 0) =
1, y > 0, -1, y < 0.
Рассмотрены случаи =-0.1 (гладкие граничные условия) и 0 (разрывные граничные условия). Вид точного решения в области слоя смешения показан на рис. 3.
Для расчетов использовался набор равномерных прямоугольных сеток размерности 8 х 32, 16 х 64, 32 х128, 64 х 256 и 128 х 512 ячеек. В табл. 2, 3 приведены порядки точности, полученные при сравнении решений на двух наиболее подробных сетках. Графики зависимости ¿2 - и С-норм от числа ячеек сетки в х -направлении представлены на рис. 4, 5.
Таблица 2
Порядки точности по С- и Ь2 -нормам, полученные в задаче Вигтона с гладкими граничными условиями
Метод K = 0 K = 1 K = 2 K = 3
Nc NL Nc NL Nc NL Nc Nl
Bassi & Rebay 1 0.93 0.96 1.9 1.92 2.82 2.95 3.74 3.8
Cockburn & Shu 0.93 0.95 1.91 2 2.86 2.97 3.79 3.98
Bassi & Rebay 2 | 0.48 | -0.41 | 1.92 | 2.01 | 2.88 | 2.99 | 3.8 | 4
Таблица 3
Порядки точности по С- и Ь2 -нормам, полученные в задаче Вигтона с разрывными граничными условиями
Метод К = 0 К = 1 К = 2 К = 3
Ыс Ыь Ыс Ыь Ыс Ыь Ыс Ыь
Bassi & Rebay 1 0.19 0.92 0.2 0.98 0.11 0.88 0.08 0.74
Cockburn & БИи 0.28 0.92 0.1 0.96 0.4 1.01 0.4 1.04
Bassi & Rebay 2 0.27 -0.5 0.27 0.94 0.25 0.95 0.26 0.91
Рис. 4. Зависимость С-нормы (а) и Ь2-нормы (б) от числа ячеек сетки в х-направлении в расчетах задачи Вигтона с гладкими граничными условиями (обозначения, как на рис. 1)
Рис. 5. Зависимость С-нормы (а) и ¿2-нормы (б) от числа ячеек сетки в х-направлении в расчетах задачи Вигтона с разрывными граничными условиями (обозначения, как на рис. 1)
Из расчетов задачи Вигтона с гладкими граничными условиями можно сделать следующие выводы:
1) все рассмотренные схемы дают близкие решения и уровни ошибок, за исключением Bassi & Rebay 2 при К = 0 (в этом случае аппроксимация отсутствует); при К > 0 эта схема, наоборот, дает чуть меньшие ошибки по сравнению с другими;
2) порядок точности по схеме Bassi & Rebay 1 даже при нечетных К при измельчении сетки приближается к оптимальному значению К +1, несмотря на то, что в общем случае эта схема
|
-4 -5 -6
-7
-8
a) 6) в)
■О"
Рис. 6. Распределение ошибки по ячейкам расчетной области в задаче Вигтона с разрывными граничными условиями. Метод Bassi & Rebay 2:
а — K = 1; б — K = 2; в — K = 3
имеет при нечетных K лишь K-й порядок точности [7]; вероятно, это связано с видом конвективных потоков в решаемом уравнении.
В расчетах с разрывными граничными условиями поведение ошибок другое:
1) в С-норме порядок точности всех схем близок к нулю. Ошибка уменьшается с увеличением K, но практически не зависит от шага сетки; уровень ошибки C ~ 0.2...1 — того же порядка, что и характерный масштаб задачи Uo = 1;
2) в ¿2 -норме порядок точности всех схем не превышает 1. Наиболее стабильно первый порядок выдерживается в схеме Cockburn & Shu, причем эта схема обеспечивает наименьший уровень ошибки из трех рассмотренных.
Особый интерес представляет характер распространения ошибки |unum - uex | из области
разрыва в граничных условиях. На рис. 6 изображены поля |unum - uex | (максимумы по ячейкам)
в логарифмическом масштабе, полученные в расчетах по схеме Bassi & Rebay 2 на сетке 32 х 128 ячеек с K = 1, 2, 3. Одинаково высокая ошибка в ячейках, примыкающих к точке начала смешения, затухает вниз по потоку до некоторого конечного уровня, причем скорость затухания и конечный уровень ошибки зависят от K. Затухание ошибки в основном связано с действием вязкости v (если положить v = 0, то ошибка лишь сносится вдоль «линий тока»). Увеличение скорости затухания ошибки с ростом K, по-видимому, связано с тем, что из-за увеличения K возрастают градиенты ошибки в ячейках, примыкающих к разрыву граничных условий, что приводит к увеличению диффузионных потоков и более быстрому размытию ошибки по пространству.
Чтобы выяснить, сохраняется ли высокий порядок точности в области гладкого решения, была проведена дополнительная серия расчетов, в которой нормы разности между численным и точным решениями определялись в области 0.875 < x < 1 — вдали от разрыва. В табл. 4 представлены порядки точности, вычисленные по полученным данным.
Таблица 4
Порядки точности по С- и Ь2 -нормам, полученные в задаче Вигтона с разрывными граничными условиями
в области 0.875 < х < 1
Метод K = 0 K = 1 K = 2 K = 3
NC NL NC NL NC NL NC Nl
Bassi & Rebay 1 0.9 0.93 1.73 1.93 2.95 2.99 3.61 3.77
СосИшт & БИи 0.9 0.93 1.94 2.01 2.95 3 3.97 4
Bassi & Rebay 2 -0.41 -0.63 1.98 2.04 2.99 3.01 3.9 4.01
Видно, что в области гладкого решения реализуется высокий порядок точности. Ошибка порядка 0(1) в точке разрыва оказывает влияние на решение лишь в ее окрестности.
Третья тестовая задача — модельный пограничный слой, описываемый уравнением
--Оу ] = 0, v = 4 -10-4,
ду1 у J
в расчетной области
0 < х < 1, 0 < у < 0.5
с граничными условиями
и(0, у, *) = ^(0, у), и(х, 0,*) = UeX (х, 0) = 0, и(х, 0.5,*) = и^х (х, 0.5),
О,(х, -0, 0 = О,(х, + 0, 0, О,(х, 0.5, 0 =диехд^у1 .
У У У ду у=0.5
Здесь и^(х, у) — точное стационарное решение
^х(х, у) = и (х, у, * ^да) = /(п), п = у/^ (х - х„), / (п) — решение краевой задачи
/"+1 п//' = 0, /(0) = 0, /(п) ^ 1.
2 п^да
Как и в предыдущей тестовой задаче, параметры на границе х = 1 брались изнутри расчетной области.
Задавались следующие начальные условия:
и (х, у, 0) = 1.
Рассмотрены случаи =-0.5 (гладкое решение без особенностей) и = 0 (особенность в точке х = 0, у = 0).
Вид точного решения в области пограничного слоя представлен на рис. 7 (изображение сжато по оси х в 10 раз).
Была построена серия вложенных сеток, самая мелкая из которых имеет размерность 64 х128 ячеек, минимальный шаг по продольной оси ~ 1.46 х 10-4, по поперечной оси кут1а ~ 1.31 х 10-4 . Шаг сетки вдоль каждой оси изменяется по закону геометрической прогрессии. Более грубые сетки были получены удалением каждого второго узла более мелких сеток.
ди + д д* дх
Рис. 7. Общий вид решения задачи о модельном пограничном слое:
а — гладкое решение; б — решение с особенностью
ц
20
40 60
i-tf
t^N,
Рис. 8. Зависимость С-нормы (а) и ¿2-нормы (б) от числа ячеек сетки в х-направлении в расчетах задачи о модельном пограничном слое с гладким решением (обозначения, как на рис. 1)
Рис. 9. Зависимость С-нормы (а) и ¿2-нормы (б) от числа ячеек сетки в х-направлении в расчетах задачи о модельном пограничном слое, решение с особенностью (обозначения, как на рис. 1)
Порядки точности, вычисленные по решениям на двух наиболее подробных сетках, приведены в табл. 5, 6. Графики зависимости L2 - и С-норм от числа ячеек сетки в х-направлении представлены на рис. 8, 9.
Расчеты задачи с гладким решением показывают, что
1) схема Bassi & Rebay 1 при четных K имеет порядок точности K +1, а при нечетных — порядок точности K;
2) схема Cockburn & Shu для всех рассмотренных K имеет порядок точности K +1;
3) схема Bassi & Rebay 2 при K > 0 имеет порядок точности K +1 и чуть меньший уровень ошибки, чем Cockburn & Shu.
Таблица 5
Порядки точности по С- и Ь2 -нормам, полученные в задаче о модельном пограничном слое без особенностей
Метод K = 0 K = 1 K = 2 K = 3
Nc NL Nc NL Nc NL Nc Nl
Bassi & Rebay 1 0.97 0.98 0.99 1 2.74 2.83 3.06 3.27
Cockburn & Shu 0.98 0.99 1.98 1.99 2.92 2.96 3.88 3.96
Bassi & Rebay 2 0.15 -0.03 2 2 2.97 3.02 3.92 4.02
Таблица 6
Порядки точности по С- и Ь2 -нормам, полученные в задаче о модельном пограничном слое с особенностью
в точке х = 0, у = 0
Метод K = 0 K = 1 K = 2 K = 3
Nc NL NL Nc NL Nc Nl
Bassi & Rebay 1 0.25 0.98 0.15 0.99 0.02 1.05 -0.23 0.69
Cockburn & Shu -0.28 0.98 0.35 1.83 0.31 1.15 0.29 1.13
Bassi & Rebay 2 -0.25 -0.06 0.48 1.8 0.25 1.17 0.17 1.09
Из расчетов с особенностью в начальной точке пограничного слоя следует, что
1) поведение С-нормы еще менее регулярное, чем в задаче Вигтона с разрывными граничными условиями. Ошибка меняется немонотонно не только при измельчении сетки, но и при увеличении K; уровень ошибки — C ~ 0.4.. 1 — снова того же порядка, что и характерный масштаб задачи u0 =1;
2) в графиках ¿2 -нормы (для схем с K = 2 и 3) видны два участка: на относительно грубых сетках порядок точности значительно превышает 1, на мелких сетках (при уровне ошибки ¿2 < 10-4) порядок уменьшается и приближается к 1; по-видимому, это связано с тем, что ячейка, в которой находится особенность, является мельчайшей, и на грубых сетках ошибка первого порядка O(h) в ней оказывается значительно меньшей, чем ошибки высокого порядка O(hK+1) в крупных ячейках. При измельчении сетки ошибки высокого порядка быстро уменьшаются, и на первый план выходит ошибка O(h). Можно сделать вывод, что на практике в областях с резким изменением решения имеет смысл делать сильные сгущения сетки и, возможно, локально снижать степень разложения решения K.
4. ЗАКЛЮЧЕНИЕ
Результаты расчетов по предложенному в настоящей статье численному методу позволяют сделать выводы о преимуществах и недостатках рассмотренных методов аппроксимации вязких членов в методе РМГ.
Метод Bassi & Rebay 1, несмотря на кажущуюся простоту, имеет существенные недостатки — пониженный порядок точности при нечетных K и некомпактный шаблон. Применять его не рекомендуется.
Метод Cockburn & Shu может быть лучшим выбором для расчетов задач, в которых стационарное решение не зависит от начальных условий, на структурированных сетках. Также он чуть лучше других работает при наличии в решении особенностей. Формулировка этого метода проще, чем методов Bassi & Rebay 1 и 2. Оптимальный порядок точности K +1 достигается при любых K, что позволяет без модификаций использовать его совместно с методом p-multigrid. Из не-
достатков необходимо отметить неочевидность обобщения метода на неструктурированные сетки и несимметричность формул, противоречащую представлениям о характере процесса диффузии.
Метод Bassi & Rebay 2 — оптимальный вариант на неструктурированных сетках, если не предполагается использовать численный метод с K = 0 (p-multigrid, а также локальное понижение K в окрестности разрывов решения). Схема имеет компактный и симметричный шаблон.
Установлено, что на точечных разрывах граничных условий происходит потеря аппроксимации, однако область высоких ошибок локализована в пространстве, и вдали от нее реализуется высокий порядок точности, соответствующий гладкому решению.
Работа выполнена при финансовой поддержке гранта РФФИ № 09-01-00243-а.
ЛИТЕРАТУРА
1. Wang Z. J. High-order methods for the Euler and Navier — Stokes equations on unstructured grids // Prog. Aerosp. Sci. 2007. V. 43, p. 1 — 41.
2. Волков А. В., Ляпунов С. В. Применение конечно-элементного метода Галер-кина с разрывными базисными функциями к решению уравнений Рейнольдса на неструктурированных адаптивных сетках // Ученые записки ЦАГИ. 2007. Т. XXXVIII, № 3 — 4, с. 22 — 30.
3. Волков А. В. Особенности применения метода Галеркина к решению пространственных уравнений Навье — Стокса на неструктурированных гексаэдральных сетках // Ученые записки ЦАГИ. 2009. Т. XL, № 6, с. 41 — 59.
4. A rno l d D. N., B rez z i F., Cockburn B., Mar ini L. D. Unified analysis of discontinuous Galerkin methods for elliptic problems // SIAM J. Numer. Anal. 2002. V. 39, № 5, p. 1749 — 1779.
5. Bassi F., Rebay S. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier — Stokes equations // J. Comput. Phys. 1997. V. 131, № 1, p. 267 — 279.
6. B a s s i F., C r i v e 11 i n i A., R e b ay S., S a v in i M. Discontinuous Galerkin solution of the Reynolds-averaged Navier — Stokes and k — ю turbulence model equations // Computers & Fluids. 2005. V. 34, p. 507 — 540.
7. Bassi F., Rebay S., Mariotti G., Pedinotti S., Savini M. A high-order accurate discontinuous finite element method for inviscid and viscous turbomachinery flows / Proceedings of 2nd European Conference on Turbomachinery — Fluid Dynamics and Thermodynamics. — Technologisch Instituut, Antwerpen, Belgium. 1997, p. 99 — 108.
8. Cockburn B., Shu C.-W. The local discontinuous Galerkin method for time-dependent convection diffusion system // SIAM J. Numer. Anal. 1998. V. 35, № 6, p. 2440 — 2463.
9. Wolkov A. V., Lyapunov S. V. Application of discontinuous Galerkin finite element method to the solution of partial differential equations. Part I. 2D scalar conservation laws // 16th IMACS World Congress. — Lausanne, Switzerland, August 21 — 25, 2000.
10. Босняков С. М., АкинфиевВ. О., ВласенкоВ. В., Глазков С. А., Горбушин А. Р., Кажан Е. В., Курсаков И. А., Лысенков А. В., Матяш С. В., Михайлов С. В. Использование методов вычислительной аэродинамики в экспериментальных работах ЦАГИ // Математическое моделирование. 2011. Т. 23, № 11, с. 65 — 98.
11. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. — М.: Мир, 1991.
12. F i dk o w s k i K. J., Oliver T. A., L u J., D ar m o f a l D. L. p-Multigrid solution of high-order discontinuous Galerkin discretizations of the compressible Navier — Stokes equations // J. Comput. Phys. 2005. V. 207, № 1, p. 92 — 113.
Рукопись поступила 27/III2012 г.
