1. На главную
  2. Электронные версии
  3. Известия РАН. Механика жидкости и газа
  4. № 3, 2023

Известия РАН. Механика жидкости и газа, 2023, № 3, стр. 125-136

СТРУКТУРА УДАРНОЙ ВОЛНЫ В КИСЛОРОДЕ

А. И. Ерофеев a, *, С. В. Русаков a, **

a Центральный аэрогидродинамический институт им. проф. Н.Е. Жуковского
Московская обл., Жуковский, Россия

* E-mail: alex.ivan.erofeev@gmail.com
** E-mail: dsmc1@mail.ru

Поступила в редакцию 12.09.2022
После доработки 17.01.2023
Принята к публикации 17.01.2023

Полный текст (PDF)

Аннотация

Представлены результаты численного изучения процессов релаксации в кислороде при высоких температурах. Столкновения частиц (атомов и молекул) описывается методами молекулярной динамики на основе траекторных расчетов в рамках классической механики. Дается описание комплекса программ для расчета релаксационных процессов в смесях высокотемпературных газов с участием внутренних мод, описывающих вращательное и колебательное движение в молекулах, диссоциацию молекул и рекомбинацию атомов в присутствии третьего тела. Процесс релаксации описан на примере изменения параметров в смеси атомарного и молекулярного кислорода с различными начальными температурами поступательных и внутренних мод. Приведены результаты расчетов структуры ударной волны в кислороде с максимальной поступательной температурой на фронте выше 5000–11 000 К. Дано сопоставление с экспериментальными данными.

Ключевые слова: релаксация, диссоциация молекул, рекомбинация атомов, ударная волна

В последние 10–15 лет изучение явлений в высокотемпературных газах с помощью точного решения динамической задачи взаимодействия молекул в газах на основе классических траекторных расчетов (метод молекулярной динамики) вышло на новый уровень благодаря успешному развитию квантово-механических методов расчета поверхностей потенциальной энергии (ППЭ) [1–5]. Несколько вариантов ППЭ построены для взаимодействия атомов азота (N + N2, N2 + N2) [6] и кислорода [5]. С помощью этих ППЭ проведено исследование элементарных процессов релаксации в азоте и кислороде при высоких температурах. Исследования релаксационных и физико-химических процессов в высокотемпературных газах, проводимые с помощью классических траекторных расчетов, используются для верификации приближенных моделей этих процессов [7].

Вместе с тем проведенное в работе [6] сравнение данных по скоростям диссоциации азота в реакции N + N2 → N + N + N для различных ППЭ, показывает, что данные, полученные с полуэмпирической ППЭ [8] – построенной с помощью потенциала LEPS – дают хорошие результаты до температур 10 000° К (т.е. в диапазоне температур, при которых проводились расчеты в [8]). Это говорит о том, что потенциал применения ППЭ, основанных на полуэмпирических моделях, не исчерпан и возможно рассмотрение на их основе сложных физических процессов в высокотемпературных газах, в том числе в смесях газов.

Расчеты структуры ударной волны в молекулярных газах с применением траекторных расчетов в методе прямого статистического моделирования (DSMC) проводятся с 1971 г. В работе [9] молекулы газа рассматривались как жесткие ротаторы, т.е. учитывался только обмен энергией между поступательными и вращательными степенями свободы. Следствием малого количества частиц в расчетной области (800) явилась большая величина статистической погрешности. С той же моделью и много большим числом молекул в расчетной области (100 частиц на ячейку) задача о структуре ударной волны рассматривалась в работе [10] в 1997 г. Кроме большого числа работ по структуре ударной волны, исследуемых методом DSMC с феноменологическими моделями взаимодействия молекул, конечно же продолжается поиск динамических моделей столкновения. Так, в работе [11] для изучения сильной ударной волны в кислороде использована модель импульсного удара молекул, позволяющей адекватно описывать процессы диссоциации молекул при высоких температурах и получить хорошее согласование с экспериментальными данными.

Для изучения возможности применения метода классических траекторий (МКТ) в высокотемпературных течениях разреженного газа в данной работе для описания взаимодействия атомов и молекул применяется подход, в котором используется простой вариант аддитивного потенциала, описывающего парные взаимодействия всех частиц, участвующих в столкновении. Потенциалы парного взаимодействия атомов описываются функцией Морзе с параметрами, определенными по работам [12–14]. Основным вариантом метода молекулярной динамики является прямое моделирование движения частиц в ансамбле, в котором все частицы одновременно взаимодействуют друг с другом. При исследовании течений разреженного газа более эффективными являются траекторные вычисления парных столкновений молекул (или тройных столкновений в случае рекомбинации атомов), поскольку в таком газе большую часть времени молекулы находятся в свободном полете. Преимущество использования МКТ в исследовании течений разреженного газа понятно, поскольку позволяет проводить расчеты без использования многочисленных феноменологических или конкретных динамических моделей столкновения, применимых в ограниченном диапазоне изменения параметров. Трудности также очевидны – нужны ППЭ, адекватно описывающие процессы взаимодействия, и большие вычислительные мощности. Но дело, на наш взгляд, идет в этом направлении и поэтому было решено начать работу, используя для начала простые, но уже опробованные ППЭ. Одной из наиболее простых в вычислительном плане задач является одномерная задача о структуре ударной волны, для которой имеются хорошие экспериментальные данные. Ниже эта задача решается в предположении, что молекулы кислорода находятся в основном электронном состоянии.

