CC BY
30
2020 Математика и механика № 67
УДК: 533.17; 519.62
DOI 10.17223/19988621/67/12
К.В. Костюшин
ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ НЕСТАЦИОНАРНЫХ ГАЗОДИНАМИЧЕСКИХ ПРОЦЕССОВ ПРИ СТАРТЕ ТВЕРДОТОПЛИВНЫХ РАКЕТ1
Представлены результаты разработки методики сквозного расчета нестационарных газодинамических процессов, протекающих при старте ракетоносителей (РН), в газодинамических трактах ракетных двигателей (РД) и внешних областях. Методика учитывает изменение геометрии поверхности горения твердого топлива в процессе работы РДТТ и изменение геометрии расчетной области с учетом динамики старта РН. Проведен анализ нестационарного силового воздействия сверхзвуковой струи на поверхность старта. Показано, что максимум силового воздействия на поверхность старта находится в окрестности дисков Маха невозмущенной струи. Проведены численные исследования газодинамических процессов при старте модельного твердотопливного ускорителя, в том числе с учетом вскрытия сопловой заглушки. Показана величина вклада составляющих силы тяги на этапе выхода на режим. Наличие заглушки на начальном этапе запуска РД приводит к скачкообразному изменению тяги, и незначительным колебаниям, которые демпфируются по мере повышения давления в камере сгорания.
Ключевые слова: твердотопливный ракетный двигатель, нестационарная газовая динамика, сверхзвуковая струя, математическое моделирование, подвижные расчетные сетки, метод Годунова, тяга.
Исследование процессов аэрогазодинамики старта ракетоносителей (РН), и разработка методов расчета входят в число сложных и актуальных научно-технических проблем, возникающих в процессе разработки, и эксплуатации ракетных комплексов (РК) [1]. Ряд проблем обусловлен нестационарным силовым воздействием продуктов сгорания твердых ракетных топлив на элементы конструкции РК. Это воздействие во многом определяется газодинамическими схемами РН, которые характеризуются значительным разнообразием форм, габаритов и функциональных особенностей. С другой стороны, динамика движения и нестационарная пространственная газодинамическая картина течения при старте определяется режимом работы ракетного двигателя (РД). Нестационарные газодинамические процессы, характерные для периода выхода ракетного двигателя на твердом топливе на режим («выброс» сопловой заглушки, разгар заряда твердого топлива) оказывают влияние на динамику движения, нестационарную пространственную газодинамическую картину течения и в общем случае определяют газодинамические нагрузки действующие на элементы конструкции РДТТ и стартовую площадку [2].
Проведение физического моделирования, особенно натурных экспериментов на начальном этапе проектирования, для этих целей связано с высокими материальными затратами. Поэтому один из подходов к решению данной проблемы со-
1 Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 19-3150049/19.
стоит в применении методов математического моделирования. Это позволяет уменьшить расходы на проектирование и сократить сроки отработки РДТТ путем замены значительной части экспериментов математическим моделированием на ЭВМ с использованием современных высокопроизводительных программных комплексов [3].
К настоящему времени применительно к отдельным исследуемым процессам разработан достаточно мощный математический аппарат [4-9], и на его основе созданы современные пакеты прикладных программ, в том числе на основе инженерных методик расчета. Дальнейшее развитие методов математического моделирования газодинамических процессов, возникающих в процессе эксплуатации РН, связано с разработкой комплексных методик, позволяющих в рамках единого газодинамического расчета определять внутрибаллистические и энерготяговые характеристики РДТТ и внешнебаллистические характеристики РН [1].
В данной статье представлены результаты разработки методики сквозного расчета нестационарных газодинамических процессов, протекающих при старте РН, в газодинамических трактах РД и внешних областях, включающих в себя: камеру сгорания с изменяемой в процессе разгара наполнителя геометрией; сопловой блок; сверхзвуковую струю взаимодействующую с обтекаемыми поверхностями РН и поверхностью старта.
Постановка задачи
Рассматривается задача математического моделирования газовой динамики старта РН на нестационарных режимах работы РДТТ с учетом изменения поверхности горения заряда твердого топлива (ТТ) и силового взаимодействия сверхзвуковых струй продуктов сгорания с поверхностью старта и элементами конструкции РН. Работа крупногабаритных РДТТ характеризуется неравномерным выгоранием заряда ТТ по длине канала, что в случае исследования нестационарных процессов приводит к необходимости учета изменения геометрии заряда. РН находится под действием газодинамических сил (сил тяги, аэродинамического сопротивления) и силы тяжести. Истекающая струя газообразных продуктов сгорания взаимодействует с поверхностью старта и элементами конструкции РН.
Схема расчетной области показана на рис. 1 и включает в себя: 1 - внешний корпус РН; 2 - заряд ТТ; 3 - газодинамический тракт РД; 4 - поверхность горения заряда ТТ; 5 - корпус камеры сгорания; 6 - сопловой блок РД; 7 - сопловая заглушка; 8 - ось симметрии; 9,10 - внешние границы расчётной области; 11 - непроницаемая поверхность старта; 12 - внешняя расчетная область.
Рис. 1. Схема расчетной области Fig. 1. Scheme of a computational domain
Предполагается что в начальный момент времени фронт горения моментально распространяется по всей поверхности заряда твердого топлива. Продукты сгорания представляют собой невязкий сжимаемый идеальный газ. Нестационарная система уравнений Эйлера, описывающая пространственные течения идеального газа, имеет вид [10]
¿(дор йО] + ^ • й8 = 0 «
йО] + Я( + /Л)-й8 = 0; (2)
ейО] + Я(е +Р) • й8 = ° (3)
где О - конечная область в трехмерном пространстве; йО = йхйуй2 - элемент объёма; £ - поверхность ограничивающая область О; й8 = пй£ - ориентированный элемент поверхности £ ; где п - внешняя нормаль к £, а й£ - элемент площади; / - время; р - плотность; р - давление; V = [и, V, w]т - вектор скорости движения газа; е - удельная внутренняя энергия, е = ре + р(и2 +v2 )/2 -
полная энергия единицы объема; I = diag [1,1,1] - единичный тензор размерности 3x3; а• Ь - скалярное произведение двух векторов.
Система уравнений (1) - (3) замыкается уравнением состояния идеального газа
"("к' (4) Для численного решения уравнений (1) - (4) используется схема Годунова, адаптированная к произвольному количеству граней в расчетной ячейке. Для построения разностной схемы используется метод конечных объёмов. Вся вычислительная область покрывается дискретными ячейками, состоящими из произвольных выпуклых конечных многогранников с объемами О,, где / = 1,2,... - номер такого объёма, и числом граней т = т (/). Каждая грань имеет площадь £у, где ] = 1,2,...,т(/). Тогда интегральные уравнения в каждом многограннике О, аппроксимируются следующим образом:
к+1 _ к т(1)
о +1 ( (V •)=(5)
¿11 У=1
)к+1 _(DV )к тЦ) тЦ)
О, ^)< А_(Р )г + I: ((V )( • 8;) + I р8; = 0; (6)
А ]=1 У=1
ек+1 _ ек т()
^е^ + 1(( + Р )( •8 у ) = °. (7)
¿м у=1
В результате получим схему Годунова [11] для произвольной расчетной сетки. Здесь 8= п, а А/ - шаг по времени. Нижний целый индекс / обозначает величины функций, отнесенные к центру масс / -го многогранника, а нижний целый индекс ] обозначает величины, отнесенные к центру ]-й грани дискретной ячейки. Верхний целый индекс к обозначает номер шага по времени. Большие величины Я , V , Р , Е обозначают соответственно плотность, скорость, давление и полную энергию на гранях дискретной сеточной ячейки. Эти величины определяются из решения задачи Римана [11] в направлении внешней нормали. Для расчета потоков могут быть использованы и другие методы, основанные на приближенных решениях задачи Риммана: Рое, Ошера, ИЬЬ, ИЬЬС, ИЬЬБ и др.[8]. Схема (5) - (7) легко обобщается на двумерное и осесимметричное приближения.
Для общего случая подвижной пространственной дискретной ячейки используется следующая схема [10]:
Здесь Б^ - скорость движения центра ]-й грани дискретной ячейки. Полуцелые верхние индексы к +1/2 обозначают значения соответствующих величин в момент времени / +1/2 А/. Уравнение (11) есть уравнение изменения объема дискретной ячейки О^. В расчетах использовалась декартовая система координат, связанная с РН.
Расчет шага по времени на подвижной расчетной сетке проводится по следующему алгоритму:
1. Определение нового положения вершин подвижных расчетных ячеек.
2. Определение площадей граней подвижных расчетных ячеек в моменты времени / +1/2 А/ и / + А/.
3. Определение скоростей движения центров граней подвижных ячеек.
4. Определение объемов подвижных ячеек в момент времени / + А/.
5. Определение параметров газа на гранях расчетных ячеек.
6. Определение локальных газодинамических параметров на новом временном слое с использованием расчетной схемы (8) - (11).
Изменение положения вершин расчетных ячеек в рассматриваемой задаче происходит за счет учета изменения поверхности заряда ТТ в процессе горения и движения РН относительно стартовой поверхности под действием газодинамических сил.
Моделирование выгорания заряда твердого топлива
Для моделирования процесса выгорания заряда твердого топлива используется методика [12], основанная на аппроксимации формы канала заряда совокупностью граней элементарных ячеек, расположенных в окрестности поверхности горения. Для определения нового положения поверхности горения и локального газоприхода используется степенной закон горения
Щ = PvV, (12)
где щ - локальная скорость горения в окрестности вершины подвижной грани на текущем шаге по времени; pv - локальное давление в окрестности вершины подвижной грани; u0, V - константы в законе скорости горения.
Новое положение вершин ячейки на поверхности горения вычисляется следующим образом:
xk+1 = xk + п хиу а/ , (13)
где х - координаты вершины подвижной грани ячейки; пх - единичный вектор
внутренней нормали к поверхности горения в точке х .
Расчет локального газоприхода т с грани расчетной ячейки, лежащей на поверхности горения, можно проводить с использованием следующего соотношения:
т = Р/иЛ, (14)
где ие - локальная скорость горения в центре грани расчетной ячейки; р/ - плотность ТТ; 8е - площадь грани расчетной ячейки.
Однако для выполнения закона сохранения массы продуктов сгорания, поступающих с поверхности горения, с учетом изменения объема ячейки, в работе использовалось следующее соотношение:
т = р/АОк / А/, (15)
где АОк - изменение объёма ячейки на текущем шаге по времени.
Алгоритм расчета с учетом перестроения расчетной сетки в окрестности поверхности горения на текущем шаге по времени строится следующим образом:
1. По известному значению давления в центре ячейки определяется локальное давление в окрестности вершин её подвижной грани.
2. С использованием (12) определяется локальная скорость горения в окрестности вершины подвижной грани.
3. Определяется единичный вектор внутренней нормали к поверхности горения в вершине подвижной грани как: пх = (пг + пг) / 2 , где п1 и пг - единичные векторы внешних нормалей граней расчетных ячеек, лежащих на поверхности горения и имеющих общую вершину х .
4. С использованием (13) определяются новые координаты вершин подвижных граней ячеек.
5. Для всех подвижных ячеек пересчитываются геометрические параметры (объём ячейки и площади граней).
6. С использованием (15) определяется локальный газоприход через грани расчетной сетки, лежащие на поверхности горения.
Моделирование движения ракетоносителя
Для определения характеристик движения ракетоносителя используется следующий алгоритм расчета:
1. По известному распределению газодинамических параметров на текущем шаге по времени определяется равнодействующая газодинамических сил, действующих на поверхности РН:
Р = Л рпёБ , (16)
где £ - поверхность внутренних и внешних элементов РН; п - единичный вектор нормали к элементарной площадке ёБ .
2. Масса РН на текущем шаге по времени вычисляется как: тРН = тРН - тк[,
где тРд - масса РН на предыдущем шаге по времени, тк[ - масса топлива сгоревшего на текущем шаге по времени, равная суммарному газоприходу с поверхности горения.
3. По известной равнодействующей газодинамических сил и силы тяжести определяется скорость и ускорение РН на текущем шаге по времени.
С учетом известных параметров движения РН расчетная сетка в окрестности поверхности старта перестраивается с использованием следующего алгоритма:
1. Определяется единичный вектор внутренней нормали к поверхности старта в вершине подвижной грани.
2. Новые координаты вершин подвижных граней ячеек определяются как к+1 к
х = х + пхигосД/, где игос - скорость ракетоносителя на текущем шаге по времени.
3. Проводится проверка на достижение заданного ограничения на геометрические характеристики расчетной ячейки в окрестности поверхности старта, по достижении которого ячейка дробится на две.
4. Пересчитываются геометрические параметры подвижных ячеек (объём ячейки и площади граней), определяются скорости движения граней.
Расчет энерготяговых характеристик
Важным вопросом при исследовании нестационарных режимов работы РД является определение энерготяговых характеристик. В общем случае тяга ракетного двигателя определяется как интеграл сил давления, действующих на внутренние и внешние поверхности РД, за исключением сил внешнего аэродинамического сопротивления [1]:
Р = Л рпёу (17)
Кроме того, согласно [2], тягу можно определить следующим образом:
Р2 = Рр + Рст + Рнест , (18)
где Рр - реактивная составляющая вектора тяги, Рст - статическая составляющая вектора тяги, Рнест - нестационарная составляющая вектора тяги.
Реактивная составляющая вектора тяги определяется как импульс продуктов сгорания, истекающих через срез сопла в единицу времени:
■У
Рр =Ярл*V—¥, (19)
Ра
где Ра - поверхность выходного сечения сопла; Упа - проекция вектора скорости
на элементарную площадку —Р.
В предположении равномерности распределения атмосферного давления по контуру внешней контрольной поверхности РН статическая составляющая вектора тяги определяется как
Рст =-Я (Ра - Рн ), (20)
Ра
где рн - давление невозмущенной окружающей среды.
Нестационарная составляющая вектора тяги определятся как изменение количества движения по свободному объему газодинамического тракта РД:
= | ШР^ ' (21)
где V - внутренний объем РД.
В предположении V = Укам + Vсuл, где Vкам - свободный объем камеры сгорания, а Vспл - свободный объем соплового блока можно записать
Р = Р + Р , (22)
нест нест.кам нест.спл > V '
где Рнест кам = — ру—V - нестационарная составляющая тяги по объему камеры
— V
кам
сгорания, а Рнестспл = — ¡¡¡Ру^ - нестационарная составляющая тяги по объему
Л V
спл
соплового блока.
Достаточно часто при исследовании преимущественно стационарных режимов работы РД в (18) пренебрегают нестационарной составляющей вектора тяги. В таком случае тяга РД может быть определена как
Рз = Рр + Рст . (23)
В работе проведен анализ тяги и её составляющих (17) - (23) на нестационарных режимах работы РДТТ.
Численные исследования
На базе изложенной методики разработан программный комплекс, позволяющий проводить расчет нестационарных течений продуктов сгорания в газодинамических трактах РД и внешних областях. Проведены численные исследования газодинамики старта модельного осесимметричного твердотопливного ускорителя. Длина РН £уск = 20.5 м; длина соплового блока Ьс = 4.3 м; длина сверхзвуковой части сопла Ьсв = 3.5 м; диаметр РН - 3 м; радиус входного сечения сопла - 0.6 м; радиус критического сечения Я = 0.5 м; радиус среза сопла Яа = 1.5 м. Рассмотрен характерный для стартовых ускорителей модельный цилиндрический одноша-шечный заряд ТТ внутреннего горения. Длина заряда твердого топлива £зар = 14.2 м; радиус внутренней поверхности заряда ТТ Явну = 0.6 м; радиус внешней поверхности заряда ТТ Явне = 1.4 м; масса ускорителя (с учетом заряда твердого топлива) -
187000 кг; плотность ТТ ру = 2200 кг/м3; температура продуктов сгорания 3400 К; равновесный показатель адиабаты газообразной смеси двухфазных продуктов сгорания 1.16; константы в законе скорости горения ТТ: и0 = 0.006 м/с, v = 0.3.
Сопловая заглушка располагается на расстоянии 1/3£св от критического сечения сопла. Расстояние от среза сопла до поверхности старта в начальный момент времени составляет 1м. Схема ускорителя показана на рис. 2.
А-А
а
Рис. 2. Схема модельного ускорителя Fig. 2. Schematic view of a booster rocket
В качестве начальных условий в расчетной области 12 (рис. 1) принимались условия невозмущенной окружающей среды у поверхности Земли: Ра = 101325 Па, Та = 300 К. На границах 1, 5, 6 и 11 задавались граничные условия непротекания; на границе 8 - условия симметрии; на гранях ячеек, лежащих на поверхности горения (граница 4), - локальный газоприход; на границе 9 - параметры набегающего потока, соответствующего скорости РН на текущем шаге по времени: Рш = Ра,
Тш = Та, Уш =[-игос,0,0] ; на границе 10 - мягкие граничные условия; до момента вскрытия заглушки на границе 7 ставится условие непротекания.
Расчеты проведены в осесимметричной постановке на блочно-структури-рованной расчетной сетке, построенной с использованием алгоритма, приведенного в работе [13]. Особенностью алгоритма является использование расчетных сеток с произвольным количеством граней расчетных ячеек, что позволяет проводить расчеты в областях с достаточно произвольной геометрией на структурированных и неструктурированных сетках. В начальный момент времени размер расчетной сетки в свободном объёме камеры сгорания составлял 200x10 ячеек, в сопловом блоке - 60x20 ячеек.
Численные исследования проведены для конфигурации РД при наличии сопловой заглушки и без неё. Вскрытие заглушки проводилось при силовом воздействии, эквивалентном равномерному распределению давления по внутренней поверхности заглушки 5, 10 и 15 бар. Предполагалось, что сопловая заглушка представляет собой мембрану, которая моментально вскрывается по достижении условия
Рвск ^ Рзаг =Я ^ , (24)
^заг
где Рвск - усилие вскрытия заглушки; Рзаг - усилие, действующее на заглушку в текущий момент времени; 5заг - внутренняя и внешняя поверхность заглушки.
Время расчета ограничено t = 2.9 с и соответствует выходу РД на номинальный режим работы и удалению РН от поверхности старта на расстояние равное
4£уск. В типовом варианте расчета общее число ячеек в момент времени t = 0 с составляло порядка 1104, а в момент времени t = 2.9 с - 8104. На рис. 3 приведено поле чисел Маха в расчетной области для варианта расчета без сопловой заглушки в моменты времени 1.3 и 2.1 с, где показана структура струи.
-20 -15 -10 -5 0 5 10 15 20 25 30 35 40
Рис. 3. Поле чисел Маха в моменты времени 1.3 (сверху) и 2.1 c (снизу) Fig. 3. Distribution of the Mach numbers at time instants of 1.3 (on the top) and 2.1 s (on the bottom)
Характер изменения максимального давления на поверхности старта по времени показан на рис. 4. Наблюдаются четыре локальных пика давления в моменты времени t = 0.83, 1.3, 1.6 и 1.9 с. Эти пики соответствуют промежуточному расположению поверхности старта в окрестности дисков Маха невозмущенной струи (рис. 3) на расстояниях 6.9, 16.1, 26.0 и 35.6 м от среза сопла. Максимум давления на преграде составляет 14.6 бар в момент времени t = 0.83 с., что соответствует взаимодействию сверхзвуковой струи продуктов сгорания с поверхностью старта в области скачка уплотнения на первом диске Маха. Минимальное силовое воздействие на преграду соответствует взаимодействию в сечениях 2-5 бочек невозмущенной струи на расстояниях 13.0, 22.2, 31.7 и 40.9 м. от среза сопла до поверхности преграды. Данный результата качественно хорошо согласуется с результатами работы [14] для неподвижной расчетной области.
0.83 1.16 1.30 1.51 1.60 1.80 1.90 2.06
Время, с
Рис. 4. Максимальное давление на поверхности старта (1) и расстояние от среза сопла РН до преграды (2) для конфигурации без заглушки Fig. 4. (1) Maximum pressure on the launch surface and (2) the distance from the nozzle exit PH to the obstacle for a configuration without plug
На рис. 5 показана динамика изменения максимального давления на поверхности старта до момента начала отрыва РН. Следует отметить, что наличие заглушки приводит к увеличению максимального давления на поверхности старта на этапе выхода двигателя на режим (от 0 до 0.03 с) с 3.8 бар (без заглушки) до 10 бар (при усилии срабатывания в 15 бар). Это обусловлено взаимодействием сверхзвуковой струи с преградой на расстоянии 1 м от среза сопла после прорыва мембраны. Данный эффект, связанный с вскрытием заглушки, нивелируется на этапе роста давления в камере сгорания и выходом двигателя на основной участок работы (рис. 5).
й ю
0.01
1 2
A
3 -4 \
\
/ M \д
! p\ Г.
j •; V *, /
/ : J
0.02
0.03 Время, с
0.04
0.05
0.06
Рис. 5. Максимальное давление на поверхности старта для конфигураций: 1 - без заглушки, 2 - заглушка 5 бар, 3 - заглушка 10 бар, 4 - заглушка 15 бар Fig. 5. Maximum pressure on the launch surface for configurations: 1, without plug; 2, 5-bar plug; 3, 10-bar plug; and 4, 15-bar plug
6
4
0
Проведен анализ энерготяговых характеристик РД на исследованных нестационарных режимах работы для конфигурации РДТТ без сопловой заглушки. На рис. 6 показано изменение х-компонент вектора тяги и осредненное по объёму камеры давление на временном промежутке от 0 до 2.9 с. Стоит отметить, что тяга,
Время, с
Рис. 6. Изменение тяги (1) и давления в камере сгорания (2), в диапазоне от 0 до 2.9 с Fig. 6. Variation of the (1) thrust and (2) pressure in the combustion chamber in the range from 0 to 2.9 s
вычисленная с использованием соотношения (17), отличается от тяги, вычисленной с использованием (18), не более чем на 0.8 % на всем временном интервале.
Изменения х-компонент векторов тяги на временном промежутке от 0 до 0.15 с, вычисленных с использованием соотношений (18), - кривая 1 и (23) - кривая 2, приведены на рис. 7. Вклад реактивной, статической и нестационарной составляющих в общую тягу (18) показан на рис. 8. Вклад каждой компоненты определялся как: Ркомп = Ркомп/Р2 , где Ркомп - составляющая тяги. На рис. 9 показано изменение нестационарной составляющей тяги (22) (кривая 1), и её компонент вычисленных по объему камеры сгорания (Рнест кам, кривая 2) и по объему соплового блока (Рнест.спл , кривая 3).
К
6 5 4
а3 н 2
1
0 0.02 0.04 0.06 0.08 0.10 0.12 0.14
Время, с
Рис. 7. Изменение тяги РД в диапазоне от 0 до 0.15 с: 1 - тяга вычисленная по (18), 2 - тяга вычисленная по (23) Fig. 7. Variation of the thrust in the range from 0 to 0.15 s: 1, the thrust calculated using (18) and 2, the thrust calculated using (23)
Время, с
Рис. 8. Составляющие тяги: 1 - реактивная, 2 - статическая, 3 - нестационарная Fig. 8. Thrust components: 1, reactive; 2, static; and 3, unsteady
В интервале от 0 до 0.01 с, пока продукты сгорания заполняют свободный объем камеры сгорания, а расход газа через выходное сечение сопла равен нулю,
общая тяга РД определяется исключительно нестационарной составляющей тяги (рис. 8, кривая 3). Этот период характеризуется резким ростом изменения количества движения продуктов сгорания по объему камеры (до 25 % от общей тяги) и сопловому блоку (до 75 %) (рис. 9).
Время, с
Рис. 9. Компоненты нестационарной составляющей тяги: 1 - нестационарная тяга, 2 - изменение количества движения по объему камеры сгорания, 3 - изменение количества движения по объему соплового блока
Fig. 9. Unsteady thrust components: 1, unsteady thrust; 2, unsteady thrust along the volume of the combustion chamber; and 3, unsteady thrust along the volume of the nozzle block
Далее на временном интервале от 0.01 до 0.03 с происходит распространение скачка уплотнения от минимального сечения к срезу сопла, что приводит к формированию сверхзвукового режима истечения продуктов сгорания в сопловом блоке РД (рис. 10) и определяет рост вклада реактивной составляющей тяги (рис. 8) в общую тягу РД. Резкий спад реактивной составляющей тяги на участке от 0.015 до 0.018 с обусловлен локальным торможением дозвукового потока в сверхзвуковой части сопла под действием отраженного от поверхности старта скачка уплотнения (рис. 11). Воздействие отраженной волны приводит к уменьшению интеграла сил давления по поверхности сверхзвуковой части сопла (рис. 12). Вклад нестационарной составляющей тяги на этом участке составляет до 40 %.
Рис. 10. Поля чисел Маха в момент времени 0.01 (сверху) и 0.03 с (снизу) Fig. 10. Distribution of the Mach numbers at time instants of 0.01 (on the top) and 0.03 s (on the bottom)
Рис. 11. Поля чисел Маха в момент времени 0.015 (сверху) и 0.018 с (снизу) Fig. 11. Distribution of the Mach numbers at time instants of 0.015 (on the top) and 0.018 s (on the bottom)
Д
£
1.5
1.0
0.5
-0.5
-1.0
-1.5
— I ...... 2
• - 3
0.02
0.04 0.06
Время, с
0.08
0.10
по всей поверхности по дозвуковой части
Рис. 12. Интеграл сил давления по соплу: 1 сопла, 2 - по сверхзвуковой части сопла, 3 сопла
Fig. 12. Integral of the pressure forces over the nozzle surfaces: 1, over the entire surface of the nozzle; 2, along the supersonic part of the nozzle; and 3, along the subsonic part of the nozzle
0
0
На временном интервале от 0.03 до 0.3 с происходит постепенный выход двигателя на режим и установление локальных газодинамических параметров в камере сгорания и сопловом блоке РД. Вклад нестационарной составляющей тяги уменьшается с 12 до 0.2 %, статической - с 20 до 3 %, при этом рост тяги происходит в основном за счет нарастающего действия сил давления на заднее днище камеры сгорания РД (рис. 13).
Влияние сопловой заглушки на тяговые характеристики модельного РД заметно на временном интервале от 0.01 до 0.03 с. Время срабатывания заглушки для конфигурации с давлением вскрытия 5 бар составляет 0.01 с, 10 бар - 0.015 с, и 15 бар - 0.02 с (рис. 14). Видно, что при наличии сопловой заглушки наблюдается резкий скачок тяги и дальнейший волновой характер выхода тяги на режим, обусловленный прохождением интенсивной ударной волны через сверхзвуковую часть сопла, образовавшейся вследствие вскрытия заглушки и серии отраженных от поверхности старта скачков.
н
0.05
0.1
0.15 Время, с
0.20
0.25
0.30
Рис. 13. Интеграл сил давления: 1 - по камере сгорания, 2 - по соплу, 3 - по внешнему контуру РН Fig. 13. Integral of the pressure forces: 1, along the combustion chamber; 2, along the nozzle; and 3, along the outer contour of the rocket
—
...... 2
.... з
у — 4
/
=3
u H
6
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10
Время, с
Рис. 14. Изменение тяги РД: 1 - без заглушки, 2 - заглушка 5 бар, 3 - заглушка 10 бар, 4 - заглушка 15 бар Рис. 14. Variation of the thrust: 1, without plug; 2, 5-bar plug;
3, 10-bar plug; and 4, 15-bar plug
Заключение
Разработана методика расчета газовой динамики старта РН и энерготяговых характеристик РДТТ на нестационарных режимах работы. Методика учитывает изменение геометрии поверхности горения твердого топлива в процессе работы РДТТ и изменение геометрии расчетной области с учетом динамики старта РН. Проведен анализ нестационарного силового воздействия сверхзвуковой струи на поверхность старта. Показано, что максимум силового воздействия находится в окрестности дисков Маха невозмущенной струи. Проведены численные исследования газодинамических процессов при старте модельного твердотопливного ускорителя, в том числе с учетом вскрытия сопловой заглушки. Показана величина вклада составляющих силы тяги на этапе выхода на режим. Наличие заглушки на
начальном этапе запуска РД приводит к скачкообразному изменению тяги, и незначительным колебаниям, которые демпфируются по мере повышения давления в камере сгорания.
Автор выражает благодарность Емельянову Владиславу Николаевичу и сотрудникам кафедры А9 Балтийского государственного технического университета «ВОЕНМЕХ» им. Д.Ф. Устинова за помощь при подготовке работы и объективную критику.
ЛИТЕРАТУРА
1. Алиев А.В., Амарантов Г.И., Вахрушев А.В. Внутренняя баллистика РДТТ / под ред. А.М. Липанова, Ю.М. Милёхина. М.: Машиностроение, 2007. 504 с.
2. Ерохин Б.Т. Теория внутрикамерных процессов и проектирование РДТТ: Учебник для высших технических учебных заведений. М.: Машиностроение. 1991. 560 с.
3. Волков К.Н., Емельянов В.Н. Вычислительные технологии в задачах механики жидкости и газа. М.: Физматлит, 2012. 468 с.
4. Липанов А.М., Бобрышев В.П., Алиев А.В., Спиридонов Ф.Ф., Лисица В.Д. Численный эксперимент в теории РДТТ. Екатеринбург: Наука, 1994. 300 с.
5. Липанов А.М., Дадикина С.Ю., Шумихин А.А., Королева М.Р., Карпов А.И. Численное моделирование внутрикамерных нестационарных турбулентных течений. Часть 1 // Вестник ЮУрГУ ММП. 2019. Т. 12, № 1. С. 32-43. DOI: 10.14529/mmp190103.
6. Vasenin I.M., Krainov A.Yu., Lipanov A.M., and Shrager E.R. Method for direct numerical simulation of turbulent gas flows in curvilinear coordinates // Computational Mathematics and Mathematical Physics. 2015. V. 55(5). P. 883-890. DOI: 10.1134/S0965542515050176.
7. Vasenin I.M. Shrager E.R. Glazunov A.A. Krainov A.Yu. Krainov D.A. Modeling of the Process of Motion of a Scramjet in the Atmosphere // Russ. Phys. J. 2013. V. 56. P. 908-913. DOI: 10.1007/s11182-013-0117-2.
8. Волков К.Н., Дерюгин Ю.Н., Емельянов В.Н., Козелков А.С., Карпенко А.Г., Тетерина И.В. Разностные схемы в задачах газовой динамики на неструктурированных сетках. М.: Физматлит, 2015. 416 с.
9. Minkov L.L., Shrager E.R. and Kiryushkin A.E. Two Approaches for Simulating the Burning Surface in Gas Dynamics // Key Engineering Materials. 2016. V. 685. P. 114-118. DOI: 10.4028/www.scientific.net/KEM.685.114.
10. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001. 607 c.
11. Годунов С.К., Забродин А.В., Иванов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. М.: Наука, 1976. 400 с.
12. Kostyushin K.V., Eremin I.V., Kagenov A.M., Zhiltsov K.N., Tyryshkin I.M., Sereda M.S., Chervakova A.V. Method of gas flows calculation in solid propellant rocket engines taking into account the combustion of solid fuel charge // IOP Conf. Ser.: Mater. Sci. 2019. 597:012077. DOI: 10.1088/1757-899X/597/1/012077
13. Kostyushin K.V., Kotonogov V.A., Kagenov A.M., Tyryshkin I.M., Glazunov A.A., Eremin I.V., Zhiltsov K.N. The universal algorithm for solving the gas dynamics equations on the mesh with arbitrary number of cell faces // Journal of Physics Conference Series. 2018. 1145(1):012048. DOI: 10.1088/1742-6596/1145/1/012048.
14. Глазунов А.А., Кагенов А.М., Костюшин К.В., Еремин И.В., Котоногов В.А., Алигасано-ва К.Л. Математическое моделирование взаимодействия одиночной сверхзвуковой струи с преградами // Вестник Томского государственного университета. Математика и механика. 2020. № 63. C. 87-101. DOI: 10.17222/19988621/63/8.
Статья поступила 20.07.2020 г.
142
K.B. KocmwMH
Kostyushin K.V. (2020) NUMERICAL INVESTIGATION OF UNSTEADY GASDYNAMIC PROCESSES AT THE LAUNCH OF SOLID-PROPELLANT ROCKETS. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika [Tomsk State University Journal of Mathematics and Mechanics]. 67. pp. 127-143
DOI 10.17223/19988621/67/12
Keywords: solid-propellant rocket engine, unsteady gas dynamics, supersonic jet, mathematical modeling, movable computational grids, Godunov's method, thrust.
The paper presents the results of the methodology developed for calculating unsteady gas-dynamic processes occurring at the launch of missiles, in the gas-dynamic paths of rocket engines, and in the external regions. The method accounts for the variation in the geometry of the solid-propellant charge in the course of solid-propellant rocket engine operation and in the geometry of the computational domain at the rocket launch. The analysis of the unsteady force impact of the supersonic jet on the launch surface is carried out. It is shown that the maximum force action is located in the vicinity of the Mach disks of the unperturbed jet. Numerical studies of gasdynamic processes at the launch of a model solid-propellant booster rocket are implemented including the case when the nozzle plug opening is taken into account. The contribution of the thrust force components at the stage of bootstrap operation is assessed. The presence of the plug at the initial stage of the engine start leads to an abrupt change in the thrust and minor fluctuations, which are damped as the pressure in the combustion chamber rises.
Financial support. The reported study was funded by RFBR according to the research project No. 19-31-50049/19.
Kirill V. KOSTYUSHIN (Research Institute of Applied Mathematics and Mechanics of Tomsk State University, Tomsk, Russian Federation). E-mail: kostushink@niipmm.tsu.ru
REFERENCES
1. Aliev A.V., Amarantov G.N. Vakhrushev A.V. (2007) Vnutrennyaya ballistika RDTT [Internal ballistics of solid-propellant rocket motors]. Ed. by Lipanov A.M., Milekhin Yu.M. Moscow: Mashinostoenie.
2. Erokhin B.T. (1991) Teoriya vnutrikamernykh protsessov iproektirovaniya RDTT [Theory of intrachamber processes and design of solid-propellant rocket motors]. Moscow: Mashinostroenie.
3. Volkov K.N., Emel'yanov V.N. (2012) Vychislitel'nye tekhnologii v zadachakh mekhaniki zhidkosti i gaza [Computational technologies in problems of fluid and gas mechanics]. Moscow: Fizmatlit.
4. Lipanov A.M., Bobryshev V.P., Aliev A.V., Spiridonov F.F., Lisitsa V.D. (1994) Chislennyy eksperiment v teorii RDTT [A numerical experiment in the theory of solid-propellant rocket motors]. Ekaterinburg: UIF Nauka.
5. Lipanov A.M., Dadikina S.Yu., Shumikhin A.A., Koroleva M.R., Karpov A.I. (2019) Numerical simulation intra-chamber of unsteady turbulent flows stimulate. Part 1. Bulletin SUSU MMCS. 12(1). pp. 32-43. DOI: 10.14529/mmp190103.
6. Vasenin I.M., Krainov A.Yu., Lipanov A.M., Shrager E.R. (2015) Method for direct numerical simulation of turbulent gas flows in curvilinear coordinates. Computational Mathematics and Mathematical Physics. 55(5). pp. 883-890. DOI: 10.1134/ S0965542515050176.
7. Vasenin I.M., Shrager E.R., Glazunov A.A., Krainov A.Yu., Krainov D.A. (2013) Modeling of the process of motion of a Scramjet in the atmosphere. Russian Physics Journal. 56. pp. 908-913. DOI: 10.1007/s11182-013-0117-2.
8. Volkov K.N., Deryugin Yu.N., Emel'yanov V.N., Kozelkov A.S., Karpenko A.G., Teterina I.V. (2015) Raznostnye skhemy v zadachakh gazovoy dinamiki na nestrukturirovannykh
setkakh [Difference schemes in problems of gas dynamics on unstructured grids], Moscow: Fizmatlit.
9, Minkov L.L., Shrager E.R., Kiryushkin A.E. (2016) Two approaches for simulating the burning surface in gas dynamics. Key Engineering Materials. 685. pp. 114-118. DOI: 10.4028/www.scientific.net/KEM.685.114.
10. Kulikovskii A.G., Pogorelov N.V., Semenov A.Yu. (2000) Mathematical Aspects of Numerical Solution of Hyperbolic Systems. Boca Raton: CRC Press.
11. Godunov S.K., Zabrodin A.V., Ivanov M.Ya., Krayko A.N., Prokopov G.P. (1976) Chislennoe reshenie mnogomernykh zadach gazovoy dinamiki [Numerical solving of multidimensional problems of gas dynamics]. Moscow: Nauka.
12. Kostyushin K.V., Eremin I.V., Kagenov A.M., Zhiltsov K.N., Tyryshkin I.M., Sereda M.S., Chervakova A.V. (2019) Method of gas flows calculation in solid propellant rocket engines taking into account the combustion of solid fuel charge. IOP Conference Series: Materials Science and Engineering. 597(012077). DOI: 10.1088/1757-899X/597/1/012077.
13. Kostyushin K.V., Kotonogov V.A., Kagenov A.M., Tyryshkin I.M., Glazunov A.A., Eremin I.V., Zhiltsov K.N. (2018) The universal algorithm for solving the gas dynamics equations on the mesh with arbitrary number of cell faces. Journal of Physics Conference Series. 1145(1):012048. DOI: 10.1088/1742-6596/1145/1/012048.
14. Glazunov A.A., Kagenov A.M., Kostyushin K.V., Eremin I.V., Kotonogov V.A., Aligasanova K.L. (2019) Matematicheskoe modelirovanie vzaimodeystviya odinochnoy sverkhzvukovoy strui s pregradami [Mathematical modeling of the interaction of a single supersonic jet with obstacles]. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika - Tomsk State University Journal of Mathematics and Mechanics. 63. pp. 87101. DOI: 10.17222/19988621/63/8.
Received: July 20, 2020
