激波和边界层相互干扰是高超声速飞行中的常见现象,其导致的局部过高热流是飞行器的安全隐患;此外,它导致的边界层分离会使飞行器气动性能改变.因此,激波与边界层相互干扰是飞行器设计中必须考虑的重要因素,Harvey等^[1]、Candler^[2]和Gaitonde^[3]分别对激波和边界层相互干扰现象进行了综述.

高超声速双锥流场存在复杂的激波与边界层相互干扰现象,国内外对其有大量研究^[4-20].国内如张静等^[4]研究了不同离散格式对高超声速双锥流场计算中激波分辨率与黏性分辨率计算结果的影响;石晓峰^[5]针对热化学非平衡效应对高超声速流场中激波反射、激波边界层相互作用的影响做了较为系统的研究;田浩等^[6]采用包括量热完全气体、平衡气体、化学非平衡气体在内的不同气体模型对不同工况下高马赫数绕流流场中的分离流动进行了研究;国外如Druguet等^[12]测试了不同的物理化学模型,包括化学反应速率、平衡常数、振动能弛豫以及输运系数等,然而其仅计算了一种来流条件且与实验结果的比较并不理想.Nompelis等^[13-15]针对双锥流动受到振动非平衡与化学反应的影响展开了一系列工作,发现只有当来流焓值不是很高且实验气体为氮气时的计算结果较为合理,当实验气体为空气时,双锥实验热流的计算结果与实验吻合不佳,且随着来流焓值的增大变得更差.由于来流中通常存在热化学非平衡,且其程度无法准确测量,为了排除来流条件的不确定性,美国嘉斯班大学布法罗研究中心于2013年开展了双锥风洞实验,实验中自由来流中不包含热/化学非平衡现象^[16].针对这组实验,Hao等^[17]研究了不同振动与离解耦合模型的影响,发现两个不同模型的计算结果相差不大.Kianvashrad等^[18]研究了化学反应的影响,发现无论低焓来流还是高焓来流,化学反应对表面热流与压力的计算结果无明显影响.Kieweg等^[19]和Ray等^[20]通过不确定性分析,认为计算与实验的差别有可能是来流条件的误差或空间分布引起.

本文从计算模型的选取角度来分析计算与实验差别,研究不同热化学反应与输运模型对双锥流动计算结果的影响.热化学模型包括完全气体模型、热力学非平衡结合冻结流以及化学非平衡流模型、热力学平衡化学非平衡模型等.虽然完全气体与热平衡模型并不适用于该问题,但能揭示内能和化学能的分布对计算结果的影响.此外,壁面催化效应的影响、Wilke/Blottner/Eucken输运模型以及更复杂的Gupta/SCEBD输运模型的影响也在研究范围内.

本文采用自主开发的CH3D~(chemically-reacting hypersonics3D)软件^[21],该软件求解全N-S方程,支持结构化多块网格、非结构化网格、以及混合网格等多种网格拓扑结构,支持热力学与化学非平衡的求解.

本文使用文献^[22]的5组元空气化学反应模型. 质量守恒方程为 \begin{equation} \label{eq1} \frac{\partial \rho _s }{\partial t} + \nabla \cdot \left( {\rho _s u} \right) + \nabla \cdot \left( {\rho _s v_s} \right) = \omega _s(1) \end{equation} 式中,$\rho _s$为第$s$组元的密度,$u$为混合物的对流速度矢量,$v_s$为第$s$组元的扩散速度,$\omega _s$为第$s$组元的化学反应生成率. 动量守恒方程为 \begin{equation} \label{eq2} \frac{\partial \rho u}{\partial t} + \nabla \cdot (\rho u\times u) + \nabla p + \nabla \cdot \tilde {\tau } = 0(2) \end{equation} 本文采用振动、电子能平衡假设,并相应引入振动、电子温度$T_{ ve}$. 振动、电子能守恒方程为

