与燃气射流耦合的易裂后盖开启过程数值分析

牛钰森，姜毅，史少岩，李静

（北京理工大学宇航学院，北京100081）

与燃气射流耦合的易裂后盖开启过程数值分析

牛钰森，姜毅，史少岩，李静

（北京理工大学宇航学院，北京100081）

对于采用易裂盖的箱式发射系统，后盖开启过程与燃气射流流场相互影响。为了研究易裂后盖的开启过程，使用基于光滑方法和重构方法的动网格技术，将后盖裂片运动过程与燃气射流流场耦合，对整个开盖过程进行了仿真计算。得到了起始冲击波的作用过程以及发射箱内各监测点的压强曲线，并分析了燃气射流对后盖裂片的冲击作用以及对其运动过程产生的影响。结果表明，与非耦合工况的计算结果进行对比，耦合工况的结果更接近于实验数据。随着后盖裂片转动角度的不断增大，旋转力矩逐渐减小，裂片的转动速度经历了先快速上升、后缓慢下降的过程。

兵器科学与技术；燃气射流；发射箱；易裂后盖；耦合；起始冲击波；动网格

0　引言

垂直发射中采用热发射方式的导弹，发动机点火后燃烧室内的压强迅速上升，达到破膜压力时堵片碎裂，根据文献［1］的研究此时会产生起始冲击波。利用冲击波的超压峰值可以使易裂后盖破裂，之后在燃气的作用下打开，为燃气的后续排导提供了通道。并且在此过程中起始冲击波在发射箱后盖上反射，沿着导弹与发射箱的间隙向前传播，最终将发射箱的易碎前盖击碎。

反射冲击波的强度直接关系到发射箱易碎前盖的打开过程，若强度过小则无法使前盖碎裂造成发射事故，若强度过大则对弹头处的精密元件造成不良影响。不同于文献［2-4］中所述的易碎盖打开过程，易裂盖破裂之后不会立即向外飞散，而是逐渐向外转动，对反射冲击波以及随后的燃气射流流场产生影响。同时燃气射流流场又会影响到易裂后盖的运动状态，二者相互耦合。本文使用动网格技术将发射箱后盖的运动与燃气射流流场耦合，并采用基于扩散方程的光滑方法和局部单元网格重构方法对网格节点位置进行更新。按照实物模型将易裂后盖处理成4块独立运动的刚体壁面，由计算流体力学（CFD）软件实时计算出燃气射流流场对后盖产生的作用力和力矩。使用3步4阶Adams-Moulton多步积分法［5］对后盖的运动微分方程进行求解。计算过程中在发射箱内的特定位置设立监测点，记录冲击波的传播速度及超压峰值。计算结果与实验测得的数据以及非耦合工况的结果进行了对比。

1　理论模型

1.1流场基本方程

在发射箱流场内存在多处马赫数大于0.3的区域，因此必须考虑气体的可压缩性。使用有限体积法对控制单元CV中的任意标量φ建立积分形式的三维可压缩非定常N-S方程［6］：

式中：ρ为密度；u为速度向量；Sφ为标量φ的原项；Γ为扩散系数。考虑燃气与空气的混合过程，使用组分模型将燃气与空气处理为两种组分，混合后气体的属性由理想气体混合定律定义［7］：

式中：p为混合气体压强；R为通用气体常数；T为混合气体温度；Yi为第i种组分的质量分数；Mm，i为第i种组分的摩尔质量；cp，i第i种组分的比定压热容。

1.2动网格方程

加入动网格后，网格节点也有运动速度，因此需要对控制方程进行修正。对控制体中任意一种标量φ，修正后其积分形式的守恒方程为

式中：ug为网格节点的速度向量；Γ为扩散系数；∂V表示控制体CV的边界；A为边界面的面积法向量。

式易裂后盖的开启过程是绕固定轴的旋转运动，因此选用基于扩散方程的光滑方法更为合适。使用此方法时由运动边界上网格节点的速度可得内部节点的速度

式中：扩散系数γ由网格节点距离边界的归一化距离d得到：

式中：α为扩散参数。

由网格节点的位移速度更新节点的坐标

式中：Xn为前一时间步网格节点的坐标向量；Xn+1为当前时间步的坐标向量。

