钝锥再入飞行器黏性干扰效应建模

2022-10-14江定武毛枚良

气体物理 2022年5期

王沛, 李锦,2, 江定武,2, 毛枚良,2

(1. 中国空气动力研究与发展中心计算空气动力研究所, 四川绵阳 621000;2. 中国空气动力研究与发展中心空气动力学国家重点实验室, 四川绵阳 621000)

引言

对高超声速飞行器而言, 随着飞行高度增加, Mach数增加, Reynolds数减少, 表面边界层厚度增加, 形成位移效应, 导致飞行器有效几何外形发生变化, 影响边界层外部的无黏流场, 而外部流场的变化也会影响壁面附近的流场。边界层与外部压缩无黏流场之间这种强耦合作用形成黏性干扰效应。这种效应会引起当地压力、摩擦阻力以及热传导等参量的剧烈变化, 进而影响飞行器的整体气动性能, 是高超声速飞行时不可忽视的复杂物理效应。

早期的黏性干扰研究主要采用理论分析的方法, 基于高超声速边界层相似律和高超声速无黏流的压力关联公式[1]: Bertram[2]结合摩阻与压力梯度的近似关系式, 研究了零厚度带攻角平板和三角翼的黏性干扰效应对压力、摩阻和气动力的影响; Whitfield等[3]完成了冷壁细长钝锥的黏性干扰效应理论及与实验的对比研究; Stollery[4]针对指数律外形给出了凹陷体、压缩拐角和膨胀拐角等若干实例的黏性干扰流场结果; Probstein等[5]评估了钝锥外形上横向曲率在弱黏性干扰区的影响; Yasuhara[6]研究了3/4幂次律轴对称外形上较弱至中等强度黏性干扰区中的相似流动。但是早期这些研究工作基本上只针对简单外形, 复杂外形面临更加复杂的流动机理, 这些黏性理论可能失效, 难以得到正确的结果。

20世纪70年代, 美国对航天飞机OV102进行了一系列飞行实验, 获得了大量的气动力飞行数据[7], 给出了与气动特性参量相关的第1、第2、第3黏性干扰参数, 并得到第3黏性干扰参数数据关联结果更为准确的结论[8]。基于此, 国内学者在黏性干扰效应方面取得了一些研究成果。庄逢甘等[9]针对OV102外形进行了初步研究, 提出了较为成熟的超声速飞行器黏性干扰效应的工程计算方法; 龚安龙等[10]对轨道器再入阶段的黏性干扰效应进行了研究, 考察了黏性干扰相关参数的关联特性; 毛枚良等[11]对升力体简化外形的高超声速绕流流动进行了模拟, 认为黏性干扰导致的气动力系数的改变量与黏性干扰参数在特定范围内具有良好的线性关系; 李维东等[12]结合参考温度方法提出了一种能够考虑黏性干扰效应的高超声速乘波体气动性能的工程预测方法; 陈坚强等[13]针对OV102外形以及高升阻比复杂外形黏性干扰效应, 提出基于Mach数和第3黏性干扰参数的关联参数和公式, 具有实现高超声速气动力数据天地相关的能力; 张益荣等[14]针对典型高超声速翼身组合体外形, 建立了完全气体条件下纵向气动力系数的黏性干扰模型, 并进行了不确定度量化分析; 龚安龙等[15]针对高升阻比面对称飞行器, 建立了黏性干扰增量与黏性干扰参数之间的线性关系式, 并采用Mach数效应、真实气体效应等进行修正建立了天地换算公式; 罗长童等[16]利用专业化机器学习算法, 得到了不同风洞试验结果共同遵守的不变规律, 为建立黏性干扰等物理效应天地换算方法提供了新思路。

国内外的研究现状表明, 传统黏性干扰的研究主要基于连续性假设成立的N-S方程, 但是真实飞行条件下随着速度和高度进一步变大, 气体分子远离平衡态, 稀薄效应变得不可忽略, 导致飞行器的气动特性相比于连续流有了较大差异, 而这一区域内目前仍然缺少合适的黏性干扰效应预测模型, 这种模型的建立对求解器提出了新要求。统一气体动理论(unified gas kinetic scheme, UGKS)是近些年来比较热门的一种跨流域多尺度求解方法, 在连续流区域可以得到精确的N-S解, 在稀薄区域也可以得到与DSMC一致的结果, 已经广泛应用于从低速到高速、从连续流到稀薄流的全速域全流域的流场模拟中[16]。