$$\begin{align*} &\frac{\partial E_{v}}{\partial t} + \nabla \cdot \left( {E_{v}u} \right) + \nabla \cdot \sum\limits_{s = 1}^{ns} {v_s} E_{{v},s} + \nabla \cdot q_{v}= \\ &\qquad Q_{T-v} + \sum\limits_{s = 1}^{ns} {\omega _s } e_{{v},s}(3)\end{align*} $$

式中$Q_{T-v}$是平动、转动能传递速率^[23]. 总能守恒方程为

$$\begin{align*} & {\frac{\partial E}{\partial t} + \nabla \cdot (E + p)u + \nabla \cdot \sum\limits_{s = 1}^{ns} {v_s} \rho _s h_s }+\\ &\qquad \nabla \cdot \left( {q + q_{v} } \right) + \nabla \cdot (u \cdot \tau ) = 0(4)\end{align*} $$

多组分混合气体的输运模型通常采用Wilke混合率^[24]、Blottner模型^[25]和Eucken公式^[26]这套方法. 使用该输运模型时,通常假设Lewis数或Schmidt数为常量. 本文除了Wilke/Blottner/Eucken模型,还采用了物理上更为复杂的SCEBD~(self-consistent binary diffusion model)模型^[27]来代替Lewis数或Schmidt数为常数的假设,采用Gupta混合率^[28]来代替Wilke混合率.

1.2.1 Wilke/Blottner/Eucken模型

Wilke/Blottner/Eucken模型被广泛用于计算多组分混合气体输运系数. 其中,Wilke混合率是Chapman-Enskog关系的一种简化,采用Blottner拟合计算各组分的黏性系数, 采用Eucken关系式计算热传导系数. 在计算扩散系数时,往往采用单一扩散系数,并且假设Lewis数或Schmidt数为常量. 本文假设Schmidt数等于0.7,扩散系数计算方法如下 \begin{equation} \label{eq5} D_s = D = \frac{\mu }{\rho S_{c}}(5) \end{equation} 其中,$\mu $是混合气体黏性,$S_{c}$是Schmidt数. Fick定律用于质量扩散的计算 \begin{equation} \label{eq6} J_s =-\rho D_s \Delta Y_s(6) \end{equation}

1.2.2 Gupta/SCEBD模型

Gupta等发展了另一种混合气体的输运模型,该模型比Wilke/Blottner/Eucken物理上更复杂, 考虑了混合气体各组分的碰撞积分. 模型中计入了重粒子与电子的影响,但本文计算中不考虑空气的电离,因此仅考虑重粒子的影响. 混合气体的黏性系数为 \begin{equation} \label{eq7} \mu = \sum\limits_s {\frac{m_{e} \gamma _{e} }{ \sum\limits_r {\gamma _r } \varDelta _{{e},r}^{(2)} \left( {T_{ve} } \right)}}(7) \end{equation} 混合气体的平动热传导系数为 \begin{equation} \label{eq8} \kappa _{t} = \frac{15}{4}k_{B} \sum\limits_s {\frac{\gamma _s }{ \sum\limits_r {a_{s,r} } \gamma _r \varDelta _{s,r}^{(2)} \left( {T_{{t}r} } \right)}}(8) \end{equation} 转动热传导系数为 \begin{equation} \label{eq9} \kappa _r = k_{B} \sum\limits_{s = mol} {\frac{\gamma _s }{ \sum\limits_r {\gamma _r } \varDelta _{s,r}^{(1)} \left( {T_{{t}r} } \right)}}(9) \end{equation} 振动、电子激发态热传导系数为 \begin{equation} \label{eq10} \kappa _{vel} = k_{B} \frac{C_{v,ve} }{R}\sum\limits_{s = mol} {\frac{\gamma _s }{ \sum\limits_r {\gamma _r } \varDelta _{s,r}^{(1)} \left( {T_{{t}r} } \right)}}(10) \end{equation} 重粒子的碰撞项,$\varDelta _{s,r}^{(1)} $和$\varDelta _{s,r}^{(2)} $的形式为 \begin{equation} \label{eq11} \varDelta _{s,r}^{(1)} = \frac{8}{3}\left[ {\frac{2M_s M_r }{\pi RT\left( {M_s + M_r } \right)}} \right]^{1 / 2}\pi \bar {\varOmega }_{s,r}^{(1,1)}(11) \end{equation} \begin{equation} \label{eq12} \varDelta _{s,r}^{(2)} = \frac{16}{5}\left[ {\frac{2M_s M_r }{\pi RT\left( {M_s + M_r } \right)}} \right]^{1 / 2}\pi \bar {\varOmega }_{s,r}^{(2,2)}(12) \end{equation} 其中$\pi \bar {\varOmega }_{s,r}^{(1,1)} $和$\pi \bar {\varOmega }_{s,r}^{(2,2)} $是碰撞积分,参见文献^[29,30].

