Численное моделирование рабочего цикла однотрубного импульсного детонационного двигателя Текст научной статьи по специальности «Физика»

УДК 533.6

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ РАБОЧЕГО ЦИКЛА ОДНОТРУБНОГО ИМПУЛЬСНОГО ДЕТОНАЦИОННОГО

ДВИГАТЕЛЯ

В.Ю. Гидаспов, Н.С. Северина

Описывается методика математического моделирования течения в ИДД в квазиодномерной нестационарной постановке. Используется сеточно-характеристический метод [4 - 13], позволяющий точно учитывать траектории движения и перепад параметров на ударных и детонационных волнах, а также зоны раздела между газами продувки и горючей смесью. Сформулирована система граничных условий, моделирующих работу стадий продувки, заполнения канала горючей смесью, инициирования и распространения детонации. Приводятся результаты численного моделирования рабочего процесса в ИДД, работающего на пропано-воздушной горючей смеси, при сверхзвуковом полете с числом Маха M = 5.

Ключевые слова: импульсный детонационный двигатель, численное моделирование, сеточно-характеристический метод с явным выделением разрывов, пропано-воздушная горючая смесь.

В современном представлении импульсный детонационный двигатель (ИДД) - это труба или связка труб, оборудованная системой подачи воздуха и топлива. Один конец трубы (тяговая стенка) закрыт или периодически закрыт в случае использования механического клапана. Другой конец трубы оборудован реактивным соплом. По мере заполнения трубы топливно-воздушной смесью производится инициирование детонации в смеси с помощью того или иного источника инициирования, в результате чего по смеси распространяется детонационная волна, которая, сжигая то-пливно-воздушную смесь, создает высокое давление на тяговой стенке [1 - 3]. Далее давление в трубе снижается после выхода в атмосферу детонационной волны, и процесс повторяется.

В работе рассматриваются квазиодномерные течения газа, влиянием эффектов вязкости, теплопроводности и диффузии пренебрегаем [4 - 14]. Такие течения реагирующего газа в каналах с пологими стенками в областях непрерывности течения описываются квазилинейной системой уравнений в частных производных [14]:

—р^ +—р^ = 0, (1) дt дx

^р^ + ^(рм 2 + p F = p ^, (2)

дt дх дх

+ 0,5м 2 + -д~ри (к + 0,5м 2 )р = 0, (3)

д д

-р^ + — рм&ь = ^, I = 1, к, N. (4)

дt дх

Здесь (1) - уравнение сохранения массы (неразрывности), (2) - импульса, (3) - энергии, (4) - уравнения, описывающие изменения химического состава, F = F (х) - зависимость площади канала от продольной координаты, N - число компонентов в смеси. Выражения для скорости образования г-го компонента получаются по обычным правилам химической кинетики [14]. Система уравнений (1) - (4) замыкается термическим и калорическим уравнениями состояний:

Р=Р( Р,Т, 71,к, 7 N), е = е( Р, Т, 71,к, 7 N).

Необходимо отметить, что для корректности системы (1) - (4) должно выполняться условие

N

I = о,

г=1

где тг - молекулярный вес г-го компонента.

В случае одномерных движений на поверхностях разрыва параметров течений выполняются соотношения Ренкина - Гюгонио, являющиеся следствием законов сохранения [14]:

Рд (В - иД ) = Рз (В - иЗ ) = т,

РД +РД(В - ид)2 = РЗ +РЗ(В - из)2, (5)

РД(В - ид)(ед + ^ + ) = РЗ(В - из)(ез + * + ).

Рд 2 Р З 2

Здесь нижний индекс Д соответствует состоянию газа перед разрывом, З -за ним, В - скорость распространения разрыва.

Рассматриваются два типа поверхностей разрыва: контактные (тангенциальные) разрывы и ударные волны. Для поверхностей первого типа т = 0:

В - ид = В - из = 0.

На КР из (5) получаем условия

иД = иЗ, РД = РЗ.

Важно подчеркнуть, что КР может быть границей раздела между газами с различными уравнениями состояния.

Для разрывов второго типа - ударных волн

В Ф ыд Ф из.

Соответственно выполняются соотношения (5), которые дополняются условиями неизменности концентраций

7дг =7зг, г =1,N.

Термодинамические свойства реагирующего газа описываются с помощью модели многокомпонентного совершенного газа в рамках допущения о равновесной заселенности энергетических уровней, отвечающих всем внутренним степеням свободы молекул и атомов. В указанном случае удельный термодинамический потенциал Гиббса имеет вид [15]

181

о (р, т, 7) = I тг (а,- 0(т) + ят 1п% (б)

г=1 р0

Здесь р0 = 101325 Па - стандартное давление; О0(Т) - температурная часть стандартных молярных потенциалов Гиббса отдельных компонентов. В справочной литературе приводятся полиномиальные аппроксимацион-

ные формулы для приведенного стандартного потенциала Ф°(Т), который связан с о0(Т) по формуле

О0 (Т) = А Н (Т0) - [н0 (Т0) - Н°(0)]- ТФ 0 (Т),

где Нг° (0) - стандартная энтальпия н(Т) при абсолютном нуле; ^Т) -

стандартная энтальпия н(Т) при Т0 = 298,15К. Для задания Ф0(х)

(х = 10-4 Т), в работе используются полиномы, аналогичные приведенным

в справочнике Глушко:

0 —2 —1 2 3

Ф, (х) = ф0 + ф1п 1п(х) + ф_2х + ф-1Х + лх + ф2х + фзх .

Здесь ф^,к = 1п,-2,...,3 — числовые коэффициенты, индивидуальные для каждого вещества. Другие термодинамические величины, используемые при математическом моделировании, выражаются через потенциал Гиббса и его частные производные.

Соответствующие (6) термическое и калорическое уравнения состояния имеют вид

N N п п йОг 0(Т)

Р = РЯ17,Т и н = 17,н0(Т),Н0(Т) = О0(Т)-Т. ,=1 ,=1 ЙТ В случае произвольного механизма из Nr обратимых химических

реакций вида

IVУ)М1 <=> IV(г)М1, г = 1,2, к, ^, г=1 г=1

