基于嵌套网格的超声速子母弹分离数值分析

张曼曼， 姜毅， 程李东， 杨桦， 刘琦

(1.北京理工大学 宇航学院, 北京 100081； 2.北京航天发射技术研究所， 北京 100071)

0 引言

子母弹武器系统是一种具有较高毁伤率的常规武器系统，在区域战争中被广泛应用，受到世界各国的广泛关注[1]。子母弹武器系统的作战效能由多个子弹之间的散布效果决定，而如何使多个子弹正常分离以达到预期的散布效果，是子母弹武器系统研制过程中面临的重要问题之一[2]。为了保证子弹正常安全分离、提高子母弹武器系统的作战杀伤能力，需要对子母弹分离流场进行分析并探讨其流动机理，研究子弹分离过程的气动特性。影响子母弹分离过程的主要因素，一方面是刚分离初始阶段子弹和母弹的飞行参数，如母弹飞行速度、子弹分离速度、子弹攻角、子弹排列方式等；另一方面是分离过程中多个子弹相互间复杂的气动干扰。子母弹间复杂流场的研究对子弹间的分离方式、排列方式、弹体结构设计等有着重要的指导意义[3]。

对于子母弹这一复杂的多组合体分离流场研究，常用的方法有风洞试验、数值模拟，最早研究可追溯至20世纪60年代初期，20世纪80年代得到快速的发展，此时美国已率先开展有关子母弹分离的风洞试验和数值模拟研究[4-6]。WOODEN等[7]采用数值计算流体力学与风洞试验相结合的方法对子母弹气动特性进行了研究；Perkins等[8]采用微型计算机程序精确模拟了马赫数为1.2时带空腔的多个子弹弹体分离流场；Cavallo等[9]采用非结构动网格方法，利用总变差不增加的有限体积法和并行计算技术对子母弹分离过程进行了数值模拟，数值模拟结果与试验结果符合较好。国内对子母弹分离流场的起步研究相对较晚。林靖明等[3]采用重叠网格技术生成计算区域网格，利用无波动、无自有参数、耗散的显示差分格式(NND)对三维Euler方程进行离散，研究了母弹激波对子弹气动特性及飞行姿态的影响；郭正[10]采用非结构动网格技术求解Euler方程，对高超声速飞行器助推器分离过程进行模拟；张玉东等[11]以分区拼接网格数值模拟方法为基础，模拟了子弹的分离过程以及子弹分离干扰流场的气动特性；王金龙等[12]采用非结构动网格方法，对时序抛撒方式下子弹在不同舱段分离的三维非定常流场进行数值模拟，得到了不同舱段下子母弹分离流场的干扰特性及子弹药气动参数变化曲线；陶如意等[13]采用风洞试验和数值模拟相结合的方法对超音速子母弹分离干扰流场特性进行了研究。以上试验研究均可以给出复杂流场中关键位置的具体数值参数，但是研究周期长、成本高；与试验方法相比，数值计算可以求解子母弹更加全面的分离流场及气动干扰特性，周期短、成本低。已有研究多采用非结构重构网格技术，所需非结构动网格数量多、计算量较大，且在子弹运动幅度较大时网格更新困难，计算不易收敛。现有采用嵌套动网格技术的研究中，嵌套动网格在边界运动时不需要重新生成网格，可以保证网格质量，同时减小计算量。然而现有研究中仅求解了无黏性的Euler方程，对于高速飞行的子母弹，黏性阻力的作用非常明显，因此求解黏性非定常Navier-Stokes (N-S)方程组才可获得更加真实的解。

本文采用结构和非结构混合的嵌套动网格技术，数值模拟了子母弹的非定常分离过程，研究了子母弹分离流场的气动特性，在此基础上分析子弹初始分离状态(如初始分离速度、初始攻角等)对分离过程的影响，为相关工程应用提供了理论指导。计算结果表明，该方法对此类复杂多体分离问题具有良好的适应性。

1 数值模拟方法

1.1 流场求解

采用数值计算方法研究子母弹分离过程，其理论基础是流体力学基本方程N-S方程，以有限体积流体微元为研究对象，建立非定常三维黏性可压缩性的N-S方程组[14]，用张量形式可表示为

1)质量守恒方程

(1)

式中：ρ为流体密度；t为时间；uj、xj分别为沿坐标轴j方向的速度、空间直角坐标。

2)动量守恒方程

(2)

式中：p为静态压力；u为速度向量；i、j为坐标轴方向，xi为沿i轴方向的空间直角坐标；ui为沿i轴方向的速度；应力张量

(3)

其中，μt为湍流黏性系数，δij为克罗内克算子，k为湍流动能。

3)能量守恒方程

(4)

式中：E为单位质量的内能；qj为热通量，qj=-λ∂T/∂x，λ为热传导系数，T为温度。

