Численное моделирование внутрикамерных нестационарных турбулентных течений. Часть 1

Для исследования газодинамических процессов, протекающих в различных технических устройствах и установках, в том числе в твердотопливных ракетных двигателях (РДТТ), важное значение имеет численное моделирование трехмерных внутренних турбулентных потоков. Трудности, связанные с постановкой физического эксперимента, выдвигают на ведущий план развитие методов математического моделирования таких течений. Улучшение рабочих характеристик РДТТ, снижение материальных затрат при разработке новых двигателей, сокращение сроков их отработки в значительной степени зависят от корректности численных исследований процессов, протекающих в двигателе, от реализуемой точности проводимых расчетов. Моделирование процесса выхода РДТТ на стационарный режим работы представляет задачу высокой сложности. Сложность моделирования течения продуктов сгорания в РДТТ обусловлена турбулентным характером течения, нестационарностью протекающих в тракте двигателя процессов, влиянием градиентов давления [1–6]. Значительный интерес также представляет и изучение влияния вдува продуктов сгорания в камеру двигателя на параметры турбулентного потока в каналах с геометрически сложной формой. Разработка новых методик для численного решения задач гидрогазодинамики всегда является актуальной задачей. С одной стороны, это объясняется необходимостью совершенствования известных, хорошо зарекомендовавших себя методов, в плане более эффективного использования ими увеличивающихся вычислительных возможностей компьютеров, с другой стороны, программной реализацией новых, обусловленных современными практическими потребностями более сложных постановок решаемых задач.

Конфигурации внутренних поверхностей реальных двигателей имеют достаточно сложную геометрию. Течение продуктов сгорания в РДТТ имеет существенно дозвуковые области, области высоких дозвуковых скоростей и область соплового сверхзвукового течения. Проведение численного моделирования потока газа во всех этих областях в рамках единого алгоритма является сложной задачей, но успехи современной вычислительной гидрогазодинамики и растущая вычислительная мощность компьютерной техники побуждают к разработке математических моделей и алгоритмов, обеспечивающих сквозной расчет внутренних течений в ракетных двигателях.

В данной работе для моделирования турбулентного течения в тракте модельного твердотопливного двигателя предложен вычислительный алгоритм, основанный на численном интегрировании уравнений сохранения, записанных в цилиндрической системе координат [5, 7–10]. Во всех существенно нестационарных течениях возникают пульсации параметров потока, что накладывает определенные требования к используемой методике. Предложенный алгоритм пригоден для сквозного расчета всего тракта двигателя, как дозвукового течения в камере двигателя, так и сверхзвукового в сопле, при наличии пульсаций давления без введения искусственной вязкости. Алгоритм основан на модифицированном методе расщепления векторов потоков, относящимся к классу методов, использующих подход Годунова [11–14]. Для обеспечения достоверности получаемых результатов при создании методик численного исследования, предназначенных для моделирования течений в зарядах реальных конфигураций, проводится проверка корректности математической модели и тестирование вычислительного алгоритма на модельных задачах. В работе приведены результаты исследований внутрикамерного нестационарного течения сжимаемого вязкого газа в модельном ракетном двигателе.

1.    Постановка задачи

Для моделирования трехмерного течения продуктов сгорания в камере РДТТ использовалась система уравнений гидромеханики в цилиндрической системе координат. В векторном виде система записывается следующим образом dQ + F + F + 1 dF + Fee =1 (^rL, + ^ + |r) + Lm , (1) ∂t ∂x ∂r r ∂θ r ∂x ∂r ∂θ где Q – вектор гидромеханических параметров (ГМП); Fx, F , Fθ, Fθθ – векторы конвективных потоков; Lx , L , Lθ , Lθθ – векторы диффузионных потоков; x, r, θ – продольная, радиальная и угловая координаты; t – время. Векторы ГМП, конвективных и диффузионных потоков имеют вид