易裂后盖从封闭到完全打开的过程中，转动角度达到了90°，运动幅度很大。若单纯采用光滑方法则会出现网格尺寸过大或过小以及畸变率过高的单元。因此需要使用重构方法对网格进行更新，以保证计算过程中网格的质量满足要求。对于体网格使用局部单元重构方法进行更新。同时，因为后盖运动影响了对称面上的网格单元，所以还需要使用区域面重构方法对变形边界面上的网格单元进行更新。网格单元的重构尺寸由尺寸函数和面单元最大畸变率进行控制。

2　计算域与后盖模型

2.1计算域模型

为便于划分网格对发射箱内以及弹体的结构进行适当简化，包含导弹的发射箱模型如图1所示。

使用混合网格对整个计算域进行划分，导弹主体部分与发射箱之间的间隙空间使用六面体结构化网格进行填充，发动机喷管出口到导流格栅的部分以及后盖运动区域使用四面体非结构化网格进行填充，整体网格如图2所示。

2.2易裂后盖模型

使用两种方式对发射箱易裂后盖进行处理。一种方式不考虑后盖与燃气射流耦合的打开过程，按照打开面积的大小将其简化为6个打开状态，如图3所示。同心方框代表了后盖打开面积占总面积的百分比，从内到外依次是10%、20%、30%、50%、70%和100%.计算过程中在指定的时间依次将相应的区域由壁面边界条件改为交界面边界条件，以这种方式处理的工况称为工况1.

图1　流场计算域结构图Fig.1 Structure diagram of computational domain

图2　流场计算域网格Fig.2 The mesh of flow field computational domain

图3　工况1中的后盖模型Fig.3 Back cover model of Case 1

另一种方式将发射箱易裂后盖处理成4瓣可以转动的刚体，如图4所示，计算过程中后盖的运动与燃气射流耦合。此方式处理的模型与真实模型一致，如图5所示，将此工况称为工况2.

图4　工况2中的后盖模型Fig.4 Back cover model of Case 2

图5　发射箱易裂后盖Fig.5 Back cover of launch canister

后盖裂片在燃气作用下不断转动，取朝向发射箱外转动为正方向，其运动微分方程为

式中：θ为后盖转动的角度；ω为后盖转动的角速度；MF为燃气流场在后盖上产生的旋转力矩，是关于时间t的函数；MG为裂片重力产生的力矩，是关于角度θ的函数；MK为后盖转动后产生的回复力矩，是关于角度θ的函数；J为裂片转动惯量。

设MF、MG与MK的合力矩为M（t，θ），并使用3步4阶Adams-Moulton多步积分法对微分方程（8）式进行求解，得式中：h为时间步长。

通过实验取得了后盖裂片的回复力矩MK，测得了多个转动角度下回复力矩的大小。并使用最小二乘法对离散的实验数据进行了拟合，得到了MK关于转动角度θ的函数曲线如图6、图7所示。

图6　裂片1、2的回复力矩Fig.6 Restoring moment curves of Debris 1 and 2

图7　裂片3、4的回复力矩Fig.7 Restoring moment curves of Debris 3 and 4

3　计算结果

3.1发射箱内流场

在距离发射箱底部0.575 m、2.665 m、4.565 m三处放置压力传感器，记录发射箱内流场的压强变化情况，如图8所示。

图8　发射箱内传感器位置Fig.8 Locations of sensors in launch canister

当燃烧室内的总压超过1.5 MPa时堵片破裂，在发动机喷管内形成一道起始冲击波。且冲击波的传播速度快于燃气射流的形成速度，因此同一时刻冲击波波阵面位于燃气与空气的交界面之前，如图9、图10所示。

图9　0.4 ms时压强分布云图Fig.9 Pressure contour at 0.4 ms

图10　0.4 ms时燃气分布云图Fig.10 Mass contour at 0.4 ms

冲击波向前传播，撞击到发射箱易裂后盖上，形成一片高压区域，使后盖上的压强急剧上升，如图11、图12所示。

图11　0.9 ms时压强分布云图Fig.11 Pressure contour at 0.9 ms

图12　后盖上的平均压强曲线Fig.12 Average pressure curve of back cover

