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

Вычислительные технологии

Том 16, № 4, 2011

Численное моделирование течений в детонационном двигателе

С.Н. Мартюшов1, Я. Г. Мартюшова2 1 Комитет по энергетике, Государственная дума, Федерального собрания

Российской Федерации, Москва 2 Московский авиационный институт, Россия e-mail: martyush@mail.ru

Для численного моделирования течений реагирующих смесей кислород — водород и водород — воздух использовалась упрощенная двухстадийная модель, включающая индукционный период и последующий период реакции. Газ предполагался невязким. В качестве тестовой задачи рассчитывались течения с переходом горение — детонация в каналах с сужениями специальной формы. Основной предмет исследования — численное моделирование течений горения и детонации, возникающих в спиновом пульсационном детонационном двигателе, с целью оптимизации геометрических характеристик спинового пульсационного детонационного двигателя, определяющей получение устойчивой стационарной спиновой детонационной волны.

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

Исследование течений горения и детонации в газовых смесях является одним из современных направлений численного моделирования в газовой динамике. Интерес к таким задачам вызван тем, что физическая природа процессов возникновения и развития детонации не полностью изучена. В этих условиях численный эксперимент дополняет натурный, а выбор математической модели (учет вязкости, эффектов турбулентности, диффузии, числа компонентов химических реакций) помогает определить механизмы, ответственные за возникновение детонации. Следует подчеркнуть, что режим детонации энергетически более выгоден, чем обычное сжигание топлива, что было установлено в работах Я.Б. Зельдовича еще в середине XX века. В связи с этим актуальна разработка различных конструкций детонационных двигателей. Дополнительные преимущества численного эксперимента в данном случае — возможность определения структуры течений внутри конструкций и многократное изменение их геометрических характеристик в процессе расчетов.

1. Математическая модель

Для численного моделирования течений реагирующих смесей водород — воздух в каналах с сужениями специальной формы, предложенными в [1], и в спиновом пульсационном детонационном двигателе конструкции, описанной в [2], использовалась упрощенная модель двухетадийной реакции [3]. Система уравнений газовой динамики для

идеального газа и упрощенные уравнения кинетики могут быть представлены в следующей интегральной форме:

V

пЯ^в + Ф = 0,

(1)

где Q = (р, т, ре, рв, ра) — вектор консервативных переменных, р, е, т плотность, энергия и вектор импульса соответственно; а = 1/т — параметр индукции, обратный к времени индукции, в ~ объемная плотность реагирующего компонента, Ф = (0, 0, 0^вр, рмр,р'шО) — источниковый член. Вектор потоков в нормальном к грани контрольного объема направлении может быть представлен в виде

Я = (т, (тт)/р + Р1, т(е + Р)/р, (тВ)), В = (в,а),

где Р = рЯТ, е = ЯТ/(у — 1) + V2/2 + вч ~ давление и энергия конечного объема, Я — газовая постоянная, Т — температура, 7 — отношение теплоемкостей, ч — энергия реакции. Для описания скоростей изменения кинетических переменных используется гипотеза Арреннуса

= -г- =-= —кхРещэ^—Ех/ЯТ),

<И

—к2Р2

в2 ехр

—Е2

— (1 — в )2 ехр —

Е2

Ч

ЯТ ЯТ

где к1,к2,Е1,Е2 — константы реакции, определенные, например, в [3]

(2)

2. Численный метод

Для пространственной аппроксимации векторов потока в нормальном к грани контрольного объема направлении использовались две известные разностные ТУБ-схемы, Схема Хартена [4] с незначительной модификацией, изложенной в [5], запишем в виде

Г.?+1/2 = 1/2

