面向新一代固体火箭的自主可控仿真技术

张 兵，吴万同，沈 治，尹宇辉，陈泽栋，高 波

(1. 中国运载火箭技术研究院，北京 100076；2. 北京宇航系统工程研究所，北京 100076)

0 引 言

2022年12月9日,由中国运载火箭技术研究院抓总研制的捷龙三号(SD-3)商业运载火箭在黄海海域成功发射,将14颗卫星准确送入预定轨道,圆满完成了首次飞行任务。

捷龙三号火箭立足于当前与未来的卫星发射市场,以成本为第一约条件,从立项研制到飞行试验成功仅用时20个月。在刚性商业化需求下,通过全面数字化转型推动了型号研制,实现了“高质量、高效率、高效益”的三高发展目标。

通过仿真试验替代实物试验是火箭研制数字化转型的有力抓手,型号队伍结合捷龙三号火箭特点与北京宇航系统工程研究所在自主可控工业软件方面的技术积淀,取消全箭气动试验、全箭模态试验、部分分离试验和部分静力静热试验,缩短研制周期近一年,解决了成本、效率、品质之间的矛盾[1]。

本文聚焦捷龙三号火箭工程研制中替代全箭风洞试验的仿真技术,描述仿真方法及其具体工程应用。首先简述了中国运载火箭技术研究院自主可控气动仿真软件的主要算法,之后分别阐述了气动特性仿真、对流换热仿真和大头颈比整流罩的脉动压力仿真的工程应用,最后给出总结。

1 自主可控仿真软件主要算法

本节先阐述可压缩湍流模拟的控制方程,介绍仿真所用的雷诺平均模拟与大涡模拟的湍流模型,随后阐述多刚体运动流场仿真方法,包括任意拉格朗日-欧拉框架和刚体动力学方程耦合求解方法,最后阐述边界的高保真度格式,以及面向未来应用需求的数据驱动增强模型。

1.1 控制方程

本文的流体仿真使用了两种方法:雷诺平均模拟与大涡模拟,根据流体力学研究的习惯,两者可沿用同一套变量和上下标符号。其控制方程为

(1)

(2)

(3)

(4)

(5)

式中:等效黏性系数满足μe=μ+μt,湍流黏性系数μt通过湍流模型封闭,湍流扩散系数Dt通过梯度扩散假设封闭。

1.2 湍流模型

在雷诺平均模拟中,湍流黏性系数通过剪切应力输运(Shear stress transport, SST)模型封闭[2]:

(6)

式中:α*为参数;F为混合函数;|S|为流动拉伸率张量的模。湍动能控制方程考虑可压缩性修正:

(7)

式中:α1和α2为膨胀可压缩性修正系数;β*为结构可压缩性修正系数[3-4]。

大涡模拟的湍流黏性通过壁面自适应局部涡黏(Wall-adapting local eddy-viscosity, WALE)模型封闭[5]:

(8)

式中:Δs代表网格尺度;Sd满足

(9)

(10)

1.3 任意拉格朗日-欧拉求解框架

在任意拉格朗日-欧拉框架下,有限体积网格单元的中心与界面均随时间移动,将网格单元移动速度记为uc,则网格单元热物理量的控制方程[6]为

(11)

(12)

(13)

上述控制方程的符号和变量与1.1节一致。若网格静止,任意拉格朗日-欧拉框架的方程回归1.1节的湍流控制方程形式。

1.4 刚体动力学方程的耦合求解

通过质心位移与绕质心旋转获得刚体运动规律,惯性系的刚体质心动力学方程为

(14)

式中:ar为加速度;ur为速度;fr为受力;mr为质量。惯性系的刚体质心旋转角速度方程为

(15)

式中:ωr为转动角速度;J为刚体惯性矩阵;M为外力对质心的力矩。

1.5 边界的高保真格式

运用中国运载火箭技术研究院自主可控仿真团队于2020年提出的边界隐式约重构方法,将网格边界处的物理量和物理量梯度同时作为未知数,通过隐式迭代方法求解,在满足边界条件的同时,保障边界的变量重构精度与计算域的内部相同,反映边界各向异性流动特征[7]。

与传统格心型有限体积法直接构造单元中心处的物理量梯度不同,本方法首先构造格点的物理量梯度,然后利用格点梯度通过加权平均组合出单元梯度。格点梯度的构造方法基于最小二乘思想,计算模板是与该格点相邻的所有单元格心,构造目标是模板内所有数据点上的线性重构值误差的平方和最小。从而获得一个以格点物理量值和格点物理量梯度值为未知量的线性方程组,并且对应的系数矩阵是对称矩阵,具有很高的求解效率。