后盖从受到冲击波的压力作用到破裂需要一定的时间。以静态实验的结果来看，当后盖的内外压差达到0.072 MPa时预制裂纹开始破裂，1 ms左右后裂纹完全破裂。在此时间内可将后盖看作为固壁面，冲击波在后盖上发生反射向发射箱前盖传播。根据文献［8］的研究，后盖的打开时间对反射冲击波的超压峰值和持续时间有很大的影响。工况1非耦合情况与工况2耦合情况对后盖的处理方式不同，因此两个工况中的后盖打开状况也不相同，导致反射冲击波的强度不同，两个工况下监测点压强曲线与实验数据的对比，如图13所示。

在冲击波的作用下，P1、P2、P3监测点的压强都经历了先快速上升、后缓慢下降的过程。其中P1监测点的压强曲线在0～5 ms的区间内有多个峰值，这是冲击波在导流面、格栅壁面和后盖之间发生多次反射产生的。P3监测点的压强曲线在14 ms左右出现了阶跃现象，这是因为冲击波撞击到前盖，再次发生反射形成反向传播的冲击波造成的。可以看出数值计算结果与实验数据基本吻合，曲线的变化趋势大体相同。但是由于对计算模型进行了简化，去掉了弹翼、脱插、沟槽等细小结构，由此造成的微弱扰动在数值计算中无法再现。按照（10）式计算数值结果的误差：

式中：x为数值计算结果，在此处分别是冲击波超压峰值和冲击波的传播速度；x*为实验结果。冲击波的传播速度按照（11）式进行计算：

式中：ti为第j个监测点的压强曲线波峰出现的时间；sj为第j个监测点与第j+1个监测点之间的距离。

图13　监测点压强曲线Fig.13 Pressure curves of monitoring points

在冲击波超压峰值上工况1的计算结果与实验数据的平均相对误差为8.29%，工况2为6.69%.在冲击波的平均传播速度上，工况1与实验数据的相对误差为5.15%，工况2为3.88%.由此可见，工况2耦合情况的计算结果误差更小。

3.2易裂后盖运动过程

易裂后盖破裂后，图11中所示的发射箱底部高压区域通过裂缝向排烟道内传播一道压缩波，如图14所示。

同时后盖裂片在发射箱内高压以及燃气射流的冲击作用下开始向外转动。旋转力矩MF随时间变化情况如图15所示。

图14　1.8 ms时压强分布云图Fig.14 Pressure contour at 1.8 ms

图15　MF变化曲线Fig.15 Curves of MF

由图15可以看出，在后盖破裂后2 ms的范围内，裂片受到的旋转力矩出现了几次波动，这是图14中所示的压缩波造成的。后盖刚刚破裂时压缩波的传播距离不远，主要作用区域集中在后盖裂缝周围，致使裂片背面的压强升高，降低了内外压差，从而使旋转力矩减小。压缩波远离裂片后，背面压强迅速下降，致使内外压差升高，旋转力矩增大。裂片不断向外转动，后盖打开面积逐渐增大使发射箱底部的压强迅速下降，如图16所示。并且随着张开角度的增大，燃气射流在裂片上的冲击作用面积不断减小，如图17所示。在这两个因素的作用下致使裂片受到的旋转力矩不断减小。

裂片转动的角度θ和转动的角速度ω按照（9）式所示的方法进行计算，结果如图18、图19所示。

裂片的转动角速度在后盖破裂后的8 ms时间内经历了一个快速上升的阶段，之后旋转力矩MF趋近于0，在裂片回复力矩MK的作用下转动角速度逐渐减小。与此对应裂片的转动角度经历了快速增大后趋于平稳的过程。

图16　不同时刻后盖附近的压强云图Fig.16 Pressure contours around debris at different times

图17　不同时刻燃气射流在裂片上的冲击区域Fig.17 Effected areas at different times

图18　裂片转动角度随时间变化曲线Fig.18 Rotating angle vs.time

图19　裂片转动角速度随时间变化曲线Fig.19 Rotating angular speed vs.time

4　结论

本文采用了基于光滑和重构方法的动网格技术，耦合求解了易裂后盖的运动微分方程与燃气流场方程，并计算了非耦合工况与其进行对比。计算结果表明：