基于此, 本文采用UGKS求解器, 以典型外形钝锥为研究对象。预测其在高空(70～110 km)、高Mach数(≥10)条件下的气动特性, 并借助黏性干扰效应的理论成果, 将Mach数、攻角、黏性系数等参数与UGKS有黏解和Euler无黏解的差量进行关联, 建立了气动特性增量预测模型。这样后续进行流场预测时, 仅需要计算Euler无黏解, 就可通过该增量模型快速高效得到相应黏性条件下的气动力参数, 这极大简化了气动特性的求解过程。最后选取相关性曲线与Pearson乘积矩相关系数进行了数据相关性分析, 采用部分新状态验证了模型预测结果能够满足工程应用的精度要求。

1 数值方法

1.1 无黏解算器

控制方程选用一般曲线坐标系下的三维可压缩Euler方程, 基于有限体积法进行离散, 采用隐式LUSGS方法进行迭代求解, 具体计算过程参见文献[17]。

1.2 黏性解算器

控制方程采用Shakhov模型方程, 其具体形式为

其中,f是气体在空间x, 速度u, 以及时间t的分布函数。τ为碰撞时间, 与黏性和压力有关,f+可以表示为

其中, 第1项g是Maxwell分布函数, 第2项是基于原始BGK模型方程的修正项, 目的是为了得到合理的Prandtl数,c和q分别为分子热运动速度矢量和热流矢量。

方程有如下积分解

e-t/τf(x-ut,0,u)

该积分解由初始分布函数和平衡态分布函数构成, 通过松弛时间τ将分子的碰撞和自由输运过程耦合起来, 从而实现分布函数在稀薄流和连续流区域的自动调节, 其具体表达式见文献[18]。

宏观守恒量Q以及应力P、热流q等变量均可以通过在速度空间上对f求矩得到

其中, dΞ=dudvdw为相空间的体积元。

基于上述模型构造了一套隐式三维大规模并行软件, 并进行了多个算例的验证, 详见文献[19-20]。

2 结果及建模

相较于返回舱, 钝锥适用于需要大横向机动能力、良好的机动性、低着陆误差和低减速负载的任务以及高进入速度和稀薄大气层的任务。本文采用钝锥外形进行高超声速黏性干扰模型的研究。钝锥球头半径取为600 mm, 半锥角10°, 全长L=3 600 mm, 力系数参考面积为1.0 m2, 力矩参考点位于质心, 坐标为(2 418 mm, 0 mm, 0 mm)。图1给出了钝锥的计算网格示意图, 表1给出了计算状态。

表1 计算状态Table 1 Computation states

2.1 流场及气动力特性分析

图2给出了Mach数10、攻角10°时有黏(100 km)以及无黏两个状态的流场压力云图和速度矢量图, 图中变量均为无量纲结果。100 km时, 黏性边界层很厚, 黏性干扰效应明显改变了物面压力和空间流场结构。

(a) Inviscid flow

(b) Viscid flow(100 km)图2 速度矢量及压力云图Fig. 2 Velocity vector and pressure contour

图3给出了两个高度下对称面及物面局部Knudsen数分布。图中局部Knudsen数变化范围为0.01～100, 而且随着高度进一步增加, 激波内部、壁面附近的局部Knudsen数进一步增加, 稀薄效应更加明显。这些局部区域内传统基于连续性假设的N-S方程已不再适用, 模拟的流场已不再准确。为了更好地捕捉流场细节, 预测真实流场, 需要采用跨流域求解器。

(a) H=70 km

(b) H=80 km图3 对称面及物面局部Knudsen数分布Fig. 3 Local Kn number on symmetry plane and wall

图4给出了Mach数为10、攻角10°时不同高度对称面固壁压力分布比较, 在球头与锥的交接位置(x/L=0.167)存在明显的压力分布拐点。迎风面大面积区域(锥身)压力分布随着高度的增加呈增加趋势。背风面大面积区域压力分布随高度增加呈现先增后减的趋势, 但是量值较迎风面大约要小1个量级, 这与物理实际基本一致。

(a) Windward

(b) Leeawrd图4 对称面压力比较Fig. 4 Comparison of pressure on symmetry plane

图5, 6分别给出了Mach数为10、攻角为10°和30°时壁面压力和摩擦阻力对轴向力系数的贡献随黏性干扰参数的变化。随着黏性干扰参数的增加, 轴向力系数增加, 压力、黏性的贡献也在增加, 而且相比压力项, 黏性项的作用占主导, 且对比发现攻角变化对轴向力的增量影响较小。

图5 攻角10°时轴向力系数特性Fig. 5 Axial force coefficient characteristics at α=10°

图6 攻角30°时轴向力系数特性Fig. 6 Axial force coefficient characteristics at α=30°

