Известия РАН. Механика жидкости и газа, 2021, № 5, стр. 89-102

АЭРОФИЗИКА ОБТЕКАНИЯ ЗАТУПЛЕННОГО КЛИНА КОНЕЧНЫХ РАЗМЕРОВ

С. Т. Суржиков^a, b, *

^a Институт проблем механики им. А.Ю. Ишлинского РАН
Москва, Россия

^b Всероссийский научно-исследовательский институт автоматики им. Н.Л. Духова
Москва, Россия

^* E-mail: surg@ipmnet.ru

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

DOI: 10.31857/S0568528121050121

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

Аннотация

Представлены результаты численного анализа процессов ионизации воздуха у поверхности затупленного клина с углом полураствора φ = 3° и радиусом затупления ${{R}_{n}}$ = 1.5 см, обтекаемого воздухом в широком диапазоне параметров набегающего потока. Обсуждаются эффекты физико-химической неравновесности при углах атаки набегающего потока α = 0 и 6°. Построены поля температур поступательных степеней свободы и колебательного возбуждения молекул N₂, концентраций атомарных и молекулярных компонент, а также электронных концентраций вблизи поверхности и в ближнем следе.

Ключевые слова: аэрофизика сверхзвукового полета, обтекание затупленного клина, модель неравновесной физико-химической газодинамики

Проблема ионизации воздуха при входе спускаемых космических аппаратов в плотные слои атмосферы Земли привлекла к себе внимание более полувека назад в связи с потерей радиосвязи с ними на участке интенсивного торможения. На всем протяжении развития практической космонавтики, вплоть до настоящих дней, эффект потери радиосвязи со спускаемым аппаратом приходится учитывать в системах обработки летной информации не только в атмосфере Земли [1], но и при исследовании планет солнечной системы [2].

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

В 60–70-х годах прошлого столетия были выполнены эксперименты по измерению концентрации электронов за фронтом ударной волны (ФУВ) [3–5]. Заметим, что вследствие малости времени протекания физико-химических и ионизационных процессов за ФУВ экспериментальная техника проведения таких опытов оказалась весьма сложной. Среди упомянутых выше работ особо следует отметить цикл работ группы В.А. Горелова [6–8], в которой были не только получены остающиеся до настоящего времени уникальные экспериментальные данные, но выполнен глубокий теоретический анализ этих процессов. Главным интегральным результатом этих исследований стало выяснение иерархической последовательности элементарных атомно-молекулярных процессов от диссоциации и ионизации на фоне многообразия химических превращений и релаксационных процессов возбуждения и дезактивации внутренних степеней свободы. Развитие данного научного направления привело к появлению многообразия расчетных моделей: от моделей, основанных на эффективных (“брутто”) реакциях процессов до “поуровневых” и радиационно-столкновительных моделей [9].

При постановке физических экспериментов на аэродинамических стендах и при движении высокоскоростных летательных аппаратов в атомосфере Земли также могут наблюдаться эффекты ионизации. Однако, в отличие от условий входа возвращающихся на Землю сверхорбитальных космических аппаратов (V > 9 км/с), когда кинетической энергии сталкивающихся частиц оказывается достаточно для прямой ударной ионизации частиц газа, при относительно малых гиперзвуковых скоростях (V < 7 км/с) возникает вопрос о первопричине появления в потоке газа значительного количества электронов. В ряде расчетно-теоретических исследований [6, 10] было установлено, что наиболее вероятной причиной первичной ионизации является процесс ассоциативной ионизации, который является, по крайней мере, двухстадийным:

${\text{N}} + {\text{O}} \to {\text{N}}{{{\text{O}}}^{ \bullet }} \to {\text{N}}{{{\text{O}}}^{{\text{ + }}}} + {{{\text{e}}}^{ - }},$

причем атомы N и O появляются в результате диссоциации компонент воздуха – молекулярного азота N₂ и молекулярного кислорода O₂ при их взаимном столкновении. Процесс диссоциации молекул часто обусловлен не просто их столкновением, но сопровождается возбуждением внутренних степеней свободы, в первую очередь – колебаний ядер в молекулах. Последующая ассоциация атомов приводит к образованию электронно-возбужденного оксида азота NO^• с последующим распадом на ион NO⁺ и электрон e^–. Так что налицо сложный многостадийный характер этого процесса.

