助推器头部跨音速脉动压力CFD数值仿真

沈 丹，苏 虹，徐珊姝

(北京宇航系统工程研究所，北京100076)

1 引言

大型运载火箭往往采用捆绑式布局，火箭在主动段以跨声速飞行时，在箭体外形突变处，会产生局部的非定常气动力，即分离流和附体流交替的流场以及激波振荡[1-3]。 基于这个原因，助推器头部和整流罩都会在跨音速阶段产生强烈的低频脉动压力。 其中，助推器的头锥不但面临自身“锥-柱”交界面形状突变带来的跨音速效应，同时还面临着与火箭芯级之间的气动干扰，即助推器锥角处的激波在芯级发生反射。 大型捆绑火箭由于尺寸大，导致风洞模型的缩放比例较大，因此针对边界层、激波-边界层干扰等影响脉动压力的局部物理现象难以模拟。

关于火箭脉动压力的研究，针对整流罩开展的计算较多。 任淑杰等[4]通过求解流场、求解非线性声场相结合的技术，数值求解了火箭跨声速表面锥柱肩部、船尾倒锥、裙柱部区典型位置处脉动压力。 赵瑞等[5]基于5 阶WENO 格式构造隐式大涡模拟方法，对跨声速来流条件下火箭整流罩脉动压力环境进行数值模拟，拟合出干扰区脉动压力的经验公式，并提出一种脉动压力环境的工程预测流程，在保证计算效率的同时，完善并提高工程经验公式的预测精度。 NASA 早在上世纪70 年就已摸索出一套关于火箭头部脉动压力的工程预示方法，并给出航天飞机脉动压力环境[6]。 然而，这些研究并未关注助推器肩部的跨声速脉动压力以及助推器与芯级之间干扰对脉动压力的影响。

本文使用计算流体力学(Computational Fluid Dynamics，CFD)方法数值求解大型捆绑火箭跨音速非定常流场，特别关注助推器头部气动干扰区域，针对直头锥、斜头锥2 种助推器头部形状，提取头锥附近壁面压力的时间历程，获取脉动压力特性与马赫数、头锥形状等参数之间的规律。

2 数值仿真

2.1 模型及网格

在全箭流场范围内构造结构化网格，拓扑结构共分为20 个块，每个助推器独立成块。 在全箭头部用O 型网格拓扑伸展出6 个，尾部O 型网格拓扑伸展出6 个，助推器、芯级包裹6 个块(4 个助推4 块，1 个芯级1 块，缝隙1 块)，助推器外侧2 块。网格和拓扑结构如图1 所示。 计算网格数500 万～700 万，在壁面第一层网格尺寸为0.01 mm。

2.2 数值计算方法

为了提高模拟非定常流动隐式推进过程的时间精度，采用伪时间子迭代技术(Pseudo Time Sub-Iteration， τTS，也称作双时间推进技术，Dual Time Stepping)，将隐式高斯赛德尔(Lower-Upper Symmetric Gauss-Seidel，LU-SGS)格式的时间精度提高到2 阶[7]。 根据经验并考虑计算精度和效率因素，取子迭代步数为3 步。