Q=	ρ ρu ρv ρw	, F x =	ρu pu 2 + p ρuv ρuw	, F r =	ρv ρuv pv2 + p ρvw	, F e =	ρw ρuw pvw pw2 + p	p —P , Fee— r	v uv v2 — w2 2vw	, (2)
	ρe		ρuh		ρvh		ρwh		vh

L x

τ xx τ xr τ xθ

UT xx +VT xr +WT xe

τ rx

Lr—

τ rr

Le—

τ θx τ θr

τθθ итех+vTer+wTee

τ rθ

UT _rx +VT _rr +WT r_e

	0
1	0
ee —- r	- τ θθ τ θr
	0

Здесь u, v, w – компоненты вектора скорости, p – давление, ρ – плотность, e – удельная энергия, h – удельная энтальпия, τ _xx , τ _rx , τ _θx , τ _xr , τ _rr , τ _θr , τ _xθ , τ _rθ , τ _θθ – компоненты тензора вязких напряжений.

Компоненты тензора вязких напряжений вычисляются так

T xx — 2M^du - ² M d ^iv Q, dx 3

T rr — 2M^dv - ² p div Q , dr 3

^T ee ^{— 2} M ( ^{1 d}w + v ] ^{- 2} M ^divQ ,

\ r dv r J 3

∂v ∂u

T xr — T rx — M V + V , ∂x ∂r

( 1 du dw\

^T x^{e — T e}x ^{— M}^ rdv + dxj.         (4)

f 1 dv dw w\

Tre = Ter — M rdv + -^ - —I, где m — коэффициент динамической вязкости, Q — (u; v; w) - вектор скорости. Дивергенция вектора скорости в цилиндрической системе координат определяется формулой d> q du   dv   1 dw  v

∂x   ∂r   r ∂θ   r .

Аналогичная система уравнений использовалась в работах [5, 7–10]. Температура и полная удельная энтальпия вычислялись по соотношениям

T —

-

C _v

u ² + v ² + w ² 2

,

h —

^u 2 ^{+ v} 2 ⁺ w ²

о      ⁺ C p T •

Уравнение состояния записывалось в следующем виде

u ² + v ² + w ²

p ^— P ^{( v -} 1) I^е--2------) •

Здесь C _v – удельная теплоемкость газа при постоянном объеме, C _p – удельная теплоемкость газа при постоянном давлении, T – температура, R – удельная газовая постоянная, ν – показатель адиабаты.

2.    Вычислительный алгоритм

Для расчета векторов конвективных потоков Fx , Fr , Fθ (2) использовалась модифицированная схема Стигера – Уорминга (схема расщепления векторов потоков), предложенная в работах [11, 12]. Применяя уравнение состояния (7), из выражений для векторов (2) можно исключить давление [11]. С использованием взвешивания по Фавру f = pf/P вектор Fx будет записан pu

F x

P ^{( v} - 1) ^e -

( v

-

3) u ² + (v

-

1) v ² + (v

-

puv

1) w ² ^

puw puh

Вектор F _x можно представить следующим образом

F x = J x Q,

∂Fx где Jx = ——— матрица Якоби. Используя (8), матрица Jx определяется так [15] ∂Q

Jx 0 v-3 ~2 1 v- 1 /~2 1 ~2\ 2 u2 +--—(v + w2) —uv —uw 0 1 — (v — 3) u v w h 0 — (v — 1)v u 0 0 0 — (v — 1) w 0 u 0 0 v — 1 0 0 0 .       (10) Собственные числа матрицы Jx равны [11] Ax1   Ax2   ^x3   u , Ax4 u + C , Ax5 u — c , (11) где c = (vp/p)r/ - скорость звука. Вектор Fx может быть записан в следующем виде [11]

F_x = J_XQ = G_xК_хG - 1 Q , x x     xx x

где G x – матрица правых собственных векторов Якобиана J x ; G _x ^{- 1} – матрица левых собственных векторов. Собственные числа Якобиана (11) являются диагональными элементами матрицы Л _х

	u	0	0	0	0
	0	u	0	0	0