В расчетно-экспериментальной работе [11] была показана высокая вероятность реализации указанного сценария развития ионизации у поверхности тела, летящего с умеренной гиперзвуковой скоростью. Эксперименты были выполнены на ударной трубе для скоростей ударных волн в диапазоне V = 4.7–6.7 км/c. Была обнаружена весьма высокая интенсивность излучения нагретой пробки воздуха за фронтом ударной волны в области спектра вакуумного ультрафиолета ($\lambda $ ∼ 190–220 нм), что подтверждает факт появления в потоке молекул NO в возбужденных электронных состояниях, обозначаемых в молекулярной спектроскопии ${{A}^{2}}{{\Sigma }^{ + }}$, ${{B}^{{'2}}}\Delta $, ${{B}^{2}}\Pi $, ${{C}^{2}}\Pi $, ${{D}^{2}}{{\Sigma }^{ + }}$ [12].

Расчет процессов физико-химической кинетики за фронтом ударной волны и спектральной излучательной способности воздуха показал хорошее совпадение с наблюдаемой в эксперименте интенсивностью излучения. А использованная кинетическая модель ассоциативный ионизации [13], применявшееся в работах [14, 15], показала высокую вероятность ассоциативной ионизации.

Это является обоснованием выбора модели физико-химической кинетики для изучения процессов ионизации при относительно низких гиперзвуковых скоростях. Однако следует иметь в виду, что использование этой модели при двух скоростях ${{V}_{\infty }}$ ∼ 7.4 км/с [14] и 5.6 км/с [16] указало на необходимость модификации констант скоростей ассоциативной ионизации при уменьшении скорости.

Кинетическая модель, положенная в основу изучения процессов ионизации у поверхности затупленного клина, при обтекании его потоками воздуха при относительно низких гиперзвуковых скоростях, основана на процессе ассоциативной ионизации, которая включает в себя процессы диссоциации молекул N₂ и О₂, образование молекул NO в результате химических реакций, в первую очередь, из-за столкновительной ассоциации атомов N и О, возбуждения электронных состояний NO^• c последующим распадом NO^• на NO⁺ и e^–.

Эта кинетическая модель используется в данной работе для анализа ионизации газа в окрестности затупленного клина в широком диапазоне параметров набегающего потока. В первой части работы расчеты выполнены для клиньев с углом полураствора φ = 3° и 9° для двух радиусов затупления ${{R}_{n}}$ = 1.5 и 2.0 см при их обтекании под нулевым углом атаки. Во второй части работы обсуждается особенность обтекания затупленного клина под углом атаки.

1. РАСЧЕТНАЯ МОДЕЛЬ

Двумерная задача ламинарного обтекания затупленного клина потоком вязкого, теплопроводного, сжимаемого химически и термически неравновесного воздуха без выделения фронта ударной волны формулируется на основе уравнений неразрывности, Навье–Стокса, сохранения энергии поступательного движения частиц в форме уравнения теплопроводности Фурье–Кирхгоффа, сохранения массы химических компонент (уравнений диффузии) и сохранения колебательной энергии в отдельных колебательных модах двухатомных молекул

(1)

$\frac{{\partial \rho }}{{\partial t}} + {\text{div}}\left( {\rho {\mathbf{V}}} \right) = 0$

(2)

$\frac{{\partial \rho u}}{{\partial t}} + {\text{div}}\left( {\rho u{\mathbf{V}}} \right) = - \frac{{\partial p}}{{\partial x}} - \frac{2}{3}\frac{\partial }{{\partial x}}\left( {\mu \,{\text{div}}{\mathbf{V}}} \right) + \frac{\partial }{{\partial y}}\left[ {\mu \left( {\frac{{\partial {v}}}{{\partial x}} + \frac{{\partial u}}{{\partial y}}} \right)} \right] + 2\frac{\partial }{{\partial x}}\left( {\mu \frac{{\partial u}}{{\partial x}}} \right)$

(3)

$\frac{{\partial \rho v}}{{\partial t}} + {\text{div}}\left( {\rho {v}{\mathbf{V}}} \right) = - \frac{{\partial p}}{{\partial x}} - \frac{2}{3}\frac{\partial }{{\partial x}}\left( {\mu \,{\text{div}}{\mathbf{V}}} \right) + \frac{\partial }{{\partial x}}\left[ {\mu \left( {\frac{{\partial {v}}}{{\partial x}} + \frac{{\partial u}}{{\partial y}}} \right)} \right] + 2\frac{\partial }{{\partial y}}\left( {\mu \frac{{\partial {v}}}{{\partial y}}} \right)$

(4)

$\rho {{c}_{p}}\frac{{\partial T}}{{\partial t}} + \rho {{c}_{p}}{\mathbf{V}}\operatorname{grad} T = \operatorname{div} \left( {\lambda \operatorname{grad} T} \right) + {\mathbf{V}}\operatorname{grad} p + {{\Phi }_{\mu }} + {{Q}_{V}} - {{Q}_{h}} + {{Q}_{d}}$

(5)

$\frac{{\partial {{\rho }_{i}}}}{{\partial t}} + \operatorname{div} {{\rho }_{i}}{\mathbf{V}} = - \operatorname{div} {{{\mathbf{J}}}_{i}} + {{\dot {w}}_{i}},\quad i = 1,\;2,\; \ldots ,\;{{N}_{s}}$

(6)

$\frac{{\partial {{\rho }_{{i(m)}}}{{e}_{{{\text{V,}}i{\text{(}}m)}}}}}{{\partial t}} + \operatorname{div} ({{\rho }_{{i(m)}}}{{e}_{{{\text{V,}}i{\text{(}}m)}}}{\mathbf{V}}) = - \operatorname{div} ({{e}_{{{\text{V,}}i{\text{(}}m)}}}{{{\mathbf{J}}}_{{i(m)}}}) + {{\dot {e}}_{{{\text{V,}}i{\text{(}}m)}}},\quad m = 1,2, \ldots ,{{N}_{{\text{V}}}}$

где

${{\Phi }_{\mu }} = \mu \left[ {2{{{\left( {\frac{{\partial u}}{{\partial x}}} \right)}}^{2}}} \right.\left. { + 2{{{\left( {\frac{{\partial {v}}}{{\partial y}}} \right)}}^{2}} + {{{\left( {\frac{{\partial {v}}}{{\partial x}} + \frac{{\partial u}}{{\partial y}}} \right)}}^{2}} - \frac{2}{3}{{{\left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial {v}}}{{\partial y}}} \right)}}^{2}}} \right],$