Nr N - (г) N ~ (г)

Щ = I (V (г) -V (г ^К(г) П (рт] У - К(г) П (РТ])У ),

Г =1 7=1 7=1

-(г )

где Мг - символ г-го вещества; у, - стехиометрические коэффициенты, константы скорости прямых

К(г) и обратных К(г) реакций связаны через константу равновесия. Для аппроксимации температурной зависимости констант скоростей прямых реакций используется обобщенная формула Аррениуса

- ( ЕЛ К(Т) = АТп ехр -Е

V Т

А, п, Е - некоторые постоянные величины, индивидуальные для каждой реакции. Для корректности проводимых расчетов все реакции должны быть обратимыми, при этом константы скорости обратной реакции должны рассчитываться через константу равновесия [14 - 15].

В данной работе для численного решения системы уравнений (1) - (5) используется сеточно-характеристический метод [4 - 13], позволяющий отслеживать положение и перепады параметров для всех типов разрывов. Рассматриваемая система уравнений (1) - (4) имеет три семейства характеристик: звуковые (С+, С") и траектории (С0), вдоль которых выполняются соответствующие соотношения:

Их

С , С : — = и ± а:

Иг

Их

1

Ии ±—Ир ± 1г Иг = 0, ра

1

С0 : — = и: Ик—Ир = 0, Иуг- — ЩИг = 0, I = 1

Иг

N.

Р

Р

где

г± И 1п ¥

Ь~ = иа--+ ^

Их

N Р

- еу. Рт

у

.2

1 2 Р еТ - ррт

Р^ Р реТ — Рте р

I=1р а\вт р р - ер Рт а - скорость звука, нижние индексы "р", "Т" и "у^ " обозначают частное дифференцирование термодинамической функции по соответствующему параметру.

Любой узел расчетной сетки относится к одному из следующих типов: входная и выходная границы (ВГ), жесткие стенки (ЖС), ударная волна (УВ), контактный разрыв (КР), точка на характеристике (ТХ), траектория газа (ТГ), фиксированная точка (ФТ). Вычисление параметров производится с использованием сеточно-характеристического метода. На рис. 1 представлены шаблоны расчета точек некоторых типов.

УВ

КР

-о-

' Л

' /\с-! / 'С /

-Ф^О-1-О-

ТГ

ТХ

л С -

/ //с 0

о / 6 -о—

Рис. 1. Расчетные схемы

183

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

Таблица 1

Типы взаимодействий

ВГ ЖС ТГ ФТ ТХ УВ КР

ВГ Х Х Х Х

ЖС Х Х

ТГ Х Х Х Х

ФТ Х Х Х Х

ТХ Х Х Х Х Х Х Х

УВ Х Х Х Х Х Х Х

КР Х Х Х Х Х

Как уже отмечалось, цикл работы импульсного детонационного двигателя состоит из стадий заполнения канала ствола горючей смесью, инициирования детонации, истечения продуктов детонации через сопло, продувка канала ствола воздухом. Первая стадия рабочего цикла - стадия заполнения - моделируется как распад разрыва с параметрами: слева - параметры стадии (расход, полная энтальпия, состав горючей смеси), справа - параметры на входе в канал. При стадии заполнения считаем, что на вход подается предварительно перемешанная горючая смесь. В результате решения задачи о распаде разрыва образуется контактная поверхность, отделяющая горючую смесь от воздуха.

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

Стадия опустошения запускается в тот момент, когда давление за клапаном подачи становится больше, чем давление перед ним на е (варьируемый параметр). На расстояние примерно одного калибра камеры идет заполнение воздухом, после чего процесс повторяется.

Моделирование граничных условий

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

Расчет граничной точки «Входная граница»

На левой границе на стадиях продувки канала воздухом и стадии заполнения канала горючей смесью будем считать режим втекания дозвуковым (рис. 2).

Рис. 2. Схема расчета входной границы

В рассматриваемом случае характеристики С и С0 направлены за пределы расчетной области, следовательно, необходимо задать число граничных условий равным числу соотношений, выполняемых вдоль данных характеристик. Предлагается считать известными расход, полную энтальпию и состав втекаемой в установку смеси газов:

N 2

IУ Н0(Т) + Ц- = Н 0,

I =1 2

ри = 0.

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

Расчет граничной точки «Выходная граница»

При расчете выходной границы в работе рассматриваются режимы дозвукового (рис. 3), звукового и сверхзвукового истечения, а также режим дозвукового втекания атмосферного газа в установку.

В случае сверхзвукового истечения все характеристики направлены внутрь расчетной области и в дополнительных граничных условиях нет необходимости. В дозвуковом режиме истечения дополнительно считается, что давление на выходе из канала равно атмосферному (заданная величина для данной высоты полета, используется вместо соотношения вдоль характеристики С"). В случае дозвукового втекания дополнительно задается условие равенства полной энтальпии и состава своим значениям в атмосфере (при начале втекания в выходном сечении образуется контактный разрыв). Если в процессе расчетов скорость газа превосходит скорость звука, то вместо соотношения вдоль характеристики С используется условие и = - а.

Рис. 3. Схема расчета выходной границы

Подача топлива

После стадии продувки воздухом необходимо построить математическую модель «впрыска» горючего. Будем считать, что в канал поступает полностью перемешенная стехиометрическая смесь горючего с воздухом, при этом расход воздуха известен. Соответствующая система уравнений имеет вид

р uF = plui F + р 2 u 2 F2 = Q,

2 2 2

(pu 2 + p) F = (pu2 + p ) Fi + (P2u| + P2) F2 = I,

puF

h +

u

2

2

: PlulFl

hl + a. l2

+ P 2u 2 F2

h2 + u2 22

= E,

PuFg i = PlulFlg i1 + P2u2 F2 g i 2 =

186

i = 1k N.

Здесь параметры с индексом 1 - параметры окислителя, с индексом 2 - параметры горючего, без индекса - параметры перемешанной смеси.

Концентрации перемешанной смеси могут быть легко найдены. В случае стехиометрической пропано-воздушной горючей смеси

= 5 или -2— = 5

р 2и 2 12с3 Н 8 0212 С3 Н 8

йТО 1

^ 02 = = 5011Ю2МC3H8,

где М с3н8 - молярная масса С3Н8. Задавая расход 02, тем самым обеспечиваем стехиометрическое соотношение горючего и окислителя в смеси. Концентрации смеси при этом

= Шп!0^, . = ж б

Далее, предполагая, что р2, Т2 известны, р2 находим из термического уравнения состояния, а и2 - из выражения для расхода 02. Зная параметры потоков воздуха и горючего, а следовательно 0, I, Е, можно определить параметры перемешанного потока решением системы нелинейных уравнений методом Ньютона.

Результаты вычислительного эксперимента

Приводятся результаты численного моделирования циклического рабочего процесса в ИДД с механическим клапаном при сверхзвуковом полете с числом Маха М = 5 на высоте И = 20000 м. Схема сверхзвукового ИДД представлена на рис. 4. Осесимметричный ИДД состоит из входного устройства 1, механического клапана 2, топливного коллектора 3, источников зажигания 4, кольцевого обводного канала 5, импульсно-детонационной камеры сгорания 6 с препятствиями 7 и выходного устройства 8 [1 - 3].

Рис. 4. Схема сверхзвукового ИДД

В данной работе расчет ведется в импульсно-детонационной камере сгорания с одним препятствием, длина камеры - 1.71 м, диаметр - 84 мм (рис. 5).

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

Рис. 5. Импульсно-детонационная камера сгорания

Рассматриваемый вариант двигателя работает на стехиометриче-ской пропано-воздушной смеси, продукты сгорания которой включают 7 компонентов (С3Н8, 02, СО, С02, Н2, Н20, N2). Для описания химических превращений в горючей смеси используется глобальная кинетическая схема из пяти реакций окисления горючего (табл. 2) [2].

Таблица 2

Глобальная кинетическая схема из пяти реакций окисления горючего

г Реакция Ц(г) А, (м^моль)Ц( )-1 п Е, К

с

г = 1 С3Н8 + 3,502 ® 3С0 + 4Н20 2 А1 0 -1

г = 2 СО + Н20 ® С02 + Н2 2 1,0 • 109 / р 0 20885,76

г = 3 С02 + Н2 ® С0 + Н20 2 3,1 •Ю10/р 0 24710,62

г II 4 Н2 + Н2 + 02 ® Н20 + Н20 3 7,0-107/ р 0,5 0 10568,70

г I 5 С0 + С0 + 02 ® С02 + С02 3 8,5 • 106/ р1,5 0 10568,70

В табл. 2 г - номер реакции, Ц(г) - молекулярность соответствующих элементарных реакций, т.е. число частиц, которые участвуют в элементарном акте химического взаимодействия, р = р/р0. Реакция 1 рассматривается как бимолекулярная. Для аппроксимации температурной зависимости констант скоростей прямых реакций используется обобщенная формула Аррениуса:

( Е ^

К(Т, р) = АТп ехр - - ,

V Т У

где А, п, Е - некоторые постоянные величины, индивидуальные для каждой реакции. Параметры А, Е для реакции 1 приведены в табл. 3 [2 - 3] в зависимости от температуры, состава смеси и давления.

Таблица 3

Параметры А, Е для реакции 1 в зависимости от температуры,

состава смеси и давления

А1 -1

р, атм 1 20 40 80 120

С3Н8 1,96 • 109 1,96 • 109 1,96 • 109 1,96 • 109 1,96 • 109 20130.85 Т < 775

1,73 -109 8,77 • 108 7,50 • 108 6,41 108 5,85 108 22647.21 Т > 775

Тестирование данного механизма привело к следующим результатам (рис. 6). Для начальных данных р0 = 101000 Па, Т0 = 288,2 К и начального состава (стехиометрия пропан-воздух) были построены ударная и детонационная адиабаты. Также была решена неравновесная задача о распространении стационарной детонационной волны с теми же начальными данными и различными скоростями детонационной волны. Видно, что на фазовой плоскости решения этой задачи (четыре прямые) приходят не на равновесную детонационную адиабату, а на другую кривую (пунктирная линия). Также, если увеличить в реакции 1 константу скорости в сто раз (при О = 2500 м/с, например), то прямая продлевается, т.е. выходит на другие равновесные значения (увеличенный фрагмент).

8е+006 7е+006 6е+006 5е+006

4 4е+006 &

3е+006 2е+006 1е+006

Ударная адиабата Детонационная адиабата Б = 1893 м/с Б = 2000 м/с Б = 2200 м/с Б = 2500 м/с

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

V, м3/ кг

Рис. 6. Р - У-диаграмма процесса

Скорость Ч-Ж по данному механизму 1893 м/с, по теории - 1826 м/с. На рис. 7 представлены распределения температуры за ударной волной. На график нанесены прямые линии - уровни равновесной температуры на равновесной адиабате для соответствующей скорости распространения детонационной волны.

Для устранения отмеченных недостатков в соответствии с работой [13] будем считать все реакции обратимыми, при этом первую реакцию запишем в виде С3Н8 + О2 ^ 6СО + 8Н2О + 5С3Н8, величину А1 умножим на 2/7.

Параметры атмосферы на рассматриваемой высоте составляют р = 5529 Па, Т = 216,7 К, р = 0,0889 кг/м . Расчет начинается с продувки камеры сгорания воздухом, который поступил через воздухозаборник, и на входе в камеру имеет следующие характеристики дозвукового потока с

189

числом Маха М = 0,407: р = 203070 Па, Т = 1007 К, р = 0.7026 кг/м3, и = 257 м/с (данные ИХФ РАН). Стадия продувки длится до формирования установившегося течения воздуха. Ниже представлена рассчитываемая временная развертка течения (рис. 7), а также распределения параметров установившегося потока.

0,0175

0,015

1,2 1,4

X, м

Рис. 7. Временная развертка течения

Заключение

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

Работа выполнена при выполнении государственного задания № 9.7555.2017/БЧ.

Список литературы

1. Импульсные детонационные двигатели // под ред. С.М. Фролова М.: ТОРУС ПРЕСС, 2006.

2. Басевич В.Я., Фролов С.М. Глобальные кинетические механизмы используемые при моделировании многостадийного воспламенения углеводородов в реагирующих течениях // Химическая физика. 2006. Т. 25. № 6.С. 54.

3. Иванов В.С., Фролов С.М. Математическое моделирование рабочего процесса и тяговых характеристик воздушно-реактивного импульсного детонационного двигателя в условиях сверхзвукового полета // Хим. физика. 2011. Т. 30. № 7. С. 48 - 61.

4. Гидаспов В.Ю., Пирумов У.Г., Северина Н.С. Математическое моделирование квазиодномерных нестационарных течений реагирующего газа с произвольным числом взаимодействующих разрывов // Вестник Московского авиационного института. 2008. № 5. С. 83 - 94.

5. Гидаспов В.Ю. Распад разрыва в детонирующем газе // Вестник Московского авиационного института. 2010. Т. 17. № 6. С. 72 - 79.

6. Гидаспов В.Ю., Северина Н.С. Численное исследование газовой детонации в ударной трубе // Вестник Нижегородского университета им. Н.И. Лобачевского. 2011. №4. Ч. 3. C. 714 - 716.

7. Гидаспов В.Ю., Пирумов У.Г., Северина Н.С. Тестирование методики моделирования нестационарных течений газа с ударными и детонационными волнами // Вестник Московского авиационного института. 2011. Т. 18. № 6. C. 119 - 124.

8. Гидаспов В.Ю., Северина Н.С. Численное моделирование экспериментов по определению времени задержки воспламенения за падающими ударными волнами // Физика горения и взрыва. 2013. Т. 49. № 4. С. 31 - 40.

9. Гидаспов В.Ю., Северина Н.С. Численное моделирование тонкой структуры цилиндрической детонационной волны в водородно-воздушной горючей смеси // Теплофизика высоких температур. 2015. Т. 53. № 4. С.556 - 560.

10. Гидаспов В.Ю., Северина Н.С. Элементарные модели и вычислительные алгоритмы физической газовой динамики. Ударные и детонационные волны: учебное пособие. М.: Факториал, 2016. 84 с.

11. Gidaspov V.Yu., Golubev V.K, Severina N.S. A software package for simulation of unsteady flows of the reacting gas in the channel // Вестник Южно-Уральского государственного университета. Сер. «Математическое моделирование и программирование». 2016. Т. 9. № 3. С. 94 - 104.

12. Гидаспов В.Ю., Северина Н.С. Некоторые задачи физической газовой динамики. М.: Изд-во МАИ, 2016. 196 с.

13. Гидаспов В.Ю., Северина Н.С. Численное моделирование детонации пропано-воздушной горючей смеси с учетом необратимых химических реакций // Теплофизика высоких температур. 2017. Т. 55. № 5. С.795 - 799.

14. Пирумов У.Г., Росляков Г.С. Газовая динамика сопел. М.:Наука, Физматлит, 1990. 368 с.

15. Термодинамические свойства индивидуальных веществ: справочное издание в 4 т.// Л.В. Гурвич [и др.]. М.: Наука, 1982.

191

Гидаспов Владимир Юрьевич, канд. физ.-мат. наук, доц., ведущий научный сотрудник, gidaspov@mai.ru, Россия, Москва, Московский авиационный институт (национальный исследовательский университет),

Северина Наталья Сергеевна, канд. физ.-мат. наук, доц., severinaamai.ru, Россия, Москва, Московский авиационный институт (национальный исследовательский университет)

NUMERICAL SIMULATION OF THE WORKING CYCLE OF A SINGLE PIPE PULSE

DETONA TION ENGINE

V.Yu. Gidaspov, N.S. Severina

The paper formulates a system of boundary conditions simulating the operation of purging stages, filling the channel with a combustible mixture, initiating and propagating detonation. The results of numerical modeling of the working process in the PDE operating on a propane-air combustible mixture in a supersonic flight with Mach number M = 5 are presented.

Key words: pulse detonation engine, numerical simulation, grid-characteristic method with explicit discontinuity, propane-air combustible mixture.

Gidaspov Vladimir Yurjevich, candidate of physical and mathematical sciences, do-cent, gidaspovamai.ru, Russia, Moscow, Moscow Aviation Institute (National Research University),

Severina Natalia Sergeevna, candidate of physical and mathematical sciences, do-cent, severinaa mai.ru, Russia, Moscow, Moscow Aviation Institute (National Research University)