CC BY
45
Вычислительные технологии
Том 13, № 1, 2008
Моделирование работы
_ «.» *
твердотопливных микродвигателем*
Л. К. ГУСАЧЕНКО
Институт химической кинетики и горения СО РАН, Новосибирск, Россия
e-mail: gusachen@kinetics.nsc.ru
Х. Милошевич Приштинский университет, Косовска Митровица, Сербия e-mail: mhrane@yahoo. com
А. Д. Рычков, Ю.И. Шокин Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: rych@ict.nsc.ru
Numerical modeling is used to investigate a domain of stable operation of solid-propellant rocket micro engines, which are used in various control systems of spacecrafts. It was established that this domain differs significantly from the on of the mid-flight engines.
Введение
"Неакустическая" неустойчивость работы твердотопливного ракетного двигателя (или газогенератора со сверхзвуковым истечением) исследовалась еще более 30 лет назад [1-3]. В этих работах методом малых возмущений получено и исследовано частотное уравнение, в линейном приближении определяющее границы устойчивости двигателя в пространстве входных параметров задачи. Позже были проведены некоторые расчеты и для нелинейной задачи, в том числе в области неустойчивости. К сожалению, упомянутые исследования не дают достаточно полной картины. Частотное уравнение проанализировано лишь для ограниченного набора параметров. В частности, не обсуждался следующий из частотного уравнения парадоксальный вывод о возможности стационарного режима работы двигателей предельно малых размеров для топлив с постоянной температурой поверхности при значении показателя V в зависимости скорости горения от давления щ = и0Ри, большем единицы. Громоздкие процедуры нахождения и анализа корней сложного частотного уравнения создавали опасность получения (при устранении иррациональности) дополнительных ложных корней, о чем предупреждал Б.В. Новожилов в [2]. Тем не менее, в более поздней работе [3] эта ошибка появилась.
* Работа выполнена при финансовой поддержке Программы междисциплинарных интеграционных исследований СО РАН (проект № 2006-28 и № 2006-113), Программы фундаментальных исследований ОЭММПУ РАН (проект № 4.1.1) и Президентской программы "Ведущие научные школы РФ" (грант № НШ-9886.2006.9).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.
Что касается численных расчетов нелинейной задачи, то они проводились главным образом для выяснения характера развития неустойчивости (о возможности получения из расчетов другой полезной информации см. ниже).
В случае двигателей малых размеров (микродвигателей) механизм их работы претерпевает качественные изменения. Возникают новые трудно устранимые виды неустойчивости стационарного режима, и обостряется проблема запуска (выхода на стационарный режим). От микродвигателей, как правило, требуется только создание заданного импульса без особых ограничений на его форму. В некоторых случаях для практики приемлемым может оказаться нестационарный режим выгорания заряда твердого топлива или даже его взрыв (импульсные режимы). Однако наиболее востребованными являются микродвигатели все же со стационарным режимом работы.
1. Математическая модель
При уменьшении размеров двигателя соответственно уменьшается параметр х = ¿с/Ьъ, где ¿с = М/М* — характерное время камеры; ¿ъ = к/и2 — время перестройки прогретого слоя горящего топлива; к — температуропроводность состава, к ~ 0.001[см2/с]; щ — линейная скорость горения; М = Vрд — масса газа в камере сгорания; V — ее объем; рд — плотность газа в камере; М* — секундный массовый расход продуктов через сопло камеры сгорания.
Если ¿с >> ¿ъ, то можно при переменном давлении использовать квазистационарный закон горения, но давление р придется находить из нестационарного баланса массы газа в камере. Я.Б. Зельдович показал [1], что при уменьшении параметра х возможна неустойчивость, для которой впоследствии был введен термин "¿*-нестабильность". Если времена ¿с и ¿ъ сравнимы или ¿с << ¿ъ (что характерно для микродвигателей), то и скорость горения щ приходится находить из нестационарной модели.
Такая нестационарная математическая модель процесса в двигателе включает в себя уравнение баланса массы М газа в камере двигателя
¿М
-— = М+ ¿Ь +
М* - М- (1)
и уравнение теплопроводности в прогретом слое твердого топлива
( дт дт \ д2 Т ор^ + щ^ = А—, <х < т(-^¿) = То, т(М) = ^ (2)
где о и р — теплоемкость и плотность топлива соответственно; А — коэффициент теплопроводности. При этом для прихода и расхода массы газа можно принять
М - Р, М* - и, М* - Р. (3)
Для нестационарной скорости горения согласно [1, 2] следует использовать зависимость иъ = иъ(Р,!), где f = дт(0,Ь)/дх, которую получают исключением начальной температуры т0 из "стационарной" эмпирической зависимости иъ = и°°(Р, т0), известной из эксперимента, и "стационарного" теплового баланса прогретого слоя
Аf = дт(0, ¿)/дх = ориъ(тз - то).
В частности, если
и0(Р,То) = С • Р* ехр(вТо), то связь величин иь, Р, / имеет вид
к/ = и(^(и) - СР^)) ' К = Л/(ср)"
(4)
Если же экспериментальная зависимость имеет вид
и0(Р, То) = А • Р^/(В - То),
то выполняется соотношение
иь(В - Тв(иь)) = А • Ри - к/. (5)
Зависимость Т,(щ) также предполагается известной из эксперимента и соответствует представлениям, изложенным в [2] и некоторых других работах в виде "универсальной зависимости температуры поверхности от скорости горения": щ = С ехр(-В/Т,). После введения безразмерных переменных
0 =(Т - То)/(Т0 - То), т = ¿(и0)2/к, £ = хиО/к, п = Р/Ро' V = и6/и0, ^ = ///0,
используя разложение Франк-Каменецкого, а также принятые в [1, 2] обозначения
г = дТ0/дТо, к = (Т0 - То)5(1п(и0))/5То,
приведем "универсальную зависимость" к виду
0, = 1 + (г/к)1п(У). (6)
В новых переменных система (1), (2) запишется как
х^т = у - п (7)
7^ + V— = т^, <е < 0, 0(-^,т) = О, 0(0, т) = 0,(V). (8)
50 Я0 5 20
--+ V— = —
дт де
Соотношения (4), (5) примут вид
к^ = V[к - (1 - г)1п^) + V 1п(п)], (9)
к^ = V[к - 1 + г 1п^)] + пи• (10)
Стационарное решение для обоих вариантов (и вообще для любого вида стационарного закона горения) выглядит одинаково:
0 = ехр(е), V = 1, п = 1- (11)
2. Аналитический метод исследования устойчивости
Линеаризуем систему (7), (8) вблизи стационарного решения (11). Пусть
П = 1 + Дпехр(ит), V =1 + ДУ ехр(ит), 9 = ехр(£) + Д9(£) ехр(ит).
Тогда
иД9 + в«ДV + Д9' = Д9", Д9(-го) = 0, Д9(0) = (r/k)ДV, Дп = ДV/(1 + хи), Д9'(0) = ДV + ДV(г - 1)/к + Д^/к.
Этот результат одинаков для обоих вариантов и вообще для любого вида стационарного закона горения. Решение для Д9 будем отыскивать в виде
Д9/ДV = А ехр(£) + А 2 ехр(О£).
Подставляя это выражение в уравнение для Д9 и приравнивая слева и справа множители при ехр(£), а затем при ехр(О£), получим:
иАг + 1 = 0, и + О = О2,
откуда Ах = 1/и и с учетом того, что Д9(-то) = 0, имеем:
О = 1/2 + 0/4 + и. Из условия Д9(0) = (r/k)ДV следует, что А2 = г/к + 1/и. Наконец из условия
Д9'(0) = ДV + ДV(г - 1)/к + Д'П»/к
для величины Д9/(0)/ДV = Ах + О А2 следует частотное уравнение
1 (г 1\ / 1 1 \ г - 1 V 1
--+ Т + - X + + и\ = 1 + —;— + --
и \к и) \ 2 V 4 I к к 1 + хи
которое без учета переменности температуры газа и при предположении об "универсальной зависимости температуры поверхности от скорости горения" принимает вид
^ +24и - 1 к(1 + хи) + ^ 1+24и + 1 (1 + хи) - (к + г - 1)(1 + хи) - V = 0.
Для случая постоянной температуры поверхности горения выполняется соотношение _
л/1 + 4и - 1
^-к(1 + хи) - (к - 1)(1 + хи) - V = 0. (12)
Предполагается, что малые возмущения зависят от безразмерного времени т как ехр(ит), так что в области устойчивости для всех корней уравнения (12) должно выполняться условие И,е(и) < 0.
В работе [2] отмечено, что при подстановке и = 0 в (12) из него немедленно следует, что потеря устойчивости при и = 0 может осуществляться лишь в малоинтересном для практики случае V =1. Это высказывание справедливо только частично, что становится
понятным при попытке выяснить, с какой стороны от значения V =1 находится область неустойчивости. Разлагая (12) при ш ^ 0, легко получить
-=Х-к • (13)
Поскольку соотношение получено в предположении |ш| << 1, оно выполняется:
1) при любых V, если х >> 1;
2) при любых х (не близких к к), если IV — 1| << 1.
Первый случай соответствует большим (маршевым) двигателям и действительно неинтересен, поскольку неустойчивость, т.е. ш > 0, достигается при V > 1, что хорошо известно. Гораздо интереснее второй случай, который предсказывает при х < к неустойчивость, если V < 1. Что касается области V > 1, то здесь возможна устойчивость, если нет других (неустойчивых) корней, что нуждается в дополнительной проверке.
В микродвигателях имеет место Х << 1 и выражение (12) принимает вид
2 к + V — 1 (л \ к (ЛЛ\
ш = (1 — ^^-7Г2 • (14)
V1 + 4" + 1 k ' ' \k + V - 1)2'
Согласно (14) в приближении Зельдовича (Ts = const) для устойчивости работы микродвигателя необходимо выполнение v > 1. Заметим, что при получении (14) из (12) использовано условие |"х| << 1, так что для утверждения о достаточности полученного выше неравенства нужна проверка, не дает ли (12) в случае v > 1 другой неустойчивый корень при |w| >> 1, когда упомянутое предположение не выполняется. Пусть |w| >> 1 и х << 1. Тогда (12) принимает вид
-(k - 1)(1+ х") - v = 0, w = V + k - 11, (15)
1 - k х
так что при v > 1 устойчивость обоих корней имеет место только в том случае, если еще и k > 1. Если же k < 1 (что вполне может иметь место), то этот корень неустойчив.
3. Численное исследование устойчивости
Для подтверждения полученных выше выводов с использованием метода малых возмущений было осуществлено численное решение нелинейной системы (7), (8) с нелинейным же граничным условием (6). Для решения уравнения (7) использовалась А-ус-тойчивая схема второго порядка точности [4], для решения уравнения (8) — разностная схема Кранка—Николсона с линеаризацией граничного условия (6).
п 2 -
0.2 0.4 0.6 0.8
Рис. 1. Изменение безразмерного давления во времени: х = 0.001, r = 0, k = 1.5, v = 1.2, <р(0) = 1.5
Рис. 2. Неустойчивый режим работы микродвигателя: % = 0.001, г = 0, к = 0.5, V = 1.5, <р(0) = 1.5
Рис. 3. Неустойчивый режим работы микродвигателя: % = 0.001, г = 0, к = 0.5, V = 0.5, <р(0) = 1.5
На рис. 1 показано изменение безразмерного давления во времени при значениях параметров, определяющих протекание процесса из области устойчивости, предсказываемой линейным анализом. Видно, что действительно имеет место выход на стационарный режим работы микродвигателя, причем характерное время полного установления согласно (14) равно 2п / ш = 60, при этом время "быстрой" подстройки давления и скорости горения к температурному градиенту согласно (15) составляет около 2п/ш = 0.002. При наборе параметров, соответствующих неустойчивому режиму работы микродвигателя (рис. 2 и 3), стационарный режим действительно отсутствует. С увеличением параметра х (увеличение размеров микродвигателя) расчеты подтверждают переход на "нормальный" режим работы двигателя, когда устойчивость достигается при V < 1, к < 1.
Список литературы
[1] ЗЕльдович Я.Б. Об устойчивости режима горения пороха в полузамкнутом объеме // ЖПМТФ. 1963. № 1. С. 67-76.
[2] Новожилов Б.В. Нестационарное горение твердых ракетных топлив. М.: Наука, 1973. 176 с.
[3] Гостинцев Ю.А., Суханов Л.А. К теории устойчивости горения пороха в полузамкнутом объеме // Физика горения и взрыва. 1974. Т. 10, № 6. С. 818-826.
[4] Рычков А.Д. Математическое моделирование газодинамических процессов в каналах и соплах. Новосибирск: Наука, 1988. 222 с.
Поступила в редакцию 23 января 2008 г.