диссипативная функция; t – время; x, y – ортогональные декартовы координаты; u, ${v}$ – проекции вектора скорости V на оси координат x и y; p, $\rho $, $T$ – давление, плотность и температура поступательного движения частиц; $\mu $, $\lambda $ – динамический коэффициент вязкости и коэффициент теплопроводности; ${{Y}_{i}}$ – массовая доля i-го компонента смеси газов общим числом ${{N}_{s}}$; c_p, ${{c}_{{p,i}}}$ – удельная теплоемкость при постоянном давлении смеси газов и отдельных компонент; c_p = $\sum\limits_i^{{{N}_{s}}} {{{Y}_{i}}{{c}_{{p,i}}}} $; ${{{\mathbf{J}}}_{i}} = - \rho {{D}_{i}}\operatorname{grad} {{Y}_{i}}$; ${{{\mathbf{J}}}_{i}}$ – плотность диффузионного потока и D_i – эффективный коэффициент диффузии i-го компонента.

Важной особенностью рассматриваемой задачи является учет неравновесного возбуждения внутренних степеней свободы двухатомных молекул. Уравнение (6) в общем виде описывает процессы колебательного возбуждения молекулярных компонент газовой смеси. В рассматриваемом случае учитывается колебательное возбуждение только молекул N₂ и O₂, поэтому далее будем считать, что каждая из упомянутых молекул кроме 3 поступательных и 2 вращательных степеней свободы обладает по одной колебательной моде движения (одной формой колебательного движения). Заметим, что для трехатомных молекул одна молекула может обладать несколькими колебательными степенями свободы (модами). Например, в молекуле CO₂ часто учитываются три колебательных моды [4, 9] молекул. Считается, что вращательное возбуждение молекул находится в равновесии с поступательными степенями свободы, что обосновано в [4, 17]. С учетом того, что колебательное возбуждение молекул может отличаться от равновесного, вводится понятие температуры колебательного возбуждения молекулы i-й молекулы в m-й моде ${{T}_{{V,m(i)}}}$, а удельная внутренняя энергия записывается в виде

(7)

${{e}_{i}} = \frac{3}{2}{{R}_{i}}T + {{\delta }_{i}}{{R}_{i}}T + {{\delta }_{i}}{{R}_{i}}\frac{{{{\theta }_{{i(m)}}}}}{{\exp \left( { - \frac{{{{\theta }_{{i(m)}}}}}{{{{T}_{{V,i(m)}}}}}} \right) - 1}},\quad {{R}_{i}} = \frac{{{{R}_{0}}}}{{{{M}_{i}}}},$

где ${{\delta }_{i}}$ = 1 для двухатомных молекул, ${{\delta }_{i}}$ = 0 для атомов.

Уравнение (6) позволяет определить удельную внутреннюю энергию m-й колебательной моды i-й молекулы

(8)

${{e}_{{{\text{V}},i(m)}}} = \frac{{{{R}_{{i(m)}}}{{\theta }_{{i(m)}}}}}{{\exp ({{{{\theta }_{{i(m)}}}} \mathord{\left/ {\vphantom {{{{\theta }_{{i(m)}}}} {{{T}_{{{\text{V}},i(m)}}}}}} \right. \kern-0em} {{{T}_{{{\text{V}},i(m)}}}}}) - 1}},$

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

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

${{Q}_{{\text{V}}}} = - \sum\limits_{m = 1}^{{{N}_{{\text{V}}}}} {{{{\dot {e}}}_{{{\text{V,}}i{\text{(}}m)}}}} ,\quad {{Q}_{h}} = \sum\limits_{i = 1}^{{{N}_{s}}} {{{h}_{i}}{{{\dot {w}}}_{i}}} ,$

${{Q}_{d}} = \sum\limits_{i = 1}^{{{N}_{s}}} {\rho {{c}_{{p,i}}}{{D}_{i}}} \left( {\operatorname{grad} {{Y}_{i}} \times \operatorname{grad} T} \right),$