在构造边界格点的物理量梯度时,开创性地将边界面心处的物理量值也作为未知量纳入梯度模板中,该线性系统在边界处是隐式的。通过迭代法对该隐式问题进行求解,并且在求解过程中需要考虑各类边界条件的约形式,最终可以同时获得满足边界约的边界格点梯度值和边界面心值。格心对流项梯度和面心黏性项梯度通过周围格点梯度的线性或非线性凸组合获得[8]。

该方法的优势包括:1)在重构近边界区域梯度时考虑边界约,重构的物理量充分体现了近边界区域的流动特征;2)对流项和黏性项的梯度构造方式统一,具有更好的数值鲁棒性;3)在模拟高超声速流动时获得比梯度限制器更好的非物理振荡抑制效果。

1.6 基于机器学习的数据驱动增强模型

对运载火箭飞行中面临的复杂流动现象,雷诺平均模拟具有工程设计阶段可接受的效率,但基准湍流模式的误差使得预测结果精确性较低;求解大涡模拟方程能够给出较高的预测精度,但计算效率较低无法支持研制阶段的迭代需求。因此对雷诺平均方程的湍流模型进行增强和修正,保证原有计算效率的基础上提升精度成为一种可行的手段。

中国运载火箭技术研究院自主可控仿真软件团队在机器学习湍流模型方面进行了详细研究和探索,在2020年提出了在基准湍流模式预测结果上叠加雷诺应力张量的方式湍流修正算法[9],并将算法应用于流动大分离问题中,以雷诺平均模拟的数值开销获得了与直接数值模拟相近的流动预测精度;在2022年进一步优化机器学习湍流建模框架,将雷诺应力张量表示定理嵌入预测框架之中,提升了机器学习湍流模型预测结果的光滑性、鲁棒性和在不同问题中的泛化能力[10]。

数据驱动模型增强研究需要大量的真实数据以支持模型训练,中国运载火箭技术研究院作为运载火箭的抓总研制单位已经积累了许多型号风洞试验和飞行试验遥测数据,利用这些真值数据开展数据驱动模型研究具有优势,通过已有结果修正得到的增强模型能够有效地应用于新型号的研制之中,提升预测精度、释放设计余量。

2 自主可控气动仿真软件工程应用

本节阐述自主可控仿真软件在运载火箭研制中的气动仿真工程应用,包含气动力仿真、滚转干扰气动优化和级间分离压力仿真。

2.1 典型状态气动力仿真

气动特性设计主要任务是根据给定的外形和飞行剖面,预测火箭六分量气动力、力矩以及动阻尼系数等气动特性,从而为弹道设计、载荷计算和控制系统设计提供依据,设计方法包括工程计算、风洞试验和数值仿真等。新一代固体火箭全面采用数值仿真方法进行全箭气动特性设计,其中基于自主研发的高保真可压缩湍流仿真方法完成了全飞行剖面的气动特性设计。

边界条件方面,箭体是绝热无滑移边界,计算域轮廓是压力远场,远场来流条件取为一级飞行阶段部分典型状态参数,见表1。

表1 捷龙三号火箭典型状态来流参数Table 1 Typical inflow parameters for the SD-3 rocket

求解方法采用雷诺平均模拟,典型状态的压力系数分布如图1所示。

图1 典型状态仿真结果的压力系数分布云图Fig.1 Pressure coefficient contours of typical simulation results

跨声速气流在整流罩的表面形成了λ激波,激波的位置随马赫数提高而向下游移动。通过湍流仿真结果修正了常规气动力设计参数。基于自主研发软件和某商业软件的典型工况气动特性仿真结果的对比如图2所示。经比较,一级飞行段轴向力系数CA偏差不超过3%,法向力系数CN偏差不超过4%,压心系数XCP偏差不超过±0.01。

图2 一级飞行段典型工况的自研软件与商业软件气动设计仿真结果对比Fig.2 Comparison of aerodynamic design simulation results between the self-developed software and commercial software for typical operating conditions in the first stage flight

2.2 滚转干扰气动仿真

预示火箭表面凸起物诱导的滚转气动干扰是姿控系统设计的关键环节[11]。由于滚转力矩数值比一般测量天平的量程小,并且风洞试验中难以同时复现飞行试验的马赫数与雷诺数,因此,相比纵向力矩测量结果,常规风洞测力试验获取的滚转力矩准确度较低。基于高精度数值仿真方法,通过采用真实的气动外形进行流动仿真,获得包括滚转力矩在内的气动特性,开展相应布局的气动外形优化。

