返回舱高雷诺数再入过程底部流动稳定性1)

朱德华 沈清 杨武兵

(中国航天空气动力技术研究院,北京 100074)

引言

飞船返回舱、再入弹头、高速导弹等都存在飞行稳定性问题.早在20 世纪60 年代研制阿波罗返回舱时人们就认识到底部流动稳定性是一个影响飞行稳定性的关键因素.底部流动失稳会产生非定常振荡的气动力,会造成钝体飞行器的非定常的振动,严重情形下会导致飞行失稳.高速钝体绕流底部流动失稳机理的认知将会对飞行器的稳定性控制提供基础理论支撑.

一直以来,由于传统的数值模拟方法数值耗散过大,往往不能获得底部流动失稳特征,而风洞试验方面试验模型需要支撑机构,模型支撑对底部流动影响很大,造成底部流动失真.因此人们在底部流动失稳现象观察与机理认识等方面一直存在理论认识上的不足和偏差.

国外一直在尝试通过风洞试验的方法研究高速底部流动稳定性.在1968 年Moseley 等[1]就针对阿波罗返回舱的飞行稳定性问题进行了研究,认识到了底部流动不稳定是导致飞行不稳定的主要因素.他们通过在背风面添加肋板或稳定翼的方法来消除底部流动的不稳定行为,增强了返回舱的飞行稳定性.

近期国外在钝体绕流底部附近的剪切层失稳方面取得了研究进展.美国NASA 的Danehy 等[2]在2006 年发现返回舱肩部剪切层会存在层流、转捩以及湍流的形态特征,并且会对底部流动产生影响.

2014 年,Combs 等[3]再次针对返回舱外形进行了研究,其发现返回舱剪切层最不稳定的状态是出现在较小的攻角下(0◦和12◦).在攻角0◦和12◦情形下剪切层为转捩状态,而在攻角24◦情形下剪切层是层流状态.Danehy 等[2]的研究由于试验模型底部存在支撑机构,底部流动整体情况以及失稳现象尚未获得.

针对钝体外形的底部流动的试验研究还有一些类似的代表性工作如文献[4-5],他们发展的试验测量技术使得认识底部流动失稳现象成为可能.

2009 年,美国CUBRC 的MacLean 等[6]针对返回舱外形进行了层流与湍流状态的数值计算和试验研究.试验是在LENS-I,LENS-II 激波风洞进行的,其在返回舱表面布置了大量的测压和测热传感器,希望用定量测量方式研究底部流动失稳.数值计算方法采用的是US3D 代码中的DES 方法.数值计算也证实了模型支架(圆形尾支撑、方形尾支撑以及圆形侧支撑)同样会对底部流动结构产生很大影响,无法正确模拟底部流动失稳.

2011 年,Brock 等[7]利用RANS 及RANS/LES 方法对高速返回舱的底部流动进行了数值模拟,获得了其底部流动的层流到湍流的数值结果.但是由于研究对象底部存在支撑,底部流动失稳特征被破坏.

针对钝体外形底部流动的数值模拟研究还有一些代表性的工作如文献[8-12],研究表明采用高精度的数值模拟方法可以获得底部流场的精细结构,有助于对底部流动失稳机理的认知.

2012 年,朱德华等[13]针对高速阿波罗返回舱底部流动稳定性开展了研究,研究发现了在低雷诺数的情形阿波罗返回舱底部流动结构存在不稳定性.

在钝体绕流的国内研究方面,2019 年涂佳黄等[14-15]在低速低雷诺数情形下的圆柱和方柱绕流的尾流干扰研究方面取得了一定的进展,但在高速高雷诺数情形下的钝体绕流的尾流相关研究较少.

综上所述,高速高雷诺数钝体绕流底部流动稳定性的研究无论是在风洞试验还是数值模拟方面均存在一定的困难,其底部流动失稳现象以及失稳机理仍然存在理论认识上的不足和偏差.随着高精度数值模拟方法的进步以及大规模并行计算能力的提升,使得高速高雷诺数钝体绕流底部流动现象认识和失稳机理的研究成为了可能.

