1 引言
多级子母弹分离的空气动力学问题，是一类典型的多体干扰、非定常复杂流动问题[1]。以机械式分离抛射为例，全过程大致可分为3个阶段进行。第一阶段，子母弹飞行过程中通常采用整流罩来保护部件免受气动力、热环境的影响，达到减阻的效果，到达一定高度和速度时需将整流罩抛离。子母弹有着飞行速度高、飞行高度低的特点，复杂的气动力将使头罩分离后的运动轨迹复杂，并可能与弹体相撞，导致灾难性的后果。第二阶段，导弹外层挡板在弹簧机构作用下以平冲的状态离开母弹，之后受到干扰区流场的影响，产生复杂的气动力，挡板的运动状态受到气动力的影响，需要设计合适的分离条件以保证挡板不与其他运动部件发生碰撞。第三阶段，抛撒过程是子母弹武器系统的关键环节，子弹抛撒分离存在严重的子弹与母弹、整流罩、挡板间的气动干扰，尤其是在超声速飞行条件下的子母弹分离，往往伴随着各运动部件间激波干扰，使得分离过程的流场结构异常复杂，可引起子弹及母弹气动力不确定性，进而影响到子母弹的分离及毁伤效果。子母弹抛撒的各阶段涉及到不同形式的多体分离问题同时又相互影响干扰，有必要对子母弹的多级分离全过程展开数值模拟。
国际上，美国率先开始了对于子母弹分离这一复杂流场的分析[2-4]。包括采用数值模拟与风洞试验相结合的方法对子母弹气动特性进行了研究[5];采用计算机程序精确模拟超声速空腔中多个子弹弹体分离流场[6];采用非结构动网格方法、总变差不增有限体积法和并行计算技术对子母弹分离过程进行数值模拟[7]等,均取得与试验结果符合较好的数值模拟结果。近年来，国内对子母弹分离流场特性的研究取得长足进展[8]。例如，对弹舱靠前和弹舱靠后2种子母弹结构进行了数值模拟[9]；讨论出舱条件、气动干扰和非定常效应等对分离运动的影响[10]；对时序抛撒方式下子弹在不同舱段分离的三维非定常流场进行了数值模拟[11]；对扰流板对子母弹翻转角速度的影响进行研究[12]；对大口径中心爆管式子母弹分离流场进行数值模拟[13]；对子母弹抛撒过程中气囊破裂射流流场结构进行仿真分析[14-15]等。前人大多是对多级子母弹分离单过程展开详细的数值模拟研究，默认开始子弹抛撒时之前抛离的整流罩与外层挡板已经远离干扰区。但是抛开头罩后导弹又飞行了一段时间，此时气动外形提前发生改变，气动系数改变，对子母弹原本规划的运行轨迹会有一定影响。如何在短时间间隔内高效安全实现三级分离，便是本研究方案探讨的目的。
三级多体分离相对单级多体分离而言，由于级间干扰的存在，除了很多共性问题以外，也有特性问题。子级部件大小、气动结构外形、分离时的姿态、分离初始条件、部件质心位置以及各级分离时序等因素，均会影响到整个分离过程的成功与否。综合过程不仅要保证各级分离单独可行，还要保证三级分离全过程的安全高效进行。只有根据各因素的安全边界进行适当调整设计，三级分离的方案在实际应用中方可实现。本研究针对虚拟模型的三级分离方案展开数值模拟，重点探讨来流迎角、各级分离初始速度、盖板气动干扰对综合过程的影响，并对设计方案的可行性与稳定性展开验证。
2 数值方法
本次数值模拟所采用的课题组自研软件，先前已对子母弹实验以及头罩分离等多体分离标模进行过验证[16]，对单级分离过程如四体全自由飞[17]、低空高动压整流罩分离[18]、抛撒投放分离、两体级间分离的实际数值模拟说明，建立的数值模拟计算系统可以在一定程度上解决外形复杂的多体运动的实际工程问题。图1为应用该软件对全自由飞数值模拟结果与实验数据的对比验证。图2为头罩分离过程的数值仿真研究。
图1 全自由飞前部俯仰角随时间变化
Fig.1 The change of the forward pitch angle in full free flight
图2 子弹穿越母弹激波分离数值计算与实验对比
Fig.2 Comparison between numerical calculation and the experiment for shock wave separation of a submissile passing through the mother shell
2.1 非结构直角网格技术
本文中数值计算网格，采用三维非结构粘性直角网格[19-21]重构的方法生成。其中主要包括两大组成部分：三维非结构粘性直角网格的生成和运动过程的网格重构。由于整流罩抛罩、外层挡板抛离、子母弹抛撒3个过程按一定的时序设计先后进行，干扰区内同时存在多个运动部件，流场物面外形复杂且时空间尺度差异较大，采用非结构网格可以更准确地对物体复杂的几何外形进行描述。直角网格具有生成简单、流动求解器容易实现等优点，同时采用局部网格加密技术提高边界附近流场的计算精度。运动物体小位移采用弹性网格变形技术，大位移直接进行网格重构。由此生成的三维非结构贴体直角网格，可以适应各种复杂外形的黏性流与湍流的计算。
2.2 ALE形式控制方程、6DOF运动方程耦合
多体分离问题是典型的动态问题，在发展针对此类问题的数值方法时，需要考虑如何捕捉运动边界。使用任意拉格朗日-欧拉形式控制方程对流场进行描述，由牛顿基本定律结合惯性/非惯性坐标系转换推得六自由度弹道方程。
耦合方法：基于当前步网格、流场，冻结各运动物体的位置姿态，求解流动控制方程，得到时间方向推进后的下一时间步的流场，积分得到运动物体受到的气动力和力矩；根据n时刻物体的位置姿态，冻结运动物体的受力，求解刚体运动方程，获得时间步内物体位置和姿态的改变量；根据物体新的位置姿态，采用动网格技术获得n+1时刻的计算网格；完成n时间步的计算，开始n+1时刻的计算。
2.3 时空间离散方法
有限体积法对网格类型、网格质量要求不高，适合于复杂外形、复杂运动过程绕流的计算，在多体分离问题模拟中应用较为广泛，本研究主要采用基于格心的ALE控制方程有限体积离散方法。
空间离散方面，采用有限体积法对流动控制方程进行空间半离散。
(1)
对于任意网格单元，M表示单元面的数目；W表示单元第i个面的面积向量；Fs和Fvs中的下标s表示变量取单元面上的值。这里粘性通量可采用中心差分求解，对流通量可以采用各种通量计算格式，但由于ALE流动控制方程中考虑了网格速度Vb的影响，因此通量格式中也需要相应的计及网格速度，如考虑了单元面运动的Roe格式。
计算对流通量时，迎风格式采用MUSCL线性重构技术以获得二阶精度，并采用Venkatakrishnan、Barth等限制器以抑制较大流动梯度区域出现振荡，由此可实现超声速弹体流场激波的精确捕捉。研究重点关注物体间分离运动且分离相对速度较小，湍流对此过程影响较小，考虑到计算效率问题，并未引入湍流模型。
时间离散方面，双时间步方法通过引入虚拟时间步，将因为稳定性条件对时间的限制转移到虚拟时间步内，物理时间上允许的时间步长将有条件地放宽。对流动控制方程进行空间离散，得到半离散型方程：
(2)
式中，Ri表示单元i各单元面无粘通量与粘性通量的积分。此半离散式中包含瞬态项，能完整描述流动随时间的变化。当采用双时间步方法模拟非定常流动时，需要在方程中引入虚拟时间项：
(3)
这里用τ表示虚拟时间。对虚拟时间步采用一阶欧拉隐式离散，物理时间步采用k阶隐式离散，离散后的通用表达形式为：
(4)
当采用动态网格计算时，单元体积随着时间是变化的，即
求解过程中先在虚拟时间上推进，当虚拟时间上达到收敛时，即认为在真实时间上从第n步推进到n+1步。
3 子母弹分离过程的数值模拟
3.1 子母弹模型
子母弹的几何模型如图3所示。该子母弹由弹体、整流罩、子弹以及外层挡板四大部分组成。本研究主要针对子母弹抛撒过程中各分离部件展开数值模拟，子母弹的其余组成部件暂且不在模型中细致的表现出来。
图3 子母弹网格模型及内部结构展示
Fig.3 The mesh of a submissile model and its internal structure
模型中各部件质量属性设计如下：设计子母弹总质量400 kg，长度为2.14 m，直径0.4 m，弹仓长度0.8 m。整流罩与外层挡板分别重4.7 kg和12.7 kg，质量均匀分布。子弹长度0.66 m，直径0.07 m，质量5 kg，4枚子弹等间距分布于弹舱内。多级分离时，在保证各分离部件不碰撞的情况下，整流罩及外层挡板的稳定性可以不予考虑，但子弹及子母弹整体的分离稳定性问题需要单独进行设计，也即分离过程中质心位置的选取。经过后续的数值模拟发现，子弹质心选取在距头部三分之一弹体长度处可保证多级分离全过程子弹与母弹的飞行稳定性。
由于各分离部件相对母弹而言尺度较小且成轴对称排布，在初期分离的较短时间内母弹状态变化不大，因而可以锁定母弹的自由度进行研究，由此产生的微小误差在可承受范围内。
3.2 子母弹运行工况、分离时序设计
假设子母弹飞行至作战区域后，在海拔高度4 km处以1.5马赫速度进行三级分离运动，进行来流条件设置。三级分离时不仅存在着各部件间的气动干扰，还存在着分离时序干扰问题。设计了以下一种分离方案：整流罩张开至45°时开始进行抛罩；抛罩分离进行至20 ms时，开始进行外层挡板的抛离；外层挡板开始抛离30 ms后，子弹在弹簧机构作用下以一定的初始速度抛出。
3.3 计算域设置及网格生成
多阶段分离时，由于各分离部件分离时间点不一致，会导致分离干扰涉及的空域较大。例如当第三阶段子弹开始进行抛撒时，第一阶段分离的整流罩已经抛出很远的距离。为了能够在分离初期对各个分离部件进行持续的数值模拟，计算域也需设置较大。
规定好计算域，对计算域进行初步的网格尺度划分。在表面网格的基础上按照直角网格生成与重构方法，采用三维非结构粘性直角网格技术对导弹进行计算域内的网格生成。仿真过程中，弹舱内子弹、头部整流罩以及弹舱处的外层挡板为运动部件，需要重点分析运动部件周围的流场信息，因此对这些部件周围流场网格进行加密设置。
4 多级分离方案可行性分析
针对上述模型，调整来流迎角与各级分离的初始速度，对设计方案的可行性与稳定性展开研究，探讨单因素对综合过程的影响。对单独进行子弹抛撒及三阶段连续进行2种分离状况数值模拟，研究头罩挡板对子弹的气动干扰，说明展开多级研究的必要性。
当自由来流攻角为0°时，对前述时序设计的适用性展开研究。对每一个分离阶段单独进行数值模拟发现：第一级分离的头罩在多数分离姿态下均可保证分离安全性；二级挡板以小于16 m/s初始速度平抛或三级子弹以小于12 m/s速度抛出时，在挡板及母弹的流场干扰作用下，挡板及子弹抛离后与母弹尾翼有发生碰撞的危险；当子弹初始分离速度大于挡板时，两侧子弹与先前抛离的外层挡板有发生碰撞的危险。因此，对于该时序方案而言，需要严格设计各级分离的初始分离状态，才能保证方案总体可行，其中第二级与第三级分离应是研究重点。
图4是标准设计状态下，来流迎角0°、挡板平抛初始分离速度V1=20 m/s、子弹抛撒初始分离速度V2=20 m/s时，对该多级分离方案的数值模拟。仿真结果进一步说明该设计方案的可行性，同时也符合短时间内实现三级分离的设计要求。在短时间内完成三级分离，多体干扰区内同时存在着整流罩、外层挡板、子弹以及母弹等多个运动部件，彼此间的相互干扰严重，激波在各运动部件间反射交错，流场复杂。观察流场信息图可知，整流罩在开始分离50 ms后基本飞离干扰区；外层挡板100 ms时已基本不再对后续的子弹抛撒产生作用；子弹在全过程分离开始90 ms后，即开始子弹抛撒40 ms后开始穿越母弹弹体激波层。后文将围绕该特定方案，调整其他参数，通过独立研究各分离因素对综合过程的影响，进一步确立可行方案。
图4 多级分离数值模拟结果
Fig.4 Numerical simulation results of multistage separation
4.1 分离姿态对方案可行性的影响
设计挡板初始分离速度V1=20 m/s,子弹初始分离速度V2=16 m/s，改变自由来流的迎角，研究设计方案在大迎角状态的分离可行性。分别对来流迎角为0°、5°、10°三种情况数值模拟，以下方子弹为例分析。图5所示为子弹分离 100 ms，俯仰角随分离时间变化。情况1、2、3分别代表来流迎角为0°、5°、10°的分离姿态。按照方案设计，质心位置大部分时期处于压心前方，保证了子弹抛撒后的动态稳定性。子弹穿越母弹激波过程中，压缩气流首先作用于弹体前端，导致短暂的失稳状态，在这20～30 ms期间压心逐渐移至质心之后，子弹快速恢复飞行稳定性。由三维仿真图像可知，虽然该过程子弹俯仰角产生较大变化，但仍可保证三级分离的安全进行。由图5可知，设计方案应尽量保证在来流迎角较小的情况下进行子弹抛撒，来流迎角10°时俯仰角上下波动约30°，说明大迎角可能会导致短时间的失稳状态扩大化，不利于子弹保持稳定状态下的安全分离。
图5 3种来流迎角下子弹俯仰角变化
Fig.5 Variation of the bullet pitch angle under three angles of attack
4.2 挡板及子弹分离速度对方案可行性的影响
第二级分离时挡板的平抛速度V1以及第三级分离时子弹的抛撒速度V2，是影响综合过程是否安全稳定的重要因素。图6所示为5种分离条件下，子弹受母弹与挡板激波层共同作用时，穿越激波过程中力矩的变化。图中M1代表上方子弹俯仰力矩，M2代表下方子弹俯仰力矩，M3、M4代表左右两侧子弹的偏航力矩。4种力矩均是对子弹分离后运动稳定性的考量，若能够保持稳定，则证明在设计的分离速度下该方案可行。
子弹穿越激波过程，受到母弹与外层挡板激波层的共同作用，根据图6对该过程子弹受力状况展开分析。首先母弹激波层作用于子弹前半部分，附加弹头向外偏转力矩，此时压心位置位于质心前方，子弹处于短暂的失稳状态；随着抛离过程的进行，压心不断后移，20～30 ms后激波层作用于子弹后半部分，子弹快速恢复运动稳定，同时激波附加弹头向内偏转力矩，原本外张的姿态角开始向内靠拢，直至完全穿越母弹激波层。通过子弹质心位置的设计，可使子弹抛离后总体上维持运动稳定。三维仿真结果显示，对于图6提出的五种分离方案，子弹抛撒过程虽有短暂的失稳状态出现，并伴随较大的俯仰力矩变化，但不会影响到分离方案的可行性与整体稳定性。
对比图6中(a)(b)(c) 3种情况，固定子弹抛撒初速度的情况下，适当增大挡板平抛速度V1，可以有效防止第三级抛出子弹与第二级抛出挡板间发生碰撞，但子弹抛离后穿越激波过程中所受的力矩变化较大，稳定性稍差。对比图6中(a)(d)(e)三幅图片，随着子弹初始分离速度V2增大，穿越激波过程的左右两侧子弹偏航力矩变化也在随之增大，即侧方位子弹分离后会产生更大的姿态角偏转。综上，方案设计应在安全分离的前提下保证V1与V2尽可能小，有利于分离后子弹的稳定飞行。
图6 子弹力矩变化
Fig.6 Change of bullet moment
4.3 2种情况子弹抛撒对比分析
对单独进行子弹抛撒与三级分离连续进行2种情况展开数值模拟。多阶段分离由于发生在较短时间内，前一个阶段的分离运动导致的流场变化对后续阶段的分离状态都会产生较大影响。由于挡板侧向分离，因此这种影响对于侧面两枚子弹最为明显。以XY平面内的一枚子弹为研究对象，对2种分离状况中子弹抛撒这一重点环节对比，探讨整流罩、盖板对子弹抛撒时运动姿态的具体影响和作用机理。
图7、图8所示为2种情况子弹抛撒时，侧方位子弹质心位置以及偏航角随时间变化情况。情况1代表单独进行子弹抛撒，情况2代表三级分离连续进行。由图7可知，该模型在子弹抛撒的初段时间内，情况1与情况2子弹分离速度相当，均可保证方案安全可靠。图8显示，2种分离情况的差异主要体现在姿态角的变化上，情况1子弹的偏航角略大于情况2，这也与先前数值模拟结果相吻合。结合图9力矩变化以及前文流场图像可知，之所以出现2种分离情况子弹偏航角相差较大的现象，是因为情况2中上一步分离的外层挡板仍在干扰区内，子弹在分离初始仍受到挡板激波的作用，作用于子弹的前半部分，附加弹头内向偏转的偏航力矩，抑制偏航角的持续增大。由此也可看出，外层挡板的存在减小了子弹分离时的力矩变化，有利于其保持相对稳定的分离状态。
图7 2种情况子弹抛撒质心位置随时间变化
Fig.7 The change of the position of the mass center vs.time under two situations
图8 2种情况姿态角随时间变化
Fig.8 The change of the attitude angle vs.time
图9 子弹力矩随时间变化
Fig.9 The change of the moment vs.time
5 结论
综合了头罩抛离、外壳分离、子母弹抛撒三类现有的分离物理问题，虚拟了一个模型，对该模型三级分离展开气动数值模拟，研究了三级分离过程中的一些共性问题和特性问题。重点对分离姿态、各子级初始分离速度以及挡板干扰3个影响因素展开研究，结果显示：
1) 改变三级分离姿态角，发现大迎角可能会导致短时间的失稳状态扩大化，不利于子弹分离时的稳定性。
2) 调整挡板及子弹初始分离速度，研究分离方案的可行性与稳定性，提出了几种具体可行的速度优化方案。
3) 干扰区内头罩抛离及挡板平抛对后续子弹抛撒有影响，挡板的存在减小了子弹分离时姿态角变化，有利于其保持相对稳定的分离状态。
4) 合理控制分离姿态、各子级初始分离速度以及挡板干扰，可以实现三级安全分离。
[1] 王巍,刘君,郭正.子母弹抛壳过程非定常流动的数值模拟[J].空气动力学学报,2006(03):285-288.
Wang W,Liu J,Guo Z.Numerical simulation of unsteady flow during shell casting of submunition[J].Chinese Journal of Aerodynamics,2006(03):285-288.
[2] STALMAEH CJ,ROMERO S.Dynamie tests of a submissile dispensing system in a supersonic wind tunnel[C]//Proceedings ofthe 11th AIAA Aerodynamic Testing Conference.Reno,NV,US:AIAA,1980:0452.
[3] RIZK M,JOLLY B.Aerodynamic simulation of bodies with moving components using CFD overset gird methods[C]//Proceedings of the 44th AIAA Aerospace Sciences Meeting &Exhibit 2006.Reno,NV,US:AIAA,2006:1252.
[4] BHANGE N P,SEN A,GHOSH A K.Technique to improve precision of kinetic energy projectiles through motion study[C]//Proceedings of Atmospheric Flight Mechanics Conference.Chicago,IL,US:AIAA,2009:5720.
[5] WOODEN P A,BROOKS W,JUBARAJ S B.Calibrating CFD predictions for use in multiple store separation analysis[C]//Proceedings of the 36th AIAA Aerospace Sciences Meeting &Exhibit.Reno,NV,US:AIAA,1998:0754.
[6] PERKINS S C,DILLENIUS M F E.Supersonic submunition aerodynamics during dispense[C]//Proceedings of the 26th Aerospace Sciences Meeting.Reno,NV,US:AIAA,1988:0335.
[7] CAVALLO P,LEE R,HOSANGADI A.Simulation of weapons bay separation flow fields using unstructured grids[C]//Proceedings of the 19th Applied Aerodynamics Conference.Norfolk,VA,US:AIAA,1999:33412.
[8] 张玉东,纪楚群.子母弹分离非定常气动特性数值模拟[C]//近代空气动力学研讨会论文集.北京:中国宇航出版社,2005:333-338.
Zhang Y D,JI C Q.Numerical simulation of the unsteady aerodynamic characteristics of the submunition separation[C]∥Proceedings of Modern Aerodynamics Symposium.Beijing:China Astonautic Publishing House,2005:333-338.
[9] 陈时通,于勇.超声速流场中子母弹分离与子弹姿态变化规律的数值模拟[J].弹道学报,2019,31(01):29-35.
Chen S T,YU Y.Numerical simulation of attitude motion of submunition during dispensing in supersonic flow field[J].Journal of Ballistics,2019,31(01):29-35.
[10] 袁武.子母弹抛撒流场仿真和气动效应研究[J].系统仿真学报,2016,28(07):1552-1560.
Yuan W.Numerical simulation and aerodynamic effect research for multi-warhead projection[J].Journal of System Simulation,2016,28(07):1552-1560.
[11] 王金龙,王浩,陶如意,等.子母弹不同舱段分离流场特性及运动特性研究[J].空气动力学学报,2016,34(04):490-496.
Wang J L,Wang H,Tao R Y,et al.The movement characteristics analysis ofinterference flow field on the separation of multi-bay clustermunition[J].Acta Aerodynamic Sinica,2016,34(04):490-496.
[12] 张曼曼,姜毅,程李东，等.基于嵌套网格的超声速子母弹分离数值分析[J].兵工学报,2019,40(1):79-88.
Zhang M M,Jiang Y,Cheng L D,et al.Numerical analysis of separation of supersonic submunition based on nested grid[J].Journal of Ordnance Equipment Engineering,2019,40(01):79-88.
[13] 张博孜，王珊珊,黄蓓,等.扰流板对子母弹翻转角速度的影响[J].南京理工大学学报,2015,39(01):7-14.
Zhang B Z,Wang S S,Huang P,et al.Influence of spoiler on overturning angular rate of submunition[J].Journal of the Nanjing Institute of Technology,2015,39(01):7-14.
[14] HUANG K.The jet flow field of submunition dispersal process and its influence on the projectile motion[D].Nanjing:Nanjing University of Science and Technology,2015:8-11.
[15] 王金龙,王浩,江坤，等.跨声速下子母弹囊式抛撒干扰特性[J].北京：北京理工大学学报,2017,37(04):348-353，359.
Wang J L,Wang H,Jiang K,et al.Analysis on the gasinterference characteristics during the transonic gasbag separation of cluster munition[J].Beijing:Journal of the Beijing Institute of Technology,2017,37(04):348-353,359.
[16] LI D,HE Y L,JI C Q.Numerical simulation of multibody separation and its application[C]//Proceedings of the 19th Annual Conference of Beijing Institute of Mechanics.2013:121-122.
[17] 李盾,纪楚群.多体全自由飞过程的数值模拟[J].空气动力学学报,2009,27(z1):90-94.
Li D,Ji C Q.The numerical simulation of multi-body’s free flight[J].Chinese Journal of Aerodynamics,2009,27(z1):90-94.
[18] 杨益农,李盾.高动压头罩分离气动特性研究[C]//北京力学会第18届学术年会论文集.2012:135-136.
YANG Y N,LI D.Research on aerodynamic characteristics of high dynamic pressure head cover separation[C]//Proceedings of the 18th Annual Conference of Beijing Institute of Mechanics.2012:135-136.
[19] 李盾,纪楚群,马汉东.三维非结构贴体直角网格的N-S方程数值模拟[J].空气动力学学报,2006,24(04):477-481,490.
Li D,Ji C Q,Ma H D.Numerical simulation of three-dimensional unstructured rectangular mesh with N-S equation[J].Chinese Journal of Aerodynamics,2006,24(04):477-481,490.
[20] 刘君,刘瑜,陈泽栋.非结构变形网格和离散几何守恒律[J].航空学报.2016，37(08):2395-2407.
Liu J，Liu Y,Chen Z Y.Unstructured deforming mesh and discrete geometric conservation law[J].Acta Aeronautica et Astronautica Sinica,2016,37(08):2395-2407.
[21] 李盾.复杂外形三维非结构粘性直角网格的生成[J].计算机与数字工程,2007,35(03):1-10.
Li D.A new method for generating three-dimensional unstructured viscous rectangular meshes with complex shapes[J].Computer and digital engineering,2007,35(03):l-10.