在计算扩散系数时,SCEBD模型采用了组分碰撞积分. 重粒子的扩散系数为 \begin{equation} \label{eq13} D_s = \left( {1-\frac{\omega _s }{\omega }} \right)\left( {\sum\limits_{r \ne s} {\frac{p_r }{pD_{s,r} }} } \right)^{-1}(13) \end{equation} 其中 \begin{equation} \label{eq14} D_{s,r} = \frac{k_{B} T_{{t}r} }{p\varDelta _{s,r}^{(1)} \left( {T_{{t}r} } \right)}(14) \end{equation} 最终,质量扩散流通量为 \begin{equation} \label{eq15} J_s =-\frac{pM_s D_s }{RT_{{t}r} }\nabla \left( {\frac{p_s }{p}} \right) + Y_s \sum\limits_r {\frac{pM_r D_r }{RT_{{t}r} }} \nabla \left( {\frac{p_r }{p}} \right)(15) \end{equation}

本文忽略空气的离解反应,高温空气化学反应采用文献^[22]的5组元(N$_{2}$,O$_{2}$,NO, N和O)反应模型,振动、离解耦合效应采用文献^[31]模型. 本文壁面催化模型采用完全催化壁及完全非催化壁两种.

采用有限体积方法来离散方程组. 对于对流项,采用修正的Steger-Warming流通量分裂^[32],结合迎风MUSCL TVD格式将精度扩展至三阶. 流通量与雅可比矩阵均采用计算单元表面处的流场数值,通过单元中心处的流场变量插值得到. 在迎风格式中,流通量被分解为正负两项 \begin{equation} \label{eq16} f_{i + 1 / 2} = f_{i + 1/2}^ + \left( {\phi _{i + 1/2}^ + } \right) + f_{i + \frac{1}{2}}^-\left( {\phi _{i + 1/2}^-} \right)(16) \end{equation} 其中$\phi $为流场变量向量,本文取$\phi = \{ p,c_s ,u,v,w, T_{{t}r} ,T_{ve} \}$. 三阶 MUSCL 差分格式为

$$\begin{align*} &\phi _{i + 1/2}^ + = \phi _i-\frac{1}{6}\Big\{ {\min od \left[ {\phi _{i-1}-\phi _i ,2\left( {\phi _i-\phi _{i + 1} } \right)} \right]} +\\ &\qquad 2\min od \left[ {\phi _i-\phi _{i + 1} ,2\left( {\phi _{i-1}-\phi _i } \right)} \right]\Big\}(17)\end{align*} $$

$$\begin{align*} &\phi _{i + \frac{1}{2}}^-= \phi _{i + 1}-\frac{1}{6}\Big\{ {\min od \left[ {\phi _{i + 2}-\phi _{i + 1} ,2\left( {\phi _{i + 1}-\phi _i } \right)} \right]} +\\ &\qquad 2\min od \left[ {\phi _{i + 1}-\phi _i ,2\left( {\phi _{i + 2}-\phi _{i + 1} } \right)} \right]\Big\}(18)\end{align*} $$