Л x =	0	0	й	0	0
	0	0	0	u + c	0
	0	0	0	0	u c — c c

Тогда вектор F x , с использованием собственных чисел может быть представлен [11]

Fx = — x 2v	2( v — 1)^ x1 + A x4 + A x5 2(v — 1)A xi u + A x4 (u + c) + A x5 (u — c) 2( v — 1) A x1 v + A x4 v + A x5 v 2(v — 1)A xi w + A x4 w + A x5 w	.                  (14)
	2(v — 1)A xi h + A x4 h + A x5 h

Основой схемы Стигера – Уорминга, как и других алгоритмов, использующих подход Годунова, является расщепление вектора потоков (14) на две части

F = F+ + -~ x x + x , где вектор Fx+ соотносится с вектором, сформированным из положительных собственных чисел матрицы Jx . Соответственно вектор Fx- соотносится с вектором, сформированным из отрицательных собственных чисел. Матрицу Лх можно расщепить аналогично

Л х = л + + л ^- .

Тогда векторы Fx, Fx , входящие в выражение (15), примут вид р + _ ( +Х л+ ( +Х  1  +         —-    ( -X -~ ( —^ 1

F x = IG xJ ^Л х IG xJ    Q , F x = Gxx) ^Л x IG xJ Q .          ⁽¹⁷⁾

Все собственные числа A +T , A x , являющиеся элементами матриц Л + , Л _x , определяются следующим образом [11]

+ х^-     х+   ^A xi^

^A x ^A x + ^A x ,        ^A x        2     ,

λ _x

λ _x - | λ _x |

Используя (14) – (18), векторы можно записать [11]

2(v

-

1)A + + A+ + A+

F ⁺ x

ρ

2(v

-

1)A +1 U ⁺ + A+ (u ⁺ + c + ) + A+ (u ⁺

-

2v

2(v

-

2(v

-

2(v

F x

ρ

2(v

-

2v

1)A xi V+ + A x4v+ + A x₅ v +

1) ^A +1W ⁺ + ^A +4 w ⁺ + ^A +5 w ⁺

-

.■w

.■w

^w

1)A +1 h ⁺ + A+h ⁺ + A +₅ h ⁺

2(v

-

^{1) A} x1 + ^A;

‘x4

+ A

‘x5

^{1) A} x1 u + ^A x4 ⁽ u + c ) + ^A x5 ^{( u}

2(v

-

-

2(v

-

2(v

^{7 +} )

cw )

,

1)A xi V ^- + A x4 V ^- + A x₅ v

l ^{) A} x1 w + ^A x4 w + ^A x5^ ^W

w

w

w

.

-

^{1) A} x1 h + ^A x4 h + ^A x5 h

Для численного интегрирования системы основных уравнений (1) по пространству применялся метод конечных объемов. Для расчета векторов потоков F i — 1/2 (Q) и F i+1/2 (Q) , соответственно на левой и правой гранях I i -го контрольного объема, необходимо определить векторы параметров Q i - 1/2 и Q i+1/2 . Параметры на границах I i -го контрольного объема Q ^- _{— 1} / ₂ и Q i+1/2 принимались равными в центре ячейки. На рис. 1 показан I _i -ый контрольный объем и параметры потока газа на его границах.

Рис. 1 . Параметры потока газа Q _{i— 1}/₂ и Q i+1/2 на гранях I i -ого контрольного объема

Далее используя выражения (15), (19), (20), вычисляем значения потока Fx. Аналогичные процедуры производятся для векторов Fr и Fθ . Векторы потоков Lx , Lr , Lθ определяются с использованием центрально-разностной схемы апроксимации уравнений (3) – (5). Затем рассчитываются элементы вектора переменных Q и, соответственно, определяются все параметры течения.

Турбулентность моделировалась методом крупных вихрей, являющимся одним из наиболее перспективных методов расчетов турбулентных течений. Данный метод широко используется в практике создания численных методик. Например, он применялся в работах [16–18]. Для моделирования пристеночных течений использовалась локальная модель вихревой вязкости (Wall-Adapting Local Eddy-viscosity, WALE) [17].

3.    Результаты расчетов

Расчеты проводились для модельного РДТТ, представленного на рис. 2. Ракетный двигатель имеет заряд твердого топлива зонтичного типа. Длина двигателя составляла 2,5 м. В начальный момент времени температура газа в камере составляла T ₀ = 294 K и давление p ₀ = 10 ⁶ Па. Шаг по времени задавался равным ∆t = 5, 0 · 10 ^{- 7} c . Скорость горения топлива определялась по формуле ω = ω 0 (p/10 ⁵ ) ^ϑ , где ω 0 – скорость горения при атмосферном давлении, ϑ – константа. При апробации вычислительной методики полагалось w = 0 .

Рис. 2 . Схема твердотопливного двигателя

На рис. 3 а), 3 б), 4 а) приведены графики изменения осевой скорости, радиальной скорости и давления в камере на начальном отрезке времени работы двигателя в четырех точках камеры двигателя. В начальной стадии работы двигателя на графиках хорошо заметен ударно-волновой характер процессов, протекающих в камере РДТТ, это связано с отражением волн от элементов камеры и поверхности топлива. Процессы, протекающие в двигателе на начальной стадии работы, существенно нестационарны. Величины параметров потока имеют сильно изменяющийся характер, особенно это заметно на графиках параметров в зазоре между передним днищем и зарядом. На рис. 4 б) приведено изменение давления по времени до выхода двигателя на стационарный режим работы. Время выхода на режим составило t ≈ 1, 1 с. Для этого момента времени на рис. 5, 7 представлены распределения газодинамических параметров течения продуктов сгорания по всему тракту двигателя, как в камере, так и в сопле. На рис. 5 приведены линии тока течения, из которых видно образование вихрей в зонтичной области камеры и областью за утопленной частью сопла. Распределения осевой скорости и давления приведены на рис. 6, 7, на которых наблюдаются большие градиенты параметров потока продуктов сгорания.

