Доклады Российской академии наук. Физика, технические науки, 2022, T. 502, № 1, стр. 55-62

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

В качестве примера использования разработанной компьютерной модели представлены прогностические расчетные данные по обтеканию затупленного с радиусом R_n = 0.5 см клина конечного размера с углом полураствора θ = 5° под углом атаки α = 8° гиперзвуковым потоком воздуха при числе Маха М = 16 в условиях сильно разреженной атмосферы. Рассмотрение клина конечного размера означает, что исследуется область течения не только в непосредственной близости от обтекаемой поверхности, но и в области ближнего следа.

Для численного моделирования обтекания затупленного клина конечных размеров использовалась компьютерная программа NERAT-2D [1, 2], в которой реализовано численное интегрирование системы уравнений механики вязкого теплопроводного химического реагирующего газа методом установления. На каждом итерационном шаге последовательно интегрировались система уравнений Навье−Стокса, уравнения сохранения массы химических соединений, уравнения сохранения энергии поступательных степеней свободы частиц и колебательного возбуждения двухатомных молекул:

где $\frac{{\text{d}}}{{{\text{d}}t}} = \frac{\partial }{{\partial t}} + {\mathbf{V}} \cdot \vec {\nabla }$; t − время; ${\mathbf{V}} = {\mathbf{i}}u + {\mathbf{j}}v$ − скорость; $u,\;v$ − компоненты скорости вдоль осей декартовой системы координат х, у; $p,\;\rho $ − давление и плотность; Т − температура поступательных степеней свободы (далее: поступательная температура); $\mu ,\;\lambda $ − коэффициенты вязкости и теплопроводности; c_p − удельная теплоемкость газовой смеси при постоянном давлении; ${{c}_{p}} = \sum\limits_i^{{{N}_{s}}} {{{Y}_{i}}{{c}_{{p,i}}}} $; ${{Y}_{i}} = \frac{{{{\rho }_{i}}}}{\rho }$ − относительная массовая концентрация i-й компоненты; ${{c}_{{p,i}}}$,${{h}_{i}}$ − удельная теплоемкость при постоянном давлении и удельная энтальпия i-й компоненты; ${{\rho }_{i}},\;{{{\mathbf{J}}}_{i}}$ − плотность и плотность потока i-й компоненты; ${{{\mathbf{J}}}_{i}} = - \rho {{D}_{i}}\operatorname{grad} {{Y}_{i}}$; D_i − эффективный коэффициент диффузии i-й компоненты; ${{e}_{{{\text{V,}}i{\text{(}}m)}}}$, ${{\rho }_{{i(m)}}},\,\,{{{\mathbf{J}}}_{{i(m)}}}$ – удельная внутренняя энергия m-й колебательной моды, плотность и вектор диффузионного потока i-й молекулы, обладающей этой колебательной модой; ${{N}_{s}}$ − количество компонентов смеси газов; ${{\dot {w}}_{j}},\,\,{{W}_{j}}$ − массовая и мольная скорость образования i-й компоненты: где X_i – мольная концентрация компонент; ${{k}_{{f,n}}}$, ${{k}_{{r,n}}}$ − константы скоростей прямой и обратной n-й реакции; $S_{{f,i}}^{n}$, $S_{{r,i}}^{n}$ − скорости прямой и обратной реакции, каждая из которых записывается в следующем символическом виде: [X_j] − химический символ реагентов и продуктов реакций; ${{N}_{r}}$ − количество химических реакций; ${{a}_{{j,n}}}$, ${{b}_{{j,n}}}$ − стехиометрические коэффициенты; ${{M}_{i}},\,{{M}_{\Sigma }} = {{\left( {\sum\limits_i^{{{N}_{s}}} {{{{{Y}_{i}}} \mathord{\left/ {\vphantom {{{{Y}_{i}}} {{{M}_{i}}}}} \right. \kern-0em} {{{M}_{i}}}}} } \right)}^{{ - 1}}}$ – молекулярный вес i-й компоненты и суммарный молекулярный вес; ${\mathbf{\hat {\Pi }}},\,\,{{\Phi }_{\mu }}$ − тензор напряжений и диссипативная функция выражаются через компоненты скорости, градиент давления и динамический коэффициент вязкости [2].