美国嘉斯班大学布法罗研究中心于2013年开展了双锥风洞实验,实验共分为6个工况,本文均进行了数值模拟, 6个工况来流条件见表1. 若无特殊说明,本文后续计算时均默认采用双温度热力学非平衡、5组元Park化学非平衡、 Gupta/SCEBD输运模型和壁温300 K的完全催化壁模型.

首先进行网格无关性测试,通过数值模拟工况1分别考察了310$\times $200,620$\times $400和1240$\times $800共三套网格,三套网格中壁面垂直方向的$y^ + $都取值为1.图1显示了三套网格的计算得出的壁面压力和热流分布,结果表明620$\times $400网格已经网格收敛. 图中显示,计算所得分离区大小要小于实验结果,透射激波位置靠前,除此之外与实验结果符合较好.

工况2至工况6的壁面热流与壁面压力分布的计算与实验结果比较分别见图2(a)与图2(b). 图中可以看出,计算与实验结果总体上符合较好,但是计算出的分离区大小比实验要小,该结果与文献^[17]一致.

本文将在6个工况中分别选取代表低焓和高焓情形的工况1与工况4,从热化学模型与输运模型入手,探讨影响流场计算结果的主要因素.

针对工况1所代表的低焓工况,此时忽略流场中的化学反应,考虑不同热力学模型的影响, 包括热力学非平衡的双温模型、热力学平衡的单温模型和比热比等于1.4的完全气体模型. 完全气体模型采用1.4的比热比,分子内能分布仅有平动态和转动态;热力学非平衡模型增加考虑了振动激发, 能量从平动、转动态传到振动态;在热力学平衡模型中平动、转动态与振动态达到平衡. 随着平动、转动态的能量逐渐转移到振动态,体现在宏观流场上即平动、转动温度降低,分离区大小减小.

图3给出了不同热力学模型对工况1中壁面热流和压力影响. 从计算结果来看,完全气体模型计算得到的分离区略大于实验结果,热力学非平衡模型的分离区其次, 热力学平衡模型计算的分离区最小,两模型分离区均小于实验结果. 虽然完全气体模型计算的分离区大小及峰值与实验符合最好,但是在第一锥度完全气体模型的热流最低, 与实验结果偏差最大.

图4给出了采用不同输运模型计算所得壁面热流分布.可以看出,采用不同输运模型计算所得热流分布差别不大.图5给出了工况1流场的马赫数云图与温度云图,从马赫数云图中可以清楚看到双锥拐角处的分离区,以及分离区下游的透射激波,从温度云图可以看到透射激波和弓形激波形成的高温激波层.对于工况1中的低焓来流条件,流场温度低于5500 K,在该温度范围Wilke/Blottner/Eucken模型与Gupta/SCEBD模型计算所得结果相似,相比之下, 输运模型对计算结果的影响明显小于热力学模型对计算结果的影响.

对于工况4所对应的高焓来流,由于流场中存在明显的化学反应,此时不仅需要考虑热力学模型的影响,同时需要考虑冻结/非平衡化学反应模型的影响.图6给出了不同化学反应模型对工况4流场中壁面热流与壁面压力计算结果的影响,图7则给出了不同化学反应模型所对应流场的马赫数云图.

从图6中可以看出：冻结流的计算结果偏离实验结果,分离区大小偏大.一方面,流场中化学反应使内能转移至化学能,宏观上表现为流场温度降低、分离区相应减小.另一方面,从图6(b)的压力分布曲线来看,冻结流的分离区下游的压力更高,需要更大的分离区来平衡这个逆向压差.从图7可以看出,冻结流计算的激波强度更大,导致分离区附近斜激波和弓形激波脱体距离更远,更厚的激波层从尺度上来说会导致更大的分离区.说明高焓流场中存在明显的化学反应,是否考虑化学反应对计算结果影响显著.