图7给出了部分状态的俯仰力矩特性。高度100 km(V′∞～0.33)及以下, 俯仰力矩变化缓慢, 符号为正, 压心在质心之前, 飞行器处于静不稳定状态。高度进一步增加到110 km, 攻角为10°， 20°情况下, 俯仰力矩符号变为负数, 压心后移到质心之后, 飞行器处于静稳定状态。且Mach数对俯仰力矩的影响较小。

图7 俯仰力矩系数特性Fig. 7 Pitching moment coefficient characteristics

2.2 模型建立

以往的研究成果表明, 第3黏性参数能够较好地将来流条件与气动特性参量关联起来[13], 所以本文基于第3黏性干扰参数建立黏性干扰预测模型。该参数表达式为

其中,T′表示边界层内的参考温度。从式中可以看出, Mach数一定的情况下, 高度增加, Reynolds数减小, 黏性干扰参数增加。为了更直观地表征黏性干扰参数与气动参量之间的关系, 所有的对比图均采用黏性干扰参数为横坐标。

图8给出了10°攻角对称面上几个典型流向位置处的压力改变量的变化情况。迎风面, 高度低于110 km时, 压力改变量与黏性干扰参数线性关系比较好。背风面, 基本看不到较好的线性关系。这说明, 即便是钝锥这种简单的外形, 压力改变量与黏性干扰参数之间也不能呈现线性关系。

(a) Windward

(b) Leeward图8 压力增量随黏性干扰参数的变化Fig. 8 Variation of pressure increment with the viscous interaction parameter

以轴向力和俯仰力矩系数为例, 给出黏性干扰模型的分析建立过程。可以类似地得到法向力系数预测模型。

图9给出了轴向力和俯仰力矩系数变化量随第3黏性干扰参数的变化, 这里的系数变化量表示UGKS求解器得到的有黏解与Euler方程求解得到的无黏解的差量。

(a) Windward

(b) Leeward图9 力系数增量随黏性干扰参数的变化Fig. 9 Variation of force coefficient increment with the viscous interaction parameter

从图中可以看到, 高度小于100 km时, 轴向力系数变化量与第3黏性干扰参数呈较好的线性关系, 但是随着高度的增加, 则逐渐偏离线性关系。因而, 本文采用二次多项式进行拟合, 并参照黏性干扰模型的推导过程[13], 初步给出了力系数的黏性干扰模型

并进一步假设

2.3 模型评估

采用45°相关线对结果进行初步分析。图10给出了黏性干扰模型预测数据(以model表示)与数值模拟计算数据(以UGKS表示)的相关性曲线, 不同攻角和Mach数条件下数据基本分布在相关线附近, 可见数据间相关性较好。

(a) ΔCA

(b) ΔCm图10 黏性干扰模型与数值模拟预测结果相关性Fig. 10 Correlation between viscous interaction model and numerical simulation results

根据上述相关性曲线, 采用相对正交距离考察黏性干扰模型数据拟合精准度

其中,di为数据点到相关性曲线的正交距离,xr为数据点在相关性曲线上投影对应的横坐标, 示意图见图11。图中横纵坐标分别为干扰模型与UGKS模拟得到的力/力矩系数增量。轴向力和俯仰力矩系数的相对正交距离计算结果如图12所示。对于轴向力系数, 高度大于70 km时, 模型预测的力系数变化量的相对偏差都低于3%; 70 km时, 相对偏差最大可达14%, 但是此时变化量的量值(约为0.236)仅为基准量(Euler无黏解, 约为1.18)的1/5。基准量叠加变化量得到的轴向力系数预测值与UGKS计算值偏差大约仅为3%。俯仰力矩系数变化量的相对正交距离具有类似的情况, 在变化量较小的时候, 相对正交距离甚至可达28%。但是此时的变化量量值比基准量量值(Euler无黏解)低1～2个量级。

图11 正交距离示意图Fig. 11 Illustration of orthogonal distance

(a) ΔCA

(b) ΔCm图12 力系数相对正交距离Fig. 12 Relative orthogonal distance of the force coefficient

采用统计学中Pearson乘积矩相关系数进行进一步的相关性分析

最后, 为了初步评估模型的准确性, 采用UGKS求解器和黏性干扰模型预测了Mach数20、攻角10°, 高度分别为100 km和110 km时的气动力增量, 具体结果及相对偏差见表2, 本文建立的模型预测的轴向力增量与UGKS求解器模拟结果偏差小于1%, 俯仰力矩增量的预测误差在10%以内。初步说明模型预测的气动增量与UGKS黏性求解器的结果基本一致, 采用黏性干扰模型可以较好地实现不同Mach数下的数据相关, 可以为工程实践提供一定的依据。