连续流失效对近空间飞行器气动特性的影响

黄 飞，张 亮，程晓丽，沈 清

（中国航天空气动力技术研究院，北京 100074）

0 引 言

传统航天飞行器在80km以上的高空大气环境中高速飞行，其流动表现出稀薄效应，而对于未来采用升力体或乘波体构型［1］的近空间（20km～100km）高超声速飞行器，由于其外形特征多为扁平状升力面，具有尖化前缘、复杂气动外形、长时间飞行和高升阻比等气动构型和飞行轨道特征，飞行环境介于稠密大气层和稀薄大气层之间，其局部小尺寸部件在中等高空，甚至在较低的高空就会出现稀薄气体效应。在稀薄流区，要实现飞行器的飞行、机动和控制，必须了解它们的气动力特性，同时为飞行器防、隔热设计提供气动热数据。而随着飞行高度的变化，摩擦力、热流等气体动力学特性较之连续流区发生很大的变化。这些变化对于传统的钝头体再入飞行器轨道影响较小，因此以往很少细致考虑稀薄气体效应的影响，大余量的设计方法对于解决这类外型和轨道的飞行器已足够，而对大升力面、尖化前缘复杂外形中低空飞行的稀薄流问题还有待研究。因此，中低Kn数的尖化前缘外形的气动力、气动热的正确预测对近空间空域飞行器的设计和工程应用都具有重要意义。

在低空问题研究中，由于大气密度高，常采用传统的CFD方法进行流场的模拟分析，而飞行器飞行过程中，随着飞行高度的变化，其所处的外部流动环境是一个渐变的过程，那么连续流假设在什么条件下失效？连续流计算方法对近空域飞行器气动力、热预测的偏差有多大？就成为目前急需解决的问题。国外已开展了许多关于连续流失效性判据［2－6］的研究工作，而对近空域稀薄效应对飞行器气动特性变化的定量分析还不多见［7－9］。

本文首先采用了圆柱算例对四种滑移模型的有效性进行了对比分析，最后针对未来近空域飞行的高超声速飞行器大升力面的构型特征，采用典型外形，选用适当的滑移模型和无滑移NS方法对近空间飞行器过渡流效应引起的气动特性变化进行了对比分析，定量地给出了稀薄气体效应对气动特性的影响规律。

1 计算方法

1.1 N－S方程

在笛卡儿坐标系下，采用层流三维非定常可压缩NS方程：

其中，Q为守恒变量；F、G、H为三个方向的无粘矢通量；Fv、Gv、Hv为三个方向的粘性矢通量。

1.2 滑移边界类型

本文采用了四种滑移边界模型，分别定义如下：Type1为经典的Maxwell滑移边界类型；Type2为Gokcen［10］滑移边界；Type3为 Lockerby［11］提出的壁函数修正法；HS为分子自由程采用硬球模型简化计算时的Maxwell滑移边界。

Type2（对应文献［9］中CFD（2））：

Type3（对应文献［9］中CFD（3））：

HS：

σ动量调节系数，α能量调节系数。

上述结果表明，在非高斯干扰明显的环境中，原有算法性能会急剧下降甚至失效，但基于Alpha稳定分布的新算法可以有效完成杂波抑制，为接下来的动目标检测奠定基础。

1.3 无量纲参数定义

偏差分析量定义如下：

其中，qnoslip为无滑移NS方程热流结果，qslip2为滑移模型Type2的热流计算结果。

2 滑移模型的适应性分析

本文滑移模型的算例验证采用文献［9］中圆柱的计算条件，来流气体为氩气，来流马赫数分别为Ma＝10、25，来流温度T∞＝300K，圆柱直径D＝304.8mm，基于直径D＝2R的努森数Kn为0.002、0.05、0.25。本算例分别采用了 Type1、Type2、Type3、HS、No＿slip等方法对圆柱扰流进行了计算分析，并与文献［9］中的DSMC结果进行了对比。此算例中的粘性系数计算方法见文献［9］，对应于DSMC的VHS模型。

图1 滑移与非滑移下Ar的流场结构（Ma＝10）Fig.1 Flow field with different methods（Ma＝10）

图1（a）为Kn＝0.002时的流场压力云图，从中可以看出，滑移方法与非滑移方法在近连续流区所得到的流场结构非常一致，滑移流方法能够准确捕捉激波结构和尾涡结构。图1（b）为Kn＝0.05时两种方法所得的流场压力云图和速度型，从图中可以发现由于在此努森数下存在较弱的稀薄效应，两种方法所得结果存在以下三点不同：滑移流方法所得的激波厚度比非滑移流方法稍厚；在物面处，滑移流方法所得的物面速度出现明显的滑移速度，使得物面处速度型不同于非滑移流所得结果；在尾部区域，两种方法所得的压力等值线明显存在不同。稀薄效应引起的以上三点不同将很有可能影响物面的热流、压力等分布。为了进行定量分析下面给出了不同滑移模型、非滑移方法、DSMC方法的物面气动特性分布。