通过仿真获得滚转气动干扰,优化了电缆罩布局方案。火箭各级发动机需要布置四排电缆罩。在初始方案中,方位角0和180°各布置两排,在优化方案中,方位角0°、90°、180°和270°各布置一排。边界条件方面,箭体和包含爆炸螺栓盒与电缆罩在内的凸起物均采用绝热无滑移壁面,计算域的外轮廓是压力远场,远场边界条件取如表2所示来流参数。

表2 捷龙三号滚转气动干扰仿真来流参数Table 2 Inflow parameters for simulations evaluating the aerodynamic inference on rolling moment of the SD-3 rocket

求解方面,Ma≥7.00时采用欧拉方程求解,其余状态采用雷诺平均模拟,可压缩修正在Ma≥2.00时运用。在固定马赫数和攻角下,取各方位角最大的滚转力矩作为对比依据,两种方案的最大滚转力矩系数如图3所示。

图3 最大滚转干扰力矩优化前后对比Fig.3 Comparison of maximum rolling moments before and after optimization

工况1和工况2分别代表原始两排电缆罩方案和优化后四排电缆罩方案。由图3可见,采用优化后的四排电缆罩方案后,一级和二级飞行段的滚转力矩系数显著下降,更有利于开展控制系统设计。

2.3 级间分离非定常压力仿真

固体火箭级间热分离设计的关键任务之一在于根据发动机性能和级间段壳体的布局,正确预示级间段内非定常高压环境,以及分离过程上面级发动机喷流作用于级间段内部的作用力,以便进行级间分离动力学仿真设计以及级间段壳体与内部分插支架等结构部件的设计[12]。

通过仿真进行捷龙三号火箭的级间压力环境精细化分析,一二级级间段的几何构型如图4所示。

图4 一二级级间段的几何构型示意图Fig.4 Geometry of the interstage section between the first and second stages

边界条件方面,箭体采用绝热无滑移壁面,喷管喉道状态根据内弹道给定,计算域的外轮廓是出口边界,地面和飞行状态的初始条件见表3。

表3 级间分离的仿真初始条件Table 3 Initial conditions of stage separation simulation

求解采用雷诺平均模拟,在柔性导爆索切割后运用任意拉格朗日-欧拉框架进行求解。在飞行状态下,以发动机堵片打开时刻为时间零点,级间段的压力系数变化规律如图5所示。发动机堵片打开后,喷管形成弓形激波,并在下面级前封头形成高压脱体激波,随后波系在上面级和下面级间反复传播,高压区域逐渐向整个级间段扩散。根据级间段内非定常高温高压环境的仿真结果,获得了下面级前裙和上面级后封头边缘的冲击载荷,评估了级间憋压分离的力学环境。

图5 级间分离过程压力系数演化规律Fig.5 The evolution of pressure coefficient during the stage separation

在地面状态的级间分离试验中,在级间段内多处位置布置了压力传感器,获取了各部位测点非定常压力变化过程。在上述地面状态仿真算例中,对照地面试验提取了相应测点位置的压力数据,其变化规律以及与试验结果的对比如图6所示。第一道冲击波作用于前封头产生了瞬时高压,随着时间变化测点压力逐渐降低,在二级后端框附近,流动存在滞止高压区,仿真结果与地面试验数据符合良好。

图6 地面状态级间段内分离过程非定常压力仿真与试验结果Fig.6 Simulation and test results of the unsteady pressure during the stage separation in ground state

3 对流换热仿真

本节阐述软件在运载火箭研制的对流换热方面的工程应用案例,包含凸起物气动加热仿真和尾舱对流换热仿真。

3.1 凸起物气动加热计算

火箭在飞行过程中,来流在凸起物局部形成干扰流场,从而引起局部高压高热环境,其热环境条件相较于其周围部段存在较大差异,按照传统方式进行大面积包络式防热会造成过度防热,增加产品重量。针对箭体表面凸起物开展气动加热环境精细化仿真设计工作有助于控制产品重量,并降低生产成本。

根据火箭表面热环境评估结果,获得热环境恶劣的典型状态,针对火箭表面的凸起物,开展气动加热精细化仿真分析工作。边界条件方面,箭体和凸起物均采用绝热无滑移壁面,计算域外轮廓的压力远场依据热环境恶劣的典型状态确定,边界条件如表4所示。

