Оценка начала ламинарно-турбулентного перехода на пластине при полете в атмосфере Марса

Полный текст

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

Вместе с тем по прогнозам [1] первые пилотируемые полеты на Марс начнутся уже в следующем десятилетии. Эскизные проекты обитаемых космических аппаратов (КА), предназначенных для доставки экспедиций на планету с безопасными перегрузками, представляются конструкциями с характерными размерами в десятки метров. В частности, рассматривается проект марсианского КА, близкого по конфигурации и размерам к аппаратам типа “Буран” (Space Shuttle).

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

При проектировании систем управления и теплозащиты КА необходимы расчеты устойчивости пограничного слоя для оценки положения зоны перехода к турбулентности, где возникает пик теплового потока. Для расчетов положения зоны перехода широко используется e^N-метод, основанный на представлении об экспоненциальном нарастании пространственных возмущений до определенного уровня амплитуды A_T к началу ламинарно-турбулентного перехода (ЛТП). Амплитуды для частотного спектра нарастающих возмущений рассчитываются на основе линейной теории [6, 7] в соответствии с соотношением

A_T = A₀exp(N_T),

где A_T – амплитуда в точке перехода; A₀ – амплитуда начального возмущения; N_T – так называемый N-фактор перехода, значение которого принимается в зависимости от конкретной задачи.

В расчетах необходимо использовать реальные характеристики марсианской атмосферы, около 96% которой составляет углекислый газ CO₂, физико-химическая кинетика которого существенно сложней кинетики двухатомных компонент воздуха. Несмотря на разреженность и низкую температуру атмосферы Марса высокие скорости на траектории посадки в пределах 2–6 км/с определяют высокие температуры в ПС за ударной волной, которые в верхней части траектории существенно превышают температуру полной диссоциации CO₂. При этом в пограничном слое возникает многокомпонентная термохимически неравновесная смесь атомов и колебательно возбужденных молекул, для учета свойств которой необходима адекватная физико-математическая модель.

В статье выполнены сравнительные расчеты положения начала зоны ЛТП для двух характерных режимов течения гиперзвукового пограничного слоя, соответствующих 66-й и 87-й секундам полета аппарата Pathfinder при посадке на Марс [5]. В данном случае они рассматриваются как последовательные положения на многовитковой траектории в процессе аэродинамического торможения.

В качестве наиболее опасного взято двумерное возмущение моды II, совпадающее с направлением несущего потока. Значение показателя степени нарастания возмущений выбрано N = 8. Для расчета стационарных течений использовались локально автомодельные уравнения плоского пограничного слоя на пластине с изотермической некаталитической поверхностью.

1. Модель течения и основные уравнения

В качестве исходной математической модели рассматривается система двухтемпературных уравнений сверхзвукового течения углекислого газа, учитывающая колебательную релаксацию и реакции диссоциации–рекомбинации с образованием молекул угарного газа CO и атомарного кислорода O:

$C{O}_2+A\rightleftarrows CO+O+A.$ (1.1)

Здесь частица – партнер по столкновению A – представляет собой одну из компонент смеси CO₂–CO–O. Рассматриваемая модель в сравнительных расчетах [3] аэродинамических характеристик при спуске конического КА в атмосфере Марса показала малое отличие от моделей с более сложной химической кинетикой.

Относительно молекул угарного газа CO принято предположение [8], что при диссоциации CO₂ молекула CO мгновенно приходит в термодинамическое равновесие со средой, а ее дальнейшую диссоциацию можно не рассматривать. При этом колебательная энергия молекул CO находится в состоянии равновесия и характеризуется статической (поступательной) температурой потока. Это связано с тем, что характерное время диссоциации CO значительно больше, чем для молекулы углекислого газа CO₂.

Система гидродинамических переменных включает в себя скорость потока u, плотность смеси ρ, массовые доли молекул CO₂ – Y₁, CO – Y₂ и атомов O – Y₃, статическое давление p, удельную внутреннюю энергию смеси, связанную с квазиравновесными внутренними степенями свободы e_t, удельную колебательную энергию молекул CO₂ – e_v1, поступательную (статическую) температуру T, колебательную температуру T_v.

Поскольку сумма массовых долей компонентов смеси равна единице, то массовая доля любой компоненты может быть выражена через массовые доли двух других компонент, в частности Y₃ = 1 – Y₁ – Y₂. Замкнутая система в терминах этих переменных включает следующие уравнения.

Уравнение неразрывности смеси в целом

$\frac{\partial \rho }{\partial t}+\nabla \left(\rho u\right)=0.$ (1.2)

Уравнения неразрывности компонентов смеси

$\frac{\partial \rho Y_i}{\partial t}+\nabla \left(\rho Y_{i}u\right)=-\nabla J_i+W_i,i=\mathrm{1,}2.$ (1.3)

Уравнения импульсов для смеси

$\frac{\partial \rho u}{\partial t}+\nabla \left(\rho uu\right)=-\nabla p+\nabla ⃖.$ (1.4)

Уравнение для плотности внутренней поступательно-вращательной энергии смеси

$\frac{\partial \rho e_t}{\partial t}+\nabla \cdot u\rho e_t+p\nabla u=\nabla \left[\lambda \nabla T+\left(⃖u\right)-\sum _i{c}_{pi}T{J}_i\right]-\sum _i{h}_{0i}{W}_i-R_{v1}.$ (1.5)

Уравнение для плотности колебательной энергии CO₂

$\frac{\partial \rho Y_1{e}_{v1}}{\partial t}+\nabla \left(\rho u{Y}_1{e}_{v1}\right)=\nabla \left({\lambda }_{v1}\nabla T_v-e_{v1}{J}_1\right)+R_{v1}.$ (1.6)

Уравнение состояния газовой смеси

$p=\rho RT\sum _i\left(Y_i/M_i\right).$ (1.7)

Здесь J_i – плотность диффузионного потока i-й компонента смеси; W_i – скорость производства i-й компоненты смеси в химических реакциях. Тензор вязких напряжений [9]

где η, η_b – динамическая и объемная вязкости смеси; R_v1 – источниковый член, описывающий колебательно-диссоциативное взаимодействие; c_pi = c_Vi + R/M_i – удельная теплоемкость при постоянном давлении, связанная с квази- равновесными внутренними степенями свободы; h_0i – удельная энтальпия образования i-го компонента; λ – теплопроводность смеси; λ_v1 – теплопроводность, обусловленная колебательными модами молекул CO₂; M_i – молярная масса i-й компоненты; R – универсальная газовая постоянная.

Плотность внутренней поступательно-вращательной энергии молекул смеси определяется следующим образом:

$e_t=\sum _{s=1}^3{Y}_s{e}_{t,s}=\sum _{s=1}^3{Y}_s{c}_{V,s}T,$

где $c_{Vi}=c_{Vi}^{tr}+c_{Vi}^{rot};c_{Vi}^{tr},c_{Vi}^{rot}$ – теплоемкости квазиравновесных степеней свободы i-й компоненты при постоянном объеме. Теплоемкости c_Vi равны $c_{V1}=5R/\left(2M_1\right),$ $c_{V2}=7R/\left(2M_2\right),$ $c_{V3}=3R/\left(2M_3\right).$ 

Здесь в равновесную теплоемкость угарного газа в связи со сделанным предположением добавлена равновесная теплоемкость ее колебательной моды c_V2 = R/M₂. Соответственно, теплоемкости при постоянном давлении есть c_p1 = 7R/(2M₁), c_p2 = 9R/(2M₂), c_p3 = 5R/(2M₃). Числовые значения энтальпий образования h_0i для компонент смеси CO₂–CO–O приведены в табл. 1 [10].

 

Таблица 1. Стандартные энтальпии образования

Компонент смеси	CO2	CO	O
h0i, кДж/моль	–393.51	-110.52	249.18

 

1.1. Коэффициенты переноса

Коэффициенты вязкости η и теплопроводности λ смеси вычисляются по формулам Вилке [11]

$f=\sum _j\frac{X_j{f}_j}{{\phi }_{jf}},{\phi }_{jf}=\sum _i\frac{X_i{\left[1+\sqrt{f_j/f_i}{\left(M_i/M_j\right)}^{1/4}\right]}^2}{\sqrt{8\left(1+M_j/M_i\right)}},X_j=\frac{Y_j}{M_j}{\left(\sum _i\frac{Y_i}{M_i}\right)}^{-1},$

где f_i – вязкость или теплопроводность i-й компоненты смеси. Вязкость i-й компоненты рассчитывается по формуле Блотнера [12]:

${\eta }_i\left(T\right)=0.1\exp \left[\right(A_i\ln T+B_i)\ln T+C_i],$

которая применима при высоких температурах вплоть до 10⁴ K. A_i, B_i и C_i здесь являются постоянными величинами [4], числовые значения которых для компонент смеси CO₂–CO–O приведены в табл. 2.

 

Таблица 2. Коэффициенты A_i, B_i и C_i вязкой модели Блотнера

Компонент смеси	Ai	Bi	Ci
CO2	–0.01952739	1.047818	–14.32212
CO	–0.01952739	1.013295	–13.97873
O	0.02031440	0.429440	–11.60314

 

Для вычисления теплопроводности i-й компоненты λ_i используются модифицированные соотношения Эйкена [13]:

${\lambda }_i={\lambda }_i^{tr}+{\lambda }_i^{rot}=\left(\frac{5}{2}{c}_{Vi}^{tr}+\frac{6}{5}{c}_{Vi}^{rot}\right){\eta }_i,i=\mathrm{1,}\mathrm{2,}3.$

Теплопроводности λ_vi, обусловленные колебательными модами молекул CO₂ и CO, вычисляются следующим образом:

${\lambda }_i^{vib}=\frac{6}{5}{c}_{Vi}^{vib}{\eta }_i,c_{V1}^{vib}=\frac{d{e}_{v1}}{d{T}_v},c_{V2}^{vib}=\frac{R}{M_2},i=\mathrm{1,}2.$

Подстановка значений теплоемкостей дает

${\lambda }_1=\left(\frac{99R}{20M_1}\right){\eta }_1,{\lambda }_2=\left(\frac{123R}{20M_2}\right){\eta }_2,{\lambda }_3=\left(\frac{15R}{4M_3}\right){\eta }_3.$

Плотность диффузионного потока определяется законом Фика:

$J_i=-\rho D_i\nabla Y_i,$

где D_i – эффективный коэффициент диффузии i-й компоненты. Аналогично [3, 14, 15] предполагается, что диффузия компонент смеси характеризуется с помощью постоянного числа Шмидта Sc = 0.5. При этом

$\rho D_i=\frac{{\eta }_i}{Sc}.$ (1.8)

1.2. Модель колебательно-диссоциативного взаимодействия

Колебательно-диссоциативное взаимодействие описывается источниковым членом R_v1, который входит с противоположными знаками в правые части уравнений (1.5) и (1.6) и записывается следующим образом:

$R_{v1}=R_{v1}^{vib}+R_{v1}^{chem},$ (1.9)

где первое слагаемое определяет вклад колебательной релаксации возбужденных колебательных мод молекул CO₂, а второе слагаемое – связь между химическими реакциями и колебательной релаксацией. При этом первое слагаемое в правой части (1.9) выражается формулой Ландау–Теллера

$R_{v1}^{vib}=\rho Y_1\left(\frac{e_{v1}\left(T\right)-e_{v1}\left(T_v\right)}{{\tau }_1}\right),$

где τ₁ – характерное время колебательной релаксации.

Молекула CO₂ имеет три колебательные моды – продольные симметричную и антисимметричную, а также дважды вырожденную изгибную, характеризуемые различными временами релаксации. В многотемпературном приближении им соответствуют три колебательные температуры. Однако близость характерных частот симметричной и изгибной мод [3] и рассматриваемый высокотемпературный диапазон позволяет характеризовать [16, 17], колебательную энергию молекулы СО₂ одной колебательной температурой T_v. В приближении гармонического осциллятора колебательная энергия молекулы СО₂ записывается следующим образом [4, 16]:

$e_{v1}=\frac{R}{M_1}\left(\frac{{\theta }_1}{\exp \left(\frac{{\theta }_1}{T_v}\right)-1}+\frac{2{\theta }_2}{\exp \left(\frac{{\theta }_2}{T_v}\right)-1}+\frac{{\theta }_3}{\exp \left(\frac{{\theta }_3}{T_v}\right)-1}\right),$

где θ₁ = 1903 K, θ₂ = 945 K, θ₃ = 3329 K − характеристические колебательные температуры первой, второй и третьей мод.

Для времени релаксации τ₁ в [17] предложена универсальная зависимость

$p{\tau }_1=\exp \left(36.5T^{-1/3}-17.71\right),$

где τ₁ измеряется в секундах, давление p − в атмосферах.

Второе слагаемое в правой части (1.9), определяющее связь между химическими реакциями и колебательной релаксацией, представляется в виде

$R_{v1}^{chem}=e_{v1}{W}_1.$

1.3. Химические реакции

В статье рассматриваются три реакции диссоциации–рекомбинации, различающиеся частицей-партнером по столкновению, которые записывается в виде (1.1). Скорости производства компонент смеси в системе (1.2)–(1.7) определяется следующим образом:

$W_1=-{\rho }^2{M}_1\left(k_f\frac{Y_1}{M_1}-k_b\rho \frac{Y_2}{M_2}\frac{Y_3}{M_3}\right)\left(\frac{Y_1}{M_1}+\frac{Y_2}{M_2}+\frac{Y_3}{M_3}\right),W_2=-\left(\frac{M_2}{M_1}\right)W_1,W_3=-\left(\frac{M_3}{M_1}\right)W_1.$

Используется чисто аррениусовская кинетика, определяемая формулами [4, 16]

$k_f=C_f{T}^f\exp \left[-\frac{{\theta }_f}{T}\right],k_b=C_b{T}^b\exp \left[-\frac{{\theta }_b}{T}\right].$

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

Принимается, что для рассматриваемых здесь реакций константы скоростей прямой k_f и обратной k_b реакций не зависят от вида частицы – партнера по столкновению. Такой выбор связан с тем, что нас в первую очередь интересует реакция диссоциации СО₂, для которой это имеет место [18].

В расчетах использовались следующие числовые данные: C_f = 3.7 × 10¹¹ м³/(кмоль × с), C_b = = 6.1 × 10³ м⁶/(кмоль² × с), f = 0, b = 0.75, θ_f = 52 500 K и θ_b = –10 240 K [4, 16, 18].

2. Обезразмеривание и линеаризация

В качестве характерных величин для обезразмеривания уравнений (1.2)–(1.7) выбраны текущее расстояние вдоль пластины x = L, параметры невозмущенного потока вне пограничного слоя, отмеченные индексом ∞: скорость U_∞, плотность ρ_∞ и температура T_∞. Коэффициенты переноса смеси η, λ и ее компонент η_i, λ_i обезразмериваются на значения коэффициетов переноса для CO₂ в невозмущенном потоке η_1∞ и λ_1∞, скорости прямой k_f и обратной k_b химических реакций – на k_f∞ и k_b∞, время колебательной релаксации τ₁ – на τ_1∞. Скорости производства компонент смеси W_i нормировались на комплекс ${\rho }_{\infty }^2{k}_{f\infty }/M_1,$ а энергии e, e_i, e_v1 и энтальпии образования компонент смеси h_0i – на комплекс T_∞R/M₁. Для обезразмеривания давления и времени используются комбинированные величины ${\rho }_{\infty }{U}_{\infty }^2$ и L/U_∞ соответственно.

В обезразмеренную систему уравнений кроме числа Шмидта (1.8) входят следующие безразмерные критерии $\mathrm{Re}={\eta }_{1\infty }L{U}_{\infty }/{\rho }_{\infty }$ – числа Рейнольдса и Прандтля $\mathrm{Pr}={\eta }_{1\infty }{c}_{p1}/{\lambda }_{1\infty };$ число Маха $M=U_{\infty }/\sqrt{\gamma R{T}_{\infty }/M_1};$ числа Дамкелера химических реакций $D{a}_d=t_{\infty }/{\tau }_{d\infty }$ и VT-энергообменов $D{a}_{VT}=t_{\infty }/{\tau }_{1\infty }.$ Здесь g = c_p1/c_V1 – показатель адиабаты; ${\tau }_{d\infty }=M_1/\left({\rho }_{\infty }{k}_{f\infty }\right)$ – характерное время реакции диссоциации; t_1∞ – характерное время релаксации колебательных мод CO₂.

Для рассматриваемой задачи лиейной устойчивости обезразмеренная система (1.2)–(1.7) линеаризовалась на стационарном решении уравнений пограничного слоя в локально параллельном приближении. Мгновенные значения газодинамических переменных представлялись в виде суммы стационарного решения и возмущений:

$u_x=U_s\left(y\right)+u^{\text{'}},u_y=v^{\text{'}},\rho ={\rho }_s\left(y\right)+{\rho }^{\text{'}},Y_i=Y_{is}\left(y\right)+{Y^{\text{'}}}_i,T=T_s\left(y\right)+T^{\text{'}},T_v=T_{vs}\left(y\right)+{T^{\text{'}}}_v,$

где индексом s отмечено стационарное решение, а штрихом – возмущения газодинамических переменных. Для обоих режимов рассчитывался альтернативный случай совершенного газа при Da_d = Da_VT = 0. Мгновенные значения коэффициентов переноса, скоростей производства компонент смеси и времени колебательной релаксации определялись следующим образом:

$\eta =\eta \left(T_s\right)+{\eta }^{\text{'}}={\eta }_s+\left({\left.\frac{\partial \eta }{\partial T}\right|}_{f=f_s}\right)T^{\text{'}}+\sum _{i=1}^3\left({\left.\frac{\partial \eta }{\partial Y_i}\right|}_{f=f_s}\right){Y^{\text{'}}}_i={\eta }_s+{\eta }_{Ts}{T}^{\text{'}}+\sum _{i=1}^3{\eta }_{Yis}{Y^{\text{'}}}_i,$

$\lambda =\lambda \left(T_s\right)+{\lambda }^{\text{'}}={\lambda }_s+\left({\left.\frac{\partial \lambda }{\partial T}\right|}_{f=f_s}\right)T^{\text{'}}+\sum _{i=1}^3\left({\left.\frac{\partial \lambda }{\partial Y_i}\right|}_{f=f_s}\right){Y^{\text{'}}}_i={\lambda }_s+{\lambda }_{Ts}{T}^{\text{'}}+\sum _{i=1}^3{\lambda }_{Yis}{Y^{\text{'}}}_i,$

$\begin{array}{c}{W}_j=W_j\left(T_s,{\rho }_s,{Y^{\text{'}}}_{1s},{Y^{\text{'}}}_{2s},{Y^{\text{'}}}_{3s}\right)+{W^{\text{'}}}_j=W_{js}+\left({\left.\frac{\partial W_j}{\partial T}\right|}_{f=f_s}\right)T^{\text{'}}+\left({\left.\frac{\partial W_j}{\partial \rho }\right|}_{f=f_s}\right){\rho }^{\text{'}}+\\ +\sum _{i=1}^3\left({\left.\frac{\partial W_j}{\partial Y_i}\right|}_{f=f_s}\right){Y^{\text{'}}}_i=W_{js}+W_{j,Ts}{T}^{\text{'}}+W_{j,\rho s}{\rho }^{\text{'}}+\sum _{i=1}^3{W}_{j,Yis}{Y^{\text{'}}}_i,j=\mathrm{1,}\mathrm{2,}\mathrm{3,}f=\left(\rho ,T,Y_1,Y_2,Y_3\right),\end{array}$

${\eta }_i={\eta }_i\left(T_s\right)+{{\eta }^{\text{'}}}_i={\eta }_{is}+\left({\left.\frac{d{\eta }_i}{dT}\right|}_{T=T_s}\right)T^{\text{'}}={\eta }_{is}+{\eta }_{i,Ts}{T}^{\text{'}},{\lambda }_i={\lambda }_i\left(T_s\right)+{{\lambda }^{\text{'}}}_i={\lambda }_{is}+\left({\left.\frac{d{\lambda }_i}{dT}\right|}_{T=T_s}\right)T^{\text{'}}={\lambda }_{is}+{\lambda }_{i,Ts}{T}^{\text{'}},$

${\tau }_1={\tau }_1\left(T_s\right)+{\tau }^{\text{'}}={\tau }_{1s}+\left({\left.\frac{d{\tau }_1}{dT}\right|}_{T=T_s}\right)T^{\text{'}}={\tau }_{1s}+{\tau }_{\mathrm{1,}Ts}{T}^{\text{'}}.$

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

$q^{\text{'}}(x,y,t)=q\left(y\right)\exp \left[i\right(\alpha x-\omega t\left)\right],q\left(y\right)=(u,\alpha v,\rho ,y_1,y_2,y_3,\theta ,{\theta }_v),$ (2.1)

где α = α_r + iα_i – комплексное волновое число вдоль координаты x, ω = ω_r + iω_i – комплексная частота, i – мнимая единица. Подстановка (2.1) в обезразмеренную линеаризованную систему дает систему уравнений для амплитуд двумерных возмущений, которая имеет вид

$D\rho +\alpha {\rho }_s\sigma +\alpha v\frac{d{\rho }_s}{dy}=0,$ (2.2)

$\frac{{\eta }_s}{Sc\mathrm{Re}}{\Delta }_y{y}_1-{\rho }_{s}D{y}_1-\alpha v{\rho }_s\frac{d{Y}_{1s}}{dy}+\frac{1}{Sc\mathrm{Re}}\frac{d{\eta }_s}{d\text{}y}\frac{d{y}_1}{dy}+\frac{1}{Sc\mathrm{Re}}\frac{d}{dy}\left(\eta \frac{d{Y}_{1s}}{dy}\right)+D{a}_d{W}_1=\mathrm{0,}$ (2.3)

$\frac{{\eta }_s}{Sc\mathrm{Re}}{\Delta }_y{y}_2-{\rho }_{s}D{y}_2-\alpha v{\rho }_s\frac{d{Y}_{2s}}{dy}+\frac{1}{Sc\mathrm{Re}}\frac{d{\eta }_s}{d\text{}y}\frac{d{y}_2}{dy}+\frac{1}{Sc\mathrm{Re}}\frac{d}{dy}\left(\eta \frac{d{Y}_{2s}}{dy}\right)+D{a}_d{W}_2=\mathrm{0,}$ (2.4)

$\frac{{\eta }_s}{\mathrm{Re}}{\Delta }_{y}u-{\rho }_{s}Du-\alpha v{\rho }_s\frac{d{U}_s}{dy}-i\alpha \epsilon +\frac{1}{\mathrm{Re}}\frac{d{\eta }_s}{dy}\left(\frac{du}{dy}+i{\alpha }^{2}v\right)+\frac{1}{\mathrm{Re}}\frac{d}{dy}\left(\eta \frac{d{U}_s}{dy}\right)=0,$ (2.5)

$\frac{\alpha {\eta }_s}{\mathrm{Re}}{\Delta }_{y}v-\alpha {\rho }_{s}Dv-\frac{d\epsilon }{dy}+\frac{\alpha }{\mathrm{Re}}\left(\frac{dv}{dy}-iu\right)\frac{d{\eta }_s}{dy}+\frac{i\alpha \eta }{\mathrm{Re}}\frac{d{U}_s}{dy}=0,$ (2.6)

$\begin{array}{c}\frac{7{\lambda }_s}{2\mathrm{PrRe}}{\Delta }_y\theta -{\rho }_{s}De-\alpha v{\rho }_s\frac{d{e}_s}{dy}-\alpha \gamma M^2{p}_s\sigma +\frac{\gamma M^2}{\mathrm{Re}}\left[2{\eta }_s\left(\frac{du}{dy}+i{\alpha }^{2}v\right)+\eta \frac{d{U}_s}{dy}\right]\frac{d{U}_s}{dy}+\\ +\frac{7}{2\mathrm{PrRe}}\left[\frac{d{\lambda }_s}{dy}\frac{d\theta }{dy}+\frac{d}{dy}\left(\lambda \frac{d{T}_s}{dy}\right)\right]+\frac{1}{Sc\mathrm{Re}}\sum _{j=1}^3{c}_{pj}\left\{{\eta }_s{T}_s{\Delta }_y{y}_j+\frac{d{y}_j}{dy}\left(\frac{d{\eta }_s{T}_s}{dy}\right)+\frac{d}{dy}\left[(\eta T_s+{\eta }_s\theta )\frac{d{Y}_{js}}{dy}\right]\right\}-\end{array}$ (2.7)

$\begin{array}{c}-D{a}_d\sum _{j=1}^3{h}_{0j}{W}_j-D{a}_d\left(e_{v1}{W}_{1s}+e_{v1s}{W}_1\right)-D{a}_{VT}{\rho }_s{Y}_{1s}\left(\frac{e_{v1}\left(\theta \right)-e_{v1}\left({\theta }_v\right)}{{\tau }_{1s}}\right)-\\ -D{a}_{VT}\left[\rho Y_{1s}+{\rho }_s{y}_1-\frac{{\rho }_s{Y}_{1s}}{{\tau }_{1s}}{\tau }_1\right]\left(\frac{e_{v1s}\left(T_s\right)-e_{v1s}\left(T_{vs}\right)}{{\tau }_{1s}}\right)=0,\end{array}$

$\begin{array}{c}\frac{6{\eta }_{1s}}{5\mathrm{Re}}{\Delta }_y{e}_{v1}-{\rho }_s{Y}_{1s}D{e}_{v1}-{\rho }_s{e}_{v1s}D{y}_1-\alpha v{\rho }_s\frac{d{Y}_{1s}{e}_{v1s}}{dy}+\frac{6}{5\mathrm{Re}}\frac{d{\eta }_{1s}}{dy}\frac{d{e}_{v1}}{dy}+\frac{6}{5\mathrm{Re}}\frac{d}{dy}\left({\eta }_1\frac{d{e}_{v1s}}{dy}\right)+\frac{1}{Sc\mathrm{Re}}\times \\ \times \left[{\eta }_s{e}_{v1s}{\Delta }_y{y}_1+\eta e_{v1s}\frac{d^2{Y}_{1s}}{d{y}^2}+e_{v1}\frac{d}{dy}\left({\eta }_s\frac{d{Y}_{1s}}{dy}\right)+\left(\frac{d{\eta }_s{e}_{v1s}}{dy}\right)\frac{d{y}_1}{dy}+\eta \left(\frac{d{e}_{v1s}}{dy}\frac{d{Y}_{1s}}{dy}\right)\right]+D{a}_d\left(e_{v1}{W}_{1s}+e_{v1s}{W}_1\right)+\\ +D{a}_{VT}{\rho }_s{Y}_{1s}\left[\frac{e_{v1}\left(\theta \right)-e_{v1}\left({\theta }_v\right)}{{\tau }_{1s}}\right]+D{a}_{VT}\left(\rho Y_{1s}+{\rho }_s{y}_1-\frac{{\rho }_s{Y}_{1s}}{{\tau }_{1s}}{\tau }_1\right)\left[\frac{e_{v1s}\left(T_s\right)-e_{v1s}\left(T_{vs}\right)}{{\tau }_{1s}}\right]=0.\end{array}$ (2.8)

Здесь

${\Delta }_y=\frac{d^2}{d{y}^2}-{\alpha }^2,D=i\alpha \left(U_s-\frac{\omega }{\alpha }\right),\sigma =\frac{dv}{dy}+iu,\epsilon =p-\frac{\alpha {\eta }_s}{\mathrm{Re}}\left({\mu }_b+\frac{1}{3}\right)\sigma ,$

$p=p_s\left[\frac{\rho }{{\rho }_s}+\frac{\theta }{T_s}+\frac{y_1+y_2\left(\frac{M_1}{M_2}\right)+y_3\left(\frac{M_1}{M_3}\right)}{Y_{1s}+Y_{2s}\left(\frac{M_1}{M_2}\right)+Y_{3s}\left(\frac{M_1}{M_3}\right)}\right],p_s=\frac{{\rho }_s{T}_s}{\gamma M^2}\left(Y_{1s}+Y_{2s}\left(\frac{M_1}{M_2}\right)+Y_{3s}\left(\frac{M_1}{M_3}\right)\right),$

$\begin{array}{c}{y}_3=-y_1-y_2,Y_{3s}=1-Y_{1s}-Y_{2s},\eta ={\eta }_{Ts}\theta +\sum _{j=1}^3{\eta }_{Yjs}{y}_j,\\ \lambda ={\lambda }_{Ts}\theta +\sum _{j=1}^3{\lambda }_{Yjs}{y}_j,{\eta }_j={\eta }_{j,Ts}\theta ,{\tau }_1={\tau }_{\mathrm{1,}Ts}\theta ,\end{array}$

$W_j=W_{j,Ts}\theta +W_{j,\rho s}\rho +\sum _{k=1}^3{W}_{j,Yks}{y}_k,e=e_{Ts}\theta +\sum _{j=1}^3{e}_{Yjs}{y}_i,e_{v1}=e_{v\mathrm{1,}Ts}{\theta }_v.$

3. Метод решения задачи

В качестве профилей стационарного течения использовались решения системы локально автомодельных уравнений пограничного слоя, которая представляется в виде

$({\rho }_s{\eta }_s{\phi }^{\text{'}\text{'}}{)}^{\text{'}}+\phi {\phi }^{\text{'}\text{'}}=\mathrm{0,}$ (3.1)

$({\rho }_s{\eta }_s{Y^{\text{'}}}_{1s}{)}^{\text{'}}+Sc\phi {Y^{\text{'}}}_{1s}+ScD{a}_d\left(\frac{4\xi W_{1s}}{{\rho }_s}\right)=\mathrm{0,}$ (3.2)

$({\rho }_s{\eta }_s{Y^{\text{'}}}_{2s}{)}^{\text{'}}+Sc\phi {Y^{\text{'}}}_{2s}+ScD{a}_d\left(\frac{4\xi W_{2s}}{{\rho }_s}\right)=\mathrm{0,}$ (3.3)

$\begin{array}{c}{\left({\rho }_s{\lambda }_s{T^{\text{'}}}_s\right)}^{\text{'}}+\left(\frac{2\mathrm{Pr}}{7}\right)\phi {h^{\text{'}}}_s+\left(\frac{\gamma \mathrm{Pr}{M}^2}{14}\right){\rho }_s{\eta }_s{{\phi }^{\text{'}\text{'}}}^2+\left(\frac{2\mathrm{Pr}}{7Sc}\right)\sum _{i=1}^3({\rho }_s{\eta }_s{h}_{is}{Y^{\text{'}}}_{is}{)}^{\text{'}}-\left(\frac{8\xi \mathrm{Pr}D{a}_d}{7}\right)\sum _{i=1}^3\frac{h_{0i\infty }{W}_{is}}{{\rho }_s}-\\ -\left(\frac{8\xi \mathrm{Pr}D{a}_{VT}}{7}\right)\left(\frac{Y_1{e}_{v1s}\left(T_s\right)-Y_1{e}_{v1s}\left(T_{vs}\right)}{{\tau }_{1s}}\right)-\left(\frac{8\xi \mathrm{Pr}D{a}_d}{7}\right)\frac{e_{v1s}{W}_{1s}}{{\rho }_s}=\mathrm{0,}\end{array}$ (3.4)

$\begin{array}{c}{\left({\rho }_s{\eta }_{1s}{e^{\text{'}}}_{v1s}\right)}^{\text{'}}+\frac{5}{6}\phi {\left(Y_{1s}{e}_{v1s}\right)}^{\text{'}}+\frac{5}{6Sc}{\left({\rho }_s{\eta }_s{e}_{v1s}{Y^{\text{'}}}_{1s}\right)}^{\text{'}}+\\ +\left(\frac{10\xi {\text{Da}}_{VT}}{3}\right)\left(\frac{Y_{1s}{e}_{v1s}\left(T_s\right)-Y_{1s}{e}_{v1s}\left(T_{vs}\right)}{{\tau }_{1s}}\right)+\left(\frac{10\xi {\text{Da}}_d}{3}\right)\frac{e_{v1s}{W}_{1s}}{{\rho }_s}=\mathrm{0,}\end{array}$ (3.5)

${\rho }_s{T}_s\left(Y_{1s}+Y_{2s}\left(\frac{M_1}{M_2}\right)+Y_{3s}\left(\frac{M_1}{M_3}\right)\right)=\left(Y_{1\infty }+Y_{2\infty }\left(\frac{M_1}{M_2}\right)+Y_{3\infty }\left(\frac{M_1}{M_3}\right)\right),$ (3.6)

где

$h_s=\sum _{i=1}^3{Y}_{is}{h}_{is}=\left(\frac{7}{2}{Y}_{1s}+\frac{9M_1}{2M_2}{Y}_{2s}+\frac{5M_1}{2M_3}{Y}_{3s}\right)T_s,h_{0i\infty }=\frac{h_{0i}}{R{T}_{\infty }/M_1},$

$Y_{3s}=1-Y_{1s}-Y_{2s},Y_{3\infty }=1-Y_{1\infty }-Y_{2\infty },$

$W_{1s}=-{\rho }_s^2\left(Y_{1s}f\left(T_s\right)-{\rho }_{s}b\left(T_s\right)\left[\frac{{\rho }_{\infty }{k}_{b\infty }}{k_{f\infty }{M}_1}\right]\frac{Y_{2s}{M}_1}{M_2}\frac{Y_{3s}{M}_1}{M_3}\right)\left(Y_{1s}+\frac{Y_{2s}{M}_1}{M_2}+\frac{Y_{3s}{M}_1}{M_3}\right),$

$W_{2s}=-\left(\frac{M_2}{M_1}\right)W_{1s},W_{3s}=-\left(\frac{M_3}{M_1}\right)W_{1s},f\left(T_s\right)=\exp \left[{\theta }_{f\infty }\left(1-\frac{1}{T_s}\right)\right],b\left(T_s\right)=T_s^{3/4}\exp \left[{\theta }_{b\infty }\left(1-\frac{1}{T_s}\right)\right],$

$k_{f\infty }=C_f\exp \left(-{\theta }_{f\infty }\right),k_{b\infty }=C_b{T}_{\infty }^{3/4}\exp \left(-{\theta }_{b\infty }\right),{\theta }_{f\infty }=\frac{{\theta }_f}{T_{\infty }},{\theta }_{b\infty }=\frac{{\theta }_b}{T_{\infty }}.$

В уравнениях (3.1)–(3.6) штрихи обозначают дифференцирование по автомодельной переменной [9]

$\zeta =\frac{1}{2\sqrt{\xi }}\underset{0}{\overset{y}{\int }}{\rho }_s\left(z\right)dz,\xi =x,$

вспомогательная функция φ(ζ) определена как

$\phi \left(\zeta \right)=2\underset{0}{\overset{\zeta }{\int }}{U}_s\left(z\right)\text{\hspace{0.17em}}dz.$

При подстановке в уравнения (2.2)–(2.8) профилей газодинамических переменных, полученных при решении системы (3.1)–(3.6), соответствующие естественные поперечные координаты определялись по формуле

$y=2\sqrt{\xi }\underset{0}{\overset{\zeta }{\int }}\frac{dz}{{\rho }_s\left(z\right)}.$

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

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

Числа Рейнольдса ${\mathrm{Re}}_{\delta T}=\sqrt{{\mathrm{Re}}_{xT}},$ определяющие начало зоны ЛТП, рассчитывались по схеме, которая ранее была использована авторами в [19]. Для определения зоны ЛТП рассматривались два вида спектральных задач. Вначале при вещественных волновых числах α = α_T и комплексных частотах ω = ω_r + iω_i рассчитывались собственные значения для возмущений, развивающихся во времени. Затем при комплексных α = α_r + iα_i и вещественных ω = ω_S вычислялись спектры пространственных возмущений. Принималось, что при y = 0 и на условной верхней границе пограничного слоя y = δ амплитуды всех возмущений обращаются в нуль. В основе этой схемы лежит связь, обоснованная в [20], где было показано, что для двумерных пространственных и временных возмущений при условии малости затухания последних w_i² << 1с точностью до O(w_i²) имеют место соотношения ω_S = ω_r, α_r = α_S, т.е. диапазоны частот пространственных и временных возмущений совпадают для одинаковых длин волн.

На первом этапе для заданного режима течения определяется диапазон частот эволюционирующих в пространстве возмущений. Для этого строится кривая нейтральной устойчивости для наиболее неустойчивых первых мод временных возмущений в координатах $({\mathrm{Re}}_{\delta },{\omega }_r).$ Верхний предел частотного интервала определяется горизонтальной касательной к верхней ветви кривой, ниже которой лежит область растущих возмущений. В качестве нижнего предела частот используется условный асимптотический предел при $({\mathrm{Re}}_{\delta },{\omega }_r).$ В найденном таким образом интервале с равномерным шагом Δω_r выбирается некоторое множество точек ω_rk, k = 1, 2, …, m.

На втором этапе для каждой из выбранных частот рассчитывается кривая N-фактора $N_{\omega }\left({\mathrm{Re}}_{\delta \left(x\right)}\right).$ С этой целью на прямой ω_rk = const внутри области неустойчивости, ограниченной нейтральной кривой, берется множество точек с абсциссами Re_dj, j = 1, 2, …, n. В каждой точке решается спектральная задача для наиболее неустойчивой пространственной моды возмущений. В результате получаем набор локальных инкрементов (декрементов) нарастания α_{i, j}. На их основе с использованием квадратурной формулы вычисляется интеграл

$N_{\omega }\left({\mathrm{Re}}_{\delta \left(x\right)}\right)=-2\underset{{\mathrm{Re}}_{\delta \left(x_0\right)}}{\overset{{\mathrm{Re}}_{\delta \left(x\right)}}{\int }}{\alpha }_i^{*}d{\mathrm{Re}}_{\delta },{\alpha }_i^{*}={\alpha }_i\sqrt{\frac{x\eta }{\rho U}},$ (3.7)

рассматриваемый как функция верхнего предела. Нижний предел ${\mathrm{Re}}_{\delta \left(x_0\right)}$ соответствует абсциссе точки на нейтральной кривой. Верхний предел увеличивается до тех пор, пока монотонное возрастание функции $N_{\omega }\left({\mathrm{Re}}_{\delta \left(x\right)}\right)$ не сменяется на убывание.

Рассчитанные таким образом кривые N-факторов для частот ω_r,1 < ω_r,2 < … < ω_r,k строятся на плоскости $({\mathrm{Re}}_{\delta \left(x\right)},N_{\omega })$ и находится огибающая семейства. Искомое число Рейнольдса Re_δT, определяющее начало зоны ламинарно-турбулентного перехода, находится как абсцисса огибающей при заданной величине фактора перехода N_T. В расчетах принималось N_T = 8.

4. Результаты расчетов

Положение зоны перехода рассчитывалось для двух характерных режимов течения в сверхзвуковом пограничном слое. Для сравнения были выбраны две точки в средней (режим 1) и нижней (режим 2) частях баллистической траектории посадки аппарата Pathfinder, соответствующих 66-й и 87-й секундам полета. В режиме 1, исходя из расчетов авторов [5], имеют место высокие температуры потока и развитый процесс диссоциации, влияние которого на ЛТП является предметом исследования. В то же время режим 2 характеризуется почти полным отсутствием диссоциации, и устойчивость определяется только условиями гиперзвукового полета в атмосфере углекислого газа, что также представляет определенный интерес. Для обоих режимов было выполнено сравнение со случаем совершенного газа при тех же внешних условиях.

Параметры набегающего потока в расчетных точках траектории взяты из работы [4]. В качестве граничных условий на условной верхней границе пограничного слоя использовались параметры потока за косой ударной волной на головной части аппарата в форме клина с углом β = 40°. Параметры потока за ударной волной рассчитывались на основе асимптотических формул теории косого скачка [9] при больших числах Маха.

Для режима 1 параметры потока за косым скачком уплотнения имеют следующие значения – скорость потока U_∞ = 2445.5 м/с, число Маха M = 12.6, давление p_∞ = 1549.5 Па, температура потока T_∞ = 2583 К, колебательная температура T_v∞ = 169 К. Температура поверхности спускаемого аппарата принималась равной T_wall = 2000 К.

Соответственно для режима 2 – скорость потока U_∞ = 1378 м/с, число Маха M = 6.9, давление p_∞ = 1062.85 Па, температура потока T_∞ = 946 К, колебательная температура T_v∞ = 179 К, температура поверхности спускаемого аппарата T_wall = 1100 К. В этом случае статическая температура потока существенно ниже температуры начала диссоциации углекислого газа. Вместе с тем она равна температуре возбуждения изгибной моды CO₂ [4], и газ остается термически неравновесным. Массовые концентрации компонент смеси в обоих режимах были выбраны близкими к составу равновесной атмосферы Y_1∞ = 0.9, Y_2∞ = 0.05 и Y_3∞ = 0.05. С учетом конечного времени колебательной релаксации колебательные температуры принимались равными температуре атмосферы в точках входа.

На условной границе ПС при y = δ для обоих режимов течения ставились следующие условия:

$U\left(\delta \right)=\mathrm{1,}{Y}_1\left(\delta \right)=Y_{1\infty },Y_2\left(\delta \right)=Y_{2\infty },Y_3\left(\delta \right)=Y_{3\infty },T\left(\delta \right)=\mathrm{1,}{T}_v\left(\delta \right)=T_{v\infty }/T_{\infty },$

а на поверхности спускаемого аппарата – условия прилипания и изотермической некаталитической стенки:

$U\left(0\right)=\mathrm{0,}{Y^{\text{'}}}_1\left(0\right)={Y^{\text{'}}}_2\left(0\right)={Y^{\text{'}}}_3\left(0\right)=\mathrm{0,}T\left(0\right)=T_v\left(0\right)=T_{wall}/T_{\infty }.$

Для нахождения локально автомодельных решений для рассматриваемых режимов течения система уравнений (2.2)–(2.7) приводилась к нормальной форме, для которой методом “стрельбы” с помощью процедуры Рунге–Кутты четвертого порядка на интервале [0, δ] решалась двухточечная краевая задача. Точкой “прицеливания” служила середина интервала ζ_c = δ/2, где требовалось совпадение значений вычисляемых величин с точностью до 10^–8. Расчеты проводились при следующих значениях параметров: для режима 1 – M = 12.6, $D{a}_d=0.669\times {10}^{-4},$ $D{a}_{VT}=5.32\times {10}^{-3};$ для режима 2 – M = 6.9, Da_d ≈ 0, $D{a}_{VT}=1.27\times {10}^{-3}.$ Для обоих режимов использовались значения параметров Sc = 0.5, Pr = 0.75, γ = 1.4, δ = 8 и рассчитывался альтернативный случай совершенного газа при Da_d = $D{a}_{VT}$ = 0.

На рисунке 1 для режима 1 представлены локально автомодельные профили газодинамических переменных невозмущенного пограничного слоя для термически неравновесного (Da_d = 0, Da_VT = 5.32 ×10^–3) и термохимически неравновесного газа (Da_d = 6.69 ×10^–4, Da_VT = 5.32 ×10^–3), которые сравниваются со случаем совершенного газа при прочих равных условиях.

По графикам температур на рис. 1в можно заключить, что последовательный учет энергоемких процессов колебательного возбуждения (кривая B) и диссоциации (кривая C) приводит к существенному охлаждению потока, хотя колебательная температура меняется слабо. Следует отметить, что полученные уровни абсолютных значений статической и колебательной температур по порядку величины близки к соответствующим значениям, полученным в работе [3] для того же режима на основе трехкомпонентной модели, но для баллистического спуска и с более сложной колебательной кинетикой.

 

Рис. 1. Автомодельные профили скорости U_s(ζ) (а) и плотности ρ_s(ζ) (а) смеси, массовых концентраций CO₂ Y1_s(ζ) (б) и CO Y2_s(ζ) (б), статической T_s(ζ) (в) и колебательной T_vs(ζ) (в) температур для режима 1: 1 – зависимости U_s(ζ); 2 – зависимости ρ_s(ζ); 3 – зависимости Y1_s(ζ); 4 – зависимости Y2_s(ζ); 5 – зависимости T_s(ζ); 6 – зависимости T_vs(ζ); A – совершенный газ (Da_d= Da_VT = 0); B – термически неравновесный газ (Da_d = 0, Da_VT = 5.32×10^–3); C – термохимически неравновесный газ (Da_d = 6.69×10^–4, Da_VT = 5.32×10^–3).

 

Спектральные задачи (1.11)–(1.17) для рассматриваемых режимов решались численно в среде пакета Matlab на основе метода коллокаций [21], который ранее был апробирован авторами в [19, 23].

Для определения диапазона частот эволюционирующих в пространстве возмущений рассчитывались кривые нейтральной устойчивости наиболее неустойчивых первых мод временных возмущений в координатах $({\mathrm{Re}}_{\delta },{\omega }_r).$ Полученные кривые для режимов 1 и 2 в полулогарифмических координатах представлены на рис. 2.

 

Рис. 2. Кривые нейтральной устойчивости наиболее неустойчивых первых мод временных возмущений: 1 – режим 1; 2 – режим 2; A – совершенный газ; B – термохимически неравновесный газ.

 

Проведенное здесь сравнение совершенного и термохимически неравновесных газов, а также сопоставление с нейтральными кривыми [22] для совершенного газа при близких гиперзвуковых числах Маха M ≈ 6 и M ≈ 12 показывает, что вид нейтральных кривых временных возмущений целиком определяется числом Маха ПС. В ближнем гиперзвуке [22, 23] происходит слияние областей неустойчивостей мод I и II. Уже при M ≈ 6 сохраняет индивидуальность только нижняя часть области неустойчивости первой моды. При M ≈ 12 мода I целиком поглощается модой II, что делает их неразличимыми. При этом происходит смещение частотной координаты критического числа Рейнольдса в нижнюю часть нейтральной кривой с образованием характерного заостренного носика [22].

Таким образом, для гиперзвуковых ПС получается единая кривая нейтральной устойчивости для первых наиболее неустойчивых мод. Можно заметить, что в режиме 2, где фактически имеет место только термическая неравновесность, рост критического числа Рейнольдса по сравнению с совершенным газом значительно меньше, чем было получено в [23] для близкого числа Маха. Это связано с тем, что в данном случае значение релаксационного числа Дамкелера Da_VT в определенной мере соответствовало реальным условиям полета в атмосфере Марса, в то время как в модельной задаче [23] было принято $D{a}_{VT}~O\left(1\right).$

Для оценки начала зоны перехода Re_δT было выбрано значение N-фактора N_T = 8, которое используется для течений с низким уровнем привходящих возмущений [6]. Нижние пределы Re_δ(x0) в интегралах (2.8) брались на нейтральных кривых, начиная с верхнего предела интервала частот растущих возмущений с шагом $\Delta \omega =1.25\times {10}^{-5}.$

Значения инкрементов (декрементов) пространственных возмущений α_i вычислялись на прямых ω = const с шагом ΔRe_δ = 10. Такой выбор шага позволил для многократных вычислений интегралов (2.8) использовать простую квадратурную формулу трапеций и получить достаточно гладкие кривые N-факторов как функций от Re_δ. Далее, используя связь ${\mathrm{Re}}_{\delta }=\sqrt{{\mathrm{Re}}_x}\text{},$ проводился пересчет N-факторов как функций от Re_x.

На рис. 3 для двух режимов течения представлены по несколько кривых из семейств N-факторов для наиболее неустойчивых первых мод и огибающие семейств. Кривые N(Re_x), соответствуют частотам ω от $4.25\times {10}^{-4}$ (кривые 1) до $5\times {10}^{-4}$ (кривые 7) для режима 1 и от $9.25\times {10}^{-4}$ (кривая 8) до 10^–3 (кривая 14) для режима 2. Частоты для обоих режимов возрастают с шагом $\Delta \omega =1.25\times {10}^{-5}.$

 

Рис. 3. Кривые N-факторов и положение ламинарно-турбулентного перехода: (а) – режим 1, ω = 4.25×10^–4 (1), 4.375×10^–4 (2), 4.5×10^–4 (3), 4.625×10^–4 (4), 4.75×10^–4 (5), 4.875×10^–4 (6), 5×10^–4 (7); (б) – режим 2, ω = 9.25×10^–4 (8), 9.375×10^–4 (9), 9.5×10^–4 (10), 9.625×10^–4 (11), 9.75×10^–4 (12), 9.875×10^–4 (13), 10^–3 (14); A – кривые N(Re_x) для совершенного газа; B – кривые N(Rex) для термохимически неравновесного газа; C – прямая N = 8; D₁ и D₁′ – точки перехода для режима 1; D₂ и D₂′ – точки перехода для режима 2.

 

Штрихпунктирная линия обозначает прямую N(Re_x) = 8, D₁ и ${D^{\text{'}}}_2$ – точки начала ЛТП для режима 1, D₂ и ${D^{\text{'}}}_2$ – точки начала ЛТП для режима 2 соответственно для термохимически неравновесного и совершенного газов. Числовые значения полученных критериев Re_xT даны в табл. 3.

Можно констатировать, что учет термохимической неравновесности в режиме 1 приводит к смещению начала ЛТП примерно на 9% по сравнению с совершенным газом. В режиме 2 соответствующее смещение еще меньше, особенно относительно результата, полученного в работе [23] для близкого по числу Маха режима.

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

Гиперзвуковым режимам в ПС соответствует монотонное возрастание числа Рейнольдса перехода Re_xT в зависимости от числа Маха [24]. Поэтому Re_xT в режиме 1 при M = 12.6 примерно в четыре раза больше, чем для режима 2 при M = 6.9.

Заключение

На основе e^N-метода выполнены сравнительные расчеты положения начала зоны ЛТП для точек в средней и нижней частях траектории посадки аппарата Pathfinder в атмосфере Марса. В расчетах использовалась трехкомпонентная модель термохимически неравновесной смеси CO₂–CO–O.

 

Таблица 3. Параметры точек перехода при N = 8

Точка перехода	D1	D1′	D2	D2′
RexT ×10-7	3.6378	3.9797	1.0141	1.0672
ωxT × 104	4.6455	4.4841	9.9823	9.8125

 

Было выбрано рекомендуемое для случая потоков с малым уровнем внешних возмущений значение N-фактора (N_T = 8), определяющего ЛТП. Частотный спектр пространственных возмущений находился по нейтральным кривым временных возмущений при одинаковых длинах волн.

Характерная для гиперзвуковой области форма нейтральных кривых со слиянием областей неустойчивости первой и второй мод полностью определяется числом Маха режима. На рассчитанных семействах кривых N-факторов строились огибающие, по которым при заданном N_T = 8 находились числа Рейнольдса перехода Re_xT.

В гиперзвуковом режиме 1 при M = 12.6 учет развитой термохимической неравновесности приводит к значительному снижению статической температуры газа в нижней части ПС. В результате начало зоны ЛТП сдвигается вниз по потоку примерно на 9% по сравнению со случаем совершенного газа. Для рассматриваемых гиперзвуковых режимов характерна монотонно возрастающая зависимость числа Рейнольдса перехода Re_xT от числа Маха. Поэтому Re_xT для режима 2 при M = 6.9 примерно в четыре раза меньше, чем для режима 1.

Финансирование

Исследование выполнено за счет гранта Российского научного фонда № 23-11-00027, https://rscf.ru/project/23-11-00027/.

Об авторах

Ю. Н. Григорьев

Федеральный исследовательский центр информационных и вычислительных технологий

Автор, ответственный за переписку.
Email: grigor@ict.nsc.ru
Россия, Новосибирск

И. В. Ершов

Федеральный исследовательский центр информационных и вычислительных технологий; Новосибирский государственный аграрный университет

Email: ivershov1969@gmail.com
Россия, Новосибирск; Новосибирск

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

Дополнительные файлы