本文将采用多区并行隐式大涡模拟方法针对返回舱类钝体外形,在高雷诺数再入条件下,研究其底部流动现象以及失稳机理,为返回舱类钝体外形外部扰动因素对底部流动的作用机理及其稳定性控制提供基础理论支撑.

1 控制方程和数值方法

控制方程为曲线坐标系下三维可压缩Navier-Stokes 方程[16],采用的是基于有限差分的高精度多区并行隐式大涡模拟方法求解该方程.对于空间离散,无黏项采用五阶WENO 格式[17],黏性项采用六阶中心差分格式.时间离散采用三阶TVD Runge-Kutta 方法[18].

壁面处采用等温无滑移边界条件,由于坐标变换在体轴上出现奇性轴,轴上点的流动参数由奇性轴周围点平均获得,远场为无反射边界条件.

需要注意的是奇性轴附近点的计算需要利用奇性轴上的值,而奇性轴上的通量为零,无法进行特征分裂,本文采取的处理方法为将五阶WENO 格式在轴方向降阶为二阶守恒型NND 格式[19],另外两个方向保持五阶WENO 格式.

2 数值模拟结果分析

计算条件选取参照的是Brock 等[7]的试验条件.具体的计算条件为(Re以1 m 为特征长度)

计算模型选取的是缩比的类猎户座返回舱外形,具体的尺寸如图1 所示.

图1 计算模型尺寸图(mm)Fig.1 Sketch of Orion mode(mm)

高速钝体绕流底部流动物理现象的认识可以为数值模拟研究提供正确的指导,因此本文利用已知的钝体绕流研究经验,初步给出了基于类猎户座返回舱外形高速高雷诺数绕流的基本流动形态(图2),为后续数值模拟研究提供一定的理论支撑.

图2 返回舱外形基本绕流形态Fig.2 Base flow patterns of return capsule

基于上述高速高雷诺数钝体绕流的基本流动形态,选取结构网格开展高精度数值模拟工作.类返回舱钝体外形具有轴对称特性,C 型网格划分是比较容易实现的,但在高马赫数,高雷诺数条件下其绕流流动结构复杂,需要在复杂流动现象区域进行局部网格加密,单一的C 型网格很难实现局部加密,因此本文选取C 型网格和H 型网格混合的方法来实现网格的局部加密.为了保证边界层的模拟精度,计算网格在壁面处第一层网格尺度为10µm,=0.8,针对类猎户座返回舱采用两套网格进行数值模拟,一套为粗网格网格量约为1800 万,一套为细网格网格量约为3300 万(主要在流向和法向进行加密,展向网格点数均为121 个).针对典型流动结构进行局部加密的类猎户座返回舱网格示意图如图3 所示.

图3 类猎户座网格示意图Fig.3 Sketch of grid

2.1 钝体绕流底部流动的网格效应

本节针对类猎户座返回舱外形不同网格下的大涡模拟结果进行分析,揭示其底部流动稳定性的网格效应.

图4 显示的是粗细两套网格同一时刻的0.1U∞速度云图以及瞬时流线.

图4 返回舱外形绕流流线分布图Fig.4 Skin friction line around the return capsule

图4 返回舱外形绕流流线分布图(续)Fig.4 Skin friction line around the return capsule(continued)

从图4 可以看出二者分离起始位置基本一致,肩部剪切层起始位置也基本一致,再附(鞍点S)位置随着网格的加密向返回舱底部靠近,底部物面处的半鞍点S′与再附鞍点S的连线是不稳定的,在剪切失稳的作用下已经出现了振荡特征.

图5 显示的是同一无量纲时间t=1.3 时刻数值模拟结果.其中返回舱中心对称面是用密度梯度显示的,尾迹发展区以及远尾迹区的涡结构是用速度梯度第二不变量Q值来描述的,这里Q=3000.

图5 返回舱外形绕流瞬时形态Fig.5 Instantaneous flow patterns around the capsule

