CC BY
76
Том ХЬV
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
2014
№ 1
УДК 544.452.42, 536.662
О РАЗЛИЧНЫХ СПОСОБАХ ОПРЕДЕЛЕНИЯ ТЕПЛОВОГО ЭФФЕКТА И ПОЛНОТЫ СГОРАНИЯ В ПОТОКЕ РЕАГИРУЮЩЕГО ГАЗА
В. В. ВЛАСЕНКО
Предложено несколько способов оценки тепловыделения за счет химических реакций и полноты сгорания топлива на основе детальной картины течения с горением, полученной в численном расчете на базе уравнений Эйлера или Рейнольдса. Точность определения полноты сгорания на основе предложенных способов оценки тепловыделения сопоставляется с точностью существующих приближенных методов на примере задачи о наклонной детонационной волне на клине и задачи о горении в плоском канале со сверхзвуковым течением на входе. Приведены примеры использования данных о тепловыделении при анализе физики течения.
Ключевые слова: горение, тепловой эффект, полнота сгорания, уравнения Эйлера, уравнения Рейнольдса, неравновесные реакции, молекулярная и турбулентная диффузия, численное моделирование.
Одной из самых важных физических характеристик течений газа с горением является количество теплоты, выделяющееся в единице массы газа за счет химических реакций. Например, в классической теории одномерного распространения волн горения в предварительно перемешанной горючей смеси газов все характеристики волны горения оказываются функциями лишь двух параметров — безразмерной скорости (числа Маха) волны горения относительно исходного газа и теплового эффекта химических реакций [1]. Теплота, выделяющаяся в камере сгорания ВРД, является основным источником тяги двигателя. При этом в большинстве практических приложений процесс горения протекает в условиях, когда топливо и окислитель не перемешаны или перемешаны не полностью, а течение имеет достаточно сложную пространственную структуру. В этих условиях количество выделившейся теплоты может существенно отличаться от значений, которые возникают при идеальном горении полностью перемешанных реагентов в замкнутом объеме. Приходится вводить понятие полноты сгорания (называемой иногда коэффициентом тепловыделения или эффективностью горения [2 — 4]), которое характеризует понижение
тепловыделения по сравнению с идеальным значением. Коэффициент полноты сгорания обязательно фигурирует в известных методиках расчета тягово-экономических характеристик ВРД [5, 2, 4]. Однако, как отмечено в [6], «до сих пор не существует надежного теоретического и даже численного метода предсказания неполноты сгорания». Полнота сгорания обычно определяется экспериментально для каждой конкретной геометрии камеры сгорания. Для этого существуют различные методы: газодинамический (в котором полнота сгорания определяется из квазиодномерной теории по результатам экспериментальных измерений продольных распределений давления и энтальпии торможения потока — см., например, [7]), по результатам анализа проб газовой смеси [8] и др. Следует отметить, что используемые в инженерных методиках уравнения баланса энергии, в которые входит полнота сгорания, являются приближенными и не учитывают особенностей развития процесса горения вдоль линий тока.
ВЛАСЕНКО Владимир Викторович
кандидат физико-математических наук, начальник сектора ЦАГИ
В последние десятилетия в связи с успешным развитием компьютерной техники появляется все больше работ, содержащих детальное численное моделирование течений с горением с учетом пространственной структуры течения, молекулярного и турбулентного смешения, конечной скорости химических реакций. Расчет дает приближенную, но достаточно полную картину течения, содержащую информацию о значениях всех параметров потока в каждой точке. Это открывает возможность анализа распределения химического тепловыделения по пространству, что в принципе может стать основой активного управления процессом тепловыделения. Кроме того, появляется возможность уточнить полноту сгорания и сравнить ее с существующими инженерными методиками. В настоящей работе предлагается несколько способов оценки химического тепловыделения по имеющейся детальной картине течения и на основе результатов численного моделирования различных течений сопоставляется точность разных методов определения теплового эффекта и полноты сгорания.
Для описания турбулентных течений в данной статье используется подход Рейнольдса. Система уравнений движения газа осредняется по времени. Поскольку рассматриваются существенно сжимаемые течения, то для отдельных параметров используются осреднения по Фавру [9 — 11]. Прямая черта означает осреднение по времени; а = ра/р — осреднение параметра а по Фавру, а" = а — а — пульсация по Фавру. Всюду предполагается суммирование по повторяющимся индексам, связанным с пространственными координатами.
Численный метод, который использовался в расчетах, описан в [12], а его дополнительные особенности, связанные с моделированием химических реакций, даны в [13]. Следует также указать на работы [14, 15], в которых описаны технология дробного шага по времени и граничное условие «закон стенки» с квазиравновесным теплообменом, используемые в данной работе.
1. ОПРЕДЕЛЕНИЕ КОЛИЧЕСТВА ТЕПЛОТЫ, ВЫДЕЛИВШЕГОСЯ ВДОЛЬ ТРАЕКТОРИЙ ОСРЕДНЕННОГО ТЕЧЕНИЯ
Рассматривается турбулентное течение многокомпонентного вязкого газа с неравновесными химическими реакциями. Пусть газ является смесью идеальных газов и температуры всех компонент совпадают. Тогда уравнение энергии можно записать следующим образом [11]:
_БЕ д р-+ —
Вг дх
рй1 +х]1й] +
к=1
= 0, (1)
В д _ д
где р — плотность; — =--+ и--лагранжева производная вдоль траекторий осредненного
Вг дг дх^
по времени течения; ЕЕ — полная энергия единицы массы газа; — сумма молекулярного и турбулентного потоков у-й компоненты импульса вдоль оси х1; ^ — сумма молекулярного и турбулентного потоков тепла вдоль х1, обусловленных градиентами температуры Т ; Н (Т) — статическая энтальпия единицы массы к-го компонента; — сумма молекулярного и турбулентного потоков массы к-го компонента смеси вдоль х1; N — число компонент газа.
Далее уравнение (1) будет преобразовано так, чтобы выделить в явном виде вклады от диффузии и химических реакций в изменение разных составляющих полной энергии Е вдоль средних траекторий объемов газа.
Всюду предполагается, что можно пренебречь корреляциями Т'У^ (У^ — массовая доля к-го компонента). Тогда уравнение состояния принимает вид:
^р у
р = = рДТ, (2)
к=1 тк
Ns,
где Я = Я0 ^ —к--газовая «постоянная» для данного состава смеси (Я — универсальная газо-
k=1 mk
вая постоянная; тк — молекулярный вес к-го компонента). Дифференциал полной энергии можно представить следующим образом:
dE = d [00. + k + ё| = d \U-U
Ns
sp f
-dk + CVdf + ^
k=1
RT
m,
dY,
k
(3)
где к = и" и' /2 — кинетическая энергия турбулентности; е — внутренняя энергия единицы мас-
^р _
сы смеси; Су = ^ ^кСук (Т) — удельная теплоемкость единицы массы газа при постоянном объ-
к=1
еме и при заданном составе газа.
С учетом (2) и (3) уравнение энергии (1) можно преобразовать к виду:
В (и,и, \ Эк ~ ВТ 1 1,1 '-Су -
Dt \ 2
Dt
Dt
р дх1
1 JYk =const
1 д Г— .- _
-1 т,1и + -1 +
р дх1 L
~J1"J ' ~ч
Jkfik ]■
N„
T dR ^ h T~-L hk
k=1
DYk
(4)
В левой части уравнения (4) стоят скорости изменения трех составляющих полной энергии (средней кинетической энергии, кинетической энергии турбулентности и тепловой энергии газа)
, дТ ЯТ дри1
вдоль средних траектории.
f 1 дрщ Л
Р дх,
= UR
— работа сил давления при за-
1 sYu =const
данном составе газа. Член
f ^-z л
dt k Dt k=1
дх1 р дх1
обусловлен изменением состава газа вдоль траек-
торий. В свою очередь, изменение состава газа связано с молекулярной и турбулентной диффузией массы, а также с химическими реакциями:
■DYk
I—-
Dt
J
дх1
Nr
- mk klWl '
(5)
i=1
где W^ — скорость 1-й реакции; — число реакций, а Vkl — коэффициент при к-м веществе в уравнении 1-й реакции. Если данное вещество расходуется в ходе реакции, то Vи > 0, а если образуется, то Vы < 0.
После подстановки (5) в (4) уравнение баланса энергии приводится к виду:
Dt \ 2
D f uU I Dk ~ DT
' '■ Cjv
Dt
Dt
f 1 дри ^ р дх1
р дх1 L J' J — ] р Z дх
1 JYk = const
р k=1
дх1
Nr
f dR 1 Nsp
^ + - Z hkmk Y,VklWl дt р k= U
(6)
Сначала будет рассмотрен случай стационарного течения. Тогда D /Dt = Ui д/ дх1 — производная вдоль средних линий тока. Пользуясь уравнением неразрывности и введя удельную теплоемкость единицы массы смеси при р = const и при заданном составе — Cp = Cv + R, можно привести (6) к виду:
D (utut
Dt
Die ~ DT"
--+ C p —
Dt p Dt
(1 d г _n 1 ^P- dh t Л
P j + q ]+p S Jk
ex,
1 Nsp N _
-ShimiSvkiW . (7)
P k=1 i=1
В невязком ламинарном течении с постоянным составом газа правая часть уравнения (7) равна нулю и, следовательно, равна нулю сумма скоростей изменения кинетической энергии и статической энтальпии вдоль линий тока. Это значит, что один вид энергии может переходить в другой, оставляя полную энтальпию смеси неизменной. При наличии молекулярной и турбулентной диффузии, а также химических реакций сумма членов, стоящих в левой части, может меняться вдоль линий тока, и члены правой части уравнения (7) характеризуют причины этого
изменения. Член
1 д г _] 1NP- dh Л
LT ии j+q, ]+_S jk, pdx, p k=1
dXi
содержит вклады молекулярной и турбулент-
ной диффузии импульса, энергии и массы. Последний член в уравнении (7) представляет собой вклад от химических реакций. Его можно еще переписать в виде
1 (Мф
_ |_ 1 Мг __Мф _
Р2 2 ИкткVк1 Щ = Р2 (<<Р )Щ , где (<<Р)1 = 2 ИкткVк1 — тепловой эффект /-й реакции при
Р I=1 ^ к=1 J Р I=1 к=1
постоянном давлении. Проинтегрировав этот член по времени вдоль линии тока, можно получить энергию (на единицу массы газа), выделившуюся в результате химических реакций:
Nr
S (Q p Wi
Z1) Л _dz
/хим _ I S :
J PiH
(8)
где ёЬ — элемент длины, вычисляемый вдоль линии тока. Здесь учтено, что дифференциал времени вдоль линии тока ^ = ёЬ /1V |.
Более простое для вычислений, но приближенное выражение для дхим можно получить, если пренебречь вкладом от диффузии массы в изменение состава вдоль линий тока и принять, что
m,
Nr
_ nv
SViiW ~—-prr (см. (5)). Тогда Р 1=1 Dt
N„
7(2) хим
-Sjhkdfk = jСpdf -jdh,
(9)
k=1
N„
где И = 2 ^к^к — энтальпия смеси. Поскольку тепловой эффект реакции очень слабо зависит от
k=1
Nsp Nsp
температуры и 2 И (Т )тк Vk/ ~ 2 Ик о тк V к1 , где Ик 0 — энтальпия к-й компоненты при стандарт-
к=1 к=1
*
ной температуре 70 , то (9) можно упростить:
* Величина hi 0 нередко называется энтальпией образования k-й компоненты, так как она совпадает с тепловым эффектом реакции образования единицы массы k-го вещества при f = f0, p = const из стандартных веществ (энтальпии которых при f = T0 принимаются равными нулю). Однако неверно отождествлять ее с потенциальной энергией химических связей в единице массы частиц k-го вещества. В самом деле, в hi0, помимо энергии связей, входят еще кинетическая энергия хаотического движения частиц при f = T0 и работа сил давления. Но при определении тепловыделения за счет химических реакций нет смысла рассматривать особо энергию связей, выделяющуюся в ходе реакций. Тепловой эффект реакции при p = const — это изменение энтальпии смеси за счет реакции, которое включает не только выделяющуюся энергию связей, но и изменение кинетической энергии теплового движения частиц смеси, вызванное изменением состава смеси в результате реакции, а также работу, совершаемую при этом силами давления.
7(3) «-
2ХИМ
к=1
Ё^Л • (10)
Интегрирование в (9) и (10) должно производиться вдоль линий тока.
Следует подчеркнуть, что (9) и (10) — последовательные ступени упрощения точного выражения (8). Каждая ступень упрощения связана с внесением погрешностей.
В нестационарном случае не удается преобразовать уравнение теплового баланса (6) так, чтобы выделить в явном виде вклад от химических реакций. Дело в том, что
Т — = ^Т—. А согласно (5),
д к=1 тк д
1 ~ 1 1 —- + иI —- = —
дt дх^ р
дх,
- т
I=1
/
Это означает, что диффузия массы и химические реакции дают вклад в дУк / д и в дУк / дх1.
Но определить отсюда, как вклад от члена тк 's^vklWl распределяется между дУк / дt и дУк / дх1,
I=1
невозможно. Соответственно, невозможно отделить часть члена Т дД/ дt, определяемую диффузией массы, от части, определяемой химическими реакциями.
Однако в случае, когда газовая «постоянная» I слабо изменяется в ходе реакции (например, при горении керосина), можно приближенно пренебречь членом с Т дД / д . Тогда и для нестационарного течения можно будет оценивать тепловой эффект реакций в потоке формулами (8) — (10).
2. ТЕПЛОВЫДЕЛЕНИЕ В СТАЦИОНАРНОМ НЕВЯЗКОМ ТЕЧЕНИИ
Для примера рассматривается задача об обтекании клина с углом 14° невязким потоком стехиометрической водородо-воздушной смеси с параметрами М = 3, р = 1/6 атм, Т = 850 К. При этих условиях на клине должна возникнуть наклонная детонационная волна. Результаты расчета данного течения на базе уравнений Эйлера с неравновесными химическими реакциями представлены на рис. 1. Для описания горения водорода в воздухе использовалась кинетическая схема Моретти [16], содержащая 8 реакций между 6 активными компонентами с участием инертного азота. Размер расчетной области лишь в 2.5 раза больше, чем длина задержки воспламенения за косым скачком уплотнения, возникающим на клине. Поэтому на рис. 1 видна так называемая «^-структура». Она образуется из-за того, что волны сжатия, создаваемые фронтом пламени, догоняют косой скачок уплотнения и усиливают его до тех пор, пока угол его наклона не станет равным углу наклона детонационной волны, соответствующему данному углу клина [17]. На рис. 1, а показано поле числа Маха, а на рис. 1, б — поле температуры и линии тока, выходящие из центров граней ячеек на левой границе расчетной области. Когда расчет закончился, было произведено интегрирование по формулам (8) и (9). Результирующие поля показаны на рис. 1, в
(поле ) и на рис. 1, г (поле 7(1им). Важно отметить, что для получения достоверных результатов при интегрировании по формуле (8) построение линий тока и аппроксимация интегралов должны проводиться с повышенной точностью, так как скорости химических реакций очень быстро изменяются по пространству. Поэтому для получения корректного поля т^им потребовалась специальная вычислительная процедура, которая не рассматривается в настоящей статье. При определении ¿им использовались гораздо более простые вычислительные методы: линии тока при их построении считались прямолинейными в пределах каждой ячейки, а интегрирование вдоль линий тока проводилось методом трапеций.
Поля т^им и т4им почти совпадают, поскольку для стационарного течения невязкого газа формула (8) должна дать те же результаты, что и формула (9). На рис. 2, а сравниваются распре-
Рис. 1. Невязкое стационарное течение с наклонной детонационной волной на клине: а — поле числа М; б — поле статической температуры [К] и линии тока; в — поле теплового
Рис. 2. Распределения тепловых эффектов д^им и ¿им [МДж/кг] в поперечном сечении расчетной области: - задача о наклонной детонационной волне на клине; б — течение в плоском канале со сверхзвуковым течением на входе, а ~ 1.2
деления д™ и 4хим
(2) в вертикальном сечении, показанном пунктиром на рис. 1. Небольшие отличия этих распределений связаны с тем, что значения д(им и д^им определяются с использованием двух существенно отличающихся приближенных численных процедур. Описание рис. 2, б будет дано ниже.
3. РАЗЛИЧНЫЕ ОПРЕДЕЛЕНИЯ ПОЛНОТЫ СГОРАНИЯ ТОПЛИВА
При рассмотрении течений в каналах (например, в камерах сгорания ВРД) большое практическое значение имеет количество теплоты Q(х), выделившееся на участке между начальным сечением канала и текущим сечением, имеющим координату х. Для определения Q(х) следует построить линии тока, выходящие из начального сечения канала, и определить количество теплоты, выделившееся за счет химических реакций вдоль этих линий тока на данном участке канала — -хим . Эту величину можно определить одним из описанных выше способов (8) — (10). Так как -хим отнесена к единице массы, то полное количество выделившейся теплоты Q( х) получится интегрированием -хим по расходу массы через данное поперечное сечение канала F (х):
Q(х) = { -химр V • n) dF, (11)
F (х)
где dF — элемент площади сечения; n — единичная нормаль к этому элементу; V — скорость газа, а -хим — количество теплоты в расчете на единицу массы газа, выделившееся вдоль струйки тока, проходящей через dF .
Определения (8) и (9) должны давать практически одинаковые значения Q( х) даже при наличии молекулярной и турбулентной диффузии. В самом деле, основное отличие этих формул связано с диффузионными потоками массы в поперечном направлении, а интегральный диффузионный поток массы в поперечном направлении должен быть равен нулю из-за наличия стенок канала. Небольшое отличие могут вносить диффузионные потоки массы в направлении течения, однако в каналах эти потоки обычно несущественны*. Поэтому для практических приложений можно порекомендовать использовать определение (9), которое, в отличие от (8), не создает при реализации в компьютерной программе трудностей, упомянутых в предыдущем разделе.
В качестве идеального случая можно рассмотреть горение равномерно перемешанной смеси в замкнутом объеме. Начальные массы топлива F и окислителя O равны соответственно GF и G0. Смесь может также содержать инертные вещества (на практике таким веществом можно считать азот, содержащийся в воздухе), масса которых равна Gj. Делается допущение, что химические процессы можно описать как однонаправленную химическую реакцию с массовым уравнением
F + LoO ^ (1 + Lo)P, (12)
где Р — продукт реакции, а Lo — массовый стехиометрический коэффициент, равный массе окислителя, расходуемой при сгорании единицы массы топлива. При протекании реакции (12) до конца максимальная теплота выделяется при коэффициенте избытка окислителя а = Go /(LoGf) = 1 [1]. При этом в расчете на единицу массы топлива в условиях T = To, р = const выделяется теплота:
Ни = hF(To) + Loho(To) - (1 + Lo)hp(To).
Эта величина называется теплотворной способностью топлива при стандартной температуре. В настоящей работе в качестве стандартной температуры используется To = 293.15 К (система Ваничева [18]).
* Течение в канале обычно можно считать тонкослойным: изменение параметров в основном происходит в направлениях, перпендикулярных скорости потока.
Итак, при заданном расходе топлива Ср максимальное количество теплоты, которое может выделиться в канале, примерно равно СрНи . Отношение реально выделившейся на данном участке канала теплоты к этому максимальному значению называют полнотой сгорания топлива п [2]*:
П(х) = х)/(СрНи ) . (13)
Величина Q(х) в (13) может быть определена по формуле (11) с использованием д(1им , ¿им или
дХим. Для стационарного течения самое точное значение полноты сгорания дает использование
тараметр°в ^ и .
В инженерной практике, однако, детальная картина течения в канале обычно недоступна, и приходится обходиться средними по сечениям характеристиками потока, которые могут быть определены экспериментально или путем расчета течения в канале на основе квазиодномерной теории (см. [19]) с привлечением экспериментальных данных или без них. Целесообразно рассмотреть приближенные выражения для полноты сгорания, которые могут быть получены в этом случае, и сопоставить их с описанными выше способами.
Реальное течение в канале заменяется на квазиодномерное течение, в котором все параметры зависят только от продольной координаты х. Параметры квазиодномерного течения в каждом
сечении — температура Т (х), скорость V (х), давление р(х), концентрации компонент газовой
смеси Ук (х) — подбираются так, чтобы продольные потоки импульса, энергии и массы каждой компоненты смеси совпадали с истинными. Если течение в канале стационарно, то баланс энергии в таком квазиодномерном течении имеет следующий вид:
СоН0ач + СрНрнач + сНнач - ^ (х) = (Со + Ср + с )Н(х), (14)
= V 2 р = =
где Н =--УкИк (Т) — энтальпия торможения единицы массы газа, а Jw (х) — теплота, от-
2 к=1
веденная через стенки канала между входным сечением канала и текущим сечением х (если теплота подводится, Jw (х) < 0). Н0ач , Нрнач и Н1нач — начальные энтальпии торможения окислителя, топлива и инертных веществ. Считается, что можно пренебречь вкладом кинетической
~ _ V 2 = = =Яр =
энергии турбулентности к . Представим Н =--+ к0) + АИ , где Ио) = ^ У^И^о ,
2 к=1
= = \ Т=
АИ ¥к (Ик
(Т) - Ик о ) = I СрёТ (здесь интегрирование ведется при фиксированном составе газа).
к=1 То
Величину Ио рассматривают как «химическую энергию» единицы массы газа, а АИ — как тепловую составляющую энтальпии. (Это не вполне точно — см. примечание, которое было сделано выше при определении «энтальпий образования» Ик о .) Вводится также «тепловая энтальпия тор-
* Тепловой эффект реакции слабо зависит от температуры. При изменении температуры Т от 293.15 К до 3ооо К величина (Ир(Т) + ЦоИо(Т) - (1 + ¿о)Ир(Т)) меняется на несколько процентов. Локальные условия течения также
влияют на протекание реакции. В свою очередь, величина тепловыделения влияет на структуру течения в канале. Поэтому определить максимальное количество тепла, которое могло бы выделиться в канале, практически невозможно. Вместо этого используется стандартная величина СрНи . Поэтому термин «полнота сгорания» является неточным. Более точным является термин «(безразмерный) коэффициент тепловыделения», который используется, например, в [3].
можения» АН = V2 / 2 + АЛ . В качестве оценки для теплоты, выделившейся в результате химических реакций, можно взять изменение «химической энергии» смеси:
Q(х) « -А(вИо) = (а¿о (Л0 )о + (Л0)Р + к- (Л0 )1 - (а¿0 +1 + к:)Л0(х)), (15)
где к- = О- /GF . Тогда уравнение энергии (14) можно переписать в виде:
Q(х) - .V (х) « GF ((а ¿0 +1 + к1)АН (х) - а ¿0 АНоач - АН^ач - к1 АН1нач ) .
Если из расчета в рамках квазиодномерной теории определить АН(х), то для полноты
сгорания возможна следующая оценка, которая используется во многих работах (см., например, [2, 20, 21]):
^^ (а ¿0 +1 + к-) АН (х) - а Ц АН 0ач - АН Г - к- АН -нач ^ £( х) - (х)
Ни ^Ни
Очевидный недостаток этой формулы — включение в полноту сгорания эффекта от тепловых потоков через стенки канала. Кроме того, формула (15) не учитывает детали развития процесса вдоль линий тока, которые учитывались при определении Q(х) по формуле (11). Наиболее точное значение полноты сгорания может быть получено по формулам (13) и (11) с использованием теплового эффекта д^им (8) или д(2м (9). Величина П^им), определенная по формулам (13)
и (11) с использованием теплового эффекта ¿им (10), а также величина п(АН), определенная по формуле (16), являются лишь неточными оценками истинной полноты сгорания. Тем не менее, при отсутствии детальной информации о картине течения можно воспользоваться оценкой п(АН).
Величина п(дхим) также может иметь практическое значение. Поскольку Л0 — константы, то определение (10) можно переписать в виде
Е Лк 0?к
у эр ^у эр
?(3) " '
хим
к=Г к=1
" ■ эр "' эр
-£{ а (к 0^ ) = -ЕА(Лк 01к)
= -А
V к=1
= -АЛ0 .
Здесь Аа = а(х) - анач — разность значений величины а в текущем и начальном сечениях струйки тока стационарного течения. Тепловыделение в рассматриваемой струйке тока с площадью сечения ёР (х) будет равно
dQ(iИы(х) = дхим(х)р(х)(((х) • п))х) = дхим(х)ёО = -А^(х)ёО = -А(ёО)(х).
Здесь было учтено, что в стационарном течении расход вдоль струйки тока не зависит от х. Далее производится интегрирование по всем струйкам тока. Поскольку оператор интегрирования (т. е. суммирования значений для разных струек тока) и оператор А коммутативны, получим
QX3и)м(х) = - \ АЙ0(х)ао = -а(|Л)ао)(х) = -а(ол0)(х),
Р ( х )
где введена средняя по расходу «химическая энергия» в данном сечении = 1 Г ^ =
Л0(х) = О ] Л0аО = Е ^кКк . Таким образом, при вычислении тепловыделения (11) с использо-
Р (х) к=1
ванием теплового эффекта д£им получается оценка (15). Поэтому
п( АН) = ^ С ^(х) = ) - СН
СрНи СрНи
Следовательно, (п(АН) -п^™)) — это обезразмеренное на СрНи количество теплоты, ушедшее в стенку между начальным и текущим сечениями канала.
Существуют и другие способы оценки полноты сгорания. В случае, если из эксперимента или расчета получена информация об изменении средней по сечению концентрации какой-либо компоненты смеси, можно ввести «химическую полноту сгорания», основанную на этом изменении. Термин «химическая полнота» был предложен в работе [2о].
Для объяснения понятия «химическая полнота» нужно снова рассмотреть идеальное горение в замкнутом объеме, в котором идет однонаправленная реакция (12). В таком процессе изменения концентраций реагентов А1- (г) = У7 (^) - ^наЧ (^НаЧ — начальная масса 7-й компоненты, У7 (г) — ее масса в момент времени г) однозначно связаны с количеством выделившейся в единице массы смеси теплоты АQ(t): АУр(г) = -АQ(г)/ Ни , АУо(г) = -L0АQ(t)/ Ни, АУр(г) = (1 + L0)АQ(г)/ Ни . Поэтому для такого идеального процесса полнота сгорания
п(г) = АQ(t) / (начНи ) может быть однозначно выражена через изменения концентраций:
) АУоС) ) (17)
Грнач L0YPнач (1 + L0) 7рнач
По аналогии с (17), для течения в канале можно предложить следующие варианты «химической полноты сгорания»:
П(А7р) = -(Ср + Со СС,№(х) = -(1 + aLо +к!)АГр(х), Ср
п№) = -(ср + со ОМ^ х) = -1+0/^ АУо^ (18)
Ц°р Цо
= (Ср + Оо + С,№( х) = НаШк. А?р( х).
К р (1 + 1ц)Ор 1 +1.0 '
Очевидно, в реальном течении из-за диффузии массы и наличия обратных реакций величины п(АГр), п(АГо) и п(АГр) не обязательно должны совпадать.
В следующем разделе все рассмотренные способы оценки полноты сгорания топлива будут сопоставлены на основе расчетов течения в канале.
4. ЗАДАЧА О ГОРЕНИИ В ПЛОСКОМ КАНАЛЕ СО СВЕРХЗВУКОВЫМ ТУРБУЛЕНТНЫМ ТЕЧЕНИЕМ ВЯЗКОГО ГАЗА НА ВХОДЕ
Для демонстрации описанных выше способов оценки тепловыделения рассматривается течение в плоском канале, изображенном на рис. 3. Канал включает: «изолятор» с медленным расширением вдоль продольной оси х; участок сильного расширения канала и участок постоянной ширины. На правом конце канала начинается сверхзвуковое сопло; оно не показано на рисунках. Поток воздуха с М « 2.5, Т « 800 К по оси канала и с толстыми пограничными слоями (50% и 30% ширины канала) втекает с левого конца канала. Три инжектора углеводородного топлива размещаются непосредственно перед правым концом изолятора (их положения показаны на рис. 3). Топливо впрыскивается перпендикулярно потоку (в направлении оси г); при двумерном расчете в плоскости (х; у) инжекция моделируется как локальные источники массы топлива, энергии и параметров турбулентности. Нижняя стенка канала считается теплоизолированной (тепловой поток через стенку qw = 0). Температура верхней стенки определяется по температуре
Рис. 3. Горение в плоском канале со сверхзвуковым течением на входе (а и 1.2):
а — типичное поле температуры Т [К] на этапе дросселирования; б — мгновенное поле Т [К] после стабилизации горения; в — мгновенное поле числа М; г — линии тока осредненного
по времени течения и изолиния и = 0; д — поле теплового эффекта ¿им [МДж/кг]; е — поле
теплового эффекта [МДж/кг]
газа в центрах приграничных ячеек из условия, что температура наружной поверхности стенки (с другой стороны стенки) Тои и 420 К в предположении, что поперек стенки мгновенно достигается постоянный тепловой поток. Детали используемой численной методологии описаны в работе [15]. Для описания горения углеводородного топлива используется кинетическая схема, аналогичная схеме [22]; она включает одну квазиглобальную реакцию разложения углеводородного топлива на СО и Н2О и 11 реакций в системе Н-С-0 с участием N2 в качестве инертной компоненты. Суммарное массовое уравнение реакции горения имеет вид (ср. (12)):
F + ¿002 ^ ¿1Н20 + ¿2СО2 ;
(19)
где ¿1 и ¿2 — массы паров воды и углекислого газа, образующиеся при сгорания единицы массы топлива; ¿1 + ¿2 = 1 + ¿0 .
Рис. 4. Развитие течения в плоском канале. Зависимости коэффициента интегральной полезной силы от времени для различных значений коэффициента
избытка окислителя а
Расчет показывает, что впрыскиваемое топливо не самовоспламеняется. Чтобы получить горение в канале, на его верхнюю стенку перед началом сопла помещен пневмодроссель (рис. 3, а). Он создает струю воздуха с высоким давлением, перпендикулярную стенке. Это приводит к запиранию канала и к образованию псевдоскачка — области течения, где переход от сверхзвукового к дозвуковому течению реализуется в серии косых скачков уплотнения, сильно взаимодействующих с отрывными пограничными слоями [23]. В отрывных зонах возникает горение, и устанавливается структура течения, показанная на рис. 3, а. Эта структура нестационарна; фронт пламени колеблется в продольном направлении. Затем пневмодроссель выключается. Этот метод дросселирования был неоднократно и успешно опробован в экспериментальных исследованиях, проводимых в ЦАГИ [24].
Для режимов течения с интегральным коэффициентом избытка воздуха а»0.6 -1.5 после выключения дросселирования была достигнута стабилизация пламени в канале. Это можно увидеть по графикам коэффициента интегральной полезной продольной силы Я (рис. 4). Для получения силы Я продольная компонента силы, приложенной к твердым стенкам, интегрировалась по всей поверхности канала и сопла. В случаях с а» 1.2 и 1.5 полезная сила, начиная с некоторого момента, становится практически постоянной. Типичная структура течения после стабилизации пламени для а» 1.2 показана на рис. 3, б — г. На этих рисунках видно, что после выключения пневмодросселя псевдоскачок располагается на участке быстрого расширения канала. Видна большая отрывная зона на нижней стенке канала и небольшая — на верхней стенке. Рециркуляционная зона на нижней стенке замыкается лишь на конечном участке расширения канала. Над этой зоной имеется сверхзвуковое течение, которое постепенно тормозится в серии скачков, размытых из-за недостаточного сеточного разрешения. После завершения расширения канала, течение становится полностью дозвуковым по всему сечению. Поток остается дозвуковым вплоть до начала сопла. Течение является почти стационарным, за исключением слабых колебаний в конце расширения канала (они видны на рис. 3, б, в; соответствующие малые колебания полезной силы можно обнаружить на кривой, соответствующей а »1.2, на рис. 4).
Полученные в расчете поля течения были осреднены по времени. По осредненному полю были построены линии тока (рис. 3, г). По изменению параметров вдоль линий тока были получены зависимости ¿им ( X), ^хим ( X) и ^хим
(х). Величины АН(х), А7Р, М02, Л^о , ЛГСС,2 определялись осреднением соответствующих параметров по расходу в каждом поперечном сечении канала. Затем по полученным величинам вычислялись оценки полноты сгорания (тепловые
эффекты ¿им'* подставлялись в формулы (11) и (13), АН(х) — в формулу (16), приращения
концентраций — в формулы, аналогичные (18)). Полученные при этом результаты выглядят следующим образом.
На рис. 3, д показано поле теплового эффекта д^м, а на рис. 3, е — поле теплового эффекта
дЦ^м . В случае невязкого течения поля д^ и д^м должны практически совпадать (см. рис. 1). Как видно из рис. 3, наличие диффузии массы приводит к сильному расхождению полей этих тепловых эффектов. На рис. 2, б представлены поперечные распределения д^ и д^ в сечении, показанном пунктирной линией на рис. 3. Поперечное распределение «истинного тепловыделения» д^им показывает, что у нижней стенки выделилось гораздо больше теплоты, чем у верхней.
При этом распределение величины д£им, которая включает дополнительный перенос тепла из-за диффузионных потоков массы, оказывается практически равномерным. Но, как и было предсказано в предыдущем разделе, осреднение по сечению (по расходу) обоих распределений, показанных на рис. 2, б, дает практически одинаковые результаты (отличие — менее одного процента). Таким образом, в реальности тепловыделение за счет химических реакций очень неравномерно, но диффузия массы интенсивно перераспределяет выделившуюся теплоту поперек канала. Для
анализа среднего по сечениям тепловыделения можно пользоваться величиной д^м, нахождение
которой не связано с такими вычислительными трудностями, как нахождение д^м .
Однако информация о пространственном распределении тепловыделения за счет химических реакций также может иметь практическое значение. Она может быть полезна для активного управления тепловыделением. Но в этом случае важно знать не то, сколько тепла выделилось вдоль линии тока от начального сечения до данной точки, а то, насколько быстро выделяется теплота вдоль линии тока в окрестности данной точки. Поэтому для визуализации областей наиболее интенсивного тепловыделения можно порекомендовать использовать не поля интегрального
тепловыделения д^ (х,у) [Дж/кг], а поля скорости тепловыделения ф(х,у) = ад^ /d¿ [Дж/(кг-м)], где производная берется вдоль линий тока. Из (8) можно получить
^р N
Е Лктк к1Ж1
ф( х, у) = ■к=1-^-. (20)
р | V |
Для построения полей этой величины не требуется вести интегрирование вдоль линий тока; она может определяться по имеющимся полям основных параметров течения после окончания расчета.
Распределения полученных рассмотренными выше способами оценок полноты сгорания представлены на рис. 5. Видно, что графики П^им) и П^им) практически совпадают. Использование д^им дает несколько завышенную оценку полноты сгорания. Наиболее часто используемая в инженерной практике оценка п(АН), основанная на формуле (16), в данной задаче привела к занижению полноты сгорания почти на 10%. Это связано с тем, что в оценке (16) из химического тепловыделения вычитается теплота, отводимая через стенки канала. Как было показано, величину потерь тепла из-за теплообмена на стенках можно определить по сдвигу между кривыми п(дхим) и п(АН). Поскольку нижняя стенка канала была теплоизолированной, то потери тепла в стенку становятся заметными лишь после того, как зона тепловыделения достигает верхней (охлаждаемой) стенки канала, т.е. на участке постоянного сечения.
Несмотря на указанные отличия, графики «термодинамических» полнот сгорания качественно подобны друг другу. Совсем другое поведение обнаруживают оценки «химических» полнот сгорания, основанные на приращениях концентраций реагентов АУк. Это можно объяснить тем, что в вязком потоке объемы газа не являются замкнутыми системами. Из-за молекулярной диффузии (усиленной турбулентностью) объем газа обменивается массой компонент смеси с ок-
Рис. 5. Течение в плоском канале после стабилизации горения (а и 1.2). Продольные распределения полнот сгорания, определенных различными способами: Х1 — вход в расчетную область; Х2 — сечение инжекции топлива; Х3 — начало быстрого расширения канала; Х4 — начало участка с постоянной шириной канала; Х5 — входное сечение сопла; х6 — начало косого среза сопла; Х7 — конец косого среза сопла
ружающей средой. Поэтому по приращению концентрации вещества нельзя судить о степени завершенности реакции. В некоторых областях приращение концентрации может быть больше, чем максимальное значение, предсказываемое уравнением (19). Таким образом, оценки полноты сгорания, основанные на формулах типа (18), не дают верной количественной характеристики процесса горения. Они характеризуют лишь особенности протекания реакций в канале. Так, по графикам ) видно, что в сопле протекают дополнительные химические процессы. Однако они
практически не меняют истинную полноту сгорания п(^хим) и П(3хим).
5. ПРИМЕРЫ ИСПОЛЬЗОВАНИЯ ДАННЫХ О ТЕПЛОВЫДЕЛЕНИИ ПРИ АНАЛИЗЕ ФИЗИКИ ТЕЧЕНИЯ
Выше были рассмотрены режимы горения в плоском канале с а и 1.2 и 1.5 , когда течение после стабилизации горения становится практически стационарным. Для режимов с а«0.6 и 1.0 расчеты показали, что стабилизированное пламя совершает интенсивные незатухающие периодические колебания в продольном направлении с частотой 250 — 300 Гц. Соответствующие колебания полезной силы видны на рис. 4. Физический механизм этих колебаний можно объяснить на основе анализа последовательных полей продольной скорости и для а и 0.6 (рис. 6) [15]. Следует отметить, что каждое поле на рис. 6 содержит изолинию и = 0 ; это позволяет видеть положения точек отрыва и последующего присоединения потока и оценивать размер рециркуляционной зоны.
Основную роль в возникновении колебаний пламени играет большая отрывная зона, расположенная в структуре псевдоскачка у нижней стенки канала. В момент, показанный на рис. 6, а, в рециркуляционную область затекает свежее топливо. Температура здесь порядка 2200 К, время нахождения топлива внутри этой зоны велико, и там начинается тепловыделение. В результате объем отрывной зоны начинает расти (рис. 6, б). Поперечное сечение сверхзвукового потока уменьшается, и скачки сдвигаются вверх по потоку. В результате псевдоскачок целиком сдвигается влево (рис. 6, б, в). Нижний инжектор попадает внутрь отрыва, и свежее топливо начинает перетекать вверх по потоку от инжектора в рециркуляционную зону. Это приводит к последующему сдвигу псевдоскачка вверх по потоку (рис. 6, г). На рис. 6, д можно видеть крайнее левое положение псевдоскачка. Рециркуляционная зона сдвигается вверх по потоку от инжектора, ее
Рис. 6. Колебания пламени в плоском канале при а и 0.6. Последовательные поля продольной скорости и [м/с] и изолинии и = 0
снабжение свежим топливом прекращается, и ее объем начинает уменьшаться. В результате псевдоскачок начинает смещаться обратно вправо. На рис. 6, е видно, что рециркуляционная зона исчезла, поэтому поток везде направлен вправо. Через некоторое время вниз по потоку от инжекторов образуется новый отрыв (рис. 6, ж), так как интенсивность косых скачков в этом сечении канала становится достаточной для отрыва пограничного слоя. Наконец, на рис. 6, з структура течения становится такой же, как и на рис. 6, а. Затем начинается новый период колебаний.
Как показывает рис. 4, амплитуда колебаний полезной силы (непосредственно связанная с амплитудой продольных колебаний пламени) понижается с ростом а. Для объяснения этого
эффекта целесообразно сопоставить поля «истинного тепловыделения» д^м для разных а (рис. 7). На этих рисунках видно, что при а и 0.6 и 1.0 основное тепловыделение сосредоточено внутри зоны возвратного течения у нижней стенки. Однако в режимах с аи 1.2 и 1.5 максимум тепловыделения достигается не в глубине рециркуляционной зоны, а на ее границе. Внутри этой зоны тепловыделение достаточно слабое. По-видимому, для этих режимов, где топливо находится в недостатке, оно успевает выгореть почти полностью в зоне смешения на внешней границе области возвратного течения и практически не проникает внутрь рециркуляционной зоны. Тепловыделение там почти исчезает, колебания объема рециркуляционной зоны прекращаются, и, соответственно, исчезают продольные колебания псевдоскачка. Это доказывает, что колебания
а — а ~ 0.6; б — а » 1.0; в — а » 1.2; г — а » 1.5
а)
3.0 7.5 8.0 'г> Утя- -
л
6) 7.0
л
я)
Рис. 8. Осредненные по времени поля десятичного логарифма
скорости тепловыделения ф , полученные для разных режимов течения в плоском канале: а — а « 0.6; б — а » 1.0; в — а » 1.2; г — а » 1.5
псевдоскачка не связаны с акустическими возмущениями, а вызываются изменением объема рециркуляционной зоны из-за выделения теплоты внутри нее. Чем беднее смесь, тем меньше топлива проникает в рециркуляционную зону и, соответственно, тем меньше колебания псевдоскачка.
На рис. 8 представлены поля логарифма скорости тепловыделения ^ ф (20), соответствующие полям ¿им на рис. 7. На основании рис. 8 также легко прийти к выводу, что при а » 0.6 и 1.0 теплота в основном выделяется внутри рециркуляционной зоны у нижней стенки, а при а» 1.2 и 1.5 тепловыделение сосредотачивается на внешней границе этой зоны. При этом для получения полей ^ ф не требуется строить линии тока и преодолевать вычислительные трудности, связанные с получением полей интегрального тепловыделения ¿им (х, У). Для режимов
с продольными колебаниями пламени построение мгновенных полей ^ ф позволяет еще и проследить динамику тепловыделения в рециркуляционной зоне и таким образом дать более надежное обоснование описанному выше механизму продольных колебаний.
6. КОРРЕЛЯЦИЯ МЕЖДУ ПОЛЕЗНОЙ СИЛОЙ И ВЫДЕЛЕНИЕМ ТЕПЛА
В КАНАЛЕ
При заданной массе топлива Ср максимальное количество теплоты, которое может выделиться при идеальном сгорании в замкнутом объеме, равно
бшах (а) =
Ни шш(1; а)
1 + а Ь0 + к1
На рис. 9 представлена полученная на основе расчетов течения в плоском канале зависимость безразмерного выделения теплоты Q(а) = 2(а)/бшах(1), где Q(а) — количество теплоты, выделившееся между входом в канал и последним вертикальным поперечным сечением сопла (сечение х7 на рис. 5). Q(а) определялось по формуле (11) на базе «истинного тепловыделения»
д^им . Для сравнения на рис. 9 также показана идеальная зависимость Qmаx(а) = Qmax(а)/Qmax(1). Расхождение между этими двумя кривыми связано с неполнотой сгорания. В противоположность идеальной кривой, выделение теплоты в канале монотонно уменьшается с ростом а, так как все меньшая часть топлива сгорает внутри рециркуляционной зоны (где условия течения благоприятны для горения) и все большая часть сгорает в основном потоке (где условия для горения плохие).
На том же рисунке показана зависимость коэффициента полезной силы Я от а, отмасшта-бированная по вертикали для сравнения с Q(а) . Очевидно, что поведения кривых Ск (а) и Q(а) качественно подобны. Это говорит о том, что изменение полезной силы при изменении коэффициента избытка окислителя а существенным образом определяется тем, как при изменении а меняется количество выделившейся теплоты. А это, в свою очередь, определяется истинной полнотой сгорания топлива п(д£им) и пСд^им).
Рис. 9. Зависимости безразмерного теплового эффекта и коэффициента полезной силы от коэффициента избытка воздуха а
ЗАКЛЮЧЕНИЕ
Проведенное исследование показывает, что для анализа пространственного распределения выделения теплоты за счет химических реакций можно порекомендовать пользоваться скоростью тепловыделения (формула (20)). Для определения интегрального количества теплоты, выделившегося между начальным и текущим сечениями канала, рекомендуется использовать тепловой
эффект ¿и., (хотя пространственное распределение ¿ra отличается от истинного поля тепловыделения).
Полнота сгорания топлива, вычисленная по формулам (13) и (11) на основе д^м, дает оценку истинной полноты с точностью ~1%. Классический способ оценки полноты сгорания, определяемый формулой (1б), содержит вклад от тепловых потоков в стенки канала. В рассмотренном примере это привело к ошибке ~10% при оценке полноты сгорания. Оценки полноты сгорания (1S), основанные на приращениях концентраций компонент смеси, непригодны. Они характеризуют только особенности протекания химических реакций в канале.
Автор глубоко признателен E. А. Mещерякову — автору ряда важных работ по расчетному определению характеристик высокоскоростных камер сгорания (см. [20, 25, 2б]). Автор советовался с ним в ходе описанных ниже расчетов горения в плоском канале; E. А. Mещеряков также взял на себя труд ознакомиться с данной статьей и сделал полезные замечания по ее содержанию.
ЛИТЕРАТУРА
1. Щетинков E. С. Физика горения газов. — M.: Наука, 19б5, 740 с.
2. Теория и расчет воздушно-реактивных двигателей / Под ред. С. M. Шляхтенко. — M.: Maшиностроение, 19S7, 56S с.
3. Болгарский А. В. Расчет процессов в камере сгорания и сопле жидкостного ракетного двигателя. — M.: Оборонгиз, 1957, 101 с.
4. Артемов О. А. Прямоточные воздушно-реактивные двигатели (расчет характеристик). — M.: компания «Спутник+», 2006, 374 с.
5. Жданов В. Т. Расчет прямоточных ВРД для больших сверхзвуковых скоростей полета // Труды ЦАГИ. 1966, с. 1 — 40.
6. СекундовА. Н. Вязкие течения, турбулентность и горение. — В кн.: ^UK^M 19S0 — 2000. Научный вклад в создание авиационных двигателей. В 2-х т. Т. 2. — M.: Maшиностроение. — 2000.
7. Орлик E. В., Старов E. В., Шумский В. В. Определение газодинамическим способом полноты выгорания в модели с горением // Физика горения и взрыва. 2004. Т. 40, № 4, с. 23 — 34.
S. Дапилец Л. А. Mегодикa определения полноты сгорания газообразного топлива в газотурбинных установках методом анализа пробы // Научно-технические ведомости СПбГПУ. 2010. Т. 2, вып. 100, с. Ш — m
9. Favre А. Equations des gaz turbulents compressibles I, II. // J. de Méchanique. 1965, V. 3, p. 361 — 390, V. 4, p. 391 — 421.
10. Friedrich R. R. Compressible turbulence // Proceedings of 2nd Space Course on Low Orbit Transportation, TU München. — 1993, V. 1, p. 15.1 — 15.45.
11. Лапин Ю. В., Стрелец M. Х. Внутренние течения газовых смесей. — M.: Наука, Ш9, 36S с.
12. Власепко В. В. О математическом подходе и принципах построения численных методологий для пакета прикладных программ EWT ЦАГИ / В сб.: Практические аспекты решения задач внешней аэродинамики двигателей летательных аппаратов в рамках осред-ненных по времени уравнений Навье — Стокса // Труды ЦАГИ. 2007, вып. 2671, с. 20 — S5.
13. Власепко В. В., Сабельников В. А. Численное моделирование невязких течений с горением водорода за скачками уплотнения и в детонационных волнах // Физика горения и взрыва. 1995. Т. 31, № 3, с. 11S — 133.
14. Власепко В. В. Численное исследование нестационарного распространения горения по каналу со сверхзвуковым течением вязкого газа // Химическая физика. 2011. Т. 30, № 7, с. 1 — 13.
15. Vlasenko V., ShiryaevaA. Numerical simulation of non-stationary propagation of combustion along a duct with supersonic flow of a viscid gas // Proc IMechE, Part G: J. of Aerospace Engineering. 2012. V. 227, N 3, p. 4S0 — 492.
16. Moretti G. A new technique for the numerical analysis of nonequilibrium flows // AIAA J. 1965, V. 3, N 2, p. 223 — 229.
17. Берлянд А. Т., Власенко В. В., Свищев С. В. Стационарные и нестационарные волновые структуры, возникающие при стабилизации детонации над поверхностью сжатия // Физика горения и взрыва. 2001. Т. 37, № 1, с. 94 — 111.
18. Ваничев А. П. Термодинамический расчет горения и истечения в области высоких температур // Технический отчет № 18, Изд. Бюро новой техники, НИИ № 1 МАП СССР. 1947, с. 1 — 26.
19. Бартльме Ф. Газодинамика горения. — М.: Энергоиздат, 1981, 280 с.
20. Мещеряков Е. А. Об эффекте замедления горения сверхзвуковых струй в расширяющихся каналах // В сб.: Моделирование в механике. Т. 1(18), № 6, — Новосибирск. 1987, с. 115 — 120.
21. Кудрявцев А. В., Тарарышкин М. С., Степанов В. А. Использование термодинамического потенциала Гиббса для определения полноты сгорания топлива газодинамическим методом // Ученые записки ЦАГИ. 2011. Т. XLII, № 4, с. 64 — 76.
22. Westbrook Ch. K., Dryer F. L. Chemical kinetic modeling of hydrocarbon combustion // Prog. Energy Combust. Sci. 1984. V. 10, p. 1 — 57.
23. Гуськов О. В., Копченов В. И., Липатов И. И., Острась В. Н., Старухин В. П. Процессы торможения сверхзвуковых течений в каналах. — М.: Физматлит, 2008, 168 стр.
24. Волощенко О. В., Зосимов С. А., Николаев А. А. Экспериментальное исследование процесса горения жидкого углеводородного топлива в плоском канале при сверхзвуковой скорости потока на входе / Модели и методы аэродинамики. Материалы I и II Международных школ-семинаров. — М.: МНЦМО, 2002, с. 75.
25. Мещеряков Е. А., Сабельников В. A. Горение водорода в сверхзвуковом турбулентном потоке в канале при спутной подаче горючего и окислителя // Физика горения и взрыва. 1981. Т. 2, стр. 55 — 64.
26. Мещеряков Е. А, Яшина В. В. Кусочно-одномерная модель псевдоскачка в каналах переменной геометрии. // Ученые записки ЦАГИ. 2013. Т. XLIV, № 5, c. 46 — 63.
Рукопись поступила 21/V 2013 г.