表4 捷龙三号凸起物气动加热计算的远场边界条件Table 4 Far field conditions for the SD-3 bump aerodynamic heating simulation

求解方法运用雷诺平均模拟,考虑湍流的可压缩性修正。压力系数的仿真结果如图7所示。仿真结果预示了高速气流在凸起物附近形成的干扰流场,仿真获得的凸起物热流密度相比包络设计结果显著下降。

图7 凸起物气动仿真的压力系数分布云图Fig.7 Pressure coefficient distribution contour of the bump aerodynamic simulation

以二级电缆罩的前封头为例,采用不同设计方法的热环境精细化设计结果对比如图8所示。采用基于仿真设计方法后,峰值热流降低超16%。借助以往风洞测热试验结果及其他商业仿真软件的校核,可以使上述仿真设计方法对凸起物的气动热环境预示更加准确,有利于提升火箭总体性能。

图8 电缆罩前封头气动加热仿真结果对比Fig.8 Comparison of simulation results of pneumatic heating of cable forward dome

3.2 尾舱热环境计算

在火箭飞行阶段,发动机喷流向大气膨胀并产生流向箭体底部的高温燃气回流,对底部产生对流加热,为此需要对底部结构和仪器设备进行热环境和热防护设计[13]。基于仿真可以模拟不同飞行状态下底部燃气喷流与外流耦合作用,进而获得空间各处对流热环境,可用于进行底部结构与仪器设备的热环境条件制定。

针对一级发动机的尾舱开展热环境精细化仿真分析,考虑防热裙的真实构型。边界条件方面,忽略喷管壁面传热,箭体采用绝热无滑移壁面,喷管喉道状态根据内弹道给定,计算域外轮廓是压力远场,典型状态根据弹道的动压和攻角确定,喷管姿态考虑0°摆角和最大摆角,典型状态见表5。

表5 捷龙三号尾舱对流环境仿真Table 5 Convection environment simulation conditions of the aft compartment of SD-3

在一级发动机尾舱,攻角平面的温度分布如图9所示。

由图9可见,发动机尾舱、喷管和外部流动形成封闭流动区域,内部为空腔涡流。涡流将喷管的高温气流卷入空腔,对内部结构造成显著的对流加热。

图9 一级发动机尾部的攻角对称面的温度分布云图与局部流线Fig.9 Temperature contours and local streamlines of the first-stage engine symmetry plane

在尾舱的热流密度测点位置处,仿真数据与遥测结果的对比如图10所示。

图10 尾舱对流热流仿真数据与飞行遥测结果对比Fig.10 Comparison of the aft compartment convection heat flux between the simulation and the telemetry data

由图10可见,仿真数据与实测结果一致。

4 脉动压力仿真

在跨声速段,箭体表面的压力脉动显著影响飞行器动态载荷,造成火箭整体弯曲振动、局部结构摆动等现象,冲击载荷在火箭内部形成噪声环境,量化跨声速段脉动压力是火箭总体设计的关键环节。本节阐述脉动压力数值风洞仿真的定解条件、求解方法与仿真结果。

4.1 大头颈比构型分析

整流罩直径与三级发动机直径比超过1.60,根据工程方法判断是否需要开展脉动压力与气动阻尼的精细化仿真分析。工程方法的判据为

(16)

(17)

式中:x和ξ为轴向坐标;r为ξ处半径;q为动压;积分的下限ξ0和上限ξ1分别是整流罩的球头后缘和箭体尾部的轴向坐标。判据C的意义是压强梯度,判据满足0.2≤C<1.8时控制舱的壁面易形成流动分离,需要开展脉动压力和气动阻尼的数值风洞仿真,判据沿箭体轴线的分布见图11。

图11 工程方法的判据结果Fig.11 Criterion distribution of the engineering method

由图11可见,判据满足0.2≤C<1.8,需要开展脉动压力仿真工作。

4.2 流场特性分析

针对跨声速飞行状态开展脉动压力数值风洞仿真,综合考虑脉动压力特性和箭体振型,给出制定动态载荷与噪声环境的跨声速工况。数值风洞的捷龙三号模型采用1∶25缩比,边界条件方面,模型表面采用绝热无滑移壁面,计算域外轮廓采用压力远场,远场边界条件见表6。

表6 捷龙三号跨声速工况远场边界条件Table 6 Far field parameters of the SD-3 under transonic conditions

