Исследование процессов истечения жидкости из малорасходных центробежных форсунок жидкостного ракетного двигателя малой тяги Текст научной статьи по специальности «Механика и машиностроение»

УДК 621.454.2 В01: 10.18287/2541-7533-2016-15-4-150-161

ИССЛЕДОВАНИЕ ПРОЦЕССОВ ИСТЕЧЕНИЯ ЖИДКОСТИ ИЗ МАЛОРАСХОДНЫХ ЦЕНТРОБЕЖНЫХ ФОРСУНОК ЖИДКОСТНОГО РАКЕТНОГО ДВИГАТЕЛЯ МАЛОЙ ТЯГИ

© 2016

Е. В. Сёмкин начальник группы конструкторского отдела, Научно-исследовательский институт машиностроения, г. Нижняя Салда, sewa20-68@rambler.ru, mail@niimashspace.ru

Обосновывается возможность применения результатов гидравлических испытаний для прогнозирования коэффициента расходного комплекса камеры сгорания двигателя тягой 13,34 Н. Приводятся результаты моделирования течения жидкости (воды) в воздушную среду через гидравлический тракт одиночных форсунок и смесительного элемента, выполненного в виде установленных соосно двух центробежных форсунок. Моделирование проводится с помощью системы осреднённых по Рейнольдсу уравнений Навье-Стокса, описывающих турбулентное течение двухфазной несжимаемой жидкости. В моделировании используются односкоростная и двухскоростная модели течения. Моделирование турбулентности выполняется по модели изотропной турбулентности BSL на специально сконструированной расчётной сетке, имеющей блочную структуру. Расчётная сетка состоит, в основном, из НЕХА элементов, ориентация которых совпадает с направлением течения в гидравлическом тракте форсунок, что повышает точность разрешения течения в пограничном слое у твёрдой стенки. Результаты расчётов сопоставлены с результатами гидравлических испытаний. Показано, что применение односкоростной модели течения в расчётной области одиночной форсунки и двухскоростной модели в расчётной области смесительного элемента даёт возможность получить результаты, удовлетворительно совпадающие с результатами измерений в гидравлических испытаниях.

Жидкостный ракетный двигатель малой тяги; центробежная форсунка; жидкофазное смешение; моделирование течения.

Цитирование: Сёмкин Е.В. Исследование процессов истечения жидкости из малорасходных центробежных форсунок жидкостного ракетного двигателя малой тяги // Вестник Самарского университета. Аэрокосмическая техника, технологии и машиностроение. 2016. Т. 15, № 4. С. 150-161. DOI: 10.18287/2541-7533-2016-15-4-150-161

Введение

В жидкостных ракетных двигателях малой тяги (ЖРДМТ), работающих на компонентах топлива АТ-НДМГ, в нашей стране наиболее часто в качестве смесительных элементов используются установленные соосно центробежные форсунки.

Применение самовоспламеняющихся компонентов топлива АТИН+НДМГ в этих двигателях накладывает ряд особенностей на организацию рабочего процесса преобразования топлива в продукты сгорания, включая процесс перемешивания и распыла. Важное значение в преобразовании компонентов топлива имеет перемешивание за время жидкофазной индукции.

По оценке, проведённой в [1], жидкофазное смешение в плёнках компонентов топлива суммарной толщиной производится по поверхности в 3-6 раз меньшей по сравнению с перемешиванием в виде капель средним диаметром (0,5 - 1,0) , что позволяет существенно упростить управление указанным процессом. Поверхность плёнок, по которой происходит взаимодействие, в случае центробежных форсунок имеет размеры (0,1...0,2)х(0,5...3) мм.

На рис. 1 представлена схема рабочего процесса в камере сгорания ЖРДМТ, предложенная в [1] и построенная для условия полного жидкофазного смешения.

Рис. 1. Схема рабочего процесса в камере сгорания ЖРДМТ: 1 - период жидкофазного смешения; 2 - период выделения и увеличения концентрации ГФПП (газообразная фаза продуктов преобразования) в камере; 3 - период газофазной реакции, появление очага пламени; 4 - подвод тепла от продуктов реакций подогрева ЖФПП (жидкофазная фаза продуктов преобразования); 5 - разложение ЖФПП, период накопления активных центров реакций; 6 - газофазная реакция продуктов разложения ЖФПП; 7 - окончательное разложение ЖФПП