Соотношение (6) показывает, что кроме интегрирования системы уравнений (1)–(5) на каждой глобальной итерации решалась система уравнений химической кинетики для определения скоростей образования химических компонент и текущего химического состава. Для расчета скорости образования химических компонент рассчитывались константы скоростей каждой из N_r прямой и обратной реакции, которые записывались в виде обобщенной формулы Аррениуса:

где ${{A}_{{f(r),n}}},\;{{n}_{{f(r),n}}},\;{{E}_{{f(r),n}}}$ − аппроксимирующие постоянные прямой (f) и обратной ($r$) реакции. Использованная в данной работке кинетическая модель высокотемпературного воздуха включала 11 компонент (N, O, e^–, N₂, O₂, NO, ${\text{N}}_{2}^{ + }$, ${\text{O}}_{2}^{ + }$, NO⁺, N⁺, O⁺) и 48 химических реакций [3].

Энтальпия и удельная теплоемкость при постоянном давлении каждой из компонент находились с использованием полиномиального представления термодинамического потенциала Гиббса для каждой n-й компоненты [4]. Коэффициенты вязкости и теплопроводности газовых смесей рассчитывались по приближенным комбинаторным формулам [5, 6], а указанные коэффициенты индивидуальных веществ определялись в первом приближении теории Чепмена–Энскога [6]. Подробности компьютерной реализации системы интегрируемых уравнений приведены в книге [2].

Совместно с (1)–(5) используется термическое и калорическое уравнение состояния идеального газа:

где ${{R}_{0}}$ − универсальная газовая постоянная.

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

В рассматриваемом случае учитывается колебательное возбуждение только молекул N₂ и O₂, поэтому далее будем считать, что каждая из упомянутых молекул кроме трех поступательных и двух вращательных степеней свободы обладают по одной колебательной моде движения. Полагается, что вращательное возбуждение молекул находится в равновесии с поступательными степенями свободы, что обосновано в [7]. С учетом того, что колебательное возбуждение молекул может отличаться от равновесного, вводится понятие температуры колебательного возбуждения молекулы i-й молекулы в m-й моде ${{T}_{{V,m(i)}}}$, а удельная внутренняя энергия записывается в виде

где ${{\delta }_{i}}$ = 1 для двухатомных молекул, ${{\delta }_{i}}$ = 0 для атомов. Температура колебательного возбуждения находится из (9). Предполагается использование межъядерного потенциала вида гармонического осциллятора.

Слагаемые ${{Q}_{{\text{V}}}},\,\,{{Q}_{h}},\,\,{{Q}_{d}}$ в правой части уравнения (4) − объемные мощности тепловыделения, обусловленные процессами колебательной релаксации в газовой смеси, химическими реакциями и диффузионными процессами переноса тепла:

где ${{N}_{{\text{V}}}}$ − число колебательных мод (в рассматриваемом случае ${{N}_{{\text{V}}}} = 2$: m = 1 для колебательной энергии N₂, m = 2 для колебательной энергии O₂); ${{h}_{i}},\,\,{{\dot {w}}_{i}}$ − энтальпия и массовая скорость химических превращений для i-го компонента смеси; ${{\dot {e}}_{{{\text{V,}}m}}}$ – источник колебательной энергии в m-й моде; ${{R}_{{i(m)}}} = {{{{R}_{0}}} \mathord{\left/ {\vphantom {{{{R}_{0}}} {{{M}_{{i(m)}}}}}} \right. \kern-0em} {{{M}_{{i(m)}}}}}$, ${{R}_{0}} = 8.314 \times {{10}^{7}}$ эрг/(K ⋅ моль) – универсальная газовая постоянная; ${{\rho }_{{i(m)}}}$, M_i(m), ${{D}_{{i(m)}}},{{{\mathbf{J}}}_{{i(m)}}} = - {{\rho }_{{i(m)}}}{{D}_{{i(m)}}}\operatorname{grad} {{Y}_{i}}$ − плотность, молекулярный вес, эффективный коэффициент диффузии в многокомпонентной газовой смеси, вектор плотности диффузионного потока i-го компонента газовой смеси, обладающего m-й модой колебательного движения; ${{\theta }_{m}}$ – характеристическая колебательная температура (${{\theta }_{{m = 1}}}({{{\text{N}}}_{{\text{2}}}})$ = 3396 К, ${{\theta }_{{m = 2}}}({{{\text{O}}}_{{\text{2}}}})$ = 2275 К).