基于有限体积法求解流场控制方程，空间上采用2阶迎风格式，时间上采用1阶后向差分格式，离散偏微分方程组并耦合隐式算法求解代数方程组，得到流场数值解。

由于子母弹是在超音速条件下分离的，其雷诺数远大于下临界雷诺数，因此有必要引入湍流模型，考虑湍流对分离过程的影响。选择基于布森斯涅克(Boussinesq)提出的涡黏性假设[15]的单方程Spalart-Allmaras(SA)湍流模型，其基本思想是用一个能包含广义黏性系数的涡黏性系数来建立雷诺应力与平均速度梯度之间的关系：

(5)

式中：

(6)

1.2 子母弹运动

子母弹运动的求解采用6自由度刚体运动方程组[17]，其中包括质心运动动力学方程组、绕各轴转动的动力学方程组、弹头质心运动学方程组、姿态角角速度方程组等，即：

质心平移运动动力学方程

mv=Fa+G，

(7)

弹体绕质心运动动力学方程

(8)

弹体平移运动学方程

(9)

弹体转动运动学方程

(10)

Runge-Kutta法是一种单步算法，计算精度很高，且计算过程中便于改变步长，可自启动计算微分方程组，在实践中应用很广[18]，故采用4阶Runge-Kutta法求解6自由度微分方程组，其算法如下：

考虑微分方程组

(11)

4阶Runge-Kutta法求解公式为

(12)

式中：下标i表示变量名，下标j表示时间步。

1.3 模型有效性验证

为了检验上述流场模型和算法的有效性，采用1.1节的流场模型和算法求解二维可压缩跨音速RAE2822翼型绕流，来流马赫数为0.725，攻角为2.92°，基于翼型弦长得到的雷诺数为6.5×106. 图1给出了翼型表面压力系数分布，试验数据来自文献[19]。由图1可知，数值计算结果与试验值吻合度较高，表明数值计算所选择的流场模型和算法是有效的。

2 子母弹分离过程的数值模拟

下面采用嵌套动网格技术对子母弹的相对分离运动进行非定常数值模拟研究。初始状态下，母弹的密封罩完全脱离，子弹暴露于气流中，不考虑密封罩脱离过程带来的扰动；计算结束时，各子弹之间的相对分离距离已经足够远，因此子弹体之间的气动干扰作用可以忽略；由于各个子弹之间的分离运动主要是径向的，相比之下，其滚转、偏航位移程度较小，因此子弹的滚转运动、偏航运动可以忽略。

2.1 计算模型

计算的模型为母弹底托和7个子弹的装配体，如图2所示。各子弹编号分别为1～7，编号为1的子弹位于中心，其余的沿周向均匀分布，每个子弹尾部均匀分布有5个翼。

建立两个坐标系：局部坐标系建立在各子弹的质心位置，从弹尾指向弹头为Z轴正向，YZ面通过弹体对称面，X轴由右手法则确定；全局坐标系与初始位置1号弹的局部坐标系重合。

2.2 边界条件

在数值计算中，由于各个子弹的分离是在超声速条件下，周围空气流动边界取压力远场，弹体表面使用黏性边界条件即无滑移条件，物面上切向速度为0 m/s、法向为无穿透条件。

考虑到装配体的外形包络为圆柱形，建立圆柱形的计算域。为满足远场边界条件的要求，计算域半径为1.8 m；为了给弹体留下足够的运动空间，计算域长度为7.5 m，右侧阴影部分为弹托和子弹，外侧为计算域(见图3)。

2.3 计算网格

网格质量直接影响到流场计算的精度，因此生成高质量网格是获得良好数值解的前提，对于子母弹复杂组合体，采用嵌套混合结构、非结构网格对计算区域进行网格划分。首先根据子弹外形及各个子弹之间的位置关系，生成子弹的嵌套子网格、子弹周边的加密网格以及计算域的背景网格。然后选定背景网格区域和嵌套子网格区域，创建重叠网格界面，设定插值方式为最小二乘插值，进行网格装配，装配完之后重叠网格中的绝大多数背景网格已被删去，只留下少量用于插值的网格，最小网格尺寸为0.5 mm(见图4)。

3 计算结果分析

对相同的子母弹外形采用相同的计算网格及边界条件，分别对子弹在850 m/s、1 000 m/s两种脱离母弹下的飞行速度工况以及在1 000 m/s飞行速度下0°、2°两种攻角工况进行了非定常数值模拟，表1为弹体分离运动的计算工况。

表1 弹体分离运动的计算工况

3.1 子母弹分离流场分析

