CC BY
36
________УЧЕНЫЕ ЗАПИСКИ Ц А Г И
Т о м XV 19 84
№ 1
УДК 521.375.826
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОЛЕБАТЕЛЬНОНЕРАВНОВЕСНОГО ДВУХФАЗНОГО ТЕЧЕНИЯ В ГДЛ НА СМЕШЕНИИ
С. Н. Веселов, Ю. Е. Маркачев
Описана процедура численного моделирования стационарного колебательно-неравновесного двухфазного течения в ГДЛ на смешении смеси газов С02— N2—Н20—02 и частиц А1203 в рамках квази-одномерного приближения с привлечением допущения о „мгновенном смешении” компонентов смеси в выделенном сечении сопла. Приведены примеры расчетов. Показано влияние размера и относительного расхода частиц на энергетические характеристики газодинамического лазера (ГДЛ) на смешении. Сделан вывод о необходимости учета влияния двухфазности среды при расчете параметров ГДЛ на смешении, в котором в качестве донорного газа используются продукты сгорания металлов в воздухе.
В настоящее время интенсивно исследуется газодинамический лазер с селективным термическим возбуждением, активная лазерная среда в котором образуется при смешении и последующем расширении газовых потоков существенно разной температуры и состава [1—4], так называемый ГДЛ на смешении.
Для получения больших расходов высокотемпературного донорного газа, необходимого для работы такого ГДЛ, можно использовать продукты сгорания углеводородных топлив в воздухе. Однако в этом случае нагретый азот содержит примеси газов (С02, Н20 и др.), значительно ухудшающих эффективность ГДЛ на смешении и зачастую ставящих под сомнение необходимость дополнительного ввода рабочего газа С02.
Вместе с тем известно, что при сжигании в воздухе ряда металлов (А1, ^Ag, В) можно получить довольно чистый нагретый до высоких температур (3000—5000 К) азот, „загрязненный" лишь мелкими частицами окислов металлов (А1203, А^О, В202).
Цель настоящей статьи состоит в оценке возможности непосредственного использования полученной двухфазной среды (газ— частицы окислов металлов) в качестве донорного газа в ГДЛ на смешении с помощью численного моделирования двухфазного колебательно-неравновесного стационарного течения в сопле ГДЛ
5—«Ученые записки ЦАГИ» № 1
65
на смешении в рамках квазиодномерного приближения с привлечением допущения о „мгновенном смешении“ компонентов газовой смеси в выделенном сечении сопла [5].
Математическая модель. Рассмотрим стационарное квазиодно-мерное двухфазное колебательно-неравновесное течение в сопле Лаваля высокотемпературной газовой смеси Ы2, 02, С02, Н20 с частицами А1203, полученной в результате сгорания алюминиевой пудры в воздухе.
В качестве математической модели такого течения используем следующую систему обыкновенных дифференциальных уравнений, представляющих собой уравнения движения двухжидкостной сплошной среды (колебательно-неравновесной газовой фазы и твердой фазы—частиц шарообразной формы одинакового размера):
(1)
(1х из ри йх
йи* <р1
-±-— (.-«л (2,
й. Го ср2
и кинетические уравнения
^^2 __ ^ ( 2 (*?2 4“ I)2 Г р , еа Ч . 3
<1х и ^ Зе2 (е2 + 2) + 2
\^{е^-е2) + ±- ср32]}; (4)
(1е4 _ 1 с1х а
^с01?з4 + — (в;*-е«)1, (6)
Ч 1
а также вспомогательные алгебраические уравнения, необходимые для замыкания системы:
= (7)
г г2 рв
_ 2го\и — ия\
Ке*=- ^ ; (8) <Р! с„ N4
** = з/-^р-г; о)
Ъ2 = рк„{е,(еа + 1)9 — 4(е3+1)ехр[— (0а — Зв2)/7’]}; (10)
<Р|4=/>*мЫ** 4- 1)ехр[(03 — 04)/Т] — е4(е3+ 1)}; (11)
Р = р/?ш1хГ; (12)
Я = сГТ + £ + -^++ (13)
рйА = С; (14)
РаИяЛ-ИГС. (15)
Обозначения общепринятые. Индексы 5 и eq относятся к параметрам частиц и к равновесным значениям кинетических параметров соответственно, /?mix, с™'х — газовая постоянная и теплоемкость при постоянном давлении смеси газов, Е — суммарная удельная колебательная энергия потока, Sn„ £со,— мольные концентрации N2 и С02 соответственно, г—радиус частиц, рв, Cs — плотность и теплоемкость вещества частиц соответственно, G — расход газовой фазы через сопло с площадью поперечного сечения A, W — относительный расход твердой фазы.
Уравнения движения (1) —(3) и (7)—(9) заимствованы из работы [6]. Там же перечислены исходные допущения и приведены выражения для коэффициента /Ь (формула Карлсона и Хогланда) и числа Нуссельта (формула Дрейка). В дальнейшем считалось, что показатель адиабаты х=1,4, число Прандтля Рг = 0,75, коэффициент вязкости смеси газов ц полагался равным коэффициенту вязкости азота, который рассчитывался по формуле Сазерленда [7].
Кинетические уравнения (4)—(6), (10), (11) представляют собой уравнения баланса среднего числа колебательных квантов е{ колебательных мод молекул С02 и Na, составленные в рамках трехуровневой модели нелинейной колебательной кинетики, nOiipo6HO изложенной в работе [8], где приведены также использованные в настоящих расчетах времена релаксации т2, *4, константы скоростей колебательного энергообмена &32, ku. При записи уравнений колебательной кинетики и вычислении времен и скоростей колебательной релаксации процессы гетерогенного колебательного энергообмена не учитывались, относительно молекулы 02 предполагалось, что она участвует в процессах V— V, V— V', V—Т обмена молекул N2 и С02 лишь как партнер по столкновению, характерные времена элементарных столкновительных процессов с участием молекулы 02 считались равными соответствующим временам для N2.
Процедура численного моделирования. Процедура численного моделирования ГДЛ на смешении в квазиодномерном приближении состоит из следующих этапов: I—Расчет стационарного квазиодно-мерного колебательно-неравновесного двухфазного течения в сопле Лаваля с переходом через скорость звука. II —Расчет параметров „мгновенного смешения1* в выделенном сечении сопла в сверхзвуковой области течения. III—Расчет сверхзвукового двухфазного колебательно-неравновесного течения и определение выходных энергетических характеристик соплового аппарата ГДЛ на смешении.
I. Для расчета двухфазного колебательно-неравновесного стационарного квазиодномерного течения в сопле Лаваля был применен разработанный ранее алгоритм решения прямой задачи теории сопла в квазиодномерном приближении [9]. В настоящей работе использована версия (ЬР -> ЬА, 8G) этого алгоритма, состоящая в том, что по заданному распределению давления РМ(х) по оси сопла и расходу газа через сопло G<1> ищется контур сопла АЮ(х), т. е. по существу решается обратная задача теории сопла в квазиодномерном приближении, затем с помощью разработанного в работе [9] алгоритма уточнения входных данных Р(х), G определяются поправки 8Я(1>(л:) и Ю^(х) с тем, чтобы контур сопла А^(х), найденный в результате решения обратной задачи с уточненными Р*2) (х) = PW (х) + S/W)^), G(2) = G<1) + 8G(1) был более близок к первоначально заданному контуру сопла А (д:).
Проведенные расчеты показали, что если в качестве начального приближения взять результаты расчета для изоэнтропического течения идеального газа в рассматриваемом сопле, то нескольких итераций (решений обратной задачи) оказывается достаточным, чтобы с хорошей точностью (в расчетах, выполненных в настоящей работе считались допустимыми отклонения 18Л (х)/А (х) |-< 10-2, | ЬР (x)jP(x)! •< 10~2, | Ю/G |<!10~2) приблизиться к заданному контуру А(х) даже при наличии существенных возмущений потока, вызванных обменом импульсом и энергией между газом и частицами.
Центр тяжести при выполнении расчета в рамках такого алгоритма падает на решение обратной задачи, где, очевидно, отсутствуют сложности, связанные с прохождением области критического сечения сопла, присущие существующим алгоритмам маршевого счета [10]. Вместе с тем время счета одного варианта задания по данному алгоритму не так велико в сравнении со временем счета по широко применяемым в настоящее время для решения подобных задач алгоритмов установления [11], где, по существу, обратная задача решается сотни раз.
Численное решение обратной задачи в нашем случае обладает хорошо известными трудностями, связанными с необходимостью численного интегрирования системы жестких дифференциальных уравнений, что обусловлено наличием физических процессов, протекающих с существенно разными характерными временами. В настоящее время разработан ряд численных методов, предназначенных специально для работы с жесткими системами обыкновенных дифференциальных уравнений. Одним из наиболее удачных алгоритмов является алгоритм Гира [12]. В настоящей работе использовалась более поздняя программная реализация этого алгоритма— стандартная подпрограмма STIFF [13].
II. Поскольку существует большое разнообразие схем смешения в ГДЛ на смешении (ввод рабочего газа в дозвуковой области течения, в области критического сечения сопла, в сверхзвуковой области) и процесс смешения является существенно двумерным (трехмерным) процессом, подчеркнем, что наш выбор местоположения сечения „мгновенного смешения" и расчет параметров в нем относится к моделированию течения в установке с дозвуковым комбинированным вводом рабочего газа, подробно описанной в работе [4]. Выбор обусловлен сопоставлением результатов расчета в квазиодномерной постановке с результатами эксперимента [4] и численным моделированием в двумерной постановке в приближении „узкого канала“ [14].
В дальнейшем считалось, что в сверхзвуковой области течения (число Маха основного потока М! ~ 3,0) вводится в направлении течения рабочий газ С02 с температурой Т2 и скоростью и2, равными температуре и скорости идеального газа с температурой торможения Т20, показателем адиабаты х = 1,4, изоэнтропически расширенному до числа Маха М2 = М4. Рабочий газ вводится при давлении р2, равном статическому давлению основного потока рх в количест-
v /~>вХ
ве, достаточном для получения заданной концентраци Ссо3 на выходе из сопла.
Удовлетворив законам сохранения массы компонент газовой и твердой фазы, импульса, суммарной и колебательной энергии потоков в сечении мгновенного смешения и предположив, что на разрыве не происходит изменения скорости и температуры частиц, а также увеличив сечение сопла Аг на величину Л2, достаточную 68
для пропускания расхода рабочего газа при параметрах смешения, получим начальные данные для расчета сверхзвукового выходного участка сопла ГДЛ на смешении.
III. Для расчета параметров потока на этом участке обычно применяют маршевый метод счета, заключающийся, в основном, в подборе градиента давления, удовлетворяющего течению в сопле заданного контура [в качестве такого контура теперь принимается А(х) = АХх)’(А-\ + А^!А-1\. Однако в настоящей работе с целью упрощения программы для расчета сверхзвуковой области течения был использован описанный выше алгоритм последовательного решения обратной задачи с уточнением распределения давления по оси сопла.
Результаты численного моделирования. Все комментируемые ниже расчеты выполнены для течения двухфазной среды в плоских соплах минимальной длины с угловой точкой. Сопло I рассчитано на получение на выходе из сопла равномерного потока идеального газа (х = 1,4) с числом М = 4,5, высота критического сечения Л* =1,2 мм. Сопло II рассчитано на М = 5,0, Л* = 0,4 мм. Координаты контура сверхзвуковой части сопл взяты из работы [15]. Прямолинейный участок дозвукового контура наклонен под углом 45° к оси сопла. Начальное сечение в обоих случаях располагалось на расстоянии 20 (Л*/2) вверх по течению от критического сечения сопла, сечение „мгновенного смешения" — на расстоянии 5(/г%/2) вниз по течению от критического сечения. В случае сопла I выходное сечение располагалось на расстоянии 88 (Л*/2), а в случае сопла II — на расстоянии 148 (Л*/2) от критического сечения сопла.
Расчеты, выполненные для течения в сопле I двухфазной среды при температуре в форкамере 7'0= 1500К и давлении Р0= Ю-ЮбПа без учета колебательной неравновесности, показали, что газодинамические параметры двухфазной среды на выходе из сопла практически не зависят от начальной разности скоростей газа (и0) и частиц (и°), поэтому в дальнейшем предполагалось, что и05 — и°. Исследование влияния начальной разности температур газа и частиц на течение не проводилось, всюду считалось, что температура частиц Т0 и газа Т° в начальном сечении совпадают.
Расчеты показали, что основное влияние на газодинамические параметры двухфазной среды в данных условиях оказывают следующие два параметра: относительный расход частиц Ш и радиус частиц г. Поступательная температура газа Тех и температура частиц Тех, давление газа рех, скорость частиц и\х на выходе из сопла существенно зависят от г — при уменьшении г от 100 мк до 1 мк при относительном расходе № = 0,1 Тех увеличивается на— 30%, Тех падает примерно вдвое, рех возрастает на~30%, иIх возрастает Ь ~ 6 раз. Было установлено, что течение двухфазной среды с частицами радиусом г= 100 мк и 1^ = 0,001 можно считать примером практически однофазного течения (частиц нет), типичным примером двухфазной среды в данных условиях можно считать течение с частицами радиусом г = 5 мк с относительным расходом 1^ = 0,1.
Прежде чем приступить к обсуждению влияния двухфазности среды на лазерные характеристики потока в ГДЛ на смешении проведем сравнение результатов расчета по предложенной схеме в случае ^ = 0,001, г=100мк и известными экспериментальными результатами [4]. Расчеты проведены для сопла I. Состав потока на входе в сопло: С02—-N3 — Н20 с весовыми концентрациями
Ссо, = 0,0001, Сн3о = 0,01. В сечение „мгновенного смешения" вводился С02 при температуре торможения Т20 = 300 1500 К в коли-
честве, достаточном для получения на выходе из сопла Ссоа = 0,2, давление в форкамере р0 = 10* 105Па. На рис. 1 приведены графики зависимости коэффициента усиления малого сигнала [8] от температуры Т0 в форкамере, полученные при различных Т20, пунктиром нанесена экспериментальная зависимость [4]. Видно, что если принять во внимание неконтролируемый в эксперименте нагрев
\о
0,5.
1200
то
тот0,к
1—Тх= 300 К, 2—Гао = 500 К, 3—ТЖ = 1000 К, 4-Тт = 1500 К, эксперимент [4]
Нис. 1
С02 в устройстве ввода, то согласие результатов расчета с экспериментом можно считать удовлетворительным. Заметно влияние на go условий „мгновенного смешения" (зависимость g0 от Т20, т. е. фактически от скорости и температуры струи рабочего газа в сечении смешения).
Расчеты, проведенные с учетом двухфазности среды, т. е. при г= 1-г-100 мк, ИР = 0,1 -*-0,25 при различных Т0, Т2й, Ссо,, Сн2о Для двух выбранных сопл также указали на существенную зависимость и располагаемой удельной лазерной мощности Р/й не только от г, Ш, но и от Т20, СсоПоэтому расчеты, иллюстрирующие далее влияние двухфазности среды на характеристики ГДЛ на смешении, проведены для сопла II при близких к оптимальным (с точки зрения получения высоких значений Р/С и приемлемых для эффективной работы резонатора g0) величинах Т0,= 1000 К и Ссо, = 0,15.
Используемая для оценки эффективности ГДЛ величина располагаемой удельной лазерной мощности определялась как
Р/0 = Т)к [ (зГ + 4х) - иг (71) Н- е? (7$) ] .
Здесь *|к — квантовый КПД лазерного перехода с длиной волны 10,6 мк молекулы С02 (т]к = 0,42), ее4х, ез* — удельные колебательные энергии колебательно возбужденного азота и колебательной моды (001) С02 в потоке С02 — N2 — Н20 — 02 на выходе из сопла, Т1 — значение колебательной температуры верхнего лазерного уровня, соответствующее нулевому значению g0, рассчитанному в предположении, что все остальные параметры, необходимые для расчета
gQ, берутся равными соответствующим величинам в выходном сечении сопла.
Начальные концентрации азота С&2, кислорода Со;, температура в форкамере Т0 после завершения процесса сгорания алюминия в воздухе являются функциями относительного расхода твердой
1—г — 1 мк, 2—г = 10 мк, 3—т — 100 мк, частиц нет
Рис. 3
1—г = 1 мк, 2—г = 10 мк, 3—г = 100 мк, 4—частиц нет Рис. 4
фазы (А1203) — параметра На рис. 2 приведены графики зависимостей С^ЛЙР), Со2(ИР), Т0 (Ц7), использованные в настоящих расчетах.
На рис. 3, 4 приведены графики зависимости выходных параметров g0, Р/О от относительного расхода твердой фазы, рассчитанные для случая течения частиц разной величины г= 1, 10,100 мк.
Для сравнения указаны g0 и Р/б, которые получаются в тех же условиях, но в отсутствии частиц (г=100 мк, = 0,001). Расчеты выполнены при /?о=Ю-10б Па, Снао = 0,01.
Видна существенная зависимость эффективности ГДЛ на смешении от размера частиц в потоке. Наиболее сильное ухудшение
выходных характеристик ГДЛ заметно в области высоких температур {Т0 = 2500 3000 К), больших относительных расходах твердой
фазы. Пунктирные кривые на этих фигурах показывают зависимость g0 и Р/G от W при повышенном давлении /?0 = 20-Ю5 Па в случае наличия в потоке частиц с г =10 мк. Видно, что с увеличением давления в форкамере уменьшение g0 и Р/G в области высоких W, Т0 происходит более быстро.
В диапазоне W = 0,1 -н 0,225 (численное моделирование, проведенное в настоящей работе, имело ограничение по температуре в форкамере установки Г0<ЗОООК в соответствии с рекомендациями по использованию примененной модели кинетики [8J в случае r= 1 мк наблюдается наличие оптимального относительного расхода W0Pt как по коэффициенту усиления g0, так и по удельной располагаемой мощности Р/G. Расчеты показали, что аналогичные оптимумы существуют и для других значений г.
Анализ полученных в расчетах колебательных температур верхнего и нижнего лазерного уровня показал, что в основном за уменьшение g0 и Р/G ответственно повышение колебательной температуры нижнего лазерного уровня, вызванное, в свою очередь, увеличением поступательной температуры и давления в результате обмена импульсом и энергией между частицами и газом.
Результаты проведенного численного моделирования двухфазного колебательно-неравновесного течения в ГДЛ на смешении указывают на то, что учет двухфазности среды при создании такого рода устройств необходим. Хотя и наблюдается значительное снижение выходных энергетических характеристик (уменьшение g0 и Р/G приблизительно на 30%) особенно при наличии мелких частиц окислов металла по сравнению с ГДЛ на смешении в качестве до-норного газа, в котором используется газовая смесь того же состава и температуры, но очищенная от частиц окислов, тем не менее возможно достаточно эффективное (P/G—70 кДж/кг, g0 — 0,5 м-1) непосредственное использование продуктов сгорания А1 в воздухе в качестве донорного газа в ГДЛ на смешении.
Авторы выражают благодарность Л. А. Заклязьминскому за постановку задачи и полезные обсуждения, В. Н. Комарову за указание неточности в уравнении (4) модели кинетики [8], использованной в настоящей работе.
ЛИТЕРАТУРА
1. Б а с о в Н. Г., О р а е в с к и й А. Н., Щеглов В. А. Тепловые методы возбуждения лазеров. — ЖТФ, 1967, т. 37, № 2.
2. В г о n f i п В. R., В о е d е k е г L. R., Cheyer J. P. Thermal lazer excitation by mixing in highly convective flow.—Appl. Phys. Lett.,
1970, vol. 16, N 5.
3. К p о ш к о В. Н., Солоухин Р. Н. Оптимальные режимы инверсии при тепловом возбуждении смешением в сверхзвуковом потоке.—ДАН СССР, 1973, т. 211, № 4.
4. Ефимов Б. Г., Заклязьминский Л. А. Экспериментальное исследование влияния условий смешения в сопле Лаваля на коэффициент усиления в сверхзвуковом потоке.—ФГВ, 1979, № 1.
5. Фомин Н. А., Хижняк С. М. Особенности выбора режима смешения С02—N2— Н20 ГДЛ,—ИФЖ, 1982, т. XLII, № I.
6. С т е р н и н Л. Е. Основы газодинамики двухфазных течений в соплах.—М.: Машиностроение, 1974.
7. Таблицы физических величин. Справочник. Под ред. И. К. Кикоина.—М.: Атомизцат, 1976.
8. Лосев С. А. Газодинамические лазеры. — М.: Наука, 1977.
9. Дмитриев Л. М., Маркачев Ю. Е. Расчет неравновесных квазиодномерных течений в соплах. — Численные методы механики сплошной среды, 1979,'т. 10, № 4.
10. Эмануэль Г. Астронавтика и ракетодинамика (Экспресс-инф. ВИНИТИ), 1965, т. 41, 5.
11. Anderson J. D. Time-dependent Analysis of Population Inversions in an Expanding Gas.—Phys. Fluids, 1970, vol. 13, N 9.
12. Gear C. W. Numerical initial value problems in ordinary differential equations.—N. Y., 1971.
13. 3 a x a p о в А. Ю., Турчанинов В. И. S77FF-nporpaMMa для решения жестких систем обыкновенных дифференциальных уравнений,—М.: ИПМ АН СССР, 1977.
14. Ефимов Б. Г., Заклязьминский Л. А., Марка-ч е в Ю. Е.—В сб.: Химическая физика процессов горения и взрыва. Кинетика химических реакций. (Материалы VI Всесоюзного симпозиума по горению и взрыву. 23—26 сент. 1980 г., Алма-Ата), Черноголовка, 1980.
15. Верховский В. П. Численный расчет плоских сверхзвуковых сопл с изломом контура. Таблицы координат сопл на числа М=3^-7.— Труды ЦАГИ, 1975, вып. 1680.
Рукопись поступила 30/ VI1982.