где ${{N}_{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-й модой колебательного движения; θ_m – характеристическая колебательная температура (${{\theta }_{{m = 1}}}({{{\text{N}}}_{{\text{2}}}})$ = 3396 K, ${{\theta }_{{m = 2}}}({{{\text{O}}}_{{\text{2}}}})$ = 2275 K).

При определении ${{\dot {e}}_{{{\text{V,}}i{\text{(}}m)}}}$ учитывались процессы колебательно-поступательной (VT) релаксации и изменения удельной колебательной энергии за счет химических реакций (VC-процессы)

(9)

${{\dot {e}}_{{{\text{V,}}m}}} = {{\dot {e}}_{{{\text{V,}}m}}}({\text{VT}}) + {{\dot {e}}_{{{\text{V,}}m}}}({\text{VC}}).$

Релаксационное изменение энергии в каждой колебательной моде рассчитывалось по приближенной теории Ландау–Теллера [13]:

(10)

${{\dot {e}}_{{{\text{V,}}i{\text{(}}m)}}}({\text{VT}}) = {{\rho }_{{i(m)}}}\frac{{e_{{{\text{V,}}i{\text{(}}m)}}^{0} - {{e}_{{{\text{V,}}i{\text{(}}m)}}}}}{{{{\tau }_{{i{\text{(}}m)}}}}}$

$e_{{{\text{V,}}m}}^{0} = {{e}_{{{\text{V,}}m}}}\left( {{{T}_{{\text{V}}}} = T} \right)$ – равновесная удельная энергия колебательного движения в m-й колебательной моде i-го компонента; ${{\tau }_{{i(m)}}}$ – характерное время релаксации m-й колебательной моды, которое определялось по приближенным соотношениям Милликена и Вайта [18] с поправкой Парка [13], ограничивающей величину ${{\tau }_{{i(m)}}}$ снизу

(11)

${{\tau }_{m}} = {{\tau }_{{{\text{VT,}}i{\text{(}}m)}}} + \frac{1}{{{{N}_{t}}{{\sigma }_{{V,i(m)}}}\sqrt {{{8kT} \mathord{\left/ {\vphantom {{8kT} {\left( {\pi {{M}_{{i(m)}}}} \right)}}} \right. \kern-0em} {\left( {\pi {{M}_{{i(m)}}}} \right)}}} }},\quad {{\sigma }_{{V,i(m)}}} = \sigma {{_{{V,i(m)}}^{'}}_{{}}}{{\left( {{{50{\kern 1pt} {\kern 1pt} 000} \mathord{\left/ {\vphantom {{50{\kern 1pt} {\kern 1pt} 000} T}} \right. \kern-0em} T}} \right)}^{2}},$

(12)

$p{{\tau }_{{{\text{VT,}}i(m)}}} = \exp [{{A}_{{VT,i(m)}}}({{T}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 3}} \right. \kern-0em} 3}}}} - {{B}_{{VT,i(m)}}}) - 18.42],$

(13)

${{A}_{{{\text{VT,}}i(m)}}} = 0.00116\eta _{{i(m)}}^{{0.5}}\theta _{m}^{{1.333}},\quad {{B}_{{{\text{VT,}}i(m)}}} = 0.015\eta _{m}^{{0.25}},$

${{\eta }_{m}}$ – приведенная масса двухатомной молекулы, ${{M}_{{i(m)}}}$ – масса молекулы i-го компонента; N_t – полная концентрация всех молекул при заданных температуре T и давлении р; ${{\sigma }_{{V,i(m)}}}$ – сечение столкновения колебательно-возбужденной частицы, характерное значение которого равно 3 × × 10^–17 см².

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

(14)

${{e}_{{V,m}}}({\text{VC}}) = {{e}_{{V,m}}}\frac{1}{2}({{\dot {w}}_{{i(m)}}} - \left| {{{{\dot {w}}}_{{i(m)}}}} \right|),$

где предполагалось, что уменьшение колебательной энергии в m-й моде в 1 см³ за 1 с пропорционально объемной скорости исчезновения молекул, имеющих эту колебательную моду.

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

(15)

${{e}_{{{\text{V,}}m}}}{{{\mathbf{J}}}_{{{\text{V,}}i(m)}}} = - {{e}_{{{\text{V,}}m}}}{{\rho }_{{i(m)}}}{{D}_{{{\text{V,}}i(m)}}}{\text{grad}}{{Y}_{i}},$

где ${{D}_{{{\text{V,}}i(m)}}}$ – эффективный коэффициент диффузии энергии колебательного возбуждения, который принимался равным определенному выше эффективному коэффициенту диффузии.

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

Для этих целей используются соотношения, получаемые в первом приближении теории Чепмена–Энскога [19], а также приближенные комбинаторные соотношения Манна, Брокау и Уилки [20] (соответствующие формулы приведены в [9]).

Массовая скорость химических превращений определялась при решении системы кинетических уравнений

(16)

${{\left( {\frac{{{\text{d}}{{X}_{i}}}}{{{\text{d}}t}}} \right)}_{n}} = {{k}_{{f,n}}}\left( {{{b}_{{i,n}}} - {{a}_{{i,n}}}} \right)\prod\limits_j^{{{N}_{s}}} {X_{j}^{{{{a}_{{j,n}}}}}} - {{k}_{{r,n}}}\left( {{{b}_{{i,n}}} - {{a}_{{i,n}}}} \right)\prod\limits_j^{{{N}_{s}}} {X_{j}^{{{{b}_{{j,n}}}}}} = \left( {{{b}_{{i,n}}} - {{a}_{{i,n}}}} \right)\left( {{{S}_{{f,n}}} - {{S}_{{r,n}}}} \right),$

где X_i – объемно-мольная концентрация i-й компоненты; N_r – число химических реакций; ${{S}_{{f,n}}}$, $S{{}_{{r,n}}}$ – скорости прямой и обратной реакции; ${{k}_{{f,n}}}$, ${{k}_{{r,n}}}$ – константы скоростей прямых и обратных реакций; ${{a}_{{i,n}}}$, ${{b}_{{i,n}}}$ – стехиометрические коэффициенты n-й химической реакции, определяемые с использованием канонической формы уравнений химической кинетики:

$\sum\limits_{j = 1}^{{{N}_{s}}} {{{a}_{{j,n}}}\left[ {{{X}_{j}}} \right]} = \sum\limits_{j = 1}^{{{N}_{s}}} {{{b}_{{j,n}}}\left[ {{{X}_{j}}} \right]} ,\quad n = 1,\;2,\; \ldots ,\;{{N}_{r}},$

где [X_j] – химические символы реагентов и продуктов химических реакций. Массовая скорость образования i-й компоненты в единице объема определяется по формуле

(17)

${{\dot {w}}_{i}} = {{M}_{i}}\sum\limits_{n = 1}^{{{N}_{r}}} {\left( {{{b}_{{i,n}}} - {{a}_{{i,n}}}} \right)\left( {{{S}_{{f,n}}} - {{S}_{{r,n}}}} \right)} ,\quad {\text{г/}}({\text{с}}{{{\text{м}}}^{3}} \cdot {\text{с}})~$

Константы скоростей прямой (f) и обратной (r) реакций определяются обобщенной аррениусовской зависимостью

(18)

${{k}_{{f(r),n}}} = {{A}_{{f(r),n}}}{{T}^{{{{n}_{{f(r),n}}}}}}\exp \left( { - \frac{{{{E}_{{f(r),n}}}}}{{kT}}} \right)$

где ${{A}_{{f(r),n}}}$, ${{n}_{{f(r),n}}}$, ${{E}_{{f(r),n}}}$ – аппроксимирующие коэффициенты; k – постоянная Больцмана; n – номер химической реакции в кинетической модели. Константы аппроксимации скоростей прямых реакций заимствовались из [13]. Константы скоростей обратных реакций рассчитывались с использованием констант равновесия [20], аппроксимированных в форме обобщенного закона Аррениуса

(19)

${{K}_{n}} = {{A}_{n}}{{T}^{{{{n}_{n}}}}}\exp \left( { - \frac{{{{E}_{n}}}}{{kT}}} \right).$

Аппроксимации констант равновесия в форме (19) тестировались в [9] сравнением с данными Парка [13], основанными на данных JANAF [22]. В [9] показано хорошее согласие констант равновесия при температурах Т > 1000 K.

Малость концентраций заряженных частиц позволила использовать условие квазинейтральности для нахождения мольных концентраций электронов. Результаты исследования степени влияния констант скоростей ассоциативной ионизации, которое оказалось не малым, обсуждались в [15].

Используемые в расчетах константы скоростей диссоциации ${{k}_{f}} = {{k}_{D}}$ подвергались модификации с целью учесть термическую неравновесность молекул. Исследование разных моделей неравновесной диссоциации [15] показало на заметное их влияние на результаты расчетов интенсивности ионизационных процессов. Однако в данной работе применялась лишь эвристическая модель Парка [13]. В этой модели температура ${{T}_{a}}$, используемая в (19), рассчитывалась для каждой из колебательных мод по формуле

(20)

${{T}_{{a,i(m)}}} = \sqrt {T{{T}_{{V,i(m)}}}} ,$

где T – поступательная температура, ${{T}_{{V,i(m)}}}$ – колебательная температура в m-й колебательной моде.

Таким образом, используемая расчетная модель, включающая систему уравнений в частных производных (1)–(6), систему дифференциальных уравнений первого порядка (16) для каждой компоненты газовой смеси и расчетные соотношения (7)–(15) и (17)–(20), позволяют решить задачу обтекания затупленного клина конечного размера вязким, теплопроводным и химически реагирующим газом с учетом ионизации и неравновесной диссоциации. Задача обтекания затупленного клина решалась без выделения ударных волн и каких-либо разрывов газодинамических параметров. На внешних границах расчетной области, кроме выходного сечения, задавались параметры невозмущенного потока. В выходном сечении, где течение также всегда было сверхзвуковым, задавались нулевые производные компонент скорости, плотности, давления, концентраций химических компонент и энергии колебательных мод вдоль соответствующих линий тока. Считалось, что на поверхности происходит рекомбинация заряженных частиц и массовые доли всех компонент восстанавливаются до условий в набегающем потоке. Задавались условия прилипания и радиационно-равновесная температура поверхности в окрестности затупления (в остальной области температура поверхности полагалась постоянной и равной ${{T}_{w}}$ = 297 K).

В расчетах использовался авторский компьютерный код NERAT-2D [9]. Применялась многосеточная многоблочная расчетная технология. Схема решаемой задачи с четырьмя расчетными блоками показана на рис. 1. Здесь представлена промежуточная расчетная сетка. Окончательное решение получалось после двукратного удвоения числа узлов.

Рис. 1.

Схема решаемой задачи, расчетная сетка и четырехблочная область (номера блоков указаны на верхней части рисунка) в задаче об ионизации воздуха при обтекании затупленного клина.

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

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

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