1）耦合工况的计算结果与实验数据更为接近，比非耦合工况的计算精度高。

2）发动机破膜后产生起始冲击波，撞击到发射箱后盖致使其破裂，并发生反射向发射箱前盖传播。

3）后盖破裂之后发射箱底部高压迅速下降，产生压缩波向外传播，随着打开面积的逐渐增大，后盖内外压差不断下降。

4）随着后盖裂片翻转角度的增大，燃气射流在裂片上的冲击作用面积不断减小。

5）后盖裂片的转动角速度经历了先快速增大、后逐渐减小的变化过程。

综上所述，该分析方法可为发射箱易裂后盖的设计验证和改进提供参考。

（

）

［1］徐强，李军.燃气射流起始冲击波形成机理的实验研究［J］.推进技术，2000，21（3）：16-18.XU Qiang，LI Jun.Experimental study on mechanism of initial shock wave in jet flow［J］.Journal of Propulsion Technology，2000，21（3）：16-18.（in Chinese）

［2］钱元，周光明，贺卫东，等.轻质复合材料发射箱盖结构设计与参数优化［J］.航空学报，2013，34（4）：826-832. QIAN Yuan，ZHOU Guang-ming，HE Wei-dong，et al.Structure design and parameter optimization of lightweight composite canister cover［J］.Acta Aeronautica et Astronautica Sinica，2013，34（4）：826-832.（in Chinese）

［3］孙志彬，王新峰，周光明.圆形复合材料定向分离［J］.南京航空航天大学学报，2012，44（6）：803-808. SUN Zhi-bin，WANG Xin-feng，ZHOU Guang-ming.Specified direction separation of circle composite fragile cover［J］.Journal of Nanjing University of Aeronautics&Astronautics，2012，44（6）：803-808.（in Chinese）

［4］王松超.某火箭系统易碎式密封盖研究［D］.南京：南京理工大学，2013. WANG Song-chao.A rocket system fragile type sealing cover research［D］.Nanjing：Nanjing University of Science and Technology，2013.（in Chinese）

［5］Burden R L，Faires J D.Numerical analysis［M］.US：Cengage Learning，2010.

［6］Versteeg H K，Malalasekera W.An introducion to computational fluid dynamics：the finite volume method［M］.England：Pearson Education，2007.

［7］ANSYS Inc.ANSYS FLUENT user's guide［M］.US：ANSYS Inc，2010.

［8］党海燕，贺卫东，鲁志刚，等.激波开盖中后盖开盖压力对前盖压力的影响研究［J］.战术导弹技术，2012（1）：54-57. DANG Hai-yan，HE Wei-dong，LU Zhi-gang，et al.Research on the effect for the opening-pressure of post-cover to the pressure of pre-cover in the shock wave opening［J］.Tactical Missile Technology，2012（1）：54-57.（in Chinese）

Numerical Analysis of Fragile Back Cover Opening Process Coupling with Jet Flow

NIU Yu-sen，JIANG Yi，SHI Shao-yan，LI Jing
（School of Aerospace Engineering，Beijing Institute of Technology，Beijing 100081，China）

The motion process of fragile back cover interacts with the flow field of gas jet.The dynamic mesh technology based on smoothing method and remeshing method is used to study the opening process of back cover，The opening process is simulated by coupling the motion process of back cover with the gas jet flow field.The initial shock wave is observed，and the pressure curves of monitoring points are plotted.The effect of jet flow on the motion of cover debris is analyzed.Compared to the results of uncoupled case，the results of coupled case are closer to the experimental data.The driving moment is getting small with the increase in rotation angle of cover debris.The rotation speed rises fast at first and then goes down slowly.

ordnance science and technology；gas jet；launch canister；fragile back cover；coupling；initial shock wave；dynamic mesh

TJ768.2

A

1000-1093（2015）01-0087-07

10.3969/j.issn.1000-1093.2015.01.013

2014-03-05

国防基础科研计划项目（B0420132805）

牛钰森（1987—），男，博士研究生。E-mail：shenzhou1987@aliyun.com；姜毅（1961—），男，教授，博士生导师。E-mail：jy2818@163.com