CC BY
10УДК 533.6+536.8
Комплексное моделирование движения летательных аппаратов и процессов горения в двигателях на основе анизотропных топлив
© Ю.И. Димитриенко, И.Д. Димитриенко, М.Н. Коряков МГТУ им. Н.Э. Баумана, Москва, 105005, Россия
Предложена модель совместного расчета аэробаллистических характеристик неуправляемых летательных аппаратов и характеристик горения в камере сгорания двигательных установок с анизотропными твердыми топливами. Модель позволяет осуществлять комплексное исследование взаимовлияния процессов анизотропного горения топлив и аэробаллистических характеристик аппаратов. Приведен численный пример комплексного моделирования, демонстрирующий возможности модели для исследования влияния анизотропии топлив на характеристики аппаратов.
Ключевые слова: горение, газодинамика, анизотропные твердые топлива, численное моделирование, комплексное моделирование.
Введение. Траектория полета летательных аппаратов с двигателями на основе твердых топлив, как известно [1-4], регулируется не только системой управления, но и особенностями процессов горения самих топлив. В связи с появлением твердых топлив с переменным составом, в том числе конструктивно-анизотропных топлив [5], для которых имеется возможность формирования переменного состава для обеспечения меняющейся тяги двигательной установки, возник интерес к взаимосвязанному, комплексному моделированию аэробаллистических характеристик ЛА и процессов горения топлив в камерах сгорания ДУ. Цель настоящей работы — разработка такого комплексного метода расчета движения летательных аппаратов и процессов горения в двигателях на основе анизотропных топлив.
Расчет аэробаллистических параметров ЛА. Рассмотрим систему уравнений, описывающую плоское движение ЛА с твердотопливной ДУ в системе координат Ох1 [4]:
dVx_ 1
т?2
-((Р(Н, *) - X - X —а - Укр Бта)-----бШ3,
(\Д V 1 1 р ) (К + н)-
„2
((Р(Н, *) - X - X,) — а + Ухкр соб а)--—-соб 3,
& т
&3 1
& т¥„
dn = cos а
dt R + H ' ( )
a = n + ф(0 - а,
t
где m = m0 - J m(t')dt' — массаЛА; P(H, t) = mg0Pf (1 - (1 - n(H)) — сила
0
тяги; X1 = kM2px (H)cx (H,M, a)SM — сила аэродинамического сопро-
тивления в осях, связанных с осью ЛД;, = 2М2р.(Я)с,(Я,М,а).„ —
подъемная сила в осях, связанный с осью ЛД; кр = —-— — коэф-
X1g - хи
фициент, характеризующий распределение масс ЛД; X — проекции управляющей силы на связанную систему координат; рда(Я) = = Я»р» (Я )6» (Я) = р» п( Я) — давление невозмущенного потока воз-
Ус
духа на высоте Я; М» = . » — число Маха для скорости
^УР» (Я Р» (Я )
движения ЛД.
Здесь также обозначены:
Ут — скорость движения центра масс ЛД (м/с) в скоростной системе координат; g0 — ускорение свободного падения (м/с2); т — масса ЛД (кг/м3); т(^') — удельный массовый расход топлива (1/с); Я — радиус Земли (м); Я — высота полета над поверхностью Земли (м); к — показатель адиабаты атмосферы; сх — коэффициент лобового сопротивления ЛД в связанной системе координат; су — коэффициент подъемной силы ЛД; р — удельная тяга двигателя в пустоте на поверхности Земли (с); Р^ — удельная тяга двигателя ЛД в пустоте на высоте Я, (с); ^ = Руд / Руд — относительная тяга на высоте Я (1); Р»(Я) — атмосферное давление на высоте Щ), (Па); р» = Я»р» (0)0» (0) = = 1 атм — атмосферное давление на поверхности Земли; п(Я) = = Р» (Я) / р» — табличная функция (1); рда(Я) — плотность воздуха на высоте Н(0 (кг/м3); 9да(Я) — температура воздуха на высоте Н(1) (К); Яда — газовая постоянная воздуха (Дж/(кг • К)); Мда — число Маха относительно невозмущенного потока воздуха; — площадь миделя ракеты (м2); а — угол атаки; ф — угол тангажа; $ — угол наклона вектора скорости ЛД к местному горизонту; п — местный угол положения ЛД относительно центра Земли и начальной точки старта; х1р, хь, x1g —
расстояния по оси ЛА от носка до центра давления, центра масс и точки приложения управляющей силы (м).
К системе (1) присоединяются начальные условия старта ЛА с начальным углом наклона
г = 0: ^ = 0, п = о, Н = 0, 3 = 30. (2)
Для плоской баллистической траектории движения ЛА (без бокового маневра) дальность полета определяется после решения системы (1) с помощью соотношения Lmax = п^. Функции рда(Н) и 0да(Н) характеризуют состояние атмосферы и являются табличными, приближенные аналитические зависимости для них имеют вид [6]
р»(Н) = р» ехр(-00Н), (3)
/ Н - Н Л
0<х, (Н) = + (+1 - 0<*,г )
V Нг+1 - Нг J
где Нг ^ Н ^ Нг + 1, г = 0, ..., 4.
Значения констант 0 Н. приведены в таблице, Р0 = 0,15 К-1.
Таблица
Значения констант 8 ., И., Вп
ад' V
Константы г = 0 г = 1 г = 2 г = 3 г = 4 г = 5
Н, км г' 0 11 25 46 54 80
0 , К 288 217 217 274 274 183
Для решения задачи аэробаллистики (задачи движения ЛА) (1), (2) должны быть заданы следующие параметры:
^ т(г ) , Pуд(г), Р уд^1 ^Ы, X1g, xls, V ^ су> Ф(г) а(г). (4)
Управление движением ЛА осуществляется с помощью выбора тяги РДг), массового расхода т (г'), угла тангажа ф(г) и угла атаки а(г). В данной работе управление полагалось выбранным таким образом, чтобы угол атаки не изменялся за все время движения ЛА.
Расчет коэффициентов аэродинамического сопротивления ЛА. Коэффициенты аэродинамического сопротивления сх, су определяются либо экспериментально, либо численным путем, для чего решается внешняя задача аэрогазодинамики движения ЛА на высоте Н, с числом Маха М и углом атаки а. Методика ее численного решения, использованная в данной работе, описана в [7].
В системе координат, связанной с направлением оси ЛА, аэродинамические коэффициенты сопротивления вычисляются по формулам
с., = с = —
11 рп[ёЪ, су = с2 =11 рп2 ёЪ,
Су = С2 = -1
Я Ъ
Я = 2 РХ2 ,
с. = с = —
1 |ри3ё Ъ,
(5)
где ч — скоростной напор на ЛА; 8м — площадь миделя.
Результаты расчетов аэродинамических коэффициентов по указанному способу на основе 3Б геометрических моделей обтекания ЛА классической формы, состоящей из сферического затупления, конуса и цилиндра, показаны на рис. 1.
Рис. 1. Результаты расчетов аэродинамических коэффициентов на основе аэрогазодинамических расчетов для ЛА с классической геометрической
формой
Система уравнений газодинамики продуктов горения топлива в рабочем тракте РДТТ. Для расчетов характеристик двигательной установки: удельного массового расхода топлива т(¿') и удельной тяги двигателя ЛА в пустоте Р.г используются формулы [4]
V
р _ оШ / _
7 £с
^оЫРо
т£с
т = р^с
(6)
где V
оиР г оШ
скорость и давление продуктов сгорания при выходе из сопла ДУ; 8ои( — площадь выходного сечения раструба сопла ДУ; р^ — плотность твердого топлива; Б — скорость горения топлива; £с(^) — площадь горения топлива в ДУ. Для расчета этих параметров воспользуемся методикой численного решения внутренней задачи газодинамики течения продуктов сгорания в рабочем тракте ДУ.
Рассмотрим модель мгновенного сгорания твердого топлива [7-9], в которой топливо, сгорая, сразу переходит в газовую фазу, состоящую из продуктов сгорания, образование к-фазы учитывать не будем. Движение продуктов сгорания в камере сгорания рассматривается в рамках системы уравнений динамики идеального нетеплопроводного газа, состоящей из уравнения неразрывности, уравнений движения и уравнения энергии. Записанная в бескоординатной форме, эта система уравнений в неподвижной системе координат, связанной с поверхностью горения твердого топлива, имеет следующий вид [5]:
—+ с-Ур + У-ру = 0, дХ
дру
Ы
дрЕ
+ с - У ® ру + У - (ру ® у + рЕ) = 0, + с-УрЕ + У-((рЕ + р) у ) = 0,
(7)
р = Яр0, е = с0 Е = е + |у|2/2.
Здесь р — плотность газа; Х — время; Е — полная энергия газа; е — внутренняя энергия; с — теплоемкость при постоянном объеме; 0 — температура газа; |у|2 = у'у. — квадрат модуля скорости; р — давление; Я — газовая постоянная (Я = М/д; ц — молекулярная масса газа; К — универсальная газовая постоянная); Е — метрический тензор; У — набла-оператор Гамильтона [10]; с — вектор движения поверхности горения.
Решение системы уравнений газодинамики (7) ищется в фиксированной неподвижной области У(Х), соответствующей начальной геометрии рабочего тракта ДУ (камера сгорания с топливом и сопло). Компоненты вектора с в подвижной системе координат Х\ связанной с поверхностью горения, вычисляются по формуле [5]
с =
» а)
2Л
1/2
(8)
где Ра. — компоненты якобиевой матрицы преобразования неподвижных координат X3 в подвижные X'. Скорость горения топлива 1) зависит от давления р продуктов сгорания на твердой поверхности горе-
ния: Б = Бг
( р Л
, здесь Б0, V — характеристики топлива; р0 = 105 Па.
^ 0 У
Вследствие неоднородности топлива его скорость горения является
переменной функцией в разных зонах топлива, поэтому характеристики топлива являются функциями координат: ^0(х), v(x).
Рассматриваются пять типов граничных условий для системы уравнений (7).
На границе Ер представляющей собой твердую непроницаемую стенку (стенки сопла), к системе (1) присоединяется граничное условие прилипания:
V = 0. (9)
На твердой границе Е2 горения твердого топлива, рассматривается условие массоприхода, которое дополняется условием задания температуры поверхности горения
V• п = --(р5-р)д 0 = ес, р
где р^ — плотность топлива; 9с — температура поверхности горения; п — вектор нормали.
В начальный период времени работы РДТТ происходит затекание в камеру сгорания горячих газов, образующихся при работе воспламенителя. На этой границе входа потока Е3 задаются следующие условия:
р = ре, ^п = уе, 9 = 9е, (10)
где ре, V 9е — заданные значения.
После воспламенения основного топлива, начинается рост давления в КС. Как только давление достигает некоторого предельного значения, разрушается заглушка в сопле и происходит истечение продуктов сгорания через сопло наружу. На дозвуковой границе выхода потока Е4 из сопла задается условие [7]
Р = Ре. (11)
На сверхзвуковой границе выхода потока Е4 не задаются никакие граничные условия.
На плоскости симметрии Е5 заряда твердого топлива задаются условия симметрии:
^ = 0, V • п = 0, д0 = 0. (12)
дп дп
Начальные условия к системе (1)-(5) имеют вид
I = 0: р = Р0, V = 0, Е = Сд, (13)
где р0, 90 — заданные значения.
После решения задачи (7)—(13) вычислим параметры продуктов сгоРания: voup p t\ vout = v(х|5 , t)• 1, pout = p(x|5 , Л, а также значе-
\ "out / \ out f
ние Sc(t) переменной поверхности горения топлива.
Методы численного решения задач аэробаллистики и внутренней газодинамики. Для решения системы задач (1), (2) и (7)-(13) вводилось два масштаба времени — «медленное» время t, соответствующее характерным значениям времени процесса выгорания топлива (секунды), и «быстрое» время т, соответствующее характерным значениям установления потока продуктов сгорания (миллисекунды). Задача газодинамики (7)—(13) решалась при фиксированных, меняющихся только параметрически, значениях времени t, поверхность горения при этом рассматривалась как функция Sc (t) . Для численного решения системы уравнений (3)-(7) был применен метод ленточных адаптивных сеток [7], использующий разностную схему типа предиктор-корректор 2-го порядка аппроксимации. Генерация адаптивной конечно-разностной сетки осуществлялась с помощью программы Sigma, разработанной на кафедре ФН-11 [11]. Введение подвижных адаптивных координат позволяет не пересчитывать разностную сетку по мере выгорания топлива, а использовать для этих целей одну и ту же исходную сетку.
Для численного решения системы уравнений (1) применялись численные методы решения систем нелинейных обыкновенных дифференциальных уравнений: методы Рунге—Кутты, пошаговый метод с явной и неявной разностной схемой, для решения которой применяется модифицированный метод Ньютона.
Результаты численного моделирования. Была рассмотрена геометрия твердотопливного заряда цилиндрической формы [12, 13]. Свойство анизотропии скорости горения топлива обусловлено наличием у него неоднородной структуры, так называемой структурной анизотропии, вследствие которой горение топлива происходит неравномерно: на участках с большими значениями скорости горения D образуются новые участки горения, что приводит к возрастанию давления в камере сгорания. Был рассмотрен случай, когда неоднородность топлива имела многокольцевую форму, в форме четырех концентрических колец, вследствие чего при горении топлива образовывались четыре кольцевых канала, распространяющихся по оси цилиндра. Изменение геометрической формы такого заряда показано в работе [5].
Были рассмотрены три варианта анизотропного топлива, отличающиеся длиной участка анизотропии L, выбраны безразмерные значения этого параметра l = 0,4, 0,6 и 0,8. В работе [5] приведены результаты расчетов газодинамических параметров движения продуктов сгорания
в рабочем тракте ДУ в различные моменты времени. На рис. 3 показана характерная кривая изменения безразмерного давления р = р / р0 в камере сгорания в окрестности клиновидной части топлива в различные моменты безразмерного времени, где р0 — характерная величина давления в камере сгорания. На начальной стадии горения из-за быстрого образования новых участков поверхности горения топлива давление повышается и достигает максимальных значений. Затем по мере выгорания площадь поверхности горения топлива уменьшается и давление падает. Вследствие кольцевой структуры анизотропного топлива, его геометрическая форма остается самоподобной, площадь поверхности горения стабилизируется, и также стабилизируется участок с пониженным давлением.
Увеличение значения параметра 1 приводило к увеличению продолжительности начального горения топлива и уменьшению (второго) маршевого участка движения ЛА.
На рис. 2 показаны сравнительные кривые изменения давления в ДУ для трех рассмотренных вариантов. Все расчеты в данном разделе приведены при углах старта $ = 80° и нулевом угле атаки.
Р А
О 0.15 0.3 0.45 0.6 0.75 0.9 1.05
Рис. 2. Расчетные кривые давления р(0 в камере сгорания РДТТ с анизотропным топливом при различных значениях параметра 1: 1.. .0,4, 2.. .0,6 и 3.. .0,8
На рис. 3 показаны сравнительные кривые удельного импульса тяги Р() ДУ в полете модельного ЛА с РДТТ на основе трех рассмотренных вариантов анизотропного топлива.
На рис. 4-6 можно видеть изменение числа Маха М(^), высоты полета #(0, дальности полета Щ) модельного ЛА с РДТТ на основе трех рассмотренных вариантов анизотропного топлива.
В результате проведенных исследований было установлено, что увеличение начального участка горения топлива приводит к уменьшению общего времени работы ДУ, увеличению скорости полета ракеты
Рис. 3. Изменение безразмерного удельного импульса тяги ДУ в полете модельного ЛА с РДТТ на основе анизотропных топлив при различных значениях параметра 1: 1.. .0,4, 2.. .0,6 и 3.. .0,8
Рис. 4. Изменение во времени числа Маха движения модельного ЛА с РДТТ на основе анизотропных топлив при трех различных вариантах параметра 1:
1.0,4, 2.0,6 и 3.0,8
на маршевом участке полета на 0,6 М, небольшому повышению высоты полета Н(Х) на маршевом участке и уменьшению максимальной дальности полета £шах за счет уменьшения суммарного времени полета ЛА.
Выводы. Предложен метод совместного расчета аэробаллистических характеристик неуправляемых летательных аппаратов и характеристик горения в камере сгорания двигательных установок с анизотропными твердыми топливами. Показано, что данный метод позволяет осуществлять комплексное исследование взаимовлияния процессов анизотропного горения топлив и аэробаллистических характеристик аппаратов. Приведен численный пример комплексного моде-
Рис. 5. Изменение во времени безразмерной дальности полета Ь модельного ЛА с РДТТ на основе анизотропных топлив при трех различных вариантах параметра 1: 1...0,4, 2...0,6 и 3...0,8
Н
2 *................
3 1
..........................1
О 0,3 0,6 0,9 1,2 1,5
Рис. 6. Траектория полета (зависимость безразмерной высоты Н от безразмерной дальности Ь) модельного ЛА с РДТТ на основе анизотропных топлив при трех различных вариантах параметра 1: 1.0,4, 2.0,6 и 3.0,8
лирования, демонстрирующий возможности метода для исследования влияния анизотропии топлив на характеристики движения летательных аппаратов.
Исследование выполнено при финансовой поддержке Министерства образования и науки Российской Федерации (номер НИР 1.5433.2011) и РФФИ (грант № 12-08-00998-а).
ЛИТЕРАТУРА
[1] Липанов А.М., Алиев А.В. Проектирование ракетных двигателей твердого топлива. Москва, Машиностроение, 1995.
[2] Калинин В.В., Ковалев Ю.Н., Липанов А.М. Нестационарные процессы и методы проектирования узлов РДТТ. Москва, Машиностроение, 1986.
[3] Ерохин Б.Т. Теория внутрикамерных процессов и проектирования РДТТ. Москва, Машиностроение, 1991.
[4] Синюков А.М., Волков Д.И., Львов А.И., Шишкевич А.М. Баллистическая ракета на твердом топливе. Москва, Воениздат, 1972, 512 с.
[5] Димитриенко Ю.И., Кулагин Ю.А., Ярмола А.П. Моделирование газодинамических процессов в камерах сгорания двигателей с анизотропными твердыми топливами. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. Спец. вып., 2011, с. 100-110.
[6] Димитриенко Ю.И. Механика композиционных материалов при высоких температурах, Москва, Машиностроение, 1997, 366 с.
[7] Димитриенко Ю.И., Котенев В.П., Захаров А. А. Метод ленточных адаптивных сеток для численного моделирования в газовой динамике. Москва, ФИЗ-МАТЛИТ, 2011, 280 с.
[8] Димитриенко Ю.И., Изотова С.Г., Ануфриев С.Н., Захаров А.А. Численное моделирование трехмерных газодинамических процессов в камерах сгорания РДДТ на основе метода геометрически-адаптивных сеток. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2005, № 3, с. 139-146.
[9] Димитриенко Ю.И., Димитриенко И. Д. Численное моделирование процессов горения смесевых твердых топлив. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2001, № 2.
[10] Димитриенко Ю.И. Тензорное исчисление. Москва, Высшая школа, 2001, 560 с.
[11] Димитриенко Ю.И., Захаров А.А. Автоматизированная система для моделирования газовых потоков методом ленточных адаптивных сеток. Информационные технологии, 2009, № 6, с. 12-16.
[12] Смирнов Н.Н., Димитриенко И.Д. Исследование конвективного горения в сжимаемом твердом топливе с продольными каналами. Физика горения и взрыва, 1990, № 4, с.14-22.
[13] Рычков А. Д. Математическое моделирование газодинамических процессов в каналах и соплах. Новосибирск, Наука, 1988.
Статья поступила в редакцию 27.06.2013
Ссылку на эту статью просим оформлять следующим образом:
Димитриенко Ю.И., Димитриенко И.Д., Коряков М.Н. Комплексное моделирование движения летательных аппаратов и процессов горения в двигателях на основе анизотропных топлив. Инженерный журнал: наука и инновации, 2013, вып. 9. URL: http://engjournal.ru/catalog/mathmodel/technic/1120. html
Димитриенко Юрий Иванович родился в 1962 г., окончил МГУ им. М.В. Ломоносова в 1984 г. Д-р физ.-мат. наук, профессор, заведующий кафедрой «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана, действительный член академии инженерных наук. Автор более 250 научных работ в области вычислительной механики, нелинейного тензорного анализа, термомеханики композитов, математического моделирования в материаловедении. e-mail: dimit.bmstu@gmail.com
Димитриенко Ирина Донатовна окончила МГУ им. М.В. Ломоносова в 1984 г. Канд. физ.-мат. наук, ведущий научный сотрудник. Научно-образовательного центра
«Суперкомпьютерное инженерное моделирование и разработка программных комплексов» МГТУ им. Н.Э. Баумана. Автор 40 научных работ по газодинамике горения твердых топлив, вычислительной механике. e-mai1:irma.dimit@gmai1.com
Коряков Михаил Николаевич родился в 1987 г., окончил МГТУ им. Н.Э. Баумана в 2010 г. Аспирант кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Автор 10 работ в области вычислительной газодинамики.