图2分别给出了马赫数等于10时，Kn＝0.002、0.05、0.25时不同滑移模型与DSMC方法所得的热流系数分布曲线，其中图2（a）为本文计算结果，图2（b）为文献［9］计算结果。从图中可以看出，在较弱稀薄效应下，即努森数Kn＝0.002时，不同滑移模型所得热流系数分布与DSMC结果吻合较好，当努森数Kn增加到0.05、0.25时不同模型所得结果之间的差距开始明显显现；随着努森数的增加，无滑移NS方法所得热流分布都一致的高于滑移方法及DSMC方法所得结果，并且从图中可以发现，类型2的滑移模型即使在较大努森数Kn＝0.05、0.25时也表现出了与DSMC相一致的结果，其它模型所得结果在较大努森数下都一致高于DSMC方法所得结果，由此可见，类型2的滑移模型对较大努森数下具有更好的适应性，其它模型的结果都将高估热流而使得防热设计趋于保守。通过与文献结果对比可以发现，本文与文献结果热流分布吻合很好。

图3分别给出了Kn＝0.002、0.05、0.25时本文与文献［9］采用不同滑移模型与DSMC方法所得的压力系数分布曲线，从图中可以看出，当努森数Kn＝0.002、0.05时，不同滑移模型所得表面压力系数分布与DSMC所得结果非常一致，当努森数增加至0.25时，各种模型所得压力分布都一致高于DSMC方法所得结果，这主要是因为强烈的稀薄效应使得NS方法所得激波结构明显薄于DSMC方法所得激波结构，NS方法产生的激波压缩性更强，使得波后压力较高。从图2热流系数分布和图3压力系数分布对比结果可以看出，相对于热流而言，压力对稀薄效应敏感性较弱。通过与文献结果对比可以发现，本文与文献结果压力分布同样吻合很好。

图4、图5分别给出了马赫数等于25时，Kn＝0.002、0.05、0.25时本文采用不同滑移模型所得的热流系数与压力系数分布曲线与文献［9］DSMC方法的结果比较，从图中可以看出，马赫数等于25时的热流、压力分布规律与马赫数等于10时的分布规律不同之处仅在于量值上的不同，其它规律与马赫数等于10时的结果相似，类型2的滑移模型在较大马赫数下也同样具有较好的适应性。

综合以上计算分析可以得出以下结论：通过本文计算结果与文献［9］所得结果的对比可以发现，本文结果与文献结果吻合较好，本文所发展的代码具有一定的可靠性；滑移模型2的适应范围较其它模型更广，它能够在较大努森数下取得相对较为满意的热流结果；在较大努森数下，滑移模型所得热流、压力结果都一致高于DSMC结果，使得防热设计趋于保守；相对于热流而言，压力对稀薄效应的敏感性更弱；经过数值模拟还可以发现，模型2虽然适应性广，但相对于其它模型而言，计算时较为耗时，收敛性较差，对网格要求较高，如果在飞行器防热设计条件允许下，可以考虑采用其它几种滑移模型。

图2 不同滑移模型在不同努森数下的物面热流系数分布（Ma＝10）Fig.2 Heat transfer rate on the surface at different Knnumbers（Ma＝10）

图3 不同滑移模型在不同努森数下的物面压力系数分布（Ma＝10）Fig.3 Pressure coefficient on the surface at different Knnumbers（Ma＝10）

图4 不同滑移模型在不同努森数下的物面热流系数分布Fig.4 Heat transfer rate on the surface at different Knnumbers（Ma＝25）

图5 不同滑移模型在不同努森数下的物面压力系数分布Fig.5 Pressure coefficient on the surface at different Knnumbers（Ma＝25）

3 梯形翼近空间气动特性的计算分析

由于滑移模型2在稀薄气动特性预测方面的适应范围更广，预测结果更为接近DSMC结果，故本文针对马赫数15，飞行高度h＝50km～80km的梯形翼主要采用了模型2进行了稀薄气动特性的计算分析。飞行攻角为10°，壁温为500K，翼前缘直径为30mm，翼根弦长约1.96m，翼尖弦长约0.63m，翼展长为0.5m，前缘后掠角约20°。

图6为本文所采用的计算网格示意图，图7给出了梯形翼不同位置处的截面示意图，其中Z＝5mm为梯形翼对称面部位，Z＝300mm为翼中部位。为了进行定量分析，以下给出了截面位置处不同方法所得结果的差异性比较分析。

图8、图9分别为攻角10°下截面Z＝5mm、Z＝300mm处不同方法所得的表面热流分布及偏差曲线结果，从上图热流分布结果可以发现，在相对较低的飞行高度，滑移方法与无滑移方法所得结果较为一致，随着飞行高度的增加，两种方法所得结果的差距逐渐增大。从图9偏差曲线图的定量结果可以明显看出，由于迎风面气体压缩，背风面气体膨胀，使得滑移的影响首先出现在背风面，故无滑移方法与滑移方法之间的差距在迎风面明显低于背风面，并且随着飞行高度的增加，稀薄效应越明显，滑移方法与非滑移方法所得的热流结果之间的偏差逐渐增大。从图中还可发现，此偏差最大值主要发生在翼前缘背风面膨胀区附近，当此种飞行器从50km至80km飞行时，在飞行器翼面大面积上的偏差约在5%至15%之间变化。

