УДК 532.517.4 Б01: 10.15350/17270529.2020.2.18
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТЕЧЕНИЯ ВЯЗКОГО СЖИМАЕМОГО ГАЗА В РДТТ С ЦЕНТРАЛЬНЫМ ТЕЛОМ
ШУМИХИН А. А., ДАДИКИНА С. Ю.
Удмуртский федеральный исследовательский центр Уральского отделения РАН, 426067, г. Ижевск, ул. Т. Барамзиной, 34
АННОТАЦИЯ. В работе представлена методика моделирования внутреннего течения продуктов сгорания в тракте управляемого твёрдотопливного ракетного двигателя с центральным телом. Приведена система определяющих уравнений, записанная в цилиндрической системе координат, описывающая трехмерный поток сжимаемого вязкого газа. Предложен вычислительный алгоритм, разработанный на основе модифицированной схемы Стигера-Уорминга. Алгоритм пригоден для сквозного расчета всего тракта двигателя, как дозвуковых зон течения, так и сверхзвуковых. Проведены исследования потока газа в РДТТ при различных положениях центрального тела, соответственно при различных значениях площади критического сечения сопла, и при различных величинах характеристик твердого топлива, определяющих его скорость горения. Получены зависимости внутрикамерных параметров потока и расходных (тяговых) характеристик ракетного двигателя от величины площади критического сечения и показателя степени в законе горения твердого топлива.
КЛЮЧЕВЫЕ СЛОВА: внутрикамерное течение газа, центральное тело, цилиндрическая система координат, вычислительная гидрогазодинамика.
ВВЕДЕНИЕ
Для исследования течения продуктов сгорания в твердотопливных ракетных двигателях (РДТТ) важное значение имеет численное моделирование потоков сжимаемого вязкого газа. Динамичное развитие методов математического моделирования процессов протекающих в РДТТ обусловлено трудностями, связанными с постановкой физического эксперимента. Актуальность разработки новых методик для численного моделирования задач гидрогазодинамики также объясняется необходимостью программной реализации новых, более сложных постановок решаемых задач, обусловленных современными практическими потребностями. Алгоритмы, используемые для исследования течения газа в РДТТ, должны учитывать турбулентный характер и существенную нестационарность потока, влияние градиентов давления на внутрикамерные процессы [1 - 5]. Для твердотопливного ракетного двигателя характерны большие разбросы тяговых (расходных) характеристик в зависимости от температуры топливного заряда и разбросов скорости горения твердого топлива. Вследствие этого вопрос управления расходной характеристикой двигателя может возникать при разработке РДТТ практически любого класса и назначения [2 - 4]. Также управление расходной характеристикой двигателя необходимо для реализации заданной циклограммы изменения тяги РДТТ в процессе полета. Управление площадью критического сечения сопла является наиболее простым и достаточно эффективным способом изменения расходных характеристик РДТТ. В данном случае система управления посредством специальных исполнительных элементов (центрального тела) активно воздействует на внутрикамерные процессы, протекающие в двигателе. Для описания трехмерного течения газа в тракте твердотопливного двигателя с центральным телом использовалась система уравнений сохранения, записанных в цилиндрической системе координат [1, 5 - 18]. Численное моделирование внутрикамерных процессов производилось с применением алгоритма, основанного на модифицированной схеме расщепления векторов потоков (схеме Стигера-Уорминга) [19, 20]. Алгоритм позволяет проводить сквозной расчет всего тракта двигателя, как дозвуковых зон течения, так и сверхзвуковых. Расчеты проводились при различных положениях центрального тела на входе в сопло, соответственно при разных площадях критического сечения, и при различных значениях показателя степени в законе горения твердого топлива.
СИСТЕМА ОСНОВНЫХ УРАВНЕНИЙ
Для моделирования трёхмерного течения продуктов сгорания в камере РДТТ использовалась система уравнений гидромеханики в цилиндрической системе координат. В векторном виде система записывается следующим образом [5 - 8, 14]
дО дЕх 1 д¥е
+ -1 ( д
д дх дг г дв
+= - —гЬ, + —гЬ + —Ь + Ь
г I дх
где О — вектор переменных; 1, К, ^, 1в — векторы конвективных потоков; , Ь, Ьв, Ьвв — векторы диффузионных потоков; х, г,в — продольная, радиальная и угловая координаты; I — время. Векторы переменных, конвективных (2) [6, 14] и диффузионных потоков (3) [7, 9, 14] имеют вид
Р ри рм рм
ри ри2 + р рт рим
рм , 1 = рт , Кг = рv2 + р , Кв = р™
рм рим рлм рм1 + р
_Ре _ р^ рvh рwh
иТх + VТхг + мТхв
иТгх + VIгг + мт,
твх твг твв
ПТо, + VI, + мт,
& 0 "
Ь О) О) = 1& — твв
здесь и, V, и — осевая, радиальная и угловая компоненты вектора скорости, р — давление, Тхх, I, твх, тхг,тгг,твг,тхв,тгв,твв — компоненты тензора вязких напряжений, р — плотность, е — удельная энергия, к — удельная энтальпия.
Компоненты тензора вязких напряжений вычисляются так [7 — 18]
о ди 2 ^
Тхх = Л-дх 3
( дv ди Л
т =т = и\\--1--I,
хг гх ^ ^ Р
I дх дг )
о дv 2
Тгг = -дг 3
твв = 2л\\
■ + г дв г) 3
Тхв = Тв = Л
Тгв = тв = Л\\
где л — коэффициент динамической вязкости, 0 = (и; V; м) — вектор скорости. Дивергенция вектора скорости в цилиндрической системе координат определяется формулой
„ ^ ди дv 1 дм V
аыО. = — + — +--+ - .
дх дг г дв г
Полная удельная энергия, полная удельная энтальпия и давление [19] вычислялись по следующим соотношениям
и + V + м и + V + м
е =-+ СТ, h =-+ СТ,
Р = Р(У-1) е 2 . 2 , 2 Л и + V + м
здесь С - удельная теплоемкость газа при постоянном объёме, С - удельная теплоемкость
газа при постоянном давлении, Т - температура, Я - удельная газовая постоянная, у - показатель адиабаты.
Векторы, входящие в систему уравнений (1), формируются с использованием уравнений сохранения, записанных в цилиндрической системе координат. Выражения для инвариантов (дивергенции, градиента, ротора) в цилиндрических координатах имеют следующий вид [21]
дЛ дЛп
дх дг г дв г
^^ф), =
дф дх &
дф дг
{^гайф)в =
/а 1 (дАг дЛ7 (гОА)х = 3 2
дг дв
(гоЛ)г = 1 дЛ1 дЛ3
г дв дх
(гоЛ)в = дЛ2 дЛдх дг
где ф - произвольный скаляр, Л = (Л; Л2; Л3) - произвольный вектор.
Уравнение неразрывности в векторном виде записывается [7, 11 - 14, 16]
+= 0.
С использованием выражения (7), уравнение (10) в цилиндрической системе координат примет следующий вид [5 - 7, 14, 17]
др дри др 1 др р _ ^ д1 дх дг г дв г В векторном виде уравнения сохранения импульса записываются так [7, 11 - 14, 18]
Здесь[12, 18]
р-— + р(—• V)— = -gгad р + т .
т = 2цф--^кну— i,
где Ф - тензор скоростей деформаций, I - единичный тензор.
Запишем уравнение (12) через инварианты для определения его проекций на оси цилиндрической системы координат. Известно [12], что
(— • V)— = го— х — + gгad Вектор Быт с учетом (13) примет вид
в1ут = цв1у{2ф)- 2 fлgгad (d/vQ).
В (15) второе слагаемое инвариантно относительно системы координат. Вектор 01у(2Ф) с использованием инвариантов описывается выражением
ПУ(2Ф) = div(2Фx ) + div(lФy )) + div(2Фz )к,
где Фх = (фх ;фху Ф ), Фу = (фух ;фуу Ф ), Ф2 = (Фх; Фу; Фгг) - векторы, составленные из
компонентов тензора скоростей деформаций Ф; г, ], к - орты осей декартовой системы координат. Компоненты тензора скоростей определяются
(2Фх)^2и + — divQ., div(2Ф 2у + — divQ, div(2Фг) = V2и + — divQ.. (17)
дх у ду дг
Тогда (16) с учетом (17) примет вид
БУ(2Ф) = V2— + gгad(div—. (18)
Используя формулу V2 — = gгad(Оу—)- го^го—) [21], запишем выражение (18)
следующим образом
БУ(2Ф) = 2gгad(diу—) - т(го—).
С учетом (19), уравнение (15) преобразуется к виду
Б^т = 4 |лgгad(ё^—)- цго1{го1—).
С использованием выражений (14) и (20), уравнение движения (12) примет форму инвариантную относительно системы координат
рд—+р
го— х—+ gгad
V 2 УУ
= -gгadp + 4 lлgгad(^^^-^го^о—).
Проецируя Б^т на оси цилиндрической системы координат с учетом (4), (5), (9), получим
(Р™т)х = 1
г V дх
дг дв
(Р^т)г = 1
г V дх
дг дв
(В^тХ= 1 (дгт^ дт дт,
+ Т
г V дх дг дв
Используя (8) - (9), проекции левой части уравнения (21) преобразуются к уравнениям
/ ^ и ди ди ду ди
\\гоИ1хИ)х =---и--V--+ V— ,
г дв дх дх дг
ду ди ди и2 и ду
\\гоИ1хИ)г = и--и--и----1---,
дх дг дг г г дв&
ди уи у ду и ди ди \\гоИ1хИ)в = V--1-----—---— + и
gгad gгad gгad
V 2 У х ^
V 2 У г ^
V 2 Ув
дг ди
г г дв г дв
ду ди
= и--+ у--+ и — ,
дх дх дх
ди ду ди
= и--+ у--+ и — ,
дг дг дг
и ди у ду и ди = г дв г дв г дв
Применяя выражения (22), (23), (24) и (8), получим проекции векторного уравнения движения (21) на оси цилиндрической системы координат Охтв:
(ди ди ди м ди Л др 1 (дгт дгт дтГЙ р\\ — + и — + V— +--1 = ——+ -|-— +-— +-—
V дг дх дг г дв) дх г V дх дг дв
(ду дv д\\ м дv м2 Л др 1 (дгт дгт дт„Й
■ — + и — + V— +----- 1 1 гх 1 гг 1 гв
дг дх дг г дв г
дг г I дх
дг дв
дм дм дм м дм 1 др 1 (дгт дгт^ дт(
р\\ — + и — + V— +--+ — 1 =---—+ -|
V дг дх дг г дв г ) г дв г&
+ - + вв
дг дв
(26) (27)
Выражения (25) — (27) сложим с уравнением неразрывности (11), предварительно умножая его на и, V, и соответственно. Таким образом, в цилиндрической системе координат уравнения движения в дивергентной форме окончательно примут вид
дри д(ри2 + р) дриг 1 дрим риг _ 1 (дгт^ дгтхг дт дг дх дг г дв г г I дх дг дв
др ^ дриг ^ д(р2 + р) 1 дрм
дг г дв г
■ + —-+ ^-^ + —:-+ ^-(v2 — м2)
дг дх - &
■ гх I • ~ гг I " гв _
■ +—---1——— т
г I дх дг
дрм + дрим | рм +1 д(р2 + р) | 2рм = 1 (дгтв , дгтвг , дтвв , г
дг дх дг г дв г г V дх дг дв в
Уравнение энергии в векторном виде записывается [11, 12, 14]
■+агг(ре+р)о = а^пт) .
Используя (7), имеем
сИу(ре + рр=д(р + р)и | д(р + р^ +1 д(р + Р)м | (р + РУ
дх дг г дв г
а^ре+р )о =
дрuh дрvh 1 дрwh рvh
■ + —-+ + дх дг г дв г Для правой части (31) с учетом (13) запишем
а™(р.т) = аг^ 2лоф—2 ¡о- а^О^,
где ОФ в формуле (34) может быть выражено через инварианты [12]
ОФ = gгad
—1 Ох гогО. 2
Тогда
а^От)=а™
— ¡О х гогО — — ¡¡О - О/гО
ду дп м дп дм
(^х гоИ1)х = V--V-----ъ м— ,
дх дг г дО дх
/_ _ч дм м2 м ду ду дп
(Ох гоИ1)г = м--1------п--ъ п — ,
дг г г дО дх дг
п ды дм дм ум м ду
(Ох гохО)в =
■-п--удх дг г г дО
С учетом (24), (37), компоненты вектора 2^*гай
-^ОхгоО- — уО1-(35)
на оси цилиндрической системы координат Охгв описываются формулами
х: т + т + мТхо,
Г. Шгх + утгг + мТО ,
О : Шох + уто + мТоо .
Тогда для у(О.т) с учетом (7), (36), (38) получим [14]:
сИу(Пт)= 1 ( дг(пТх + уТхг + мТхо) I дг(пТгх + уТгг + мТгО) , д(пТОх + уТОг + ™ ТОо)Л (39)
г V дх дг дО )
С использованием (33), (39) уравнение энергии в цилиндрической системе координат (31) окончательно примет вид [14]:
дре дрпк друк 1 дрмк рук
дХ дх
дг г дО
дг (птхх + утхг + Т ) , дг (птгх + утгг + Т ) , д (пТОх + уТОг + Т )
Таким образом, с использованием выражений (11), (28) - (30), (40), формируются векторы (2) - (3), входящие в систему уравнений (1).
ВЫЧИСЛИТЕЛЬНЫЙ АЛГОРИТМ
В разработанном вычислительном алгоритме использовалась схема Стигера-Уорминга (метод расщепления векторов потоков), предложенная в работах [19, 20]. Вектор ^ можно представить следующим образом
рх = Ах<2 , (41)
где А = —~ - матрица Якоби. Собственные числа матрицы Ах равны [19]
Лх1 = Лх 2 = Лх3 = 11 ■
Лх4 = п + с,
Лх5 = п - С >
где с = (ур / р)12 - скорость звука.
С использованием (41, 42) вектор ^ может быть представлен
Рх =р 2(У - + Лх^ + ^ 2У
Основой схемы Стигера-Уорминга является расщепление вектора Кх на две части
К = К! х К-,
где вектор Кх+ соотносится с вектором, сформированным из положительных собственных чисел матрицы А. Соответственно вектор Кж~ соотносится с вектором, сформированным из отрицательного собственных чисел. Аналогично можно расщепить собственные числа матрицы Якоби
Соответственно, используя (43) - (45). можно записать
Для дискретизации системы основных уравнений (1) по пространству применялся метод конечных объёмов. По времени дискретизация производилось с использованием одношаговой явной двухслойной схемы. В используемом алгоритме параметры потока газа на границах объёмов задавались согласно рис. 1.
! Й-3/2 = Э-1 Эх1/2 = Э-1 & Э-1/2 = Э Эх1/2 = Э & Э-1/2 Й&х1 Эгх 3/2 = Эх 1 1
; к-з/2 (Э-1) Кь2 (э-1) ! К-1/2 (э) к х1/2 (Э ) ; К-1/2 (йх 1) Кх3/2 (йх1)
Рис. 1. Параметры потока газа 0-1/2, Эг-1/2 и QiХl/2, QiХl/2 на гранях ^ -го контрольного объёма
Для расчёта векторов потоков К-^Й-ш) и Fiх^Й^/г), соответственно на левой и правой гранях Ц -го контрольного объёма, необходимо определить векторы параметров и Эх1/2. В используемом алгоритме векторы параметров и на границах Ц -го
объёма принимались равными векторам параметров в центре объема Э . Таким образом значения векторов на гранях Ц -го контрольного объёма будут равны Эгх1/2 = . = Э
и 0+1/2= <, 0-1/2= <1. Собственные числа на гранях также определяются с использованием значений векторов < в центре Ц-го контрольного объёма Л+^2(<^1), Л-_1/2(<) и Л++1/2(<), Л-+1/2(<г+1). Тогда, учитывая (44, 45), собственные числа и потоки на гранях I -й ячейки определяются следующим образом
Ъ-1/2 ==2 (<-1Л+1/2 (<-1 )) + +2 < , Л--1/2 (0 )) , Л-1/2 (0,-1 ) + |Л-1/2 (0,-1 )|
Л+1/2 (<-1 ) = ■
Л-1/2 (0 ) =
л 1/2 (о, )-\\л 1/2 (о, )|
Ъ +1/2 = Ъ+1/2 (<, & ЛЪ+1/2 (<, )) + Ъ+1/2 (<,+1& К+1/2 (<,+1 )) & (49)
Л1/2 () + |Л 1/2 (<2. )|
Л++1/2 (0 ) = Л1/2 (<+1 ) =
л+1/2 (а+1) |л+1/2 (<,+1 )|
Данная процедура выполняется для всех ячеек расчётной сетки. В результате, с использованием выражений (44) - ( 47), будут определены значения потоков , для всей расчетной области. Аналогичные операции производятся для векторов и . Значения вектора Ъвв также вычисляются с использованием значений вектора параметров в центрах контрольных объемов. Диффузионные потоки Д, Д, Д, аппроксимировались методом центральных разностей, с использованием выражений (3) - (5). Алгоритм также описан в работе [5]. Дискретизация по времени выполнялась с использованием явной схемы 1 -го порядка точности (схема Эйлера, одношаговая явная двухслойная схема). Конвективные и диффузионные потоки на текущем временном слое , , , Ъдв и Ьх, Д, Ь", Ь"вв рассчитывались с использованием значений вектора переменных также на текущем временном слое <". Далее с использованием вычисленных потоков определялись элементы вектора переменных на следующем временном слое <"+1. В итоге, по известным значениям <"+1 соответственно рассчитывались все параметры газа на " +1 временном слое.
Турбулентность моделировалась методом крупных вихрей, являющимся одним из наиболее перспективных методов расчётов турбулентных течений. Данный метод в частности применялся в работах [22, 23].
РЕЗУЛЬТАТЫ РАСЧЕТОВ
Длина двигателя составляла Ь = 2.2 м. В начальный момент времени температура газа в камере задавалась Т0 = 294К и давление р0 = 10бПа. Количество ячеек в области составляло 14700. Шаг по времени задавался равным АХ = 1.0-10 7е. Скорость горения топлива определялась по формуле с = с0(р/105у, где с0 - скорость горения при
атмосферном давлении, У - константа. Угловая скорость полагалось м = 0, т.е. задача решалась в осесимметричной постановке. Расчетная сетка представлена на рис. 2.
Расчёты производились при различных положениях центрального тела.
Рис. 2. Расчетная сетка
На рис. 3 приведены распределения давления в тракте двигателя при показателе степени в законе горения топлива 3 = 0.7 и трёх значениях площади критического сечения (минимальном £„„ = 0.0181 м2, промежуточном £„„ = 0.02475 м2
максимальном
$кр = 0.0314 м2).
р, Па /• 1 ЛЛ!-! | л X"
— 6.170Е+06
— 3.500Е+06
— 3.000Е+06
— 2.500Е+06
— 1.600Е+06
— 1.000Е+06
■ 1.000Е+05
р, Па
- 1.954Е+06
- 1.930Е+06
— 1.600Е+06
- 1.400Е+06
- 1.000Е+06
- 4.000Е+05
— 1.000Е+05
р, Па
=1 1.012Е+06
_ 1.011Е+06
— 7.000Е+05
э.шиь+оз
— 3.000Е+05
— 2.000Е+05
Рис. 3. Распределение давления в двигателе при 3 = 0.7 : а) - минимальная площадь £ ;
б) - промежуточная площадь Б ; в) - максимальная площадь Б
На рис. 4, рис. 5 показаны соответственно распределения осевой скорости и плотности при минимальной площади критического сечения. Полученные зависимости давления в камере и тяги двигателя от площади критического сечения при различных значениях 3 представлены на рис. 6.
и, м/с
— 1800
Рис. 4. Распределение осевой скорости в двигателе при минимальной площади £ и 3 = 0.7
р, кг/м3
-щ 4.500
■ 1.166 0.500
Рис. 5. Распределение плотности в двигателе при минимальной площади £ и 3 = 0.7
р, Па 6.00Е+06
Рис. 6. График зависимости (1 — максимальная £ , 2 — промежуточная £ , 3 — минимальная £ ):
кр кр кр
а) — давления от показателя степени в законе горения; б) — тяги от показателя степени в законе горения
ЗАКЛЮЧЕНИЕ
Предложен алгоритм и разработана программа для моделирования течений сжимаемого вязкого газа в тракте твёрдотопливного ракетного двигателя с центральным телом.
Проведены тестовые расчёты двигателя при различных положениях центрального тела и при различных величинах показателя степени в законе горения топлива S .
Представлены распределения параметров потока продуктов сгорания в РД^, полученные в результате моделирования.
Определены зависимости давления в камере и тяги двигателя от площади критического сечения при различных значениях S . Полученная в результате расчетов глубина регулирования тяги менялась от 2.52:1 при S = G.5 до 5.G3 : 1 при S = G.?.
СПИСОК ЛИТЕРАТУРЫ
Numerical simulation of a compressible viscous gas flow in solid-fuel rocket engine with a central body
Shumikhin A. A., Dadikina S. Yu.
Udmurt Federal Research Center, Ural Branch of the Russian Academy of Science, Izhevsk, Russia
SUMMARY. The paper presents a method for modeling the internal flow of combustion products in the path of a controlled solid-fuel rocket engine with a central body. A system of defining equations written in a cylindrical coordinate system describing the three-dimensional flow of a compressible gas is given. A computational algorithm based on a modified Steger-Warming scheme is proposed and implemented. The algorithm is suitable for end-to-end calculation of the entire engine path, both subsonic and supersonic flow zones. Numerical studies of the gas flow in the solid-fuel rocket engine at different positions of the central body, respectively, at different values of the critical area of the nozzle cross section, and at different values of the characteristics of solid fuel that determine its burning rate. The dependences of the gas flow parameters inside the chamber and the flow (traction) characteristics of the rocket engine on the size of the critical cross-section area and the degree indicator in the law combustion of solid fuel are obtained.
REFERENCES
Шумихин Андрей Александрович, кандидат физико-математических наук, научный сотрудник, Институт механики УдмФИЦ УрО РАН, тел. (3412)20-34-76, e-mail: shumihin@udman.ru
Дадикина Светлана Юрьевна, младший научный сотрудник, Институт механики УдмФИЦ УрО РАН, e-mail: sveta-dadikina@yyandex. ru