¥3 + . + £К.+1/2 (1/2Ф(а$+1/2)(у5 + ^5+1) —

г=1

— С(а

. +1/2

+

а

. +1/2

д. = ИшИет(а 5+1/2 ,а.-1/2), Ф(г) = О(г) — Хг2,

г Г ф(а5+1/2)(д.+1— д.)/а5+l/2, а5+1/2 = 0

'5+1/2 1 0, а!.|1/0 = 0,

5+1/2

0(г) =

|г|, |г| > 8

(г2 + 82)/28, |г| <8.

(3)

где а.+1/2 = ~Ц'5+1/2^5+1 — QJ) — компоненты вектора характеристических переменных в дельта-форме, а.+1/2, К.+1/2, Ь.+1/2 — собственные значения, правые и левые собственные векторы матрицы А = Переход к характеристическим переменным

является общим элементом ТУБ (а также ЕМО) схем, обеспечивающим расщепление системы уравнений газовой динамики (1) на отдельные уравнения.

Второй использованной была схема Чакраварти — Ошера [6], в соответствии с которой вектор потока в нормальном к грани контрольного объема направлении представляется в виде

^.+1/2 = Н,.+1/2 + X; ((1 + ¿)5?Ь+1/2 + (1 - А5+1/2К}+1/2/4-

г=1

6

- X ((1 + 5)^+1/2 + (1 - 6)5^+1/2^+1/2^+1/2/4, (4)

г=1

где

Н

■¿+1/2

^¿+1) + г(д.) - X (Л+ - л*-) «2,^+1/2^^+1/2

г=1

а 1,^+1/2 = ^¿+1 /2- Яз), а 2,3+1 /2 = ^¿+1 /2-

а3,3+1/2 = Ц+1/2(Я3+2 - Я3+1),

«1,3+1/2 = ттто(1(а1,.+1 /2, Ьа2,3+1/2), 52.+1/2 = ттто(1(а2.+1/2,

«2,3+1/2 = ттто(1(а2,.+1/2,Ь а3,j+l/2), 5з.+1/2 = ттто(1(аЗ^^ Ьа2,j+l/2),

где Ь = (3-6)/(1-6) — параметр сжатия схемы, 6 = -1 и 1/3 — параметр интерполяции, определяющий порядок схемы, Лг+, Лг- — собственные значения матрицы А =

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

дп+2 = ¿(2д^)дп; £(2Дг) = (дг)^(дг). (5)

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

$ = 30 +

$ = (330 + + Д*ВД1))/4, (6)

= (30 + 232 + дщдг2))/з.

Такая временная дискретизация производилась с применением схемы Чакраварти — Ошера, что обеспечивало третий порядок точности алгоритма по всем переменным.

При решении задач горения и детонации оператор шага по времени Ь(О) обычно разбивается на газодинамический и кинетический операторы шага. В общем случае точность и устойчивость алгоритма зависит от порядка решения газодинамического и кинетического этапов оператора шага по времени [7]. В использованном алгоритме явный оператор шага по времени Ь^) в направлении организован как единый в следующей последовательности:

\ тли+1/2 -пп+1/2

1) вычисляются потоки г^+1/2, г^_1/2 через грани ячеики в координатном направлении ]]

2) вычисляются новые значения концентраций реагирующих компонентов в ячейке при помощи интегрировании системы обыкновенных дифференциальных уравнений кинетики методом Гира и соответствующее изменение полной энергии в источниковом члене Ф;

3) новые значения газодинамических переменных находятся с учетом вычисленных

т^п+1/2 т^п+1/2

конвективных потоков г ^+1/2 , г ^_1/2, концентрации реагирующих компонентов и изменения полной энергии, обусловленных химическими реакциями.

Структурированные эллиптические сетки для разностной задачи строились по алгоритму томпсоновского типа [8], основанному на решении системы уравнений Пуассона [9].

3. Моделирование течений в каналах сложной формы

В качестве тестовой задачи рассчитывались течения с переходом ударная волна — детонация в осесимметричных каналах с одним и двумя сужениями специального вида [1]. Визуализация течений производилась картиной изопикт (серые линии) и изолиний плотности водорода (черные линии) (рис. 1). Сужения при размерах канала 0 < х < 8, у < 0.5 задаются образующими, определяемыми по формулам

1) х = 1.6 - 2.965у2, 2) х = 1.354у2 + 1.3, 3) х = 2.6 - 3.658у2. (7)

В начальный момент времени горючая смесь заполняла весь канал за исключением камеры в левой части канала, заполненной не реагирующим газом под высоким давлением. Этот газ при Р1 = 26 атм и Т1 = 794 К в начальный момент времени прорывает диафрагму и образует ударную волну высокой интенсивности, движущуюся направо. Ударная волна инициирует быстрое горение газа, которое для конфигурации рис. 1, б переходит в детонационную волну. По мнению авторов [1], детонационная волна появляется в "детонационном мешке" между двумя сужениями специальной формы и выходит на стационарный режим в цилиндрическом сужении (см. рис. 1,6). Возникшая таким

Рис. 1. Изопикты (серые линии) и изолинии плотности горючего компонента (черные линии) для канала с одним (а) и двумя (б) сужениями