Процесс преобразования компонентов топлива идёт в три стадии. Начальная стадия протекает при температуре 27°С за время г1 «10— с. Вторая стадия - за время г2 «10 2 с при температуре от 100 до 200°С. Третья стадия становится заметной при температуре от 250 до 420°С.

Анализ причин низкой экономичности процессов преобразования компонентов топлива в продукты сгорания в камерах сгорания ЖРДМТ, проведённый рядом авторов [2-4], показывает, что одной из основных причин является несовершенство процессов перемешивания в период жидкофазного перемешивания.

В частности, в [2] расчёт потери удельного импульса тяги ЖРДМТ показывает, что несовершенство процессов перемешивания компонентов топлива в жидкой фазе приводит к неравномерному распределению соотношения компонентов топлива в поперечном сечении камеры сгорания. Из [3] следует, что потери расходного комплекса, обусловленные низкой эффективностью смешения компонентов топлива, составляют 25-32%. Эпюра соотношения компонентов в основном формируется у форсунок и слабо меняется вплоть до среза сопла [1; 4].

В [1] на основе сопоставления результатов испытаний на натурных компонентах топлива АТИН+НДМГ с результатами гидравлических испытаний составлена зависимость расходного комплекса камеры сгорания ЖРДМТ от гидравлических характеристик центробежно-струйной форсунки. Его максимальное значение, полученное при огневых испытаниях, соответствует определённому значению соотношения кинетиче-

ских энергий компонентов топлива

т V2 '"о* о

туГ

= 1,1 -1,2 суммарному углу конуса распыла

2аЕ = 108 -110 °, определённому в процессе гидравлических испытаний.

Учитывая большое влияние эффективности процессов перемешивания компонентов за время жидкофазной индукции на расходный комплекс камеры сгорания и выявленную зависимость расходного комплекса от параметров форсунок в гидравлических

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

£

К пер =~м, (1)

где Кпер - коэффициент перемешивания; 5см - толщина слоя перемешанной жидкости

(воды) в рассматриваемой точке при гидравлическом испытании; 5Е - суммарная толщина пелён жидкости (воды) в рассматриваемой точке при гидравлическом испытании.

Параметры 5см и 5Е зависят от геометрических размеров гидравлического тракта

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

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

В процессе моделирования рассматривается течение воды от входа в тангенциальные каналы подвода к форсункам до выхода из сопла форсунок и от среза сопла форсунки до места разрушения конуса распыла на капли.

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

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

Представляя по аналогии с [5] двухфазное течение как течение непрерывных и взаимодействующих между собой частей двухфазного потока, состоящих из жидкости (воды) и газа (воздуха), в процессе решения будем отслеживать положение границы раздела между жидкой и газообразной фазами течения, учитывать взаимодействие фаз и распад водяной плёнки на капли за пределами гидравлического тракта центробежной форсунки под действием сил тяжести, поверхностного натяжения и инерционных сил потока двухфазной жидкости в свободном объёме, прилегающем к форсунке. Расчётные области с геометрическими размерами форсунок приведены на рис. 2 - 4.

Расчетная область

Рис. 2. Геометрия расчётной области форсунки «Г»

Рис. 3. Геометрия расчётной области форсунки «О»

Рис. 4. Геометрия расчётной области смесительного элемента

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

Для случая, когда скорость и давление воды и воздуха считаются общими в расчётной области (односкоростная модель), система уравнений имеет вид:

Л

- + ) = 0, (2)

А dt'

2

P = Z raPa , (4)

:(pu)+v.(pu ®u)-v.(^(vu+(vu )т )) = SM _buoy + Fap-Vp, (3)

a=1 2

^ = Z Wa , (5)

a=1

Z Га= 1, (6)

pa= p для всех a = 1 или 2. (7)