从图5 可以看出,细网格肩部剪切层失稳靠前,细网格发卡涡结构已经出现在了尾迹发展区,而粗网格发卡涡结构出现在远尾迹区.发卡涡结构的出现和破碎行为预示着转捩的发生,可以预测继续加密网格会使发卡涡结构的出现略微提前但不会移动到分离区域,因此继续加密网格来捕捉尾迹发展区的转捩行为已非必要.当然远尾迹区的湍流行为加密网格还是十分必要的,再压缩波会受到远尾迹区的湍流影响,从图5(a)和图5(b)的再压缩波角度的对比中可以看出由于细网格转捩更早,流动更早的向湍流发展因此其再压缩波的角度更小.

最后给出类猎户座返回舱外形不同网格分布下头部对称线上无量纲压力(p=¯p/ρu2∞)的分布对比曲线图,如图6 所示.从图中可以看出细网格(grid2)相比粗网格(grid1)激波捕捉能力略有提升,驻点处的无量纲压力值稍有增加,细网格下驻点压力为0.935,无黏球头绕流驻点压力为0.93[20](其中M∞=6.41状态为高精度差值结果).可以看出本文提出的奇性轴处理方法在驻点上的计算结果是精确的.

图6 头部对称线上的无量纲压力分布Fig.6 Pressure distribution on the symmetrical axis

本文的研究只限于底部以及尾迹发展区的流动稳定性,出于计算量考虑针对返回舱外形绕流研究采用的是细网格.后续将针对远尾迹区继续加密网格开展研究.

2.2 类猎户座返回舱外形底部流动失稳特征

前文借鉴高速钝体绕流底部流动研究经验初步给出了类猎户座返回舱外形基本流动形态,下面分析数值模拟结果与其进行相互验证.

头部区域主要出现的物理现象是激波和驻点.图7 中返回舱表面显示的是压力云图,可以看出,物面压力分布呈现轴对称形态,并且随着计算时间推进头部对称性保持很好,证实了本文采用的奇性轴处理方式的合理性.

图7 返回舱外形绕流瞬时流动结构及监测点分布Fig.7 Instantaneous flow patterns and monitor point around the return capsule

底部区域主要出现的物理现象是流动分离,在雷诺数超过临界值后本质上是结构不稳定,其在扰动的存在下极易出现非对称结构失稳,出现周期性以及非周期性的结构振荡特征,并且流场结构失稳与流动失稳存在相互影响.

肩部区域主要出现的物理现象是剪切层.图7中的中心对称面是采用密度梯度来显示的,可以看出,高马赫数气流经过弓形激波后变为亚声速气流,亚声速气流经过返回舱肩部由于其外形的收缩特征导致气流出现膨胀加速,即出现膨胀波,气流经过膨胀波后达到超声速状态,此时返回舱底部由于压力梯度的存在出现流动分离,流动分离内部气流速度很低为亚声速状态.这种流动形态的存在导致返回舱肩部出现明显的剪切区域即剪切层,剪切层本质上是不稳定的,在扰动的作用下其会出现失稳特征.从图中可以看出,剪切层经过短暂的平直段后开始出现基频涡,此时剪切层对流马赫数Mc=(u1−u2)/(a1+a2)近似为0.8 左右,存在明显的压缩性效应.

尾迹发展区以及远尾迹区主要出现的物理现象是转捩、湍流以及类卡门涡街.从图7 中速度梯度第二不变量Q值的分布可以看出肩部剪切层失稳直接导致了在尾迹发展区出现了发卡涡结构并随着流向距离的发展迅速破碎,尾迹流动开始出现转捩并向湍流化发展.

综上所述,数值模拟获得的类猎户座返回舱绕流符合高速高雷诺数钝体绕流的基本特征.图7 中黑色圆圈示意位置为压力信号监测点,在下文中用于对其基本特征形成机制开展详细分析.

2.2.1 底部流动结构失稳

类猎户座返回舱外形的底部流动存在大面积的流动分离,其流动拓扑结构不同时刻具有不同特征.中心对称线附近流动的拓扑结构主要是由流动分离形成的,底部物面处的半鞍点S′与再附鞍点S的连线本质上是不稳定的,在水平和垂直扰动的作用下其会出现拉伸和振荡特征,导致局部的拓扑结构发生变化.