образом детонационная волна догоняет фронт ударной волны. Для конфигурации рис, 1, а возникающий при тех же начальных параметрах падающей ударной волны фронт быстрого горения не переходит в детонацию, В работе |1| были проведены численные расчеты, согласно которым дня инициации детонации в каналах с одним сужением (см, рис, 1,а) требуется существенно более высокое давление в камере предварительного смешения горючей смеси, чем дня капана с двумя сужениями (см, рис, 1,6), Картины изолиний показывают, что дня капана с одним сужением фронт горения отстает от головной ударной волны, а дня канала с двумя сужениями переходит в детонационную волну, которая догоняет головную ударную волну.

4. Моделирование течений в спиновом детонационном двигателе

Спиновый, или ротационный, детонационный двигатель был предложен Б.Ю. Войце-ховским |10|. В настоящее время такие двигатели конструируются и исследуются рядом авторов, в том число численно |11|. В |2| предложена компоновка спинового двигателя с раздольной подачей подготовленного воздуха и горючего компонента с образованием ротационной детонационной волны не на цилиндрической поверхности, а на конусе. Конструкция |2| (рис, 2) состоит из кольцевого тороидального канала 1, но которому воздух с определенными значениями температуры и давления распространяется в кольцевую коническую камеру сжигания 2. Водород впрыскивается через кольцевые сверхзвуковые сопла, расположенные либо в точке сопряжения кольцевого тороидального канала и камеры сгорания 5, либо в начало конической поверхности камеры сгорания 6. Процесс спиновой детонации на конической поверхности инициируется зажиганием алюминиевой фольги. В экспериментах |2| была получена устойчивая стационарная импульсная и спиновая детонация водородо-воздушной смеси, порожденная волной разрежения на границе тороидального канала и конической камеры сгорания. При конструировании двигателей данной конструкции должны быть оптимизированы, в частности, такие показатели как смешение топлива и подаваемого сжатого воздуха, временные и пространственные параметры инициации детонационной волны, давление и плотность смеси, концентрация горючего компонента на поверхности конуса, тяга двигателя, причем последний должен существенно увеличиться при переходе от цилиндрической конструкции к конической. Варьирование формы тороидальной полости, изображенной па рис, 2, может существенно повлиять на указанные параметры. Дня численного исследования оптимизации конструкции были выбраны сужения специальной формы (7), предложенные в |1|, При этом капан в форме кругового цилиндра заменяется на тороидальную область. Знание механизма возникновения горения

1 5

<N / „ /

и " / >п.

/

6 2

Рис. 2. Схема ешшовшх) импульеншх) детонационншх) двигателя [2|

и перехода к детонации в "детонационном мешке", образованном сужениями (7), может быть полезным при конструировании рассматриваемого типа двигателя. Цель оптимизации состояла в следующем:

1) улучшение перемешивания горючего компонента с потоком воздуха;

2) ускорение инициации горения и детонации смеси;

3) повышение концентрации горючего компонента смеси на конической поверхности;

4) обеспечение нераспространения горения и детонации вверх по потоку.

Расчет производился в несколько этапов. На первом этапе на установление рассчитывалось протекание не реагирующего газа высоких давления и температуры через тороидальную конфигурацию слева направо. В начальный момент времени воздух с высокими давлением и температурой Р1/Рте = 15, Т/Тте = 2.5 поступал с левого конца расчетной области, заполненной воздухом с атмосферными параметрами. Расчет продолжался до установления стационарного распределения в расчетной области. На втором этапе осуществлялся "вбрызг" под различными углами к основному потоку водорода из коллектора через сверхзвуковое сопло, расположенное в "детонационном мешке". При этом параметры сопла задавались таким образом, чтобы давление на срезе сопла превышало давление воздуха на левом конце рассматриваемой области. Краевые условия на кольцевом сопле определялись аналогично заданию краевых условий в спиновом цилиндрическом детонационном двигателе [1]. Обозначая: Р0 _ давление в расчетной области на границе кольцевого сопла, Рт, рт, Тт — давление, плотность и температура водорода в коллекторе, ркг , икг, итах — критическая плотность, критическая скорость и максимальная скорость сопла, Р1, Р2 — расчетные режимы дозвукового и сверхзвукового истечения сопла, удовлетворяющие уравнению

(Р/Рт)1/7(1 - (Р/Рт)<7-1)/7)1/2 = (2/(7 + 1}}1/(7-1)((7 - 1)/(7 + 1))1/2(Яг/5).

Тогда на срезе кольцевого сопла задаются следующие краевые условия:

( и = 0, Ро >Рт,

\ Р = Ро, Р1 <Ро <Рт, (8)

Р = Р1, Ро > Р1,

Рис. 3. Изолинии плотности смеси (серые линии) и горючего компонента (черные линии) для трех компоновок двигателя; цифрами 13 показаны сужения, описываемые формулами (7)

Рис. 4. а — расчетная сетка для геометрии рис. 3, в (изображена каждая десятая координатная линия), сверху — внутренняя коническая поверхность, на заднем плане — цилиндрическая наружная поверхность области, справа — выходной срез расчетной области; б — изолинии плотности смеси (серые линии) и горючего компонента (черные линии) на конической поверхности, переменные по осям абсцисс и ординат соответствуют номерам у, к узлов сетки на конической поверхности

ju = umax (1 - (P/Pm -1)/Y )1/2, P2 < Po < Pm, , .

\ p = Pm(P/Pm)1/Y, 1 j

( puS = PkrUkrSkr, Po < P1, , 4

l7P/(7 - 1)p + u2/2 = umax/2. (ш>

На втором этапе вычислялась тяга двигателя по формуле

P = 2п J(Pwall - Po)|y|nxdS.

S

Первые два этапа рассчитывались в осесимметрической постановке для трех компоновок тороидального канала, "детонационного мешка" для смешения воздуха и водорода и конической камеры сгорания. На рис. 3 приведены картины изолиний плотности смеси и горючего компонента для трех компоновок геометрии двигателя. Сравнение полученной в расчетах величины тяги показало явное преимущество компоновки, представленной на рис. 3, в.

На третьем этапе в конической части запускалась спиновая детонационная волна. Расчет проводился в трехмерной геометрии. На рис. 4, а представлена расчетная сетка размерности 60 х 240 х 120 узлов (для сравнения проводились расчеты на более грубой сетке 30 х 150 х 120 узлов, показавшие качественное совпадение результатов с первым случаем). На рис. 4, б приведена картина изолиний плотности смеси и водорода на поверхности конуса в камере сгорания (координаты соответствуют номерам j, k узлов сетки).

Таким образом, результаты расчетов показали преимущество компоновки, приведенной на рис. 3, е. В планы авторов входит продолжение численных исследований в трехмерной постановке с целью получения устойчивой стационарной спиновой детона-

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

[1] Семенов И.В., Уткин П.С., Марков В.В. Численное моделирование двумерных детонационных течений на многопроцессорной вычислительной технике // Вычисл. методы и программирование. 2008. Т. 9. С. 119-128.

[2] Bykovskii F.A., Zhdan S.A., Vedernikov F.F. Continuous spin and pulse detonation of hvdrogen-air mixtures in supersonic flow generated by a detonation wave // Proc. of 22nd ICDERS. Minsk, Belarus, 2009.

[31 Korobeinikov V.P., Levin V.A., Markov V.V. et al. Propagation of blast waves in combastible gas // Astronáutica Acta. 1972. Vol. 17. P. 529-537.

[4] Y ее H.C., Warming R.E., Harten A. Implicit total variation diminishing (TVD) schemes for steady-state calculation //J. Сотр. Phvs. 1985. Vol. 57. P. 327-361.

[5] Мартюшов C.H. Расчет двух нестационарных задач дифракции на основе явной TVD-схемы Хартена // Вычисл. технологии. 1996. Т. 1, № 2. С. 82-89.

[6] Chakravarthy S.R., Osher S. A new class of high-accuracy TVD-schemes for hyperbolic conservation laws // AIAA Pap. 1985. 85-0363.

[7] Mahmoudi Y., Mazaheri K. Operator splitting in simulation of detonation structure // Proc. of 22nd ICDERS. Minsk, Belarus, 2009.

[81 Thompson J.F., Warsi Z.U.A., Mastín C.W. Numerical Grid Generation. N.Y.: North Holland, 1985. 685 p.

[9] Мартюшов C.H. Numerical grid generation in computational field simulation // Proc. of 6th Intern. Conf. Greenwich Great Britain, 1998. P. 249-256.

[10] Войцеховский Б.Ю. Стационарная детонация // Докл. АН СССР 1959. Т. 129(6). С. 1251-1256.

[11] Zhdan S.A., Bykovskii F.A., Vedernikov F.F. Mathematical modeling of a rotating detonation wave in a hydrogen-oxvgen mixture // Combustion. Explosion and Shock Waves. 2007. Vol. 43(4). P. 449-459.

Поступила в редакцию 17 июня 2011 г., с доработки — 23 июня 2011 г.