Здесь a - составляющие части двухфазного потока: вода (1) и воздух (2); ra -объёмная доля контрольного объёма расчётной области, занятая фазой a ; SM buoy - источник импульса силы тяжести в расчётной области:

SM _buoy =(P-Pref )g , (8)

где p - текущая плотность в расчётной области; pref - плотность жидкости в расчётной области в начале расчёта, задаваемая пользователем;

Fa/3 - сила поверхностного натяжения, действующая на границе раздела фаз:

Faf3 = fap^aj3 , (9)

fap=-^iapnap+V , (10)

5aP=\VraP\, (11)

где а - сила поверхностного натяжения воды при контакте с воздухом на границе раздела фаз; nap - вектор нормали на поверхности контакта фаз, направленный от первичной жидкости к вторичной; Vs - оператор градиента по поверхности контакта фаз; kap - поверхностное искривление поверхности раздела фаз:

kap = V.nap. (12)

При проведении решения принимается, что обмен массами между фазами протекает равновесно.

Для случая, когда у воды и воздуха скорости в каждой рассматриваемой точке расчётной области не равны, а давление одно и то же, система уравнений имеет вид:

Я

~{rapa)+v*{rapaua) = 0; (13)

Gt

д_ Gt

(14)

Gt(DPaUa) + V*( (paUa ®Ua)) = SM buoy + Faß -raVpa + • (raJua (vua + (VUa))) +Ma;

P = t ГаРа ; (15)

a=1

2

J = t raJa ; (16)

a=1

t d = 1; (17)

pa= p для всех a = 1 или 2. (18)

Дополнительно к уже рассмотренным компонентам в уравнении (14) появился дополнительный член Ma, который описывает передачу импульса в процессе взаимодействия фаз потока на границе их раздела. При расчёте Ma представляется как

Ma= MDß . (19)

Здесь MDß - учитывает силу сопротивления при взаимодействии двух фаз:

MDß= CDPaßAaß \Uß-Ua\Uß-Ua), (20)

где CD - коэффициент сопротивления принят в расчёте равным 0,47; pDß - плотность на границе раздела фаз, определяемая из выражения

Paß= raPa+ rßPß , (21)

Aaß=|Vra| . (22)

В системах уравнений (2)-(12) и (13)-(22) вязкость представляется по уравнению

J = Jk + JT , (23)

где турбулентная вязкость определяется по изотропной модели турбулентности BSL. В расчёте турбулентность считалась гомогенной.

Расчётная область делится на части центральными секторами с осью, совпадающей с осью форсунки, составляющими одну четвёртую часть всей моделируемой области, показанной на рис. 2 - 4. В качестве упрощения принимается, что течение рабочего тела внутри одного такого центрального сектора периодически повторяет течение внутри остальных секторов.

Сектор расчётной области показан на рис. 5. На этом же рисунке даны типы граничных условий на отдельных поверхностях расчётных областей.

В табл. 1 суммируются граничные условия для каждой выделенной поверхности расчётной области.

Поверхность быходной гртицы

_ Поериность переодических /

и граничит ¡¡слабни

Б у&Еличено

Поберхность периодических граничных ус/юбий

Поберхность периодических | граничных услодий

В увеличено

Поверхность бходнай границы

Поберхность тбердои стенки

^ Поверхность быходнои границы

ТПоберхность периодических граничных ус/юбий

Поверхность бходной гршицы

Рис. 5. Расчётная область с граничными условиями

Таблица 1. Граничные условия расчётной области

Тип выделенной поверхности Граничное условие

Поверхность входной границы Избыточное полное давление на входе в тангенциальные каналы Р .кс вх' 2 см соответствует давлению на входе в форсунки при гидравлических испытаниях. Объёмная доля воды на входе в расчётную область составляет 1,0; объёмная доля воздуха - 0

Поверхность твёрдой стенки Гладкая или шероховатая твёрдая стенка с условием непротекания на ней. Угол смачивания стенки водой равен 40°. Заданная в расчётах шероховатость стенки составляет 0,4 мкм

Периодические граничные условия Условия непрерывности течения через границу периодических граничных условий: для ск°р°сти УШутр = УШр , для давления Рвнутр = Рнар , для объёмной доли ^";нутр = ЯСр, , Я>£утр = Рщр,