В экспериментальных исследованиях закономерностей ионизации в ударных волнах [3–8, 11] было установлено, что по мере уменьшения давления в набегающем потоке газа состояние газа за ударной волной все более становится неравновесным. Несомненно также определяющее влияние скорости набегающего потока на ионизацию газового потока. Хорошо известно, что в неравновесной физико-химической газовой динамике и аэрофизике ионизационных процессов существует множество комбинаций исходных данных, кардинальным образом изменяющих картину течения. Поэтому важно отметить основные тенденции и закономерности, например, при обтекании неравновесным газом удлиненных аэродинамических конфигураций простейшей конфигурации. В данной работе в исследованных расчетных случаях последовательно уменьшаются давление и плотность набегающего потока, что позволяет проследить за постепенным увеличением роли неравновесных процессов. Скорость потока положена постоянной скорости ${{V}_{\infty }}$ = 4 × 10⁵ см/с.

Матрица рассчитанных вариантов дана в табл. 1. Для справки в правой колонке приведены теоретические значения давления торможения ${{p}_{*}} = {{\rho }_{\infty }}V_{\infty }^{2}$.

Таблица 1.

Исходные данные расчетов

№ расчетной серии	${{p}_{\infty }}$, эрг/см³	${{\rho }_{\infty }}$, г/см³	${{p}_{*}}$, эрг/см³
1	55 300	0.889 × 10^–4	1.42 × 10⁷
2	2870	0.400 × 10^–5	6.40 × 10⁵
3	798	0.103 × 10^–5	1.70 × 10⁵

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

На рис. 2–4 представлены результаты расчетов обтекания клина с углом полураствора φ = 3°. Для условий вариантов 1 и 2 из табл. 1 отметим характерные особенности течения газа в указанных условиях.

Рис. 2.

Поля поступательной температуры (а, сверху) и колебательной температуры N₂ (а, снизу), продольной скорости $Vx = {u \mathord{\left/ {\vphantom {u {{{V}_{\infty }}}}} \right. \kern-0em} {{{V}_{\infty }}}}$ (б, сверху) и числа Маха (б, снизу), концентрации электронов (в, сверху) и ионов NO⁺ (в, снизу) при ${{p}_{\infty }}$ = 55 300 эрг/см³, ${{\rho }_{\infty }}$ = 0.889 × 10^–4 г/см³, ${{T}_{\infty }}$ = 217 K и ${{V}_{\infty }}$ = 4 × 10⁵ см/с. Здесь и далее цена деления цветовых шкал обозначается в экспоненциальной форме (например, 2.20E+03 означает 2200).

Рис. 3.

Распределения вдоль координаты y поступательной и колебательной температур N₂ (а), а также концентрации электронов в 6 сечениях вдоль оси x при ${{p}_{\infty }}$ = 2870 эрг/см³, ${{\rho }_{\infty }}$ = 0.400 × 10^–5 г/см³, ${{T}_{\infty }}$ = 250 K и ${{V}_{\infty }}$ = 4 × 10⁵ см/с.

Рис. 4.

Поля числовых концентраций атомов N (а, сверху) и окиси азота NО (а, снизу), атомов кислорода (б, сверху) и молекул О₂ (б, снизу), электронов (в, сверху) и ионов NO⁺ (в, снизу) при ${{p}_{\infty }}$ = 2870 эрг/см³, ${{\rho }_{\infty }}$ = 0.400 × × 10^–5 г/см³, ${{T}_{\infty }}$ = 250 K и ${{V}_{\infty }}$ = 4 × 10⁵ см/с.

На рис. 2а показаны поля поступательной температуры Т и колебательной температуры молекулы азота ${{T}_{{V1}}}$. Наибольшие значения поступательной температуры достигаются в окрестности лобового затупления клина. Здесь температура близка к температуре торможения T ∼ 4500 K. Колебательная температура N₂ заметно отличается от поступательной в зоне ближнего следа, где реализуется отрывное течение разрежения (см. рис. 2б). На рис. 2в приведены поля концентраций электронов и ионов окиси азота NO⁺.

В рассматриваемом на рис. 2в случае электронная концентрация в сжатом слое над поверхностью клина достигает величины n_e ∼ 10¹¹ cм^–3. При этом распределение концентраций электронов и ионов NO⁺ по толщине сжатого слоя показано на рис. 2в, откуда хорошо видно, что толщина ионизованной области достигает почти 10 см, а степень ионизации падает незначительно вдоль образующей клина. Распределения поступательной и колебательной температур N₂ в сжатом слое показывают хорошую термализацию внутренних степеней свободы.

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

Однако уменьшение давления и плотности в набегающем потоке привело к значительной неравновесности в потоке. На рис. 3а показаны распределения поступательной температуры и колебательной температуры N₂ по высоте сжатого слоя в разных сечениях вдоль оси х, откуда хорошо видно, что в максимумах это различие достигает ∼1000 K. Такое отличие температур связано с замедленной релаксацией колебательно возбужденных состояний молекул в течении у боковой поверхности клина аналогично тому, как это наблюдается в отрывных течениях, где температура поступательных степеней свободы падает гораздо быстрее.

Следствием указанной неравновесности является различие в электронных концентрациях, рассчитанных в приближении локального термодинамического равновесия поступательных и колебательных степеней свободы и по модели неравновесной диссоциации N₂ и O₂. Распределения концентраций электронов в этом случае поперек сжатого слоя в разных сечениях вдоль продольной координаты показаны на рис. 3б (соответственно E_q(N_e) и Ne_q(N_e)).