在湍流计算方面，采用基于剪切应力输运模型(Menter's Shear Stress Transport， SST)的雷诺平均/大涡模拟( Reynolds Averaged Navier-Stocks/Large Eddy Simulation，RANS/LES)混合方法，在高频、小尺度运动占主导地位的近壁区域使用高效可靠的湍流模式，而在低频大尺度运动占优的非定常分离流动区域则使用LES 类方法[8-10]。 模拟跨声速火箭助推器头部的非定常激波/边界层干扰流场，包括流动分离、激波振荡等现象。

图1 计算网格和拓扑结构Fig.1 The mesh and topology structure

在并行计算方面，采用基于并行计算(Message Passing Interface，MPI)的数值算法，可以获得更大规模的网格分辨率，以期得到更符合实际的数值预测结果。

2.3 工况及初边值

按照运载火箭实际飞行弹道规律，挑选如表1 所示的工况作为主要的计算工况，随马赫数变化，攻角和雷诺数按照典型弹道得到。

表1 计算工况Table 1 Calculation cases

同时对直头锥、斜头锥2 种头部构型的助推器进行对比。 其中，直头锥助推器广泛应用于我国现役火箭，而斜头锥则应用于我国新一代大型运载火箭。 本文中，直头锥助推器的半锥角为15°，斜头锥助推器的最大半锥角为30°。 2种几何构型的示意图如图2 所示。 工况1～6 的来流参数均是根据火箭飞行中的典型时刻确定的。

图2 剖面示意Fig.2 Section plane position

图2同时给出了A-J 共10 个压力数据提取剖面的位置，A、B、C、D、G、H 为助推器上的6 个半圆，E、F 为助推表面的1/4 圆，I、J 为芯级表面的1/4 圆。 半圆上按顺时针方向分别在0°、45°、90°、135°、180°5 个点布置测压点，1/4 圆上按顺时针方向分别在0°、45°、90° 3个点布置测压点。

3 计算结果

3.1 流场特征

直头、斜头锥助推器的时均压力云图如图3所示。 由图可见，在跨音速时期助推器头部的自由来流马赫数小于1，气流经过头锥压缩增速，在拐角处接近音速，经过拐角后膨胀加速，在肩前部产生1 道斜激波，在肩后部产生1 道正激波的结尾激波，2 道激波前后跳动振荡。 所不同的是，直头锥的拐角角度小，气流经过拐角时产生一系列波系与芯级壁面相互作用。 而对于斜头锥，外侧拐角角度大(是斜头锥的2 倍)，气流在拐角处膨胀严重，形成较大的分离，流动在超声速附体流和亚声速分离流之间反复交替[11]。

图3 Ma=0.8 时均压力云图Fig.3 Pressure contour(Ma=0.8)

直头、斜头锥助推器的Q ＝-0.0001 等值面压力云图如图4 所示。 其中Q 为速度梯度张量第2 不变量，其等值面可以用来表现三维流场的涡结构。 由图可见，直头锥的肩部涡系对芯级表面影响很大，二者之间的干扰严重。 相对而言，斜头锥除背风面以外，助推与芯级之间少有涡结构出现，二者之间的干扰轻微。

3.2 均方根压力系数

均方根脉动压力数据也是火箭脉动压力最重要的设计指标，将提供给载荷计算所用。 均方根脉动压力系数按式(1)计算：

式中，PLU为局部非定常压力，PLS为局部非定常压力的时间平均值，T 为取样时间，q 为来流动压。 则针对上述仿真结果进行统计，得到各个工况下均方根脉动压力系数的中位数值、峰值及其所在剖面，如表2 所示。

图4 Ma=0.8 等Q 面压力云图Fig.4 Pressure contour of Q(Ma=0.8)

表2 各工况均方根脉动压力系数统计Table 2 CP-rmsstatistics

由此可见，工况1(Ma ＝0.7、a ＝4.0°)时，直头锥助推器首先迎来脉动特征。 此时虽然来流马赫数仅有0.7，但助推器肩部局部已经进入跨音速区域，跨音速激波不稳，在助推肩部振荡。 此时各剖面的脉动压力分布规律如图5 所示，图中脉动压力的高值都出现在助推与芯级间的狭缝处，是助推直头锥肩部波系与芯级之间的干扰所致。

工况2(Ma ＝0.8、a ＝6.0°)时，直头锥和斜头锥均反映出一定的脉动压力特性，且分布比较随机。

图5 工况1 直头锥均方根脉动压力分布规律Fig.5 Case 1，Prmsdistribution on positive nosecone

工况3(Ma ＝0.9、a ＝8.0°)开始，直头助推器脉动压力显著减小，而斜头助推器在此时迎来脉动压力极值，最高值达0.19821，位于D 剖面的135°位置上，是助推器头锥外侧。 显然是由助推器头锥肩部折角处产生激波振荡所主导的脉动，计算中各截面的激波振荡频率在15～50 Hz 之间。 此时各剖面的脉动压力分布如图6所示。 极值出现的截面D-135°测点上，出现了图7 所示的无量纲功率谱密度，存在33 Hz 的主频。

图6 工况3 斜头锥均方根脉动压力分布规律Fig.6 Case 3，Prmsdistribution on oblique nose cone

图7 工况3D-135°测点无量纲功率谱密度Fig.7 Case 3，Amplitude spectrum on D-135°

随着马赫数的进一步增加，2 种构型的跨音速效应均明显减少，脉动压力的均方根随之减小。对于工况5 和工况6，无论直头锥还是斜头锥，所引起的脉动压力均可忽略。

这一结果与文献[6]中的结论是相符的，即：在一般的运载火箭上升飞行过程中，对于15°锥角的锥柱交界面，在Ma ＝0.74 左右发生由激波振荡引发脉动压力极值，而对于30°锥角的锥柱交界面，则在Ma ＝0.90 左右发生由激波振动引发脉动压力极值。

从本文的计算可见，直头锥助推器先于斜头锥发生激波振荡，按照一般的运载火箭轨迹，在Ma ＝(0.7，1.0)这一飞行时段，动压的变化范围并不大，均在20 000～25 000 N/m2量级左右，即来流的能量相差不大。 但15°锥角脉动压力系数极值为0.119，30°锥角脉动压力系数极值为0.198，二者相差显著。 可见在飞行中，大的外折角出现脉动压力极值的时间靠后，但强度更强。

除膨胀转折角度的差异以外，直头锥和斜头锥还有1 个显著不同，即与芯级的干扰关系。显然，直头助推器的折角靠近芯级，而斜头助推器的平直段靠近芯级，前者与芯级表面的干扰更为显著。 在来流速度不断增加、全场完全进入超音速流动后，直头助推器头激波在芯级表面上发生反射，工况4～5 中，芯级表面(IJ 剖面)的脉动压力显著回升。 而对于斜头锥，此时助推器肩部出现脉动极值，而芯级表面却仍然十分平静。

3.3 与试验值的比对

符合工况2 和4 的斜头锥外形风洞试验在中国空气动力研究与发展中心的FL-24 风洞进行，试验使用压阻式脉动压力传感器，传感器布置在上述各截面处，采样频率为5 kHz。 将数值计算得到的均方根脉动压力系数与试验数据进行对比，如图8 所示。

对于工况2(Ma ＝0.8)，计算与试验值在量级和分布规律上比较接近。 而对于工况4，计算得到某些大峰值，这些峰值出现在背风助推器或助推器的背风面，而试验值却并没有体现。 当Ma在0.9～1.0 时，背风面出现大的激波震荡，这种震荡由激波和边界层干扰，对来流马赫数和表面边界层的细微变化十分敏感，因此在试验中未必能够捕捉到。

图8 均方根脉动压力系数计算值与试验值比较Fig.8 Comparison of calculated and experimental CP-rms values

4 结论

1)针对直头锥和斜头锥2 种不同的助推器外形，由于外折的膨胀角度不同，迎来最大均方根脉动压力系数的先后和量级不同。 直头锥在Ma ＝0.7左右率先达到再附条件并实现激波振荡，斜头锥在Ma ＝0.9 左右才达到脉动压力极值，但脉动压力系数极值更大；