При определении ${{\dot {e}}_{{{\text{V,}}i{\text{(}}m)}}}$ учитывались процессы колебательно-поступательной (VT) релаксации по приближенной теории Ландау–Теллера [7]:

и изменения удельной колебательной энергии за счет химических реакций (VC-процессы): где $e_{{{\text{V,}}m}}^{0} = {{e}_{{{\text{V,}}m}}}\left( {{{T}_{{\text{V}}}} = T} \right)$ − равновесная удельная энергия колебательного движения в m-й колебательной моде i-го компонента; ${{\tau }_{{i(m)}}}$ – характерное время релаксации m-й колебательной моды, которое определялось по приближенным соотношениям Милликена и Вайта [7] с поправкой Парка [3], ограничивающей величину ${{\tau }_{{i(m)}}}$ снизу: ${{\eta }_{m}}$ – приведенная масса двухатомной молекулы, ${{M}_{{i(m)}}}$ – масса молекулы i-го компонента; ${{N}_{t}}$ – полная концентрация всех молекул при заданных температуре T и давлении р; ${{\sigma }_{{V,i(m)}}}$ – сечение столкновения колебательно-возбужденной частицы, характерное значение которого равно 3 × 10^–17 см². Диффузионный перенос энергии колебательного возбуждения молекул учитывался за счет их диффузии.

Таким образом, сформулированная вычислительная модель включает в себя систему уравнений в частных производных (1)–(5), систему дифференциальных уравнений первого порядка (6) для каждой компоненты газовой смеси, расчетные соотношения (7)–(15) и кинетическую модель [3]. Эта модель позволяет решить поставленную задачу.

Расчет обтекания затупленного клина проводился сквозным по расчетной области способом без выделения ударных волн. На внешних границах, кроме выходного сечения, задавались условия невозмущенного потока. В выходном сечении, где течение всегда было сверхзвуковым, задавались нулевые производные компонент скорости, плотности, давления, концентраций химических компонент и энергии колебательных мод вдоль соответствующих линий тока. На поверхности обтекаемого клина задавались условия рекомбинация заряженных частиц и условия восстановления массовых долей всех компонент до условий в набегающем потоке, а также условия прилипания и радиационно-равновесная температура поверхности в окрестности затупления (в остальной области температура поверхности полагалась постоянной и равной ${{T}_{w}}$ = 300 К).

Расчеты выполнены при следующих параметрах набегающего потока: ${{p}_{\infty }}$ = 798 эрг/см³, ${{\rho }_{\infty }}$ = = 0.103 × 10^–5 г/см³, ${{T}_{\infty }}$ = 271 К, число Маха М = 16, угол атаки α = 8°. Результаты расчетов показаны на рис. 1–3. Выбор значения скорости обтекания обусловлен наличием экспериментальных данных [8, 9], на которые можно ориентироваться при проведении исследовательских расчетов.

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

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

Наглядным свидетельством отсутствия термодинамического равновесия в потоке являются распределения температур поступательных и колебательных степеней свободы, показанные на рис. 2. У наветренной стороны клина формируется относительно тонкий температурный пограничный слой, структура которого плохо различима на указанных рисунках. В окрестности критической линии тока максимальная температура поступательных степеней свободы составляет 8290 К. Температура колебательного возбуждения растет от фронта ударной волны по направлению к поверхности и вблизи внешней пограничного слоя достигает своего максимального значения 3800 К, после чего опять снижается до температуры поверхности, заданной в данных расчетах равной 300 К. По мере отхода от критической линии тока вдоль наветренной и подветренной поверхностей температура поступательных степеней свободы быстро падает. На осевом расстоянии х = 5 см (см. рис. 2) максимальная температура вблизи нижней поверхности составляет 4200 К, а вблизи верхней – 3500 К. На расстоянии 10 см это уже 3300 и 2600 К. Вблизи кормовой кромки клина максимальные поступательные температуры у наветренной и подветренной сторон клина составляют соответственно 2300 и 1700 К. Таким образом, наблюдается закономерное снижение температуры поступательных степеней свободы вдоль всей обтекаемой поверхности, что также видно на рис. 2.

