CC BY
51
УДК 519.63+532.5-1/-9 Модификация 2В численного кода для газодинамических течений в полярных координатах
Е. А. Филистов
Кафедра Физики Московский государственный строительный университет Россия, 129337, Москва, Ярославское шоссе, 26
Модифицирован двумерный численный код для моделирования совершенного газового потока. Код построен на основе эффективного экономичного метода расчета нестационарных течений идеального газа с использованием полностью консервативной разностной схемы строго дивергентных газодинамических уравнений в эйлеровых переменных в полярных координатах.
Ключевые слова: газодинамика, численное моделирование.
1. Введение
В общей структуре течения как различных газодинамических задач, так и в прикладных задачах астрофизики, как правило, присутствуют разрывные решения, которые с математической точки зрения представляются наиболее интересной особенностью систем нелинейных гиперболических уравнений. Использование математических особенностей гиперболических уравнений позволяет строить эффективные численные методы, несмотря на то, что получающиеся решения должны быть почти разрывными. Общее направление развития численных методов для решения уравнений газовой динамики лежит на пути уменьшения численной вязкости при сохранении монотонности конечно-разностной схемы расчёта течений с разрывами. При этом одним из наиболее важных свойств разностной схемы является требование их полной консервативности, то есть выполнение дискретных аналогов законов сохранения в дифференциальной форме. Известно [1—4], что схемы, не обладающие этим свойством, могут давать решения, весьма далёкие от истинного (в частности, ударные волны, движущиеся с неправильными скоростями). Для обеспечения консервативности схемы естественно использовать запись схемы в потоковом виде, то есть когда искомые решения — функции получаются в результате разностного дифференцирования функции соответствующего потока. Однако в настоящее время, как правило, при моделировании совершенного газового потока в рамках осесимметричного приближения используется полудивергентная форма газодинамических уравнений в эйлеровых переменных (см., например, [5]).
В данной работе представлена существенно дивергентная форма системы газодинамических уравнений в эйлеровых переменных в полярных координатах. Для численных расчетов сложных нестационарных течений применен эффективный экономичный метод с использованием полностью консервативной разностной схемы аппроксимации дифференциальных уравнений. Для моделирования совершенного идеального газового потока построен двумерный численный код решения системы уравнений газовой динамики.
2. Основные уравнения 2.1. Уравнения газодинамики в декартовой системе координат
Стандартные уравнения газодинамики для нерелятивистской среды (уравнения Эйлера) дают возможность определения временной эволюции макрохарактеристик динамической системы. Основными физическими параметрами движений
Статья поступила в редакцию 15 ноября 2012 г.
идеального газа в ньютоновской теории являются вектор массовой скорости движения газа v (r,t), плотность газа р (r,t), газодинамическое давление P (r,t), где г — радиус-вектор «жидкой частицы» вещества, t — время. Удобно иметь дело с дивергентной формой уравнений классической газодинамики, которые в тензорных обозначениях и без учёта массовых сил и диссипативных процессов могут быть записаны в виде [6]:
для плотности потока массы газа:
dtp = -д, (pit); (1)
для плотности потока импульса:
dt (pvz ) = -д3 (pvV + P g*); (2)
для плотности потока суммарной энергии:
9t {py+ре) = О9?+ре+р)); (3)
где е = e(P,p) — удельная внутренняя энергия на единицу массы газа, §ij — метрический тензор. Уравнения (1)—(3) записаны в ковариантной форме, чтобы их можно было использовать в произвольной криволинейной системе координат.
Для замыкания системы уравнений (1)—(3) к ним надо присоединить уравнение состояния
p = p(p, T) (4)
(например, уравнение состояния идеального газа P = (7 — 1)ер), связывающее между собой давление, плотность и температуру среды, а также уравнение непрерывности для энтропии
dt(Ps ) = —д, (pSv?), (5)
где произведение pSv3 представляет собой плотность потока энтропии; энтропия единицы массы идеального газа определяется выражением S = Cv ln(P/p7), 7 = Cp/Cv > 1 — показатель адиабаты, Ср и Cv — удельные теплоёмкости газа при постоянном давлении и объёме. В случае адиабатического движения газа (при отсутствии поступления энергии извне) dS/dt = 0. В случае изотермического потока уравнение энергии (3) исключается (вследствие радиоактивного излучения энергия не сохраняется), а давление определяется соотношением P = с^р, с^ = — T,
р
Ст — изотермическая скорость звука, р — молекулярный вес, — — универсальная газовая константа. На практике, впрочем, вместо изотермических уравнений можно использовать газодинамические уравнения с уравнением состояния идеального газа. Отношение теплоёмкостей при этом берётся достаточно близким к единице (скажем, 7 =1, 001), что, согласно известной формуле из статистической физики 7 = (а + 2)/а, соответствует очень большому а — числу степеней свободы частиц, из которых состоит газ, так что если газ нагревается (например, ударной волной), то нагрев распространяется по всем степеням свободы и оказывает лишь незначительное влияние на давление и температуру.
Уравнения (1)—(5) представляют собой систему нелинейных гиперболических уравнений и, очевидно, что для большинства задач газодинамики не представляется возможным получить аналитические их решения — возможен лишь численный расчёт. Большим преимуществом при таких расчётах обладают полностью консервативные (дивергентные) разностные схемы, для которых справедливы законы сохранения в разностном виде. Для случая декартовых координат однородная система двумерных уравнений газовой динамики в консервативной (дивергентной) форме может быть записана как
dtq + dxF + дуG = 0 .
(6)
Векторы консервативных переменных q и потоков Е, С определяются следующими выражениями:
q
( р\
ри
pv \РЕ)
(
F
ри
ри \
2 +P
puv puh
(
G
pv puv
\
pv
2 + P
\ pvh /
где
V = {и ,v}
tf2
tf2
PP
E =— + £ h = — + £ +- , £ = --—
2 2 p (7- l)p
Описывающая движение газа газодинамическая система уравнений (6) представляет собой систему квазилинейных дифференциальных уравнений в частных производных гиперболического типа. Гиперболичность системы уравнений газовой динамики означает, что матрицы Якоби дЕ/дq, дС/дq имеют действительные собственные числа и могут быть приведены к диагональному виду, то есть имеют полный набор линейно независимых собственных векторов. Свойство гиперболичности очень важно для построения конечно-разностных схем для решения газодинамических уравнений. Характеристические свойства гиперболической системы (6) определяют общее направление развития численных методов для решения этих уравнений и используются для уменьшения численной вязкости при сохранении монотонности схемы.
В газодинамике принято различать два вида течений — потенциальные и вихревые (частным, но распространённым случаем последних являются турбулентные движения). Обычная газодинамика изучает в основном потенциальные течения, чаще всего имеющие определённую симметрию движения (одномерные, плоские, центрально-симметричные). В реальной газодинамической среде такие правильные движения встречаются сравнительно редко; чаще всего движение газа имеет вихревой и турбулентный характер.
2.2. Дивергентная форма уравнений газодинамики в полярных
координатах
В полярных координатах (г, tf) уравнения Эйлера (6) также можно записать в дивергентной полностью консервативной форме. Связывая старые и новые векторы скоростей и потоков соотношениями
q = rq, F = r(F ■ cos tf + G ■ sin tf), G = (—F ■ sin tf + G ■ cos tf), находим
dtq + drF + d#G = 0 ,
(7)
где компоненты вектора консервативных переменных q и составляющие Е и (С векторов потоков имеют вид
Í Р\ рм pv \РЕ)
( pvr \ pvr u + Pai2 pvr v + Pai3 \ phvr J
G
Í pv# \ pv$u + P&12 pv#v + P613 V Phv^ )
Здесь, как и прежде, р = р(г,$,£) — плотность; V = V = (иг — век-
тор скорости; Р = Р (г, ¿) — давление, а для сокращения записи используем
г
г
следующие обозначения:
р = 7 - 1, Р = 2 -а12 = cos ft, ai3 = sin ft, bi2 = — sin ft , bi3 = cos ft,
u = vrai2 + v$bi2 , v = vrai3 + v$bv¿ , v = {vr ,v#}T .
F д G
Матрицы гиперболичности л = -7-г и В = для рассматриваемой системы
aq aq
уравнений (7) имеют вид
Л
с такими значениями элементов:
Í 0 ai2 ai3 0 0 bi2 bi3 0
(121 (122 423 pai2 , в = 1Г Ь21 b22 b23 pbi2
(131 а>32 033 pai3 b3i Ь32 b33 pbi3
\Й41 Ü42 Ü43 JVr j \b4i Ь42 b43 1Щ j
v?
а,21 = — uvr + Ру а\2, v2
a^i = —vvr + ру а13,
(ь
а41 = — IV I h — ру I
г2
b2i = — UV# + Ру bi2,
v2
Ьз1 = —vv-в + ру 6i3,
641 = — ^ ( h — py
(h — ^т)
022 = vr + /3 uai2, аз2 = щ + P uai3, a42 = hai2 — pvr u,
b22 = ^^ + P ubi2,
Ьз2 = —Vr + /9u&i3, 642 = hbi2 — pv#u,
a23 = — v$ + ^va^, азз = ^ + P vai3, a43 = hai3 — pvr v;
&23 = ^ + P vbi2,
b33 = v# + P vbi3,
b43 = hbi3 — pv#v .
Вычисления приводят к стандартным выражениям для действительных собственных значений матриц гиперболичности:
Лд = { vr — с; vr; vr; vr + с } , Лв =
v$ — с v$ v$ + с ; ; ;
где с = л/7Р/р — адиабатическая скорость звука (в изотермическом случае используется изотермическая скорость звука ст). Заметим, что для построенных
матриц л и В тождественно удовлетворяются соотношения Е = л - 4 и С! = в ■ ц. Кроме того, как матрицы простой структуры с действительными собственными числами, матрицы л и В обладают полным набором левых \т и правых гт собственных векторов.
Общий случай системы нелинейных гиперболических уравнений содержит информацию о распространении и взаимодействии нелинейных волн. Как показано в работах [7-9], такую информацию можно получить если точное решение нелинейной задачи заменить её приближенным (линейным) решением с помощью специального выбора газодинамических комплексов и ¿¿я между соседними ячейками в разностной схеме.
Выбирая матрицы Д(ц^) и в(цв) такими, чтобы удовлетворялись соотношения Ея -Еь = л(ц*а)^ (¿я -Сь) и ССя - ССь = в(ц*в )■ (Ся -¿Ы, и решая следующие из них уравнения (здесь и далее введено обозначение д = гр)
a(qvr )л + v*rAv#aa(q)A — v#aa(qvr )л — v*aa(qv#)A = 0,
a(qv2r)A + vfAa(q)A - 2v*raa(qvr)л = 0 , a(qhvr)л + hAv^AA(q)A - hAA(qvr~)a - v*AA(qh)A = 0, для матрицы и
a(pvrv#)B + v*.Bv*rba(p)B - v*^ba(pvr)b - v*Ba(pv$)b = 0 ,
a(pvj)b + vlbA(p)b - 2v*bA(pv#)b = 0,
a(phv$)b + h*bv**bA(p)b - h*ba(pv$)b - v**ba(pK)b = 0,
для матрицы b(q*b) A(*) = (*)r - (*)l — обозначение разности соответствующих параметров (индексы л или B характеризуют принадлежность вычисляемой величины потоку F или G соответственно), находим соотношения для средних значений газодинамических комплексов q^ и q^:
* ,--л/qar var + VqAL vAL
qA = v<laR • <laL , VA = - v
hA =
л/OaR + л/OaL
\¡ QaRhaR + V QaL haL jOaR + jOaL '
* = I Уча! • cAL + /qjr • c\R p q*A / ^
Л V J^1 + 2 (JqJE + JqIR)r \ A> '
* ,--л/pbRvbR + y/pbLvbL
pb = vpbr • pbl , vb = ■ v
hr =
x,fPbR + t/PBL PbRhbR + V PbLhbL jPbR + /PBL '
c* = /JPB¿ • c2bl + /1ВД • 4д + p p*B /A ^ y
B V /Abl + JPbR 2 (/pBx + jpbr)2 v B/
= 0,5($вя + $вь). Такие значения ¿с^ и ^ при условии К = Ь приводят к
- ^Е ас
выполнению естественных соотношений = и о(с{*в) = .
Приращения характеристических переменных авА = ) • А(ссд) и а,^ =
!™(Чв) • А^в) для системы уравнений газовой динамики (7) имеют вид
А^д4 = ^ (А(г • Р)л Т ч\с\а(ьг)а) , А^3 = тщ2 (сл • аш ± 9лс*аа(ь$)а - а(г • Р)л) , А%'4 = 212 (А(г • Р)в Т Р*вс*ва(г • ад)в),
аЛ'3 = 1 2
(c*b2a(q)b т р*вс*вa(r • vr)в - a(r • P)b)
"В — 2г*2
Промежуточные значения скорости звука с*д и с* вычислены по формуле Эйн-фельдта [10], не содержат разностей больших величин и, как следствие, не приводят к потере точности.
3. Метод аппроксимации уравнений газодинамики
Результаты сравнения различных методов численного моделирования астрофизических задач (заметим, что впервые такое сравнение было проведено в работе [11]), в которых ударные волны и контактные разрывы играют существенную роль, показывают, что наиболее корректные решения дают TVD методы (Total Variation Diminishing), предложенные в работах Роу [7-9], а также в работах [12] и [13]. Однако многочисленные тестовые расчёты различных газодинамических задач (см., например, [11]) по схеме Роу показывают, что в ряде случаев эта схема может давать решение, соответствующее нефизическим разрывам, и, с другой стороны, сильно размазанное решение около сильных и слабых разрывов. Причиной возникновения неэволюционных разрывов в модели Роу является низкая численная вязкость схемы.
Наиболее простой и эффективный метод предотвращения возникновения неэволюционных стационарных ударных волн путём увеличения численной вязкости схемы предложен в работе [10]. Идея заключается в замене в схеме Роу действительных собственных значений матриц гиперболичности Ax(q*) на min(Ax(q*), А1) для и > 0 (для исключения нефизических скачков необходимо взять min от
(А^А^^)), что увеличивает |АХ| и, как следствие, увеличивает численную вязкость только для неэволюционных волн, в то время как на эволюционные волны эта операция влияния не оказывает. Для и < 0 необходимо произвести замену A3(q*) на max(A3(q*),A3+1) (для исключения нефизических скачков надо взять max от (^3+1, . В результате получается схема Роу-Эйнфельдта (постро-
енная на базе схемы Лакса-Фридрихса [14,15]) на равномерной сетке.
Следуя работам [7-10], строим 2D численный код для для однородной системы уравнений газовой динамики в меридиональном сечении
q"+1 - qt,3 + ^1+1/2,3 - Fг-1/2,3 + Gг,3+1/2 - Gг,0-1/2 _ о
г + Дг + а§ ,
где т, Дг, Д$ — шаги сетки. Здесь радиальный поток строится следующим обра-
зом:
Fг+1/2,3 _ г'° 2 1,3 - 1 |^(q*+1/2,j)|AsI+1/2,jr^(q*+1/2,;7) ,
д^,'^1/2,3 _ -((г*+1 f*+1,3 - т 4**+1/2,3c*+1/2,3(vrz+1,3 - Vr^,3)) ,
%+1 /2,3
а 2,3 Двг+1/2,з
2-(С*г+1/2,з(Ь+1,3 - Чг,;) ± С*+1 /2,;Я* К г+1,; - V#%,3)-(Гг+1рг+1;-ГгРг,;))
+1/2,3 4 '
vдг+1,з Уг+1,з + Уг,з
Q1+1/2,3 _ Vh+1,3 <1г,3 , v*+1/2,3 _
л/Чг+1,3 +
,3 hi+1,3 + V%~3 h h*+1/2,3 _ + ^
/
_ Cj,3 + VV+73 cj+1,3 + У У^г,3 дг+1,3 _ _)
"г+1/2,3 + v^ + 2 (^ + V4+3)2 ( +1,3 ,3)
Набор правых г^ (4*) собственных векторов матрицы Л содержится в матрице
/
и* - с*а12
V* - с*а13
Н* - с*ь*
V*2 2
V*2 2
Н* + с*у*
\
Аналогично строится выражение для меридионального потока с набором правых ) собственных векторов матрицы В
(
1
1
1
1
\
и* — Г*Ь* с о12 и* - с*&13 и* + с*Ь*13 и* + С*&?2
(4*) = V* — Г*Ь* С °13 V* + с*Ь*12 V* - с*Ъ*12 V* + С*Ь*13
V*2 * * — - С V* V*2 1 * * — + С V*
Н* Н* * * - с V** Н* + сху
На основе схемы Роу-Эйнфельдта с помощью некоторых ограничивающих процедур обычно строятся схемы повышенного порядка аппроксимации (схемы 2-го - 3-го порядка аппроксимации) на гладких решениях [12,16]. Такие схемы остаются монотонными, то есть не дающими нефизических осцилляций около разрывов. Дополняя потоки схемы Роу-Эйнфельдта их антидиффузионными добавками можно получить схему Роу-Эйнфельдта-Ошера [13,17].
С учётом условия монотонности схемы в линейном приближении шаг по времени в схеме Роу-Эйнфельдта-Ошера для уравнений газовой динамики на каждом временном слое выбирается в виде
Т = (Г Ш1П
{|АГ|}
4
|/ 5 - <р +(1 + у)?
где Н — минимальный пространственный шаг, 0 < а ^ 1 — параметр схемы, и /3 — свободные параметры, причём последний должен удовлетворять соотношению
3 — <р
1 < Р ^ /9Шах, где /9тах = -. Схема имеет третий порядок аппроксимации,
1 - у
если <р = 1/3 [11,13].
Схема Роу-Эйнфельдта-Ошера относится к ТУБ методам решения систем нелинейных гиперболических уравнений газовой динамики. Свойства ТУБ схем подробно исследованы и хорошо математически обоснованы. Преимущество этих схем и в относительной простоте программирования, что также является весьма важным моментом. Двумерные тестовые расчёты по схеме Роу-Эйнфельдта-Ошера [18] показывают, что приближенное решение передаётся с большой точностью и, практически, не даёт нефизических осцилляций за ударными волнами. Автором на основе описанного метода были проведены численные расчёты сверхзвуковых течений газа, развивающихся под влиянием внешнего спирального гравитационного потенциала [19,20]. Расчёты показали, что метод достаточно хорошо передаёт контактные разрывы и почти не размазывает ударные волны.
1
1
1
1
*
*
*
*
*
*
*
*
4. Заключение
Для численных расчётов сложных нестационарных течений идеального газа представлена существенно дивергентная форма системы газодинамических уравнений в осесимметричном приближении в полярных координатах, которая используется для построения явной полностью консервативной разностной схемы третьего порядка аппроксимации по пространственным координатам в областях гладкости решения.
На основе описанного метода автором разработан двумерный численный код для моделирования газодинамических течений на эйлеровой цилиндрической сетке, с помощью которого были проведены численные расчёты идеального сверхзвукового потока газа, развивающегося под влиянием внешнего спирального гравитационного потенциала [19,20]. Расчёты показали, что метод обеспечивает хорошую пространственную локализацию разрывов в нестационарном потоке газовой среды.
Литература
1. Тихонов А. Н., Самарский А. А. О сходимости разностных схем в классе разрывных коэффициентов // ДАН СССР. — 1959. — Т. 124. — С. 529-532. [Tikhonov A. N., Samarskii A. A. On the Convergence of Difference Schemes in the Class of Dischontinuous Coefficients // Dokl. Akad. Nauk SSSR, Moscow, 1959. — Vol. 124. — P. 529 — 532. ]
2. Самарский А. А, Попов Ю. П. Разностные методы решения задач газовой динамики. — М.: Наука, 1980. — С. 352. [Samarskii A. A., Popov Yu. P. Difference methods for solving problems of gas dynamics. (Raznostnye metody resheniya zadach gazovoi dinamiki), Moscow: Nauka 1980. — P. 352. ]
3. LeVeque R. J. Numerical Methods for Conservation Laws. — Basel, Birkhauser-Verlag, 1990. — P. 228.
4. Kimoto P. A., Chernoff D. F. Convergence Properties of Finite-Difference Hydrodynamics Schemes in the Presence of Shocks // Astrophys. J. Suppl. Ser. — 1995. — Vol. 96. — Pp. 627-641.
5. Studies of Equilibrium Configurations for a Gaseous Cloud Near a Gravitating Center / M. V. Abakumov, S. I. Mukhin, Y. P. Popov, V. M. Chechetkin // Astronomy Reports. — 1996. — Vol. 40, No 3. — Pp. 366-377.
6. Филистов Е. А. МГД ударные волны в самогравитирующем газе // Вестник РУДН, Серия «Физика». — 2001. — № 9(1). — С. 71-73. [Filistov E. A. MHD shocks in flow of gravitating gas. Bulletin of Peoples' Friendship University of Russia, Moscow, 2001. — No 9(1). — P. 71-73 ]
7. Roe P. L. The Use of the Riemann Problem in Finite-Difference Scheme // Lect. Notes Phys. — 1980. — Vol. 141. — Pp. 354-359.
8. Roe P. L. Approximate Riemann Solvers, Parameter Vectors and Difference Schemes // J. Comp. Phys. — 1981. — Vol. 43. — Pp. 357-372.
9. Roe P. L. Characteristic-Based Schemes for the Euler Equations // Ann. Rev. Fluid Mech. — 1986. — Vol. 18. — Pp. 337-365.
10. Einfeldt B. On Godunov-Type Methods for Gas Dynamics // SIAM J. Numer. Anal. — 1988. — Vol. 25, No 2. — Pp. 294-318.
11. Кузнецов О. А. Численное исследование схемы Роу с модификацией Эйн-фельдта для уравнений газовой динамики // Препринт ИПМ им. М. В. Келдыша РАН. — 1998. — № 43. — С. 44. [Kuznetsov O. A. Numerical Analysis of the Roe's Scheme with Einfeldt's Modification for the Equations of Gasdynamics. — 1998. — No 43. — P. 44. — URL: http://library.keldysh.ru/preprint.asp? id=1998-43 ]
12. Harten A. High Resolution Schemes for Hyperbolic Conservation Laws //J. Comp. Phys. — 1983. — Vol. 49. — Pp. 357-393.
13. Chakravarthy S. R., Osher S. A. A New Class of High Accuracy TVD Schemes for Hyperbolic Conservation Laws // AIAA Pap. — 1985. — No 85-0363. — Pp. 1-11.
14. Lax P. D. Weak Solutions of Nonlinear Hyperbolic Equations and Their Numerical Computations // Commun. Pure Appl. Math. — 1954. — Vol. 7, No 1. — Pp. 159197.
15. Friedrichs K. O. Symmetric Hyperbolic Linear Differential Equations // Commun. Pure Appl. Math. — 1954. — Vol. 7, No 2. — Pp. 345-392.
16. Sweby P. K. High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws // SIAM J. Numer. Anal. — 1984. — Vol. 21, No 5. — Pp. 995-1011.
17. Вязников К. В., Тишкин В. Ф., Фаворский А. П. Построение монотонных разностных схем повышенного порядка аппроксимации для систем уравнений гиперболического типа // Мат. Моделирование. — 1989. — Т. 1, № 5. — С. 95-120. [Vyaznikov K. V., Tishkin V. F., Favorskii A. P. Construction of Monotone Difference Schemes of an Increased Order of Approximation for Systems of Equations of Hyperbolic Type // Mat. Model. — 1989. — Vol. 1. — No 5. — P. 95-120 ]
18. Woodward P. R., Colella P. The Numerical Simulations of Two-Dimensional Fluid Flow with Strong Shocks // J. Comp. Phys. — 1984. — Vol. 54, No 1. — Pp. 115173.
19. Filistov E. A. Polygonal Structure of Spiral Galaxies // Astronomy Reports. — 2012. — Vol. 56, No 1. — Pp. 9-15.
20. Filistov E. A. Correlation of a Bar and Spirals at Formation of Galaxies // Astronomical and Astrophysical Transactions. — 2012. — Vol. 27, No 2. — Pp. 215-220.
UDC 519.63+532.5-1/-9
Modification of the Two-Dimensional Numerical Code for Gas-Dynamical Flows in Polar Coordinates
E. A. Filistov
Department of Physics Moscow State University of Civil Engineering 26, Jaroslavskoje shosse, Moscow, 129337, Russia
The numerical method for solution of the gas-dynamical equations in strict divergent form has been modified. The two-dimensional numerical code for perfect non-stationary gas-dynamical flows simulation on the polar grid is constructed. This code is based on the explicit quasimonotonic high resolution TVD-scheme.
Key words and phrases: gas-dynamics, numerical simulation.