Выходная граница Избыточное давление Рых = 0,0^2. см Объёмная доля воздуха на выходной границе расчётной области (извне в расчётную область) - 1,0. Объёмная доля воды на выходной границе расчётной области (извне в расчётную область) - 0

Гидравлический тракт форсунки и присоединённый свободный объём на начало расчёта Давление окружающей среды Рокр = 1.0^. см Температура окружающей среды и воды поступающей в расчётную область неизменна на протяжении всего расчёта и равна Т = 20°С

Расчётная сетка

Для проведения расчёта была сконструирована специальная расчётная сетка, которая строилась с помощью программного средства ANSYS CFD. Результат приведён на рис. 6.

Сетка имеет блочную структуру и состоит, в основном, из НЕХА элементов. Структура блоков позволяет проводить сгущение узлов там, где находится конус распыла в гидравлическом тракте форсунок и свободном объёме, и отслеживать место положения поверхности раздела фаз.

Структура расчётной сетки показана на рис. 7.

Её элементы в расчётной области сориентированы, как правило, вдоль линий тока воды, что повышает точность расчёта и уменьшает количество узлов, необходимое для точного разрешения рассчитываемых величин.

Ориентация узлов расчётной сетки вдоль линий тока обеспечивается выбором угла а корневого угла конуса распыла из результатов предварительных гидравлических испытаний или расчётов на сетке без ориентации расчётных узлов. В расчётах количество узлов сетки изменялось от 5000000 до 15000000.

Результаты расчёта

На рис. 8 - 10 приведены факелы распыла объёмной фракции воды в ходе моделирования течений в расчётных областях (рис. 2 - 4). При этом использовалась одно-скоростная модель описания течения в форсунках «Г» и «О» (рис. 2, 3) и двухскорост-ная модель описания течения в смесительном элементе (рис. 4). Для сравнения на рис. 8 - 10 показаны результаты фотографирования факелов распыла этих форсунок в процессе гидравлических испытаний.

В табл. 2 приведены секундные массовые расходы воды через входные границы расчётной области и результаты их измерений при гидравлических испытаниях.

На рис. 11 - 13 показаны результаты моделирования течения в гидравлических трактах форсунок смесительного элемента (рис. 2 - 4).

На рис. 11 - 13 также обозначены размеры радиусов выходных сечений сопел форсунок и толщин пелён воды в выходных сечениях сопел форсунок и смесителя. Расчёт толщины пелены воды в выходных сечениях сопел форсунок «Г» и «О» по методике [7] даёт величину 0,038 мм.

Таблица 2. Секундные массовые расходы

Форсунка «Г» (рис. 2) Форсунка «О» (рис. 3) Смеситель

Расчёт 2,0 г/с 2,63 г/с Полость «О» - 2,46 г/с Полость «Г» - 1,91 г/с

Эксперимент 1,85 г/с 2,54 г/с Не измерялось

Рис. 8. Сравнение моделирования факела распыла форсунки «Г» в виде изоповерхности для объёмной доли воды 0,007 (рис. 2) с результатами фотографирования при гидравлических испытаниях

Рис. 11. Результаты моделирования течения воды в гидравлическом тракте форсунки «О»

эксперимент

Рис. 9. Сравнение моделирования факела распыла форсунки «О» в виде изоповерхности для объёмной доли воды 0,007 (рис. 3) с результатами фотографирования при гидравлических испытаниях

Рис. 10. Сравнение моделирования факела распыла форсунки в виде изоповерхности для объёмной доли воды 0,2 с результатами фотографирования при гидравлических испытаниях

Рис. 12. Результаты моделирования течения воды в гидравлическом тракте форсунки «Г»

Рис. 13. Результаты моделирования течения воды в гидравлическом тракте смесителя

Анализ результатов