2)直头助推器与芯级的干扰效应更显著，将在芯级表面由于激波反射引发脉动压力；

3)随着来流速度的提高，跨音速效应明显减少，脉动压力均方根随之减小；

4)计算和试验的均方根脉动压力数值符合良好，验证了本文中数值算法在工程应用中的可行性。

参考文献(References)

[1]NASA.Buffeting during atmosphere ascent[R].NASA SP-8001，1964.

[2]鲍贤栋.对两种头部外形气动噪声的预测分析[J].强度与环境， 1980(2)：31-42.Bao X D.Prediction of aerodynamic noise on two different nose cone shape [J].Structure ＆ Environment Engineering，1980(2)：31-42.(in Chinese)

[3]沈丹，吴彦森，岑拯.芯级与助推器头部气动干扰流场数值模拟[J].导弹与航天运载技术， 2013(6)：42-46.Shen D，Wu Y S，Cen Z.Numerical simulation of aerodynamic interaction characteristics between the rocket and a booster's nosecone [J].Missile and Space Vehicles， 2013(6)：42-46.(in Chinese)

[4]任淑杰，张收运，闫桂荣.基于RANS/NLAS 的火箭跨声速脉动压力环境预示[J].固体火箭技术，2011，42(4)：418-422.Ren S J，Zhang S Y， Min G R.A prediction of fluctuation pressure conditions with transonic rocket by RANS /NLAS method [J].Journal of Solid Rocket Technology，2011，42(4)： 418-422.(in Chinese)

[5]赵瑞荣， 吉利， 王超，等.跨声速条件下火箭整流罩表面脉动压力环境分析与预测[C]/ /中国力学大会， 上海，2015：158.Zhao R R， Ji L， Wang C H，et al.Analysis and prediction of fluctuating pressure environment on the surface of rocket fairing under transonic condition [C]/ /CCTAM， Shanghai，2015：158.(in Chinese)

[6]Kenneth J P， Robertson J E.Prediction of space shuttle fluctuating pressure environments， including rocket plume effects[R].NASA-CR-124347.

[7]赵慧勇，乐嘉陵.双时间步方法的应用分析[J].计算物理， 2008， 25(3)：253-258.Zhao H Y， Le J L.Application analysis on dual-time stepping[J].Chinese Journal of Computational Physics， 2008， 25(3)：253-258.(in Chinese)

[8]肖志祥.湍流模型在复杂流场数值模拟中的应用[J].计算物理， 2003，20(4)：335-340.Xiao Z X.Applications of turbulence models in simulation of complex flows [J].Chinese Journal of Computation Physics，2003， 20(4)：335-340.(in Chinese)

[9]肖志祥，陈海昕.采用RANS/LES 混合方法研究分离流动[J].空气动力学学报， 2006，24(6)：218-222.Xiao Z X， Chen H X.Simulation of separation flows with RANS/LES hybrid methods [J].Acta Aerodynamica Sinica，2006， 24(6)：218-222.(in Chinese)

[10]肖志祥， 陈海昕.激波/边界层干扰流动数值模拟格式的应用研究[J].西北工业大学学报，2004， 22(6)：800-806.Xiao Z X， Chen H X.The simulation of shock/boundary interaction flows with several temporal and spatial schemes [J].Journal of Northwestern Poly Technical University， 2004， 22(6)：800-806.(in Chinese)

[11]操小龙.锥-柱体外形脉动压力及抖振载荷响应研究[J].战术导弹技术，2010(1)：67-72.Cao X L.The Research of fluctuating pressure and buffeting load on the cone-column body[J].Tactical Missile Technology， 2010(1)：67-72.(in Chinese)