1. ОСНОВНЫЕ ПОЛОЖЕНИЯ

1. Рассматриваются процессы релаксации в смеси атомов и молекул, в котором происходят столкновения различных типов:

– столкновение двух молекул с обменом поступательной, вращательной и колебательной энергиями; это столкновение может привести к диссоциации одной или обеих молекул;

– столкновение атома с молекулой, которое также может закончится диссоциацией молекулы;

– столкновение двух атомов:

– упругое;

– столкновение в присутствии третьего тела, в котором может произойти рекомбинация атомов. Третьим телом может быть как атом, так и молекула.

2. Столкновения частиц рассчитываются методом молекулярной динамики – решаются системы дифференциальных уравнений, описывающих траектории частиц в классическом приближении. Уравнения движения (уравнения Гамильтона) записываются в системе координат, связанной с центром масс. Подробный вывод уравнений движения, описание различных алгоритмов даны в [15–18].

3. Используется модель парного взаимодействия всех атомов (единичных или в составе молекул), участвующих в столкновении. Каждое взаимодействие описывается потенциалом Морзе

(1.1)
$U(R) = D\left[ {\exp ( - 2\alpha (R - \operatorname{Re} )) - 2\exp ( - \alpha (R - \operatorname{Re} ))} \right]$
параметры которого (D, α, Re) определяются для каждой пары, а $R$ – расстояние между частицами.

Применение потенциала Морзе обусловлено тем, что он и его модификации довольно часто использовались при изучении процессов релаксации; параметры потенциала определены [12–14] и прошли проверку временем; он входит, как составная часть, в более сложный потенциал LEPS; для ангармонического осциллятора Морзе имеется точное решение уравнения Шредингера и определены колебательные квантовые уровни; прост в реализации.

4. Программы расчета уравнений движения сталкивающихся частиц разработаны для общего случая – все атомы могут быть различными, т.е. пригодны для описания смеси газов.

Пример столкновения молекул приведен на рис. 1. Взаимодействие между атомами в молекуле обозначено пружинками. Расстояния между частицами ${{R}_{{14}}}$ и ${{R}_{{23}}}$ не показаны, чтобы не перегружать изображение.

Рис. 1.

Схема взаимодействия молекул.

5. Разработаны программы, которые независимо решают описанные выше типы столкновений, они объединены также в единый комплекс. Это объединение дает возможность изучать и чисто релаксационные процессы, и может быть включено в программу расчета течений разреженного газа методом Монте-Карло (метод DSMC) как “блок столкновений”.

6. Начальные данные определяются так: поступательные скорости атомов и молекул определяются из максвелловского распределения при заданной температуре газа Ttr0, внутренняя энергия молекулы (вращательная и колебательная) определяются больцмановскими распределениями с температурами Trot0 и Tvib0; начальное расстояние R между центрами масс сталкивающихся частиц определяется из условия |U(RAB0)| ≪ Etr (Etr – начальная поступательная энергия относительного движения частиц).

7. Решение уравнений движения проводится методом Рунге–Кутты (первые три шага) и четырехточечным методом Адамса с постоянным шагом по времени ht вдоль всей траектории. Оба метода являются методами 4-го порядка. Начальный шаг определялся равным ht = 0.0025tхар. От величины ht зависит точность расчета, которая контролировалась выполнением закона сохранения энергии с относительной погрешностью $\left| {\delta H{\text{/}}H} \right| < \varepsilon $ в конце траектории (как правило, величина $\varepsilon \leqslant 0.001$, H – гамильтониан системы). В случае невыполнения этого условия расчет повторялся с теми же начальными данными и величиной шага по времени ht/2 и т.д. Допускалось 6-кратное повторение расчета. Из-за сильного увеличения времени расчета траектории дальнейшее уменьшение шага по времени не проводилось и выбирались новые начальные условия.

8. Для приведения системы уравнений движения к безразмерному виду вводятся характеристические величины, удобные при решении задач молекулярной физики:

(1.2)
${{L}_{{{\text{хар}}}}} = {{10}^{{--8}}}\,{\text{см = }}1,\quad {{t}_{{{\text{хар}}}}} = {{10}^{{--14}}}\,{\text{c}},\quad {{E}_{{{\text{хар}}}}} = 1.66 \times {{10}^{{--12}}}\,{\text{эрг}}$

Они приняты для удобства согласования с начальными данными: массы атомов и молекул задаются в атомных единицах массы (1.66 × 10–24 г), линейные параметры потенциалов – в ангстремах. Характерное время tхар = 10–14 c имеет порядок периода колебаний атомов в молекуле.

9. В качестве состава исходной смеси в программе принята укороченная версия атмосферы Земли, а именно:

– атомы азота N и кислорода О;

– молекулы азота N2, кислорода О2, окиси азота NО;

и задается процентное содержание компонентов. Молекулы окиси азота NО, если не заданы в начальной смеси, могут образовываться только в результате рекомбинации атомов N и О, т.е. прямое, непосредственное из свободных элементов соединение азота с кислородом в результате обменных реакций:

${\text{О}} + {{{\text{N}}}_{2}} \to {\text{NО}} + {\text{N}},\quad {{{\text{N}}}_{{\text{2}}}} + {{{\text{О}}}_{2}} \to 2{\text{NО}}$
в текущей версии пакета программ не предусмотрено.

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

10. Принятые в работе значения параметров потенциалов Морзе (1.1) приведены в табл. 1.

Таблица 1
Система D, эВ α, Å–1 Re, Å
N–N 9.904 2.6895 1.094
O–O 5.21 2.78 1.207
N–O 6.623 2.83 1.15
N–N2, O2, NO 0.007892 2.045 4.0
O–N2, O2, NO 0.007892 2.045 4.0

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

Атом-атомное взаимодействие в случае, если тройного столкновения не происходит, описывается просто как упругое столкновение. Если тройное столкновение состоялось, то взаимодействие атомов, образующих диатом, описывается потенциалом по первым трем строкам табл. 1, а их взаимодействие с третьим телом (атомом или молекулой) по двум нижним строкам табл. 1.

11. Одним из наиболее трудных вопросов в методе классических траекторий является формулирование начальных условий для внутренних степеней свободы молекул. Дело в том, что в классической механике Гамильтониан свободно движущейся двухатомной молекулы имеет вид:

(1.3)
$H = \frac{{\mu {{{\dot {r}}}^{2}}}}{2} + \frac{{{{M}^{2}}}}{{2\mu {{r}^{2}}}} + U(r) = {{E}_{{int}}} = {{E}_{{rot}}} + {{E}_{{vib}}}$

Здесь $\mu = {{m}_{1}}{{m}_{2}}{\text{/}}({{m}_{1}} + {{m}_{2}})$ – приведенная масса атомов в молекуле, r – расстояние между атомами, ${\mathbf{M}} = \left[ {{\mathbf{r}},{\mathbf{p}}} \right]$ – угловой момент, p – обобщенный импульс, $U(r)$ – потенциальная энергия взаимодействия атомов; ${{E}_{{rot}}},\;{{E}_{{vib}}},\;{{E}_{{int}}}$ – вращательная, колебательная и полная внутренняя энергия молекулы.

Угловой момент является интегралом движения и поэтому в уравнении (1.3): ${{M}^{2}} = {\text{const}}$

Запишем гамильтониан в виде:

(1.4)
$H = \frac{\mu }{2}{{\dot {r}}^{2}} + {{U}_{{eff}}}(r),\quad {{U}_{{eff}}}(r) = \frac{{{{M}^{2}}}}{{2{{\mu }}{{r}^{2}}}} + U(r),$

Рассматриваемая задача сводится к движению частицы массы ${{\mu }}$ в сферически симметричном поле с эффективным (центробежным) потенциалом ${{U}_{{eff}}}(r)$. Записывая угловой момент как

(1.5)
${{M}^{2}} = {{\hbar }^{2}}j(j + 1)$
где $\hbar $ – постоянная Планка, j – вращательное квантовое число, можно построить зависимость ${{U}_{{eff}}}(r,j)$, приведенную на рис. 2. В данном случае эта формула используется как удобное определение квадрата момента при построении зависимости ${{U}_{{eff}}}(r,M)$ в виде функции от j – ${{U}_{{eff}}}(r,j)$.

Рис. 2.

Эффективные потенциальные кривые для кислорода: потенциал Морзе + центробежный потенциал.

Из рис. 2 видно, что диссоциация молекулы имеет место тогда, когда ее внутренняя энергия ${{E}_{{\operatorname{int} }}} > {{U}_{{eff}}}({{R}_{{j{\text{max}}}}})$, а ${{R}_{{j\max }}}$ есть положение максимума ${{U}_{{eff}}}(r,j)$ при заданной величине углового момента (1.5). Для параметров потенциала взаимодействия атомов в молекуле О2, приведенных в таблице, максимальное число j, при котором еще существует связанное состояние молекулы ${{j}_{{{\text{max}}}}} = 226.5$. Величина jmax используется в качестве одного из маркеров при анализе состояния молекулы в конце траектории: если j > jmax, то молекула считается диссоциированной.

Из уравнения (1.3) видно, что разделение внутренней энергии молекулы на вращательную и колебательную в общем случае невозможно, так как величина $r$ является переменной. Для модели жесткого ротатора расстояние между атомами

(1.6)
${{r}_{e}} = {\text{const}},\quad {{E}_{{rot}}} = \frac{{{{M}^{2}}}}{{2{{\mu }}r_{e}^{2}}},\quad {{E}_{{vib}}} = 0$

По аналогии с моделью жесткого ротатора за величину мгновенной вращательной энергии ${{E}_{{rot}}}$ можно принять значение, даваемое в каждый момент времени соотношением Erot = ${{M}^{2}}{\text{/}}2{{\mu }}{{r}^{2}}$, в качестве характеристики энергии вращательного движения молекулы взять среднее значение ${{\left\langle {{{E}_{{rot}}}} \right\rangle }_{\tau }}$ за период колебания $\tau $. Поскольку величина ${{E}_{{int}}}$ также является интегралом движения, то характеристическим значением колебательной энергии будет