针对工况4,图8给出了3种不同热力学模型的计算结果,完全气体模型偏离实验结果较大,与工况1相比,热平衡与热非平衡两模型之间的计算结果差别变小.图9为热非平衡模型针对工况1与工况4计算得到的平动、转动温度与转动、振动温度差值的云图.可以看出,两工况均在第一个锥体附近存在较大的温度差值,即存在明显的热非平衡效应,但此处两工况振动温度低,因此热非平衡影响较小.但在弓形激波层内,由于更高的来流密度和更强的激波,此处工况4的热非平衡要远小于工况1,因此热非平衡对工况4的结果影响较小.

这里计算热平衡和非平衡模型时,均采用了非平衡化学反应.比较图8中完全气体结果与图6中的化学冻结流结果,发现两者差别不大,因此认为图8中完全气体模型结果的巨大差别是由化学反应引起.总的来说,针对工况4,化学反应模型是影响结果的主要因素.

图10给出了工况4中壁面催化对壁面热流计算结果的影响.可以看出,壁面催化对第一锥面影响很小,其影响在第二锥面更为显著,这是由于第一锥面处温度低,化学反应程度低,而第二锥面附近的弓形激波层内温度高,化学反应显著.与实验结果相比,考虑壁面催化效应的计算结果与实验结果吻合的更好.高焓情形,必须考虑壁面催化的影响.

图11给出了不同输运模型对热流计算结果的影响.与图4中低焓工况结果相比,对于高焓工况,Gupta/SCEBD模型计算的热流及其峰值要明显高于Wilke/Blottner/Eucken模型.针对工况4,流场中的最高温度高于10 000K,在该温度范围内,两输运模型计算结果之间的差别显著,且Gupta/SCEBD模型与实验结果符合更好.

图12为使用Gupta/SCEBD模型时流场的Schmidt数云图.可以看出,流场中的Schmidt数并非常量,而是介于0.57$\sim $0.75之间.Gupta/SCEBD模型与Wilke/Blottner/Eucken模型的主要区别之一在于后者假设Schmidt数为常量,很明显,当流场温度较高时,该假设并不成立,此时应采用Gupta/SCEBD模型.

本文针对双锥风洞实验中的激波和边界层相互干扰中的热化学和输运模型影响进行数值分析.从低焓至高焓状态,本文模拟了所有6个实验工况.计算结果所示分离区较实验结果偏小,除此之外,壁面压力和温度的计算结果与实验符合较好.

针对工况1所代表的低焓状态,通过比较不同热力学模型的计算结果发现,分子内能的分布对分离区大小以及壁面热流峰值和压力的峰值的计算结果影响显著.随着平动、转动态的能量逐渐转移到振动态,体现在宏观流场上就是平动、转动温度降低,分离区大小减小;低焓状态,Gupta/SCEBD模型与Wilke/Blottner/Eucken模型的计算结果相近,热力学模型对计算结果影响显著.

对于工况4的高焓状态,化学反应对计算结果的影响显著.当考虑化学反应时,内能转移至化学能,计算所得分离区显著变小.工况4的热平衡与热非平衡所得计算结果之间的差别较工况1要小,这是由于工况4密度更高,弓形激波层内非平衡效应小.高焓状态下壁面催化的影响也更为显著,考虑壁面催化时的计算结果与实验结果吻合的更好.对于不同的输运模型,针对第二锥段的壁面热流与热流峰值,Gupta/SCEBD模型的计算结果要明显高于Wilke/Blottner/Eucken模型,且与实验吻合的更好.根据Gupta/SCEBD模型所得流场中Schmidt数在0.57~0.75之间变化,常用的Schmidt数为常数的假设并不合理,计算时应采用变Schmidt数模型.