图8 显示的是时间间隔t=1.4 ∼1.7 内其局部拓扑结构的变化规律.可以看出图8(a)中拓扑结构主要是SNS连线结构,这是中心对称线处SS连线结构受扰动作用形成的一种拓扑结构.图8(b)中拓扑结构主要是SS连线结构,其在垂直扰动的作用下已经偏离中心对称线.图8(c)中拓扑结构又重新变为了SNS连线结构,至此形成了一个周期的振荡过程,其无量纲振荡频率f≈1/0.3 ≈3.3.

图8 底部流线的瞬时拓扑形态Fig.8 Topological structure patterns of base flow

为了确认上述流动结构失稳现象,任意选取中心对称面上底部物面半鞍点S′与再附鞍点S的连线附近一点作为压力信号监测点(参见图7 监测点1),其无量纲压力随任选无量纲时间(t=¯t/(L/u∞),其中L为特征长度为1 m)间隔内的分布如图9 所示.

图9 监测点1 任选时间间隔内的压力分布Fig.9 Pressure distribution of monitor point 1

针对此压力信号进行傅里叶分析获得其频谱特征,从图10 中可以看出其存在无量纲频率为3.3左右的振荡频率(1/t),这与拓扑结构振荡频率基本吻合.

图10 基于监测点1 压力信号的频谱特征Fig.10 Spectrum on pressure at the monitor point 1

2.2.2 肩部流动失稳

肩部剪切是高速钝体绕流普遍存在的物理现象,其在高雷诺数情形下极易出现失稳.任意选取肩部剪切层开始出现失稳以及失稳过程中位置作为压力信号监测点(参见图7 监测点2 和3),其压力随时间的分布如图11 和图12 所示.

图11 监测点2 任选时间间隔内的压力分布Fig.11 Pressure distribution of monitor point 2

图12 监测点3 任选时间间隔内的压力分布Fig.12 Pressure distribution of monitor point 3

针对上述两点的压力信号进行傅里叶分析获得其频谱特征,从图13 和图14 中可以看出二者均存在一个无量纲频率为3.3 左右的拓扑结构振荡频率,这是底部分离造成的剪切层的整体振荡频率,而另一占优频率约为10 左右,这是剪切层自身失稳产生的,对应着剪切层的基频.随着剪切层的空间发展,更多更高频率的激发代表着基频涡演化出了更多尺度的涡,剪切层失稳加剧.

从底部流动结构失稳和肩部剪切失稳的分析中可以看出,二者在各自主导区域内均具有各自的特征,相互作为扰动源推动各自的失稳历程.

图13 基于监测点2 压力信号的频谱特征Fig.13 Spectrum on pressure at the monitor point 2

图14 基于监测点3 压力信号的频谱特征Fig.14 Spectrum on pressure at the monitor point 3

2.2.3 尾迹区耦合失稳

在尾迹发展区可以预见上述两种失稳模式将会产生复杂的相互作用,出现耦合失稳,导致尾迹区出现转捩、湍流以及类卡门涡街等复杂物理现象.任意选取尾迹发展区一点作为压力信号监测点(参见图7监测点4)分析其频谱特征.

从图15 可以看出,底部结构失稳与肩部剪切失稳相互作用后各自的占优频率仍然存在,整个频谱呈现出宽频特征,更多高频能量的激发意味着转捩的发生.另外需要注意的是两种失稳模式作用产生了约为0.95 的无量纲占优频率,这可能对应着类卡门涡街失稳频率.

为了验证存在类卡门涡街物理现象,在远尾迹区任意选取再附激波附近一点作为压力信号监测点(参见图7 监测点5)分析其频谱特征(图16).

图15 基于监测点4 压力信号的频谱特征Fig.15 Spectrum on pressure at the monitor point 4

图16 基于监测点5 压力信号的频谱特征Fig.16 Spectrum on pressure at the monitor point 5