(1.7)
${{\left\langle {{{E}_{{vib}}}} \right\rangle }_{\tau }} = {{E}_{{int}}} - {{\left\langle {{{E}_{{rot}}}} \right\rangle }_{\tau }}.$

Однако не ясно, будут ли вычисленные таким образом энергии внутренних мод определять равновесные значения соответствующих температур ${{T}_{{rot}}}$ и ${{T}_{{vib}}}$. Сложности определения внутренних энергий молекулы имеются и при квантовомеханическом рассмотрении [19–22]. В [23, 24] при определении конечного состояния молекул для вращательного квантового числа используют соотношение типа (1.5), а колебательное квантовое число n определяют из квазиклассического правила квантования [25]

$\oint {pdq = 2\pi \hbar \left( {n + \frac{1}{2}} \right)} $
где интеграл берется по полному периоду классического движения частицы, а q и p – обобщенные координаты и импульс.

В данной работе принят последовательный классический подход, в котором нет необходимости в разделении внутренней энергии на вращательную и колебательную. Начальные условия при столкновениях определяются знанием полной внутренней энергии и ее углового момента. Эти же величины определяются и в конце траекторного расчета (если, конечно, молекула не диссоциировала). Далее используется 2-х температурное приближение [26, 27]. Полагая, что время вращательной релаксации мало, вращательная температура ансамбля молекул принимается равной поступательной температуре, а колебательная температура определится из соотношения (1.7).

При формировании начальных условий для траекторного расчета применялись два метода:

1) – если выбранная молекула еще не участвовала в столкновении, то ее вращательная ${{E}_{{rot}}}$ и колебательная ${{E}_{{vib}}}$ энергии определялись из больцмановского распределения при начальной температуре ансамбля. Вращательное число j вычислялось из соотношения ${{E}_{{rot}}} = {{B}_{e}}j(j + 1)$. Далее по величинам ${{E}_{{int}}} = {{E}_{{rot}}} + {{E}_{{vib}}}$ и j определялись точки поворота для центробежного потенциала и величины начальных значений обобщенного импульса.

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

Детальное изложение определения начальных данных приведено в работах [15–18].

12. В данной работе полное сечение столкновения молекул определяется с помощью данных по зависимости коэффициента вязкости от температуры. Исходным является выражение для коэффициента вязкости [28]:

(1.8)

Здесь M – молекулярный вес, Т – температура [К]. $T{\kern 1pt} * = kT{\text{/}}\varepsilon $ – приведенная температура; $\sigma $ – диаметр столкновений, Å; $\varepsilon {\text{/}}k$ – параметр потенциальной функции межмолекулярного взаимодействия, К; ${{\Omega }^{{(2.2)\diamondsuit }}}(T{\kern 1pt} *)$ – приведенный интеграл столкновений. Для потенциала Леннарда-Джонса величины ${{\Omega }^{{(2.2)\diamondsuit }}}(T{\kern 1pt} *)$ затабулированы в диапазоне температур $0.3 \leqslant T{\text{/}}T* \leqslant 400$ [28]. Для азота $\sigma = 3.681\,{\AA}$ и $\varepsilon {\text{/}}k = 91.5K$, для кислорода $\sigma = 3.433\,{\AA}$ и $\varepsilon {\text{/}}k = 113K$. Сечения столкновений в траекторных расчетах определялись как:

(1.9)
${{S}_{{cr}}} = \pi R_{{cr}}^{2},\quad {{R}_{{cr}}} = {{A}_{r}}{{R}_{{coll}}} = {{A}_{r}} \times \sigma \times {{[{{\Omega }^{{(2.2)\diamondsuit }}}]}^{{1{\text{/}}2}}}$
где ${{A}_{r}}$ – варьируемый параметр порядка единицы.

Расчетные данные для ${{R}_{{coll}}}$ аппроксимировались приближенными зависимостями от температуры $T$.

Для кислорода применялись следующие аппроксимации для разных диапазонов температур