Результаты моделирования свободного факела распыла форсунок, показанные на рис. 8, 9, свидетельствуют о том, что с помощью односкоростной модели описания течений двухфазной среды в форсунках «О» и «Г» можно получать параметры распыла моделируемых форсунок с удовлетворительной точностью. Результаты расчёта корневого угла распыла этих форсунок и формы свободного факела распыла хорошо согласуются с результатами гидравлических испытаний. Массовые расходы через форсунки в расчёте, по сравнению с результатами гидравлических испытаний, отличаются на 8% для форсунки «Г» и на 3,5% для форсунки «О». В расчёте с удовлетворительной точностью определяется место распада пелены факела распыла и перехода к течению воды в виде капель в воздухе. Расчёт течения в выходном сечении сопла форсунки позволяет определить толщину пелены воды с учётом гидравлических потерь. Расчётные толщины пелены воды в 1,5...1,8 раза больше тех, что были получены по теории течения в идеальной форсунке Г.Н. Абрамовича [7], что свидетельствует о большом влиянии гидравлических потерь в форсунке на её параметры. Вместе с тем, размер капель после перехода от плёночного течения к капельному в форсунках моделируется с большой погрешностью. Сравнение результатов расчёта с результатами гидравлических испытаний показывает, что для форсунки «Г» расчётные размеры капель значительно превышают размеры капель в реальном распыле. В случае форсунки «О» капельный распыл вообще не моделируется. Вместо этого ниже по течению от поверхности пелены, где происходит распад плёнки на капли, снова формируется сплошная пелена жидкости. Одной из возможных причин низкой точности является недостаточная разрешающая способность расчётной сетки в местах распада пелён воды на капли.

Моделирование течения в гидравлическом тракте смесителя двигателя тягой 13Н с помощью двухскоростной модели расчёта позволяет с высокой точностью определять суммарный корневой угол факела распыла и характер капельного течения после столкновения пелён. Кроме того, можно рассчитать место и угол встречи пелён внешней и внутренней форсунок и, тем самым, найти параметры перемешивания плёнок жидкости (воды) форсунок окислителя и горючего.

Таким образом, имеет место удовлетворительная точность расчёта параметров перемешивания воды по формуле (1). Сравнивая данные огневых и гидравлических испытаний и также расчёты при одних и тех же параметрах форсунок, можно сопоставлять коэффициенты расходного комплекса и перемешивания Кпер (1), разрабатывая рекомендации по увеличению эффективности преобразования компонентов топлива в продукты сгорания в ЖРДМТ тягой 13Н уже на этапе проектного расчёта.

Заключение

Результаты проведённой работы свидетельствуют о возможности создания методики определения параметров центробежных форсунок на основе численных методов однородной модели турбулентности BSL и односкоростной (для одиночной форсунки) или двухскоростной (для смесителя двигателя тягой 13,34 Н) системы уравнений движения двухфазной среды в виде осреднённых по Рейнольдсу уравнений Навье-Стокса, которая с удовлетворительной точностью позволяет рассчитать параметры распыла форсунок.

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

компонентов топлива в продукты сгорания в камере двигателя уже на этапе проектных расчётов.

Библиографический список

1. Андреев Ю.З. Исследование зависимостей характеристик ЖРДМТ тягой 50...400 Н на топливе АТ+НДМГ от основных параметров двухкомпонентной соосной центробежной форсунки и струйных форсунок завесы: дис. ... канд. техн. наук. Нижняя Салда, 2004. 181 с.

2. Градов В.Н. Исследование внутрикамерных рабочих процессов в ЖРДМТ с целью повышения их экономичности: дис. ... канд. техн. наук. Куйбышев, 1974.

3. Царапкин В.С. Особенности процесса преобразования жидкого ракетного топлива НДМГ+№04 в камерах ЖРДМТ: Автореферат дис. ... канд. техн. наук. Ленинград, 1982. 18 с.

4. Неустойчивость горения в ЖРД / под ред. Д.Т. Харрье, Ф.Г. Рирдона. М.: Мир, 1975. 816 с.

5. Дейч М.Е., Филиппов Г.А. Газодинамика двухфазных сред. М.: Энергоатомиз-дат, 1981. 472 с.

6. ANSYS CFX-Solver, Release 10.0: Theory. ANSYS Europe Ltd, 2005. 266 p.