а)

б)

а)

Рис. 4 . Изменение давления по времени: а) 1 – в зазоре между передним днищем и зарядом; 2 – на оси симметрии переднего днища; 3 – за утопленной частью сопла; 4 – в зонтичной части; б) изменение давления на оси симметрии переднего днища до выхода на стационарный режим работы

Рис. 3 . Изменение осевой (а) и радиальной (б) скоростей по времени: 1 – в зазоре между передним днищем и зарядом; 2 – на оси симметрии переднего днища; 3 – за утопленной частью сопла; 4 – в зонтичной части

б)

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

Заключение

На основе модифицированного метода расщепления потоков предложен алгоритм и разработана программа, позволяющая проводить сквозное моделирование внутренних нестационарных турбулентных течений сжимаемого вязкого газа в тракте РДТТ. Проведены тестовые численные исследования выхода модельного твердотопливного ракетного двигателя на стационарный режим работы. Представлены результаты мо-

Рис. 5 . Линии тока течения газа

Рис. 6 . Распределение осевой скорости

и, м/с

I 2.45Е+03

I 2.42Е+03

I 2.30Е+03

I 2.00Е+03

■ 1.60Е+03 — 1.20Е+03 — 9.00Е+02

- 6.00Е+02

- 3.00Е+02 О.ООЕ+ОО

р, Па

Рис. 7 . Распределение давления

7.3Е+06 7.2Е+06 7.1Е+06

6.7Е+06 6.1Е+06 З.ЗЕ+06 9.6Е+05 4.0Е+05

2.0Е+05 О.ОЕ+ОО

делирования внутрикамерных процессов, протекающих в тракте двигателя. В результате расчетов определено время выхода двигателя на стационарный режим.