与尾迹发展区的频谱特征曲线对比可以看出,无量纲频率约为0.95 的占优频率仍然存在且能量更高.为了直观的认识此类卡门涡街振荡现象,将不同时刻远尾迹区中心对称面上的涡量分布显示在图17 中.从图中可以看出,高频的小涡被低频流动结构带动形成类卡门涡街形式的振荡流动现象,振荡的频率约为1.0 左右与频谱分析结果基本一致.

卡门涡街的振荡频率公式为f=S t×(U∞/D),无量纲参数S t在高雷诺数情形下约为0.27,U∞为来流速度,D为绕流物体直径.将本文相关参数无量纲化代入可得f=1.06,这与分析获得的频率基本吻合,可以认为远尾迹区的振荡流动现象为类卡门涡街振荡现象.

2.3 钝体绕流底部流动失稳机制

图17 远尾迹湍流区涡结构振荡特征Fig.17 Oscillation characteristics of vortex structure in the far wake turbulent region

综上所述,高雷诺数状态下的类猎户座返回舱底部流动与低雷诺数状态下的阿波罗返回舱底部流动一样均存在流动不稳定问题,但二者在底部流动失稳模式和流动状态方面存在明显差异.高速高雷诺数状态下的返回舱类钝体绕流的底部流动稳定性主要存在两类失稳模式,一类是肩部剪切区域具有明显压缩性特征的剪切失稳模式,另一类是底部分离区域具有振荡特征的流动结构失稳模式.这两种模式会在尾迹发展区以及远尾迹区出现耦合失稳,促使尾迹发展区出现转捩行为,远尾迹区出现湍流行为以及类卡门涡街振荡行为.而高速低雷诺数状态下的返回舱类钝体绕流底部流动稳定性主要存在的是流动结构失稳模式,远尾迹区是层流振荡行为.

困扰人们的钝体绕流肩部高热流问题,底部阻力问题,钝体飞行时的振动问题等均严重依赖于以上两类失稳模式及其耦合效应的正确模拟.

本文研究指出,在模拟高速高雷诺数状态下的钝体绕流问题时,底部分离的起始位置的正确捕获十分重要,这决定了钝体外形肩部剪切层的形成.另一方面是肩部剪切层失稳行为必须正确捕捉,其决定了分离涡再附点位置.尾迹发展区的转捩行为以及远尾迹区的湍流行为也需要准确模拟,其影响的是再压缩波的产生以及强度.

需要注意的是,本文的研究中采用的是隐式大涡模拟方法,针对返回舱类钝体外形的绕流数值模拟并未加入任何人为的扰动信息,此种模拟方法可以用来模拟静音风洞或者真实飞行条件下的钝体绕流问题.对于噪声风洞的钝体绕流问题,需要在此方法的基础上考虑扰动信息的引入,后续将开展此项研究工作.

3 结论

本文采用大涡模拟方法,针对高雷诺数返回舱类钝体外形再入过程的底部流动稳定性开展了研究.细致的刻画了类猎户座返回舱外形的底部流动形态以及稳定性特征.从肩部剪切失稳、底部流动结构失稳,尾迹发展区以及远尾迹区的耦合失稳等多个角度分析了类猎户座返回舱外形的底部流动失稳机制.研究发现了其底部流动稳定性主要存在两类失稳模式即肩部剪切失稳模式以及底部流动结构失稳模式,二种模式存在耦合效应;同时在远尾迹湍流区存在类卡门涡街形式的振荡行为.这些结论为理解返回舱外部扰动因素对底部流动的作用机理及返回舱稳定性控制提供了基础理论支撑.

具体结论如下:

(1)高速高雷诺数状态下的钝体绕流底部流动稳定性均主要存在两类失稳模式,一类是肩部剪切区域具有明显压缩性特征的剪切失稳模式,另一类是底部区域具有振荡特征的流动结构失稳模式.

(2)两类失稳模式存在耦合效应,远尾迹区的湍流流动在耦合效应的作用下会出现类卡门涡街振荡行为.

(3)底部流动形态以及失稳特征的正确模拟有助于解决钝体绕流肩部高热流,底部阻力,钝体飞行时的振动等问题.