Поля числовых концентраций атомарных и молекулярных компонент, а также электронов и ионов NO⁺ для рассматриваемого случая неравновесных условий показаны на рис. 4. Отметим хорошую корреляцию в распределениях концентраций частиц с температурным распределением. Числовая концентрация молекул О₂ возрастает в окрестности ударной волны.

При уменьшении скорости ${{V}_{\infty }}$ для двух исследованных вариантов до величины ∼3 × 10⁵ см/с наблюдается ожидаемое падение поступательной и колебательной температур в сжатом слое. Заметно уменьшается концентрации электронов в сжатом слое и ионов NO⁺. Наибольшая концентрация электронов не превышает n_e ∼ 10¹¹ см^–3. Это означает, что при дальнейшем уменьшении скорости движения затупленного клина концентрация электронов будет быстро уменьшаться.

При дальнейшем уменьшении давления и плотности набегающего потока течение становится в большей степени неравновесным. Если на относительно больших скоростях колебательная температура ${{T}_{V}}$ молекул N₂ в сжатом слое у боковой поверхности превышает поступательную температуру Т (как показано на рис. 3), то при низких скоростях наблюдается обратная ситуация: ${{T}_{V}} < T$. Это объясняется тем, что при относительно малых скоростях на больших высотах колебательная температура ${{T}_{V}}$ не успевает достичь поступательной температуры Т.

Картина течения значительно усложняется при ненулевом угле атаки. На рис. 5–9 показаны распределения газодинамических функций и концентраций электронов при обтекании затупленного клина под углом атаки α = 6°. При еще меньшей плотности набегающего потока (вариант 3 из табл. 1), а значит – при более высокой степени неравновесности.

Рис. 5.

Поле чисел Маха при обтекании затупленного клина под углом атаки α = 6° при ${{p}_{\infty }}$ = 798 эрг/см³, ${{\rho }_{\infty }}$ = = 0.103 × 10^–5 г/см³, ${{T}_{\infty }}$ = 271 K и ${{V}_{\infty }}$ = 4 × 10⁵ см/с.

Рис. 6.

Поле давления (Pres = ${p \mathord{\left/ {\vphantom {p {{{p}_{*}}}}} \right. \kern-0em} {{{p}_{*}}}}$) при обтекании затупленного клина под углом атаки α = 6° при ${{p}_{\infty }}$ = = 798 эрг/см³, ${{\rho }_{\infty }}$ = 0.103 × 10^–5 г/см³, ${{T}_{\infty }}$ = 271 K и ${{V}_{\infty }}$ = 4 × 10⁵ см/с.

Рис. 7.

Поле температуры поступательных степеней свободы (в K) при обтекании затупленного клина под углом атаки α = 6° при ${{p}_{\infty }}$ = 798 эрг/см³, ${{\rho }_{\infty }}$ = 0.103 × 10^–5 г/см³, ${{T}_{\infty }}$ = 271 K и ${{V}_{\infty }}$ = 4 × 10⁵ см/с.

Рис. 8.

Поле температуры колебательного возбуждения молекул N₂ (в K) при обтекании затупленного клина под углом атаки α = 6° при ${{p}_{\infty }}$ = 798 эрг/см³, ${{\rho }_{\infty }}$ = 0.103 × 10^–5 г/см³, ${{T}_{\infty }}$ = 271 K и ${{V}_{\infty }}$ = 4 × 10⁵ см/с.

Рис. 9.

Числовая концентрация электронов (в см^–3) при обтекании затупленного клина под углом атаки α = 6° при ${{p}_{\infty }}$ = 798 эрг/см³, ${{\rho }_{\infty }}$ = 0.103 × 10^–5 г/см³, ${{T}_{\infty }}$ = 271 K и ${{V}_{\infty }}$ = 4 × 10⁵ см/с.

Распределение чисел Маха (рис. 5) и давления (рис. 6), как и ожидалось, свидетельствует о появлении зоны разрежения потока над верхней гранью клина. Вблизи нижней грани создается область повышенного давления, которое, однако, более чем на порядок меньше давления торможения, достигаемого в окрестности затупленной кромки.

Распределение температуры поступательных (рис. 7) и колебательных степеней свободы молекулы N₂ (рис. 8) указывает на сжатие потока у нижней грани клина. И, наоборот, вблизи верхней грани наблюдается расширение потока, что является прямым следствием обтекания клина под углом атаки. В следе температура поступательных степеней свободы имеет ряд характерных особенностей: имеются заметный скачок температуры поперек нижней ударной волны и закономерное значительное повышение температуры там, где поток затормаживается с возрастанием давления от сходящихся потоков.

Примечательным является то, что конфигурация области повышенной числовой концентрации электронов (рис. 9) хорошо коррелирует с конфигурацией области повышенных колебательных температур (рис. 8). Это еще раз подтверждает важность корректного решения задачи о неравновесной диссоциации и релаксации колебательного возбуждения, что оказывает значительное влияние на ионизационные процессы. Ранее, при решении задачи неравновесного излучения фронта ударной волны [23] также было показано, что эффективная температура электронного газа оказывается весьма близкой колебательной температуре двухатомных молекул воздуха в силу эффективности резонансных процессов электронно-колебательного взаимодействия.