(1.10)
${{R}_{{coll}}}({\AA}) = \left\{ \begin{gathered} 13{{T}^{{ - 0.23}}},\quad T < 283\,{\text{K}} \hfill \\ 5.42414{{T}^{{ - 0.077554}}},\quad 283 \leqslant T \leqslant 45{\kern 1pt} {\kern 1pt} 000\,{\text{K}} \hfill \\ \end{gathered} \right.$

Для азота

(1.11)
${{R}_{{coll}}}({\AA}) = \left\{ \begin{gathered} 13.02{{T}^{{ - 0.2258}}},\quad T < 230\,{\text{K}} \hfill \\ 6.1652{{T}^{{ - 0.08643}}},\quad 230 \leqslant T \leqslant 35{\kern 1pt} {\kern 1pt} 000\,{\text{K}} \hfill \\ \end{gathered} \right..$

2. РЕЛАКСАЦИЯ ПРИ ВЫСОКОЙ НАЧАЛЬНОЙ ПОСТУПАТЕЛЬНОЙ ТЕМПЕРАТУРЕ МОЛЕКУЛ

В этом разделе рассматриваются релаксационные явления в газах, связанные с поступательным – вращательным – колебательным обменом энергии молекул и атомов при их столкновениях с возможными процессами диссоциации молекул и рекомбинации атомов. Это рассмотрение проводится с ансамблем молекул и атомов размерности D0, т.е. без учета пространственного смещения частиц, что позволяет рассматривать процессы на достаточно больших временных интервалах, добиваясь (если возможно) прихода системы в равновесное состояние. Вероятность столкновения пары частиц принималась равной ${{p}_{{coll}}} \approx g\sigma (g) \approx {{g}^{{0.6}}}$. Сечение столкновений определялось по формулам (1.10) и (1.11). Предполагается, что начальный состав газа молекулы, число которых ${{N}_{{M0}}}$. Начальные температуры: поступательная ${{T}_{{tr0}}}$, вращательная ${{T}_{{rot0}}}$, колебательная ${{T}_{{vib0}}}$ определялись так, что ${{T}_{{tr0}}} = {{T}_{{rot0}}} \gg {{T}_{{vib0}}}$. Временной параметр процесса ${{K}_{{step}}}$ определен так, что один шаг соответствует примерно $\Delta t = 0.5\tau $, где $\tau $– среднее время свободного пробега частиц.

Основной задачей при проведении расчетов являлось выяснение роли рекомбинации атомов в установлении равновесного состояния системы. Вероятность тройного столкновения, т.е. столкновения двух атомов в присутствии третьего тела, которым может быть атом или молекула, определяется как [29]:

${{P}_{3}} = \frac{4}{3}\pi R_{0}^{3} \times {{n}_{3}}$
где ${{R}_{0}}$ – размер зоны взаимодействия двух атомов, возможность рекомбинации которых изучается, ${{n}_{3}}$ – числовая плотность третьих частиц. Величина R0 в общем случае может зависеть от относительной скорости атомов, а значит, от температуры газа. Как будет показано ниже, в процессе релаксации при указанных выше начальных температурах вначале происходит сильное падение поступательной температуры ансамбля, обусловленное диссоциацией молекул. На этом этапе атомов еще немного и их рекомбинация слабо влияет на релаксационный процесс. В дальнейшем падение температуры замедляется, и релаксация происходит при не очень больших изменениях температуры. Это позволяет предположить, что вероятность тройного столкновения можно рассматривать как постоянный параметр задачи и изучать процесс релаксации при его различных значениях. Можно отметить, что такая же картина имеет место за фронтом сильной ударной волны.

На рис. 3 показаны результаты расчета процесса релаксации при начальных температурах, моделирующих условия эксперимента [30–33] сразу за фронтом ударной волны. Расчеты проводились для различных значений вероятности тройного столкновения ${{P}_{3}}$. При ${{P}_{3}} = 1$, когда каждое столкновение атомов сопровождается наличием в зоне их взаимодействия третьего тела (атома или молекулы), вероятность рекомбинации атомов велика. Тогда равновесие в системе достигается достаточно быстро – примерно за 4000 шагов по времени или за $t = 2000\tau $, равновесная температура (осреднение проводилось на интервале ${{K}_{{step}}} = 6001{\kern 1pt} - {\kern 1pt} 20000$ шагов):

$\left\langle {T_{{trM}}^{{eq}}} \right\rangle = 5465 \pm 63\,{\text{K}},\quad \left\langle {T_{{trА}}^{{eq}}} \right\rangle = 5491 \pm 125\,{\text{K}},\quad \left\langle {T_{{\operatorname{int} }}^{{eq}}} \right\rangle = 5582 \pm 57\,{\text{K}}$
Рис. 3.

Процесс релаксации в кислороде при TtrM0 = Trot0 = 10 820 K и Tvib0 = 500 K(a), Tvib0 = 300 K(б–г): (a) P3 = 1, (б) P3 = 0.1, (в) P3 = 0.01, (г) P3 = 0.001.

Как видно, с учетом статистической погрешности, равновесная температура различных мод очень хорошо совпадает, хотя $\left\langle {T_{{\operatorname{int} }}^{{eq}}} \right\rangle $ превышает $\left\langle {T_{{trM}}^{{eq}}} \right\rangle $ примерно на 2%. Следует отметить, что большая погрешность при вычислении $\left\langle {T_{{trА}}^{{eq}}} \right\rangle $ определяется меньшим числом атомов $\left\langle {{{N}_{A}}} \right\rangle = 1153$ ± ± 13 по сравнению с числом молекул $\left\langle {{{N}_{М}}} \right\rangle = 3423 \pm 6.5$ на том же интервале числа шагов.

С уменьшением вероятности тройного столкновения ${{P}_{3}}$ до 0.01, количество рекомбинаций атомов уменьшается, равновесные температуры также становятся меньше, величина релаксационной зоны увеличивается до $t \approx 10{\kern 1pt} {\kern 1pt} 000\tau $ (${{K}_{{step}}} \approx 20{\kern 1pt} {\kern 1pt} 000$) при ${{P}_{3}} = 0.1$ и до $t \approx 50{\kern 1pt} {\kern 1pt} 000\tau $ (${{K}_{{step}}} \approx $ 100 000) при ${{P}_{3}} = 0.01$. Наконец, при ${{P}_{3}} = 0.001$ количество рекомбинаций атомов уменьшилось настолько, что при продолжающейся диссоциации молекул равновесное состояние системы не достигается и при ${{K}_{{step}}} = 250{\kern 1pt} {\kern 1pt} 000$. Температуры различных мод стали меньше 4000 К, а диссоциировало примерно 25% молекул.

3. СТРУКТУРА УДАРНОЙ ВОЛНЫ В КИСЛОРОДЕ

В этом разделе приводятся результаты расчета структуры ударной волны в кислороде для условий, соответствующих серии экспериментов, описанных в работах [30–33]. “В экспериментах на ударной трубе получены профили поглощения света в кислороде в интервале длин волн 200–260 нм в диапазоне температур 4000–10 800 K. С помощью этих данных измерены профили колебательной температуры молекул кислорода за фронтом ударной волны. Метод определения колебательной температуры кислорода основан на сравнении результатов измерений поглощения и детального расчета спектров поглощения кислорода в системе Шумана-Рунге” [30]. Погрешности измерения температур различны на разных участках фронта ударной волны. Погрешность определения колебательной температуры составляет ±20 и ±10% на восходящей и нисходящей ветви ее изменения соответственно. На приводимых ниже рисунках экспериментальные погрешности отмечены вертикальными отрезками прямых линий.

Расчеты проводились для задачи движения поршня в покоящемся молекулярном газе (плотность ${{n}_{\infty }}$, температуры всех мод ${{T}_{\infty }} = 300$ K) со скоростью ${{V}_{\infty }}$ в системе координат, связанной с поршнем. По этой причине выбиралось такое число Маха, при котором максимальные поступательные температуры газа на фронте ударной волны ${{T}_{{tr,\max }}}$ соответствовали значениям, полученным в эксперименте и обозначенным как ${{T}_{0}}$. Параметры газа до ударной волны указаны в подписях к рисункам, а температуры отнесены к температуре газа перед ударной волной. Шаг расчетной сетки (размер ячейки) до ударной волны принимался равным $\Delta {{x}_{1}} = {{\lambda }_{\infty }}{\text{/}}4$, где ${{\lambda }_{\infty }}$ – средняя длина свободного пробега молекул в покоящемся газе. За фронтом ударной волны шаг – $\Delta {{x}_{2}} = {{\lambda }_{\infty }}{\text{/}}10$. При уменьшении шага сетки $\Delta {{x}_{2}}$ в два раза результаты расчетов практически не изменялись. Количество молекул в ячейке до ударной волны принималось равным ${{N}_{\infty }}$ = 20–40 (это изменение N∞ не влияло на результаты расчета). Общее число частиц в расчетной области, как правило, не превышало 500 000. Расчет проводится методом установления. Поэтому шаг по времени $\Delta t = \Delta {{x}_{1}}{\text{/}}{{V}_{\infty }}$. Поскольку расстояние между фронтом ударной волны и поршнем увеличивается с течением времени, то построение профилей плотности и температур проводилось с привязкой к некоторому значению плотности молекул ${{n}_{c}}$ на фронте ударной волны – обычно ${{n}_{c}} = 3{{n}_{\infty }}$ (т.е. близкому к значению срединной плотности на фронте волны) – и усреднением по 100–200 шагам. Поэтому приводимые на графиках данные не представляются гладкими линиями, а имеют статистический разброс. На приведенных ниже рисунках представлены: поступательные температуры молекул ${{T}_{{trM}}}$ и атомов ${{T}_{{trA}}}$; величина ${{T}_{{int}}} = {{E}_{{int}}}{\text{/}}2k$, характеризующая энергию внутренних степеней свободы молекул (1.3); колебательная температура ${{T}_{{vib}}} = 2{{T}_{{int}}} - {{T}_{{trM}}}$, определенная в двухтемпературном приближении, когда вращательная температура ${{T}_{{rot}}} = {{T}_{{trM}}}$; плотности молекул ${{n}_{M}}$ и атомов ${{n}_{A}}$. Температуры и плотности отнесены к ${{T}_{\infty }}$ и ${{n}_{\infty }}$ соответственно.

На рис. 4а приведены результаты расчетов, полученные для случая, когда вероятность тройного столкновения ${{P}_{3}} = 1$ (т.е. когда при каждом столкновении двух атомов в их окрестности имеется третья частица). Видно, что процесс установления квазиравновесной температуры происходит достаточно быстро за фронтом ударной волны – уже на расстоянии $x{\text{/}}{{\lambda }_{\infty }} \cong 25$ все температуры сравнялись и далее от фронта они согласованно изменяются (уменьшаются). Это уменьшение вызвано тем, что за фронтом волны начинается процесс диссоциации молекул (см. изменение плотности атомов ${{n}_{A}}$), для которого характерное время релаксации велико. Как видно, даже на расстояниях $x{\text{/}}{{\lambda }_{\infty }} \cong 225$ равновесное состояние еще не наступило.

Рис. 4.

Структура ударной волны в кислороде. Расчет: (a) М∞ = 11, N∞ = 40, P3 = 1; (б) М∞ = 11.5, N∞ = 20, P3 = 0.01, ${{T}_{{tr,{\text{max}}}}} = 10{\kern 1pt} {\kern 1pt} 740$ K. Эксперимент: □ ⚫ – ${{T}_{0}} = 10{\kern 1pt} {\kern 1pt} 820$ K, Р∞ = 0.8 тор.

На рис. 4б расчет структуры ударной волны проведен при тех же значениях параметров газа, но вероятности тройных столкновений ${{P}_{3}} = 0.01$ – на два порядка меньше, чем в предыдущем расчете. Уменьшение количества рекомбинаций атомов привело к уменьшению всех температур за фронтом волны, т.е. равновесная температура за фронтом стала меньше, чем в предыдущем случае с более высоким уровнем рекомбинации атомов. Как видно из приведенных данных, диссоциация молекул начинается практически сразу на фронте ударной волны, о чем свидетельствует образование атомов ${{n}_{A}}$, причем поступательная температура атомов ${{T}_{{trA}}} \cong {{T}_{{trM}}}$. На этом же рисунке приведены результаты экспериментов, полученных в [31]. Они показывают, что расчетные профили температур не только качественно, но достаточно хорошо и количественно (с учетом погрешности), согласуются с экспериментальными данными, показанными маркерами. Из этого сравнения следует, что вероятность тройных столкновений ${{P}_{3}} \leqslant 0.01$.

Далее на рис. 5 дано сопоставление расчетных данных с экспериментальными, полученными при меньших числах Маха и при меньших поступательных температурах на фронте волны. Из этих данных следует, что при уменьшении температуры релаксационная зона за фронтом ударной волны увеличивается, степень диссоциации молекул уменьшается и при ${{T}_{0}} = 5300$ К атомарный газ присутствует в очень малом количестве. И в этих случаях можно отметить хорошее количественное согласие расчетных и экспериментальных данных.

Рис. 5.

Расчет: (a) М∞ = 9.9, N∞ = 30, P3 = 0.01, ${{T}_{{tr,{\text{max}}}}} = 8290$ К, (б) М∞ = 7.62, N∞ = 40, P3 = 0.01, ${{T}_{{tr,{\text{max}}}}} = 5370$ K. Эксперимент: ${{T}_{0}} = 8620$ К (a), Р∞ = 1 Topp; ${{T}_{0}} = 5300$ K (б), Р∞ = 2 Торр.

ЗАКЛЮЧЕНИЕ

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

Столкновение частиц описывается как парное взаимодействие всех атомов друг с другом на основе классических уравнений движения, в качестве потенциалов парного взаимодействия используется потенциал Морзе. В общем случае все атомы могут быть различными, но в реализованном варианте рассматриваются варианты, в которых участвуют атомы азота N, кислорода О, молекулы N2, О2 и NО. Образование молекулы NО в результате обменных реакций не предусмотрено, но молекула NО может быть образована при рекомбинации атомов N и О.

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

При больших числах Маха (и высоких температурах на фронте) исследована структура ударной волны в кислороде. Полученные данные хорошо согласуются с данными экспериментов в кислороде при температурах газа на фронте волны от 5600 до 10 800 К.

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

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

  1. Jaffe R., Schwenke D.W., Chaban G. Vibrational and Rotational Excitation and Relaxation of Nitrogen from Accurate Theoretical Calculations // AIAA 2008-1208. P. 14.

  2. Jaffe R., Schwenke D.W., Chaban G. Vibrational and Rotational Excitation and Dissociation in N2–N2 Collisions from Accurate Theoretical Calculations // AIAA 2010-4517. P. 13.

  3. Varga Z., Paukku Y., Truhlar D.G. Potential energy surfaces for O + O2 collision // J. Chem. Phys. 2017. V. 147. 154312.

  4. Paukku Y., Varga Z., Truhlar D.G. Potential energy surface of triplet O4 // J. Chem. Phys. 2018. V. 148, 124314.

  5. Grover M.S., Torres E., Schwartzentruber T.E. Direct molecular simulation of internal energy relaxation and dissociation in oxygen // Phys. Fluids. 2019. V. 31. 076107.

  6. Jaffe R.L., Schwenke D.W., Grover M., Valentini P., Schwartzentruber T.E., Venturi S., Panesi M. Comparison of quantum mechanical and empirical potential energy surfaces and computed rate coefficients for N2 dissociation // AIAA 2016-0503. P. 25.

  7. Погосбекян М.Ю., Сергиевская А.Л. Моделирование реакции диссоциации кислорода в термически неравновесных условиях: модели, траекторные расчеты, эксперимент // Химическая физика. 2018. Т. 37. № 4. С. 20–31.

  8. Esposito F., Armenise I., Capitelli M. N–N2 state to state vibrational relaxation and dissociation rates based on quasiclassical calculations // Chem. Phys. 2006. V. 331. № 1. P. 1–8.

  9. Macpherson A.K. Rotational temperature profiles of shock waves in diatomic gases // J. Fluid Mech. 1971. V. 49. № 2. P. 337–351.

  10. Koura K. Monte Carlo direct simulation of rotational relaxation of diatomic molecules using classical trajectory calculations: Nitrogen shock wave // Physics of Fluids. 1997. V. 9. № 11. P. 3543–3549. https://doi.org/10.1063/1.869462

  11. Luo H., Alexeenko A.A., Macheret S.O. Development of an impulsive model of dissociation in direct simulation Monte Carlo // Phys. Fluids. 2019. Vol. 31. 087105.https://doi.org/10.1063/1.5110162

  12. Konowalow D.D., Hirschfelder J.O. Intermolecular potential functions for nonpolar molecules // Phys. Fluids. 1961. V. 4. № 5. P. 629–636.

  13. Konowalow D.D., Hirschfelder J.O. Morse potential parameters for O–O, N–N, and N–O interaction // Phys. Fluids. 1961. V. 4. № 5. P. 637–642.

  14. Гордеев О.А., Калинин А.П., Комов А.Л., Люстерник В.Е., Самуйлов Е.В., Соколова И.А., Фокин Л.Р. Потенциалы взаимодействия, упругие сечения, интегралы столкновений компонентов воздуха для температур до 20 000 К. Обзоры по теплофизическим свойствам веществ. ТФЦ. – М.: ИВТАН. № 5 (55). 1985. 100 с.

  15. Ерофеев А.И., Русаков С.В. Применение классических траекторных расчетов столкновения молекул для вычисления коэффициентов переноса и изучения истечения разреженного газа в вакуум // Ученые записки ЦАГИ. 2020. Т. LI. № 5. С. 13–28.

  16. Karplus M., Porter R.N., Sharma R.D. Exchange reactions with activation energy. I. Simple barrier potential for (H, H2) // J. Chern. Phys. 1965. V.43. № 9. P. 3259–3287.

  17. Lordi J.A., Mates R.E. Rotational relaxation in nonpolar diatomic gases // Phys. Fluids. 1970. V. 13. № 2. P. 291–308. https://doi.org/10.1063/1.1692920

  18. Полак Л.С., Гольденберг М.Я., Левицкий А.А. Вычислительные методы в химической кинетике. М.: Наука, 1984. 280 с.

  19. Ступоченко Е.В., Лосев С.А., Осипов А.И. Релаксационные процессы в ударных вонах. М.: Наука, 1965. 484 с.

  20. Никитин Е.У. Теория элементарных атомно-молекулярных процессов в газах. М.: Химия, 1970. 456 с.

  21. Jaffe R.L. The Calculation of High-Temperature Equilibrium and Nonequilibriunl Specific Heat Data for N2, O2 and NO // AIAA-87-1633.

  22. Capitelli M., Colonna G., Giordano D., Maraffa L., Casavola F., Minelli P., Pagano D., Pietanza L.D., Taccogna F. Tables of Internal Partition Functions and Thermodynamic Properties of High-Temperature Mars-Atmosphere Species from 50 K to 50000 K. ESA STR-246, ESA Publications Division. ESTEC, Noordwijk, The Netherlands, 2005, 267 p.

  23. Jaffe R.L., Schwenker D.W., Panesi M. First principles calculation of heavy particle rate coefficients. Hypersonic nonequilibrium flows: Fundamentals and recent advances / Ed. E. Josynla// AIAA. 2015.

  24. Bender J.D., Valentini P., Nompelis I. et al. An improved potential energy surface and multi-temperature quasiclassical trajectory calculations of N2 + N2 dissociation reactions // J. Chem. Phys. V. 143, 054304 (2015). https://doi.org/10.1063/1.4927571

  25. Ландау Л.Д., Лифшиц Е.М. Квантовая механика. М.: Физматгиз, 1963, 704 с.

  26. Park C., Assessment of a two-temperature model for dissociating and weakly ionizing nitrogen // J. Thermophysics and Heat Transfer. 1988. V. 2. № 1. P. 8–16.

  27. Losev S.A., Makarov V.N., Pogosbekyan M.J., Shatalov O.P., Nikol’sky V.S. Thermochemical nonequilibrium kinetic models in strong shock waves on air. AIAA-94-1990.

  28. Гиршфельдер Дж., Кертисс Ч., Берд Р. Молекулярная теория газов и жидкостей. М.: ИЛ, 1961. 930 с.

  29. Зельдович Я.Б., Райзер Ю.П. Физика ударных волн и высокотемпературных физических явлений. М.: Наука, 1966, 688 с.

  30. Забелинский И.Е., Ибрагимова Л.Б., Шаталов О.П. Измерение колебательной температуры кислорода за фронтом ударной волны в условиях термической и химической неравновесности // Изв. РАН. МЖГ. 2010. № 3. С. 159–167.

  31. Ibraguimova L.B., Sergievskaya A.L., Shatalov O.P. Dissociation Rate Constants for Oxygen at Temperatures up to 11000 K // Fluid Dynamics. 2013. V. 48. № 4. P. 550–555.

  32. Ibraguimova L.B., Sergievskaya A.L., Levashov V.Yu., Shatalov O.P., Tunik Yu.V., Zabelinskii I. E. Investigation of oxygen dissociation and vibrational relaxation at temperatures 4000–10800 K // J. Chem. Phys. 2013. V. 139. 034317. https://doi.org/10.1063/1.4813070

  33. Ибрагимова Л.Б., Левашов В.Ю., Сергиевская А.Л., Шаталов О.П. Моделирование колебательно-диссоциационной кинетики кислорода при температурах 4000–11 000 К // Изв. РАН. МЖГ. 2014. № 1. С. 131–139.

Дополнительные материалы отсутствуют.

Инструменты

предыдущая статья выпуска содержание выпуска

Известия РАН. Механика жидкости и газа

Архивы выпусков Информация о журнале Отправить рукопись в журнал