表1为不同高度下滑移方法和无滑移方法所得气动力特性的结果比较，从中可以看出，随着飞行高度的增加摩阻系数、轴向力系数、法向力系数都有所增加，其中摩阻增加较大，故而可以发现，升阻比随着高度的增加逐渐减小。从不同方法的偏差结果可以看出，随着飞行高度的增加，气体稀薄度增加，稀薄效应增强，传统无滑移假设下的NS方法已不能正确描述壁面边界条件，故而两种方法所得结果的偏差量逐渐增大。

图6 计算网格示意图Fig.6 Grid in the computation region

图7 不同位置处的截面示意图Fig.7 Slice at different stations

图8 不同位置处的热流分布结果Fig.8 Heat flux at different stations

图9 偏差分析结果Fig.9 Error at different stations

表2为不同方法所得峰值热流随高度的变化，从中可明显看出，随着高度的增加，气体密度降低，飞行器前缘的峰值热流降低较快，而采用滑移与无滑移两种方法所得的结果偏差则随着高度的增加而增大，这主要是由于高度的增加使得稀薄度增大，传统连续流假设下的无滑移方法不再有效。

表1 不同方法下的气动力特性结果Table 1 Aerodynamics with different methods

表2 不同方法下的峰值热流结果Table 2 Peak heat flux with different methods

4 结 论

高超声速飞行器飞行过程中经历了不同大气环境密度的变化，本文针对中低空域大升力面飞行的高超声速飞行器的结构特征，采用典型外形研究了连续流失效后不同滑移模型对表面气动特性的预测能力，最后给出了大升力面梯形翼近空间气动特性的变化规律。结果表明：

（1）本文所得规律与文献结果吻合很好，所建立的模拟手段具有一定的可靠性。

（2）滑移模型2比其它模型适应范围更广，能够在较大努森数下取得相对较为满意的结果。

（3）在较大努森数下，除了滑移模型2以外，其它模型所得热流、压力结果都一致高于DSMC结果，使得防热设计趋于保守。

（4）相对于热流而言，压力对稀薄效应的敏感性较弱。

（5）在本文所研究的尺寸下，当大升力面梯形翼的飞行高度从50km增加至80km时，滑移方法与无滑移方法在翼面大面积上所得的热流偏差量约在5%至15%之间变化；而驻点处峰值热流的偏差量约在1.6%至14.5%之间变化；相应的升阻比偏差结果约在0.06%至14.38%之间变化。

［1］李怡勇，沈怀荣.发展近空间飞行器系统的关键技术［J］.装备指挥技术学院学报，2006，17（5）.

［2］BIRD G A.Breakdown of translational and rotational equilibrium in gaseous expansions［J］.AIAAJournal，1970，8（11）：1998－2003.

［3］TIWARI S.Coupling of the boltzmann and Euler equations with automatic domain decomposition［J］.Journal ofComputationalPhysics，1998，144：710－726.

［4］CAMBEROS J A，SCHROCK C R，McMULLAN R J.Development of continuum onset criteria with direct simulation Monte－Carlo using boltzmann＇s H－theorem：review and vision［R］.Proceedings of the 9th AIAA／ASME Joint Thermophysics and Heat Transfer Conference，San Francisco，California，June 2006.

［5］BOYD I D，CHEN G，XANDLER G V.Predicting failure of the continuum fluid equations in transitional hypersonic flows［J］.PhysicsofFluids，1995，7（1）：210－219.

［6］WANG W L.A hybrid particle／continuum approach for nonequilibrium hypersonic flows［D］.［Ph.D.Thesis］.The University of Michigan，2004.

［7］BOYD I D，PADILLA J F.Simulation of sharp leading edge aerothermodynamics［R］.AIAA Paper 2003－7062.

［8］LOFTHOUSE A J，BOYD I D，etc.Effects of continuum breakdown on hypersonic aerothermodynamics［R］.AIAA Paper 2006－993，2006

［9］lOFTHOUSE A J，SCALABRIN L C，BOYD I D.Velocity slip and temperature jump in hypersonic aerothermodynamics［R］.AIAA Paper 2007－0208，2007.

［10］GOKCEN T，MACCORMACH R W.Nonequilibrium effects for hypersonic transitional flows using continuum approach［R］.AIAA Paper 1989－0461，1989.

［11］LOCKERBY D A，REESE J M，GALLIS M A.Capturing the knudsen layer in continuum－fluid models of nonequilibrium gas flows［J］.AIAAJournal，2005，43（6）：1391－1393.