Температура колебательных степеней свободы ведет себя иначе. Как уже отмечалось, вблизи критической линии тока температура колебаний N₂ достигает в сжатом слое максимальной величины 3800 К. Затем, по мере движения вдоль нижней и верхней поверхностей указанная колебательная температура тоже начинает снижаться, но гораздо медленнее: на расстоянии 5 см – 3500 и 2900 К, 10 см – 3000 и 2300 К, а в вблизи кормовой кромки – 2000 и 2300 К. Примечательно, что у наветренной поверхности, где наибольшая плотность и давление, наблюдается весьма эффективная термализация колебательных степеней свободы – колебательная температура приближается к поступательной. Над подветренной поверхностью наблюдается сильное разрежение потока и, как следствие, ухудшаются условия для термализации колебательных степеней свободы, поэтому колебательная температура изменяется заметно медленно, а на значительном расстоянии вдоль поверхности вообще прекращает изменяться, в результате чего колебательная температура начинает превосходить поступательную температуру в значительном объеме возмущенного течением воздуха (сравните поля температур на рис. 2). Обратим внимание на то, что в области ближнего следа, где сходятся потоки от нижней и верхней поверхностей и давление и плотность возрастают, создаются условия для эффективной термализации колебательных степеней свободы. В результате не только подрастает поступательная температура, но и колебательная температура приближается к поступательной – газовый поток приходит в состояние термодинамического равновесия.

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

Анализ экспериментальных данных по ионизации в окрестности гиперзвукового экспериментального аппарата RAMC [8, 10] показал, что наиболее вероятным основным механизмом ионизации в условиях, близких к рассматриваемым в данной работе, является механизм ассоциативной ионизации, который может быть представлен в следующем символическом виде:

т.е. заряженные частицы NO⁺ и e^– возникают в результате распада возбужденных молекул ${\text{NO}}{\kern 1pt} *$, являющихся, в свою очередь, продуктами ассоциации атомов, образовавшихся при диссоциации молекул ${{{\text{N}}}_{{\text{2}}}}{\text{,}}\,\,{{{\text{O}}}_{{\text{2}}}}$. Заметим, что несмотря на убедительные свидетельства в пользу такого механизма, в силу его многостадийности, в литературе имеются большие разногласия в задании определяющих его кинетических постоянных. В данной работе использовалась кинетическая модель [3]. Числовая концентрация электронов, полученная в расчетах, показана на рис. 3. Обратим внимание на очевидные корреляции с распределением газодинамических функций, обсуждавшихся выше. Наибольшая концентрация наблюдается в окрестности критической линии тока, ${{N}_{e}}$ ~ 5.6 × 10¹³ см^–3. Вдоль обтекаемой поверхности электронная концентрация падает по-разному. Вблизи нижней наветренной стороны кормовой кромки ${{N}_{e}}$ = 9.3 × × 10¹⁰ см^–3, а вблизи верхней подветренной стороны – 3.7 × 10¹⁰ см^–3. За кормовой кромкой, в области отрывного течения, концентрация электронов резко падает, а в области следа, где сходятся потоки газа с нижней и верхней поверхностей, снова наблюдается рост ${{N}_{e}}$. Анализ прямой и обратной скоростей реакции ассоциативной ионизации скоростей показывает, что основным источником заряженных частиц является область вблизи критической линии тока. Практически сразу за областью затупления лидирующей кромки клина начинает превалировать область рекомбинации заряженных частиц, которая над нижней и верхней поверхностями клина идет с разной скоростью из-за сильно различающихся условий в потоке.

Отмеченные закономерности ионизации в окрестности затупленного клина конечных размеров остаются подобными в разных условиях обтекания. Наиболее сильными факторами, влияющими на уровень электронных концентраций, являются скорость набегающего потока и такие термодинамические условия в нем, как давление и плотность. При одних и тех же давлении и плотности набегающего потока с уменьшением скорости резко падает интенсивность процессов ионизации. Так, например, при тех же условиях в набегающем потоке, что и выше, но при М = 13 характерными значениями электронных концентраций у поверхности клина и в ближнем следе становятся ${{N}_{e}}$ ~ 3.0 × 10⁹ см^–3, а при М = 10 – ${{N}_{e}}$ ~ ~ 5.0 × 10⁴ см^–3. При увеличении давления в набегающем потоке примерно в 15 раз при числе М = 10 уровень электронной плотности приближается к ${{N}_{e}}$ ~ 10⁷ см^–3.

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

однако эта задача требует более детального исследования.

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