在超声速初始分离条件下，子弹头部产生较强的斜激波，斜激波相交并反射，反射激波与弹体表面相互作用，随后激波从弹体表面反射相交而产生干扰，如此弹体之间的斜激波反复产生反射相交的运动；在靠近弹体尾翼的过程中，各个子弹之间的激波强度逐渐减弱，致使在弹翼前缘处形成低压区域,如图5所示。对比弹体表面整体压强分布，弹体尾翼处绝对压强数值较低；由于母弹底托的存在，7个子弹尾翼与底托之间形成漩涡区域,如图6所示，气流速度梯度较大，使得该区域的压力梯度较高，由此产生迫使各个子弹散开分离的气动力和气动力矩。

初始分离状态下，1号子弹周围的压强分布基本均匀，1号子弹不会产生明显的偏航运动和俯仰运动，如图7所示，t1～t4为分离过程中4个典型时刻。由于1号子弹与各个外围子弹体之间的斜激波相交干涉作用，2号～7号子弹体尾翼外围压强高于内侧压强，由此将产生迫使各个子弹的弹尖散开的气动力矩。

在气动力和气动力矩的综合作用下，各个子弹逐渐开始分离运动，分离初始阶段弹体表面仍存在较强的斜激波。在子弹装配体分离过程中，各个弹体表面及其周围的压强不断在变化，尤其是2号～7号子弹体尾翼两侧的压强，一面处于高压，另一面则处于低压，即总有一个迫使弹体散开和旋转的气动力矩存在。从图7中可以看出，1号子弹体各个尾翼两侧的压强变化始终处于相对均匀状态。

工况2下子母弹分离流场的变化规律与工况1大致相同，不再重复叙述。

工况3下的流场变化规律与工况1对比出现不一致性，不同点在于：工况3下子弹体头部的斜激波强度增强，使得子弹体之间的气动耦合作用增强，母弹底托表面的压强增大，底托受到较大的作用力，进一步加快了其与各子弹体之间的相对分离运动；工况3下中间子弹体与其上方、下方各个子弹体之间的压强数值大小分布不对称，致使各个子弹体之间的气动耦合作用出现差异，尤其是中间子弹的上方、下方子弹之间，即2号、3号、7号子弹体相对1号子弹之间的径向运动幅度要大于4号、5号、6号子弹，如图8所示，外圈各个子弹相对中间1号子弹向外分离运动的幅度大小明显不同。由此可知，攻角的存在不仅加剧了各个子弹之间的气动耦合效应，还使得流场不对称性进一步加强，从而使得2号～7号子弹体相对1号子弹径向分离的空间位移变化出现不对称现象。

3.2 子弹分离过程的气动参数对比分析

3.2.1 子母弹初始分离时子弹飞行速度对多个子弹分离过程的影响分析

由3.1节工况1下的流场变化规律分析可知：2号～7号子弹的运动变化规律基本一致，而1号子弹的运动变化与外围弹体有较大差异，因此选取1号和2号子弹为代表，分别对比分析两个子弹体在不同初始飞行速度下的数值模拟气动参数，进而分析不同的初始分离速度对各个子弹分离过程的具体影响。1号、2号子弹能否安全快速分离主要取决于弹体径向与轴向的气动力和力矩。

表2给出了1号、2号子弹在工况1、工况2下分离过程中Y方向的速度及作用力(全局坐标系)与X方向的力矩(局部坐标系)变化曲线。从表2中可以看出：在同一种工况相同时间内，气动力、角速度等气动参数表现出明显的周期性振荡，这是因为弹体内外压强交替变换所致(见图7)，由此说明1号、2号子弹体的分离运动是一种往复摆动过程。由气动力与力矩曲线可知，2号子弹所受气动力与力矩大于1号子弹；从3.1节的流场分析可知，1号子弹不会产生明显的偏航运动和俯仰运动，因此说明外圈子弹径向运动的幅度均大于中间1号子弹。

通过工况1、工况2下各个子弹分离过程的仿真气动参数结果分析可知：子母弹初始分离时子弹飞行速度的变化对中间子弹分离运动姿态变化的影响大于周围弹体，具体表现在工况2下1号子弹气动参数数值大小围绕工况1下的气动参数曲线上下波动的幅度较大，这主要是因为弹体结构周向不对称和周围的子弹弹翼分布不对称；子弹飞行速度的增大加快了各个子弹体之间的相互分离运动，具体表现在工况2下2号子弹的气动参数变化波峰波谷出现的时间要比工况1提前，1号、2号子弹体之间的气动耦合作用时间缩短(见表2)；子弹飞行速度的增大加大了2号子弹体所受的气动力与力矩，2号子弹相对1号子弹的径向分离运动幅度均有所增大。

3.2.2 子母弹初始分离时子弹飞行攻角对多个子弹分离过程的影响分析