求解采用大涡模拟,贴体网格的最大壁面距离不超过20 μm,空间离散与时间积分均为二阶精度。仿真结果的瞬态压力系数和涡判据等值面如图12所示。

图12 脉动压力典型状态仿真结果的压力系数云图和涡等值面Fig.12 Pressure coefficient contours and vortex iso-surface of typical conditions with pressure fluctuation

由图12可见,在Ma=0.80的亚声速算例中,整流罩外表面存在显著压力脉动;在Ma=1.10的算例中,整流罩的压力脉动被向下游延伸的膨胀波抑制,由于全箭一阶振型节点在三级发动机后缘,压力脉动显著影响动态载荷的状态在亚声速段。

4.3 脉动压力分析

针对动态载荷较为严酷的速度范围,开展细化的脉动压力系数分布和功率谱密度分析。细化的脉动压力数值风洞仿真的状态,远场边界条件见表7。

表7 细化捷龙三号脉动压力数值风洞仿真的远场边界条件Table 7 Far field conditions of the refined SD-3 pressure fluctuation simulation

在箭体的迎风和背风母线上,均方根脉动压力系数(cprms)分布如图13所示。

由图13可见,Ma=0.80,α=0°算例的脉动压力系数最大,约17%。根据遥测数据,整流罩正锥的脉动压力测点位置处的最大脉动压力系数约18%,与仿真预示结果一致。根据飞行工况的速度和动压与地面状态的比拟关系,对脉动压力数值风洞仿真的脉动压力功率谱密度数据做转化,并将结果与遥测数据对比,对比结果如图14所示。在高频段,受限于传感器的采样频率,无法进行数据的有效对比,在100 Hz以下仿真与遥测的脉动压力功率谱密度较为一致。

图13 箭体均方根脉动压力系数分布Fig.13 Root-mean-square pressure coefficient distribution on the rocket surface

5 结 论

本文针对新一代固体运载火箭捷龙三号的研制特点,结合中国运载火箭技术研究院在仿真设计经验和自主可控工业仿真软件研发方面的技术积累,阐述了捷龙三号火箭全面取消全箭风洞试验后开展的气动仿真分析与设计工作。

中国运载火箭技术研究院工业仿真软件团队运用并发展了多项行业先进的仿真技术,包括针对高速流动进行可压缩性修正的湍流模型,保持边界重构精度与各向异性流动特征的边界格式以及能够基于机器学习提高湍流模拟精度与效率的数据驱动增强型湍流模型,形成了具有较高模型保真度和数值精度的自主可控仿真工具,并在型号研制中的多个环节开展了仿真试验与分析,发挥了重要的支撑作用。

在气动力仿真方面:1)参考典型状态开展了全箭湍流模拟,获得的气动力和气动力矩系数结果被用于验证并优化常规气动力设计参数;2)面向凸起物进行精细化仿真设计,分析由电缆罩布局方案造成的滚转干扰,通过优化电缆罩布局方案显著减小了一级和二级飞行段的滚转力矩;3)面向一二级热分离开展仿真工作,获得级间憋压环境下二级发动机尾喷管弓形激波与相关波系的演化规律,仿真结果被用于级间段静力分析与冲击载荷精细化设计。

在对流换热精细化仿真方面:1)选取典型工况开展凸起物的热环境精细化仿真,基于热流密度比拟关系,获得凸起物随弹道的气动加热规律;2)参考弹道,选取大动压工况开展发动机尾舱对流换热的精细化仿真设计,根据空腔涡流造成的对流加热制定发动机尾舱热防护条件。

在脉动压力仿真方面:1)采用工程算法,给出火箭轮廓属于不稳定分离构型的结论,指导脉动压力与气动阻尼仿真工作;2)依据典型工况,开展脉动压力数值风洞仿真,结合全箭一阶振型,辨识动态载荷的设计工况,并获得控制舱仪器噪声环境的设计工况;3)基于整流罩和控制舱的脉动压力仿真结果获得相应的无量纲脉动压力功率谱密度,制定动态载荷设计的输入条件。

捷龙三号火箭的工程研制首次确立了“将仿真试验当成实物试验来对待”的创新研制路线。对于传统设计存在差异化的项目,加大仿真比重;对于成本巨大或无法进行地面试验的项目,应用仿真技术管控风险。仿真试验大幅替代了实物试验,同时弥补了传统地面试验中无法有效模拟真实尺寸和飞行环境等方面的局限性,充分体现了自主可控仿真软件在未来运载火箭研制中的价值与优势,为今后商业火箭研发模式沉淀了成功经验。