CC BY
63
ПРИКЛАДНАЯ МАТЕМАТИКА И МЕТОДЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ
УДК 519.6
Ю. И. Димитриенко, М. Н. Коряков, А. А. Захаров, Е. К. Сыздыков
РАЗВИТИЕ МЕТОДА ЛЕНТОЧНО-АДАПТИВНЫХ
СЕТОК НА ОСНОВЕ СХЕМ TVD
ДЛЯ РЕШЕНИЯ ЗАДАЧ ГАЗОВОЙ ДИНАМИКИ
Предложен вариант метода ленточно-адаптивных сеток, основанный на схеме TVD 2-го порядка точности. Представлены результаты тестирования численного метода на основе решения задач распада разрыва, обтекания клина и ступеньки. Показано, что разработанный метод на основе схемы TVD обеспечивает более высокую точность численного моделирования, чем метод ленточно-адаптивных сеток на основе схемы Мак-Кормака, однако время расчета этими методами примерно одинаково.
E-mail: dimit.bmstu@gmail.com
Ключевые слова: газовая динамика, численные методы, схемы TVD,
ленточно-адаптивные сетки.
Схемы TVD в последнее время стали достаточно популярны при решении задач газовой динамики благодаря повышенной точности численного решения, получаемой с их помощью, особенно для задач с наличием сильных разрывов [1, 2]. Однако их, как правило, используют в ортогональных системах координат, поэтому опыт применения этих схем для задач с областями сложной формы пока не велик. Цель настоящей работы — развитие алгоритма построения схем TVD по методу Хартена для численного решения системы уравнений газовой динамики в произвольных неортогональных системах координат, которые применяются при использовании метода ленточно-адаптивных сеток [1].
Система уравнений газовой динамики. Рассмотрим систему уравнений динамики идеального газа, которая в адаптивных криволинейных координатах Xг может быть записана в виде [1]
OU дVl _ ~dt + OX- = 0'
(1)
где
и = Vg
p
pQ)v
pQ2vj pQfv3 pE
V = Vg
( pvl
Q)(fjvlv3 + pglj) Q2{pvlvj + pglj) Q3j(pvtvj + pglj) vl(pE + p)
(2)
I |2
— координатные столбцы; р — плотность газа; E = cv9 + —--полная
2
энергия газа; cv — теплоемкость при постоянном объеме; 9 — тем-
I |2 г
пература газа; |v| = угуг — квадрат модуля скорости; p = р—9 —
р
давление, R — универсальная газовая постоянная; р — молекулярная масса газа; д = det gij — детерминант метрической матрицы gij = = гг • rj = QkQjSki; V — компоненты вектора скорости в базисе гг ада-
dx3
птивной системы координат Xг; ^ = . — компоненты якобиевой
дХг
матрицы преобразования декартовых координат в адаптивные; E — метрический тензор; V = ггтр|1 — набла-оператор в криволинейных координатах Xг; 7 = о.р/си — коэффициент Пуассона.
Кроме системы (1), имеющей дивергентный вид, рассмотрим также систему уравнений газовой динамики недивергентного вида [1]
д U д V3
--b Pl-= 0,
dt + 3 dXi 0,
(3)
где
U =
( Р \
pv1 pv2
— 3
pv \pE )
V1 =
í pvг \ pv'V1 + pS%l
pv%V2 + pSi2
pvf'v3 + pSl3
V vi (pE + p) )
(4)
а V1 — компоненты в декартовой системе координат.
Схема ТVD для дивергентной системы. Рассмотрим способ получения ТУБ-схемы 2-го порядка аппроксимации, основанный на методе Хартена [3] (методе модифицированного потока), который заключается в применении ТУБ-схемы 1-го порядка аппроксимации к соответствующим образом модифицированным потокам функций:
/ = (/ + — д), где А = ——, а д — модифицирующая добавка (будет V А '
введена ниже). Тогда конечно-разностные соотношения, аппроксимирующие уравнения Эйлера (1), можно представить в виде явной схемы
un+1 = и - a^V" - vk-ni/2)-
- a2(vkvni/2 - vi") - a3(v;
fc+1/2
- V3-1/2) =
3
= U - E[Aí(Vk+i/2 - Vk-1/2)], (5)
i= 1
где иП = И(Хк, ¿п) — значения функций в к-м узле ленточной адаптивной сетки (нумеруются одним индексом) в момент времени ¿п;
А* = 2А1/(Хк+ - Хк_),
( Е, г = 1; ( Вк, г = 1;
к+ = } Як, г = 2; к- = < Ьк, г = 2;
{ ик, г = 3; [ Бк, г = з.
Здесь Ек,Як,ик,Ьк,Бк,Вк — индексы в адаптивной системе координат для узлов, находящихся "спереди", "справа", "сверху", "слева", "снизу", "сзади" по отношению к узлу к, Vк+1/2 — столбцы численного потока в г-м направлении, которые определяются как
V
k+1/2
= 1 (vr + Vk+ + 4 Rl
2 у 'к 1 - к+ 1 А* "к+1/2 • Фк+1/2^ . (6)
В формуле (6) УГ = У*(Хк,Г); 1/2 = (К£ + ИГ)/2; ИГ = = Иг(Хк, ¿п); Иг — матрица правых собственных векторов матрицы Сг = дУ*/дИ, имеющей следующий вид:
Gl =
0
+ (7 - 1)P1|v|2
2
(Y - 1)P2 |v|2
2
(Y - 1)P3 |v|2
—vlv2 +
—vlv3 +
€М '<Y-M_ H
Pi
Vl - (y - 2)€1P1l
pilv2 - (y - i)v1p2
P-lv3 - (y - i)v1P3 P[H - (y - i)vlv1
P P2
P
P3
P2v1 - (y - 1)v2P1 P3v1 - (y - 1)v3P1 (y - 1)P1 vl - (y - 2)v2P2 p3v2 - (y - i)v3P2 (y - 1)P2
p2v3 - (y - 1)v2P3 vl - (y - 2)v3P3 (y - 1)P3
P2H - (y - 1)vlv2 P3H - (y - 1)vlv3
|v|
где H =
Yv
+ —--интеграл Бернулли; а2 — квадрат скорости
Y - 1 2
звука. Матрица Rk правых собственных векторов имеет следующий
2
0
вид:
Ri
v1 — asi
— 2
v2 — as2
v3 — as3
v1
v2
v3
0
s2 s1 0
0
—s3 0 s1
v1 + as1 v2 + as2 v3 + as3
2
vi | v |2 _1 _2 _3 _1 v
H — a— —— v1s2 — v2s1 v3s1 — v1s3 H + a
Г
V/gii.
(7)
если s1 = 0;
Ri
1 1
v1 — as1 v1
v2 — as2 v2
v3 — as3 v3
H— vi a—= v/gii |v|2 2
0
s2 —s1 0
0 0
s3 s2
v1 + as1 v2 + as2 v3 + as3
- v1s2 — v2s1 v2s3 — v3s2 H + a
v/gii.
если s2 = 0;
(8)
Ri
1 1
v1 — as1 v1
v2 — as2 v2
v3 — as3 v3
H— vi a .— v/gii |v|2 2
0
—s3 0 s1
0 0
s3 —s2
v1 + as1 v2 + as2 v3 + as3
—— v3s1 — v1s3 v2s3 — v3s2 H + a
V/gii.
если s3 = 0. Здесь Sj = Pj
(9)
Столбец Ф ненты:
k+1/2 численной диссипации имеет следующие компо-
Ф
k+1/2
n
=
fc+1/2
i %,n
5,fc+1/2
)T;
(10)
li,n i,n . i,n / /\i i,n . i,n \ i,n
Ф>+1/2 = j + gj,fe+ — aj,fc+1/2 + 7j,n+1/2)aj,fe+1/2;
i,n i,n iy> • n~i,n I ~i,n i,n м
j = sj',fe+1/2 maxt0, mm(|^.i';n+1/2|,gj;fe+1/2 sj',fe+1/2)];
sj;fe+1/2 = sign(
1
1
1
1
1
~i,n 1 г / / \i i,n \ /\г i,n \2i i,n
g,j',fe+1/2 = aj,fe+1/2) - (A ai,fc+1/2) Jajf,fc+1/2;
Y
(gj,n+ - gj,n)/а
j,k+1/2
jf,k+1/2'
0,
а а
= 0;
г,п f,k+1/2 г,п = о.
j,k+1/2
Величины а*'" = а* , ¿п) — это значения собственных чисел а* матриц Сг в точке к в момент времени ¿п, которые определяются по формулам
аг2 = уг; а3 = а4 = V; аг5 = V + (11)
представляющие собой перепады вектора газодинамических переменных И" в характеристической форме, находятся из соотношений
а\ = v1 - a\/~äf%; Величины а*'П+1/2
а
г,п "fc+1/2
»i,n
z' i, = (аl,
г,п fe+1/2,
,аг5,Пк+1/2)Т = (R+1/2)-1(Un+ - un). (12)
Здесь (К£ 1/г)-1 = ((К'+ )-1 + (КГ)-1)/2, а (КГ)-1 = (Кг)-1(Х,Г) — матрицы левых собственных векторов матрицы Сг, которые имеют следующий вид [4]:
(4 )-1 =
(y - 1)|v|2 + 2a-V (1 - y)^1 - asi (1 - y)^2 - as2 (1 - y)v3 - as3 Y - 1
4a2 2a2 - (y - 1)|v |2 2a2 (y - 1)г1 2a2 (y - 1)v2 2a2 (y - 1)г3 2a2 Y - 1
4a2, v2--S2 v9" a2 s2 -S3 a2 s2 -1 a2 s2 s3 a2 0 0
. si —S3 - г3 vs" s1 s2s3 s1 1 - s3
si - s1 s1
(Y - 1)|v|2 - 2a-v 7ЦЯ (1 - y)^1 + asi (1 - y)^2 + as2 (1 - y)v3 + as3 Y - 1
4a2 2a2 2a2 2a2 2a2 J
если s1 = 0;
К )-1 =
(Y - 1) | v |2 +2a -V (1 - y)^1 - as1 (1 - y)V2 - as2 (1 - y)^3 - as3 Y - 1
4a2 2a2 2a2 2a2 2a2
2a2 - (y - 1)|v |2 (y - 1)v1 (y - 1)v2 (y - 1)г3 Y - 1
-i 4a2 —'т= s1 - v1 V s" a2 1 - s1 a2 -s1 a2 s1s3 a2 0
s2 s2 - s2
-3 7Z г3--^ s3 V s" s1s3 s3 s3 - 1 0
s2 s2 s1
(Y - 1)|v|2 - 2a-V (1 - y)^1 + as1 (1 - y)v2 + as2 (1 - y)^3 + as3 Y - 1
(13)
(14)
4a2 2a2 2a2 2a2 2a2 J
если s2 = 0;
(R )-1 =
(y - 1)|v|2 + 2a-= V g (1 - y)w1 - as1 (1 - y)V2 - as2 (1 - y)v3 - as3 y - 1
4a2 2a2 2a2 2a2 2a2
2a2 - (y - 1)|v|2 (y -1)«1 (y - 1)v2 (Y - 1)v3 y - 1
4a2 a2 a2 a2 a2
w1--7= si V g" s1 - 1 S1S2 -S1 0
S3 S3 S3
—= S2 - v3 yg" S1S2 1 - s2 -S2 0
S3 S3 S1
(Y - 1)|v|2 - 2a-= у g" (1 - y)«1 + as1 (1 - y)v2 + as2 (1 - y)v3 + as3 Y - 1
4a2
2a2
2a2
2a2
2a2
(15)
если s3 = 0.
Функция численной вязкости ф определяется по формуле
, |z|, z > е;
фЫ = ^
V ) ^ (z2 + е2)/2е, z < е,
(16)
где е — параметр диссипации. При е = 0 диссипация будет наименьшей; чем больше е, тем более диссипативной является разностная схема.
Для недивергентной системы (3) конечно-разностные соотношения схемы ТУБ 2-го порядка, примененной для ленточно-адаптивной сетки, имеют вид
ип+1 = ип
V
V
V
1 ' n k+1/2 -V1 'n , v k —1/2 V1 '« , Vk+1/2 _ V1 'n , v k—1/2 V1 '« , v k+1/2 _ "V1 'n v k—1
xf' n - Xk1 'n xr n - Xk2 'n -1^3 ' n -1^3 ' n - Xk
2 n k+1/2 - V k-1/2 V2 'n, Vk+1/2 - "k—1/2 V2 'n, Vk+1/2 _ "V2 'n v k—1
X1'n - Xk1 'n xr' n - Xk2 'n v3 ' n v3 ' n - Xk
■3' n k+1/2 - "k-1/2 V3 '«, Vk+1/2 - "k—1/2 "V3 'n, Vk+1/2 _ "V3 'n v k—1
v1 ' n - Xk1 '« xR'n - Xk2 'n -1^3 ' n -1^3 ' n Xk
- At
- At
- At
где численные потоки представляются в виде
1 /Т-Т1„П т—,-¿,«,4 JIM 1
2 \ k ■ k+/_ ¿,k 1
T/i>n _ -1" /Л7"V" , у¿'n^ Р-Jn | ' Р1 . « Лг'п
+ Vfc+)Pi'fc + Kfe+i/2 • фк+1/2;
k+1/2
(17)
(18)
Р/кГ = Р* (Хк, £г) — компоненты обратной матрицы Якоби перехода из декартовых координат к адаптивным. Матрица К'+1/2 правых собственных векторов получается из (7)-(9) заменой Р/ на , ьг на аг, дгг на 1 для всех г.
Столбец Фк'+1/2 определяется по формулам (10), (11), в кото-
гг _-г г _-г г
рых следует провести замену а/к+1/2 ^ а/к+1/2, где а/к+1/2 = = а/' Г+1/2Р/к+1/2, а а/" = а/(Хк, £г) — собственные значения матрицы Сг, в декартовой системы координат имеющие вид
а1 = аг - а; а2 = -аг; аг3 = -аг; а4 = аг; аг5 = аг + а. (19)
Матрицы левых собственных векторов (Кк+1/г) также получаются из соотношений (13)-(15) заменой Р/ на , ьг на -уг, дгг = 1 для всех г.
Была проведена серия численных тестов для апробации разработанного метода ТУБ. В качестве сравнения использовались результаты задач, полученные по методу ленточно-адаптивных сеток на основе схем типа Мак-Кормака [1], а также численные расчеты на основе схемы ККОО, взятые из [5].
Задача о распаде произвольного разрыва. В этой классической задаче [1] в качестве начальных данных были заданы скачки безразмерной плотности и давления: при х ^ 0,5 р =1, р =1; при х > 0,5 р = 0,1, р = 0,125 (все параметры здесь и далее, в том числе геометрические размеры и время, — безразмерные). На рис. 1 показаны распределения плотности, давления и скорости газа в момент времени £ = 0,25. Решение по методу ТУО имеет значительно более высокую точность решения, в том числе для скачка плотности и давления на ударных волнах и на контактном разрыве. "Размазывание" решения на скачках для метода ТУО происходит на одной-двух разностных ячейках, а для метода Мак-Кормака — на трех-четырех. Шаг расчетной сетки Ах = 0,0025.
Течение газа в канале со ступенькой. Сверхзвуковой поток идеального газа (число Маха М = 3, коэффициент Пуассона 7 = 7/5) затекает в прямоугольный канал размером 3 х 1, на расстоянии 0,6 от левой грани которого располагается ступенька высотой 0,2. Начальные условия: р = 0,05, р = 0,5, ух = ьу = 0; граничные условия входа потока: р =1/7, р = 1, ух = 3, ьу = 0. На рис. 2 показаны результаты решения этой задачи методом Мак-Кормака и по разработанной схеме ТУО с разными значениями параметра диссипации е, а также численные результаты, полученные в ИПМ им. М.В. Келдыша [5]. Решение по методу ТУО при е = 0,01 обладает более высокой четкостью линий уровня, чем по схеме Мак-Кормака, — "размазывание" ударных волн
Рис. 1. Графики плотности (а), скорости (б) и давления (в), рассчитанные аналитическим методом (1) и численными методами TVD (2) и Мак-Кормака (5)
происходит на двух ячейках, в то время как для метода Мак-Кормака на трех-четырех ячейках. Для схемы ТУБ с е = 0,02 "размазывание" происходит на трех ячейках. Шаг расчетной сетки Ах = Ау = 0,0125.
Обтекание клина. Клин размером 3,2 х 2,2 с углом 30° обтекается сверхзвуковым потоком. Начальные условия: р = 1, р = 1,4,
г
Рис. 2. Результаты численного решения (линии уровня для безразмерной плотности), полученного методами Мак-Кормака (а), ТVD при е = 0,01 (б), е = 0,02 (в) и в ИПМ им. М.В. Келдыша (г) для задачи об обтекании сверхзвуковым потоком газа ступеньки в канале
Vх = ^у = 0; граничные условия: р =117, р = 8, ух = 8,25, г>у = 0. На рис. 3 показаны линии уровня плотности в момент времени £ = 0,2 с. Как и в случае задачи об обтекании ступеньки, метод ТУЭ обеспечивает более высокую четкость картины течения с меньшим "размазыванием" решения на ударных волнах. Шаг расчетной сетки по адаптивным координатам АХ1 = АХ2 = 0,03.
Метод на основе схемы ТУЭ обеспечивает более высокую точность численного моделирования, чем метод ленточно-адаптивных сеток на основе схемы Мак-Кормака, но время расчета у этих методов оказалось примерно одинаковым вследствие значительно более сложного алгоритма вычислений у схемы ТУЭ, несмотря на более крупный шаг по времени.
Выводы. Предложен алгоритм численного решения задач газовой динамики на основе схемы ТУЭ 2-го порядка и метода ленточных адаптивных сеток. Проведенное тестирование для ряда классических
Рис. 3. Результаты численного решения (линии уровня для безразмерной плотности), полученного методами Мак-Кормака (а), ^VD при е = 0,1 (б) и в ИПМ им. М.В. Келдыша (в) [5] для задачи об обтекании клина сверхзвуковым потоком газа
задач показало, что алгоритм схемы TVD обеспечивает более высокую точность численного моделирования, чем метод ленточно-адаптивных сеток на основе схемы Мак-Кормака, однако машинное время расчета задач с помощью этих методов примерно одинаково вследствие существенно более сложного алгоритма вычислений для схемы TVD, несмотря на более крупный шаг по времени.
Работа выполнена при частичной поддержке гранта Президента РФ МК-4234.2010.8.
СПИСОК ЛИТЕРАТУРЫ
1. Димитриенко Ю. И., Захаров А. А. Метод ленточных адаптивных сеток в газовой динамике. - М.: Изд-во НТЦ "Университетский", 2008. - 180 с.
2. ГильмановА. Н. Методы адаптивных сеток в задачах газовой динамики. - М.: Наука. Физматлит, 2000. - 248 с.
3. H a r t e n A. High resolution schemes for hyperbolic conservation laws. - 1982. -72 p.
4. R o h d e A. A. Computational study of flow around a rotating disc in flight. - PhD thesis, Florida, 2000. - 169 p.
5. Галанин М. П., Грищенко Е. В., Савенков Е. Б., Токарева С. А. Применение RKDG метода для численного решения задач газовой динамики. - Препринт ИПМ им. М.В. Келдыша РАН. № 52. - М., 2006. - 24 с.
Статья поступила в редакцию 20.01.11
Юрий Иванович Димитриенко родился в 1962 г., окончил в 1984 г. МГУ им. М.В. Ломоносова. Д-р физ.-мат. наук, профессор, заведующий кафедрой "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана, действительный член Академии инженерных наук. Автор более 160 научных работ в области вычислительной механики, нелинейного тензорного анализа, термомеханики композитов, математического моделирования в материаловедении.
Yu.I. Dimitrienko (b.1962) graduated from the Lomonosov Moscow State University in 1984. D. Sc. (Phys.-Math.), professor, head of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Full member of the Russian Academy of Engineering Sciences. Author of more than 160 publications in the field of computational mechanics, nonlinear tensor analysis, thermomechanics of composite materials, mathematical simulation in science of materials.
Михаил Николаевич Коряков родился в 1987 г., в 2010 г. окончил МГТУ им. Н.Э. Баумана. Аспирант кафедры "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана. Автор ряда работ в области вычислительной газодинамики.
M.N. Koryakov (b. 1987) graduated from the Bauman Moscow State Technical University in 2010. Post-graduate of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Author of several publications in the field of computational gas dynamics.
Андрей Алексеевич Захаров родился в 1983 г., в 2005 г. окончил МГТУ им. Н.Э. Баумана. Канд. техн. наук, доцент кафедры "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана. Автор 20 научных работ в области вычислительной газодинамики.
A.A. Zakharov (b. 1983) graduated from the Bauman Moscow State Technical University in 2005. Ph. D. (Phys.-Math.), assoc. professor of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Author of 20 publications in the field of computational gas dynamics.
Елтуган Кимашевич Сыздыков родился в 1956г., в 1979 г. окончил КГТУ им. А.Н. Туполева. Канд. техн. наук, заместитель генерального директора ОАО ГосМКБ "Радуга" им. А.Я. Березняка. Автор более 50 научных работ в области механики аэрокосмических конструкций.
Ye.K. Syzdykov (b. 1956) graduated from the Kazan State Technical University n.a. A.N. Tupolev in 1979. Ph. D. (Eng.), deputy general director of JSC GosMKB "Raduga". Author of more than 50 publications in the field of mechanics of airspace structures.