ЗАКЛЮЧЕНИЕ

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

Даны примеры условий, при которых с хорошей точностью соблюдается локальное термодинамическое равновесие. Выбранные для иллюстративных расчетов условия обтекания являются, в известной мере, пограничными – при уменьшении скорости потока относительно рассмотренных, интенсивность ионизационных процессов резко падает. Заметим, что согласно [1, 2], при числовых концентрациях электронов ${{N}_{e}}$ ∼ 10¹⁰–10¹¹ см^–3 обычно наблюдается блокировка радиосвязи (зависит от частоты радиосигналов).

Работа выполнена по теме государственного задания (№ госрегистрации АААА-А20-120011690135-5) и при поддержке гранта РФФИ 19-01-00515.

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

Мартин Дж. Вход в атмосферу. Введение в теорию и практику. М.: Мир, 1969. 320 с.
Суржиков С.Т. Расчетный анализ ионизации сжатого слоя при входе космического аппарата SCHIAPARELLI в плотные слои атмосферы Марса // Изв. РАН. МЖГ. 2020. № 3. С. 80–93.
Wilson J. Ionization rate in air behind high speed shock waves // Phys. Fluids. 1966. V. 9. № 10. P.1913–1921.
Ступоченко Е.В., Лосев С.А., Осипов А.И. Релаксационные процессы в ударных волнах. М.: Изд-во “Наука”, Главная редакция физико-математической литературы, 1965. 484 с.
Залогин Г.Н., Лунев В.В., Пластинин Ю.А. Ионизация и неравновесное излучение за ударными волнами в воздухе // Известия АН ССР. МЖГ. 1980. № 1. С. 105–112.
Горелов В.А., Кильдюшова Л.А., Чернышев В.Л. Об измерении ионизации воздуха за сильными ударными волнами // ТВТ. 1983. Т. 21. № 3. С. 449–453.
Горелов В.А., Киреев А.Ю., Шиленков С.В. Неравновесные ионизационные процессы за сильными ударными волнами при высоких скоростях распространения в воздухе // Физико-химическая кинетика в газовой динамике, 2007. Т. 5. 12 с. http://chemphys.edu.ru/issues/2007-5/articles/55/
Горелов В.А., Кильдюшова Л.А. Особенности процессов ионизации и излучения за сильными ударными волнами в воздухе // ПМТФ. 1987. № 6. С. 23–28.
Суржиков С.Т. Компьютерная аэрофизика спускаемых космических аппаратов. Двухмерные модели. М.: Наука, 2018. 543 с.
Park C. Nonequilibrium Hypersonic Aerothermodynamics. N.Y. Wiley-Intern. Publ., 1990. P. 358.
Kozlov P.V., Surzhikov S.T. Nonequilibrium radiation NO in shocked air // AIAA paper 2017-0157. 16 p.
Ельяшевич М.А. Атомная и молекулярная спектроскопия. М.: Эдиториал УРСС, 2000. 896 с.
Park C. Review of Chemical-Kinetics Problems of Future NASA Missions, I: Earth Entries // Journal of Thermophysics and Heat Transfer. 1993. V. 7. № 3. P. 385–398.
Candler G.V., MacCormack R.W. 1991 Computation of Weakly Ionized Hypersonic Flows in Thermochemical Nonequilibrium. Journal of Thermophysics. 1991. V. 5. № 3. P. 266–273.
Суржиков С.Т. Двумерный численный анализ ионизации потока в летном эксперименте RAM-C-II // Химическая физика. 2015. Т. 34. № 2. С. 24–42.
Surzhikov S.T. Numerical Rebuilding of Shock Layer Ionization for Two Flight Tests // AIAA paper № 2016-4305. 2016. 32 p.
Базаров И.П., Геворкян Э.В., Николаев П.Н. Термодинамика и статистическая физика. М.: Изд-во Московского ун-та, 1986. 310 с.
Millikan R.C., White D.R. 1963 Systematic of Vibrational Relaxation // J. Chemical Physics. V. 39. № 12. P. 3209–3212.
Гиршфельдер Дж., Кертис Ч., Берд Р. Молекулярная теория газов и жидкостей. М.: Изд-во иностранной литературы, 1961. 929 с.
Берд Р., Стьюарт В., Лайтфут Е. Явления переноса. М.: Изд-во “Химия”, 1974. 687 с.
Гурвич Л.В., Вейц И.В., Медведев В.А. и др. Термодинамические свойства индивидуальных веществ. М.: Наука, 1978. 495 с.
JANAF Thermochemical Tables. 3rd Ed. Chase M.W., Jr., Davies C.A., Downey J.R., Jr., Fririp D.J., Mc Donald R.A., Syverud A.N. // J. Phys. Chem. Ref. Data. 1985. V. 14. Suppl. 1.
Суржиков С.Т. Расчет неравновесного излучения ударных волн в воздухе с использованием двух моделей // Изв. РАН. МЖГ. 2019. № 1. С. 99–114.

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

Инструменты

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

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

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