由3.1节工况3下的流场变化规律分析可知：2号、3号、7号子弹的运动变化规律基本一致，4号、5号、6号子弹的运动变化规律基本一致，而1号子弹的运动变化与其有较大差异，因此选取1号、2号和5号子弹为代表，分别对比分析3个子弹体在其不同初始分离攻角下的数值模拟气动参数，进而分析不同初始分离攻角对各个子弹分离过程的具体影响。弹体间能否安全快速分离，主要取决于与弹体沿径向、切向运动相关的作用力及力矩。以下着重分析与1号、2号、5号子弹径向、切向运动相关的作用力与力矩等气动参数。

表2 1号、2号子弹在工况1、工况2下分离过程的气动参数变化曲线

表3给出了1号、2号、5号子弹在工况2、工况3下分离过程中X轴方向、Y轴方向气动力与X轴方向力矩变化的曲线，其中气动力为全局坐标系，气动力矩为弹体局部坐标系。从表3中可以看出：1号子弹在工况3下的气动力与力矩数值围绕工况2下的气动参数曲线上下来回波动，与前面3.2.1节对比分析1号子弹在工况1、工况2下的气动参数的变化趋势一致；工况2下2号、5号子弹的气动参数数值变化曲线基本吻合，这是因为2号、5号子弹体对称分布于1号子弹两侧，在无攻角条件下两子弹体与周围各子弹体之间的气动作用力大小相等。

通过对工况2、工况3下各个子弹分离过程的对比分析两种工况下2号、5号子弹Y轴方向的位移曲线可知，工况2下2号、5号子弹相对1号子弹的切向分离运动的幅值变化基本相同，即Y轴方向的运动位移数值曲线基本吻合；对比工况2，工况3下2号子弹的径向位移变大，5号子弹径向位移变小，即2号子弹向外散开，5号子弹向内靠拢，严重不利于各个子弹之间的安全分离。由此可知，攻角的增大加大了各个子弹之间气动力的不稳定性，使得各个子弹之间受力不均匀，具体如图9、图10所示。

表3 1号、2号、5号子弹在工况2、工况3下分离过程的气动参数变化曲线

仿真气动参数结果分析可知：初始攻角的增大同样也会加快各个子弹体之间的相互分离运动，具体表现在工况3下1号、2号、5号子弹的气动参数变化波峰波谷出现的时间要比工况2提前，即各个子弹体之间的气动耦合作用时间变短，原因与3.1节对工况1下的流场分析相同；初始攻角的增大不仅加强了子弹之间的气动耦合效应，还加大了各个子弹气动受力的不稳定性，各个子弹所受气动参数大小出现明显差异，具体表现在工况3下2号、5号子弹的气动参数数值比工况2大，且工况3下2号子弹所受的气动力与力矩数值均大于5号子弹，2号子弹气动参数变化曲线的周期性振荡幅度大于5号子弹，由此出现外圈各个子弹相对中间1号子弹向外分离运动的幅度大小不一致现象，尤其是中间子弹体上方的子弹，这与前面3.1节描述的攻角增大的分离流场特点相一致。

4 结论

本文采用嵌套网格技术对子母弹在超声速条件下的分离运动作了研究和分析，得出主要结论如下：

1)采用嵌套动网格技术，利用有限体积法和4阶Runge-Kutta法求解采用SA湍流模型的雷诺平均N-S方程和6自由度刚体运动方程，可以很好地模拟超声速下子母弹相互分离的复杂干扰流场。

2)在超声速初始分离条件下，各个子弹头部形成较强的斜激波，斜激波相交产生干扰；激波反射到子弹弹体上，且在相邻两个子弹之间发生反射，产生复杂气动干扰效应，如此反复反射相交；从弹头到弹尾，激波强度逐渐减弱，并在弹翼前缘处形成低压区域；由于母弹底托的存在，7个子弹尾翼与底托之间形成漩涡区域，气流速度梯度较大，导致该区域的压力梯度较高，形成各个子弹分离的初始气动力和气动力矩。

3)初始分离时子弹飞行速度越大，各个子弹之间的气动耦合效应越强，外圈子弹与中间子弹分离运动所需时间缩短，中间子弹与周边子弹体之间的径向分离运动幅度加大，因此可以通过提高各个子弹体的初始飞行速度，保证子弹体之间的安全快速分离。

4)随着子弹初始分离飞行攻角的增大，各个子弹之间的气动耦合耦合作用增强、子弹分离时间加快的同时，子弹体之间的分离运动状态变化不稳定性增强，气动参数大小出现明显差异，外圈子弹相对中间子弹向外分离运动的幅度大小不同，即中间子弹体上方子弹向外分离运动的幅度变大，而中间弹体下方子弹向外分离运动的幅度减小，且有向中间子弹靠拢的趋势，严重不利于各个子弹之间的安全分离。因此在各个子弹初始分离时刻，为了保证各个子弹体之间可以安全快速地分离，应尽量减小子弹的初始攻角。