7. Абрамович Г.Н. Теория центробежной форсунки // В сб.: «Промышленная аэрогидродинамика». Вып. 3. М.: БНИ ЦаГИ, 1944. С. 84-88.

RESEARCH OF LIQUID OUTFLOW FROM LOW-EMISSION SWIRL INJECTORS OF A LOW-THRUST LIQUID ROCKET ENGINE

© 2016

Ye. V. Semkin Head of Group, Design Engineering Department, Research & Development Institute of Mechanical Engineering, Nizhnyaya Salda, Russian Federation, sewa20-68@rambler.ru

The paper justifies the possibility of using cold-flow test data to predict the combustion efficiency for a 13.34 N thruster combustion chamber. The results of modeling working liquid (water) flowing through hydraulic passages of single injection elements and coaxial swirl injectors are presented. We performed the simulation using the Reynolds-averaged Navier-Stokes equations that describe a turbulent flow of a two-phase non-compressible liquid. We used both one-velocity and two-velocity flow models. We modeled the turbulence by the BSL isotopic turbulence model on a computational block-structure grid. The computational grid consisted mostly of HEXA elements. The orientation of the elements coincided with the flow direction in the injector hydraulic channels, which improved the flow resolution accuracy in the boundary layer near the solid wall. We compared the results of calculation with the cold-flow test data. The use of the one-velocity flow model in the design area of a single injector and the two-velocity model in the design area of the mixing element yielded the results that agreed satisfactorily with the cold-flow data.

Liquid-phase mixing; swirl injector; flow modeling.

Citation: Semkin Ye.V. Research of liquid outflow from low-emission swirl injectors of a low-thrust liquid rocket engine.

Vestnik of Samara University. Aerospace and Mechanical Engineering. 2016. V. 15, no. 4. P. 150-161.

DOI: 10.18287/2541-7533-2016-15-4-150-161

References

1. Andreev Yu.Z. Issledovanie zavisimostey kharakteristik ZhRDMT tyagoy 50...400 N na toplive AT+NDMG ot osnovnykhparametrov dvukhkomponentnoy soosnoy tsentrobezhnoy forsunki i struynykh forsunok zavesy. Dis. kand. tekhn. nauk [Dependence of 50...400 N NTO+UDMH liquid rocket engine characteristics on the basic parameters of a bipropellant coaxial swirl injector and jet injectors of the cooling film. Cand. Sci. (Eng.) Dissertation]. Nyzhnyaya Salda, 2004. 181 p.

2. Gradov V.N. Issledovanie vnutrikamernykh rabochikh protsessov v zhidkostnom raketnom dvigatele maloy tyagi s tsel'u povysheniya ikh ekonomichnosti. Diss. cand. techn. nauk [Research of chamber processes of liquid rocket thrusters with the aim of increasing their cost efficiency. Cand. Sci. (Eng.) Dissertation]. Kuybyshev, 1974.

3. Tsarapkin V.S. Osobennosti protsessa preobrazovaniya zhidkogo raketnogo topliva NDMG+N2O4 v kamerakh zhidkostnykh raketnykh dvigatelei maloy tyagi. Avtoreferat. diss. cand. techn. nauk [Peculiarities of UDMH+N2O4 liquid propellant conversion in combustion chambers of liquid rocket thrusters. Extended abstract of Cand. Sci. (Eng.) Dissertation]. Leningrad, 1982. 18 p.

4. Harrje D.T., Reardon F.G. Liquid propellant rocket combustion instability. Washington: National Aeronautics and Space Administration, 1972. 1314 p.

5. Deich M.E., Philippov G.A. Gazodinamika dvukhfaznykh sred [Two-phase fluid gas dynamics]. Moscow: Energoatomizdat Publ., 1981. 472 p.

6. ANSYS CFX-Solver, Release 10.0: Theory. ANSYS Europe Ltd, 2005. 266 p.

7. Abramovich G.N. Teoriya tsentrobezhnoy forsunki. Vsb.: «Promyshlennaya aerogi-drodinamika». Vypusk 3. Moscow: BNI TSAGI Publ., 1944. P. 84-88. (In Russ.)