基于计算流体力学与刚体动力学耦合的高速旋转弹丸弹道计算方法

钟阳，王良明，吴映锋,2

(1.南京理工大学 能源与动力工程学院, 江苏 南京 210094; 2.63961 部队, 北京 100012)

0 引言

先进的武器装备是国防建设的必要条件，射程、精度、威力等是武器系统的主要战术技术指标，武器系统的研制常从弹道设计[1]开始，因此弹道计算在整个武器系统研发过程中有着非常重要的地位。传统的弹道计算方法是利用基于气动力和力矩系数的弹道模型，通过龙格- 库塔(Runge-Kutta)法或阿当姆斯预报- 校正法进行数值求解。传统弹道计算方法需要建立弹箭气动力模型，对于旋成体弹丸，气动模型相对准确；而对于气动外形复杂甚至非对称的弹箭，建立准确的气动模型非常困难。因此，对于外形复杂的弹丸，如何准确获取其动态变化过程和弹道飞行特性，是弹道学所面临的一个难题。

近十多年发展起来的计算流体力学(CFD)和刚体动力学(RBD)耦合计算方法，为弹道仿真提供了一种新的技术途径。相对于传统方法，该方法能够实时模拟弹箭在大气中真实的飞行情况，即使再复杂的弹箭，也能从实时变化的流场参数中获取气动力和力矩，在无需提供气动模型情况下完成弹道计算。另外，该方法能够实时提供弹箭飞行过程中弹体表面和周围空气的压力、密度、温度等全息流场参数，为研究人员提供丰富的数据资料，有利于新型弹箭的研制。文献[2]通过欧拉角建立弹体坐标系与地面坐标系的关系，在弹体坐标系下采用任意拉格朗日- 欧拉(ALE)控制方程对流场进行求解，首次将CFD/RBD耦合技术成功应用到弹道计算中。进而以某尾翼稳定弹为对象，采用阿当姆斯预报- 校正法对流场和6自由度刚体动力学方程组进行耦合求解，并通过试验结果验证了方法的有效性。该方法被进一步应用于超音速、跨音速和亚音速弹丸模拟中，通过弹道仿真结果和试验结果对比再一次验证了方法的正确性[3-5]。文献[6]在耦合CFD/RBD技术中引入入口边界条件，研究了喷流对弹道的影响规律，验证了在耦合计算中加入进口边界来模拟脉冲修正弹弹道的可行性。文献[7-11]通过弹丸飞行动力学模型推导出气动参数估计模型，根据耦合CFD/RBD计算结果估计出弹丸的气动系数。文献[12-16]在耦合CFD/RBD技术中加入飞行控制系统(FCS)，建立耦合CFD/RBD/FCS方法并模拟了鸭舵受控尾翼弹的姿态运动，通过与实验数据对比验证了方法的正确性。文献[17]采用CFD/RBD和非结构动网格技术计算均匀流场下不同密度平板的运动轨迹，分析了碎片对飞机造成的二次伤害。文献[18]基于自主开发软件HUNS3D，结合非结构嵌套网格技术，使用阿当姆斯预估- 校正法耦合求解6自由度刚体运动方程，研究了子母弹分离过程。文献[19]采用刚性动网格和CFD技术，在规定运动下模拟了带鸭舵弹丸的运动过程，研究了锥形运动下弹丸非定常空气动力学特性。

综上所述，文献[19]采用的是规定弹丸运动的方法，不能直接用于研究弹丸的自由运动。文献[17-18]这类方法的计算域通常包含物体运动轨迹，并不适用于研究射程达几十千米的弹丸。相比较而言，文献[2-16]基于弹体系建立流动关系更为合理。然而，对于高速旋转弹丸，基于弹体系建立流动关系也存在不足之处，高速旋转弹丸会使远场产生很大的附加速度，造成远场计算不准确[20]。现有的CFD/RBD耦合方法一般采用阿当姆斯预报- 校正法，该方法较容易实现，但不能自启动[1]。在流场初始耦合时只能采用当前值来代替过去值，这必将在初始阶段损失时间精度，且不易变步长。

本文针对已有研究的不足，基于ALE形式的流动模型，建立控制体表面运动和弹轴坐标系运动的耦合模型。在此基础上提出一种基于4阶Runge-Kutta法的CFD/RBD耦合计算方法，为研究高速旋转弹丸的真实弹道和流场提供参考。

1 CFD数值计算方法

1.1 基于弹轴运动的ALE形式控制方程

弹丸飞行时姿态不断变化，与弹丸固连的计算网格具有附加速度，ALE形式的Navier-Stokes方程(简称N-S方程)能够有效地求解具有网格速度的流场。由于弹丸高速旋转，在弹体系下建立流动关系时，远场流动速度相对于壁面附近非常大，从而导致计算不准确[20]。而弹轴坐标系有效去除了弹体高速自转，克服了这一问题。基于弹轴运动建立具有ALE形式的N-S方程并对流动进行求解，其积分形式为

(1)

(2)

Fc为静态网格下的对流通量，v∂Ω为附加在控制体表面的反变速度。

当弹轴系的平动速度矢量和转动角速度矢量分别为vp和ωa时，附加在控制体表面的反变速度为

v∂Ω=n{vp+ωa×(r-rg)}，

(3)

式中：n=[nx,ny,nz]T为控制体表面的单位外法向矢量，nx、ny、nz分别为单位外法向矢量在x轴、y轴、z轴上的分量；r为控制体表面微元的矢径；rg为弹丸质心矢径。

1.2 修正简单低耗散迎风矢通量分裂格式

(1)式中对流通量的计算方法是CFD的研究热点之一。简单低耗散迎风矢通量分裂格式(SLAU)是著名迎风型矢通量分裂格式(AUSM)系列的一个发展方向，具有低耗散、格式简洁、低马赫数下无需调整参数等优点，其改进型SLAU2[21]在压力通量中加入了与来流马赫数相关的耗散项，继承了原格式优点，同时实现了比其他全速格式更简洁的形式。引入弹丸运动后，基于SLAU2的计算格式为

(4)

(5)

cij=(cl+cr)/2,

(6)

式中：cl和cr分别为控制面左右两侧音速。

压力通量为

(7)

(8)

(9)

M为法向流动速度与音速比值，M=(v·n)/cij.

1.3 边界条件、时间推进及湍流模型

(10)

(11)

式中：下标i表示控制体第i个面；nf为面的个数；vi为控制体第i个面上的流动速度矢量；ni为控制体第i个面的单位外法向矢量；ci和ΔSi分别为离散后控制体第i个面的音速和面积。

计算(1)式中黏性通量的湍流黏性时采用Spalart-Allmaras(S-A)湍流模型，该模型具有计算量小、稳定性好、计算精度满足工程要求等优点，在气动力计算和航空航天等领域得到广泛应用。S-A模型发展了较多版本，本文通过某弹丸三维外部绕流数值实验发现，标准可压S-A湍流模型[23]对于高雷诺数和低雷诺数的适应能力较强，且在尾部网格不对称情况下能够收敛到较好的结果，因此本文采用该湍流模型计算湍流黏性系数。

2 基于CFD的弹丸飞行动力学模型

CFD坐标系通常以弹头为原点，Oxc轴与弹轴重合指向弹尾为正，Oyc轴垂直Oxc轴向上为正，Ozc轴满足右手法则。其他坐标系无特殊说明均参考文献[1]，弹轴坐标系无特殊说明均指第1弹轴坐标系。旋转稳定弹丸质心运动的动力学方程组通常采用弹道坐标系Otxtytzt下的力系形式[1]。为了使CFD计算的空气动力和空气动力矩能够耦合到弹道模型中，需建立转换关系。首先将CFD坐标系下的力系转换到第1弹轴坐标系下，再将第1弹轴坐标系下力系转到第2弹轴坐标系下，最后将第2弹轴坐标系下的力系转到弹道坐标系下。CFD坐标系到第1弹轴坐标系的转换矩阵为

Lac=diag[-1 1 -1].

(12)

第1弹轴坐标系到弹道坐标系的转换矩阵为

(13)

(14)

式中：vt为弹丸质心运动速度；θa为速度高低角；ψd为速度方向角；φd为弹轴方位角；m为弹丸质量；g为重力加速度；Fcξ、Fcη、Fcζ分别为空气动力在CFD坐标系Oxc轴、Oyc轴、Ozc轴上的分量；Lc为系数矩阵，

Lc=diag[1/m1/mvtcosψd1/mvt].

(15)

弹丸质心运动的运动学方程组[1]为

(16)

式中：vpx、vpy、vpz分别为弹丸质心速度矢量vp在地面坐标系Oexeyeze中Oexe轴、Oeye轴、Oeze轴上的分量。由(12)式和基准坐标系到弹轴坐标系的转换矩阵[1]，得到控制体表面附加平动速度分量为

(17)

式中：vcξ、vcη、vcζ分别为控制体表面附加平动速度矢量；φa为弹轴高低角。

根据(12)式，得到弹丸总空气动力矩矢量M在弹轴坐标系下的3个分量为

(18)

式中：Mcξ、Mcη、Mcζ分别为CFD数值计算获得的弹丸总空气动力矩矢量在CFD坐标系3轴上的投影。将(18)式代入弹丸在弹轴坐标系下绕质心转动的动力学方程组[1]，可得

(19)

式中：C为极转动惯量；A为赤道转动惯量；ωξ、ωη和ωζ分别为弹丸角速度矢量在弹轴坐标系上3个轴上的投影分量。由(12)式得到壁面角速度为

(ωcwξ,ωcwη,ωcwζ)=(-ωξ,ωη,-ωζ).

(20)

弹轴坐标系角速度矢量ωa在弹轴坐标系3轴分量为(ωaξ,ωaη,ωaζ)=(ωζtanφd,ωη,ωζ)[1]，通过转换矩阵(12)式得到控制体表面附加角速度为

(ωcξ,ωcη,ωcζ)=(-ωζtanφd,ωη,-ωζ)，

(21)

对于高速旋转弹丸，由于|ωζtanφd|远小于|ωξ|，由(3)式可知，在弹轴坐标系下控制体表面附加速度远小于弹体坐标系，特别是远场。

弹道方程组中绕质心转动的运动学方程组和δa、δd、β3个角关系方程参见文献[1]，不再赘述。

3 CFD/RBD耦合计算方法

本文弹道和流场耦合计算框架如图1所示，本节主要研究图1中的耦合模块。

图1 CFD/RBD耦合计算框架Fig.1 Framework of CFD/RBD coupling calculation

(14)式、(16)式、(19)式和文献[1]中的绕心转动运动学方程组以及3个角关系方程，共同组成了由CFD数值计算的空气动力和空气动力矩驱动弹道方程。(17)式、(20)式和(21)式将弹道方程解算出的状态量反馈给CFD求解器，从而形成了耦合机制。为叙述方便，(17)式、(20)式和(21)式统称为接口模型。

耦合计算采用经典的4阶Runge-Kutta法。已知n时间层面值(tn,y1,n,y2,n,…,ym,n)，tn表示数值计算到第n步的时间，ym,n表示第n步下第m个参数，步长为h. 则求解n+1时间层面各函数值的公式为

(22)

式中：

(23)

fi表示第i个微分方程的右端表达式。由(23)式可知，每一个Runge-Kutta子步都需要气动数据的支持。有2种较简单的气动力和力矩耦合方法：第1种为定值法，借鉴系数冻结法的思想，将n时间层面CFD计算出的瞬态气动力和力矩冻结，即在Runge-Kutta 4个子步中保持不变；第2种为插值法，利用n-1和n时间层面的气动力和力矩，在4个子步中进行线性外插。本文提出一种紧耦合实现方法，使每一个子步中的气动力和力矩都为CFD计算出的瞬态结果。具体实现步骤如下：

步骤1保存n时间层面的流场参数，并积分出CFD坐标系下弹丸受到的空气动力和力矩，通过弹道方程右端子式计算出ki1.

步骤2将流场参数设置为步骤1保存的n时间层面；根据ki2右端括号中的参数形式，利用步骤1计算出的ki1和接口模型计算出控制体表面附加运动速度。将时间步长设为0.5h，采用双时间步将流场推进到tn+0.5h时刻；积分出弹丸受到的空气动力和力矩。利用弹道方程右端表达式计算出ki2.

步骤3将流动情况返回到n时间层面；根据ki3右端参数形式，利用步骤2计算出的ki2和接口模型计算出控制体表面附加运动速度；双时间步步长设置为0.5h将流场推进到tn+0.5h时刻；积分出作用在弹丸上的空气动力和力矩；利用弹道方程右端表达式计算出ki3.

步骤4将流动数设置为tn时刻的值，通过步骤3计算出的ki3来获得ki4右端参数值，由接口模型获取控制体表面附加运动情况；以h为步长，采用双时间步将流场推进到tn+h时刻并获取弹丸上受到的气动力和力矩；由弹道方程右端表达式计算出ki4.

步骤5由(22)式计算出n+1时间层面的弹道参数，并通过接口模型计算出控制体表面附加运动；将流场返回至tn时刻，以h为步长，利用双时间步将流场推进至tn+h时刻；积分出当前的空气动力和力矩。

图2 流场和弹道紧耦合计算过程Fig.2 Flow field and ballistic tight coupling calculation process

通过以上5步可完成流场和弹道方程从n时间层面到n+1时间层面完整的4阶Runge-Kutta紧耦合计算过程，计算流程示意图如图2所示，图2中j=1，2，3，4表示Runge-Kutta4个步骤。需要注意的是，步骤2～步骤4中涉及到的tn+0.5h和tn+h时刻只是中间过程，步骤5中采用双时间步推进到的tn+h时刻为n+1时间层面。

4 仿真验证与分析

4.1 M549弹丸定态运动和规定运动模拟

计算模型采用美国M549旋成体弹丸，具体结构尺寸参见文献[24]。利用Solidworks商业软件进行几何建模，如图3所示。

图3 M549弹丸模型Fig.3 M549 projectile model

将几何模型导入CFD前处理软件ICEM中进行网格划分，如图4所示。文献[25]对旋成体弹丸的网格划分进行了详细研究，并通过与试验数据对比验证了网格划分的合理性。本文在网格划分过程中部分技术参数参考此文献，如远场位置、近壁第1层网格厚度和黏性层网格延展率等。为使弹丸头部和尾部网格具有良好的正交性，进行O型剖分，网格总数约为186万。

图4 M549弹丸计算域网格Fig.4 Computational domain mesh of M549 projectile

通过使弹丸做定态飞行和规定运动获得弹丸的气动特性，将结果与Spinner数据[24]进行对比，以检验本文CFD数值方法的可靠性。

(24)

图5 俯仰力矩系数迟滞回线Fig.5 Hysteresis loop of pitching moment coefficient

表1 马赫数为2时M549弹丸气动系数Tab.1 Aerodynamic coefficients of M549 projectile for Ma=2

总之，在马赫数为2.0时，静态气动系数和动态气动系数的仿真结果相对Spinner数据的误差分别小于6%和10%，数值仿真的气动系数较可靠。

4.2 M549弹丸自由飞行仿真验证与分析

文献[27]基于传统弹道计算方法，采用旋成体弹丸气动模型和CFD仿真气动系数计算某旋成体弹丸弹道，并与靶道试验结果对比，验证了CFD数值方法结合传统弹道计算方法预测旋成体弹丸实际弹道的可行性。本文借鉴文献[27]的方法计算M549旋成体弹丸弹道，并以该弹道结果为基准，研究第3节中3种力和力矩耦合方法以及不同时间步长对耦合计算方法模拟弹道能力的影响。

弹道起始条件为：射角45°；射向0°；高低攻角和飞行马赫数分别沿用4.1节中的3°和2.0；方向攻角0°；弹轴无摆动。流场仿真条件为：双时间步子迭代最大计算步数为200；收敛条件为10-2；流场初始化采用4.1节中对应条件下的计算收敛结果。

对于传统6自由度弹道计算，时间步长通常取5 ms左右[1]。为了研究3种力和力矩耦合方法求解弹道的能力，取步长h=5 ms进行仿真，得到攻角关系曲线如图6所示。通过观察图6中的攻角变化情况发现：采用定值法效果最差，仿真攻角结果发散，计算弹道失败；插值法和紧耦合法均模拟出一般旋转稳定弹丸的双圆运动，其中紧耦合法计算结果与基准结果更加接近，即紧耦合精度最好。

图6 h=5 ms时3种耦合方法对比Fig.6 Comparison of 3 coupling methods for h=5 ms

对于非定常流场，为了准确反映流动细节，时间步长要小得多。流场与弹道耦合计算的主要矛盾是非定常气动参数的收敛程度，对某155 mm口径弹丸进行非定常气动力计算时发现，时间步长取0.5 ms和0.05 ms时滚转力矩系数差别仅为4.67%[28]，取时间步长为0.5 ms时兼具精度和效率。

图7 弹道参数仿真结果Fig.7 Simulated results of trajectory parameters

采用紧耦合方法，在时间步长分别为5 ms和0.5 ms条件下进行仿真，仿真结果与基准结果如图7所示。由图7(a)可见：总攻角δt初始为3°，在陀螺稳定力矩作用下逐渐衰减；步长5 ms下总攻角计算结果偏大，衰减明显不足；步长0.5 ms下衰减明显加快，且衰减速度能够跟随基准结果；从二者均值上看，初值均在4°左右，当时间为1.805 s时，前者衰减到1.714 8°，后者衰减到1.170 4°，小步长比大步长衰减效率提高了约24%. 从图7(b)和图7(c)中弹轴摆动曲线以及图7(e)中速度方向角曲线可以看出：当步长为5 ms时仿真结果在第1个快圆运动周期内与基准结果基本一致，但之后不能保持，虽然曲线形状基本一致，但幅值偏大；当步长为0.5 ms时，仿真结果与基准结果始终保持较高的符合度。小步长下速度高低角与基准结果基本重合，大步长下在波谷处略低，如图7(d)所示。弹道侧偏ze反映了地面坐标系Oexeyeze中弹丸偏离射击面Oexeye[1]的程度，如图7(f)所示，仿真结果与基准结果符合较好，在t=1.80 s时，0.5 ms步长和5 ms步长下的弹道侧偏量与基准结果偏差分别为2.18%和4.92%.

4.3 修正组件不转时二维弹道修正弹模拟与分析

旋转稳定二维弹道修正弹通过后体高速旋转来保持飞行稳定。控制时，修正组件滚转角固定在某个位置，对弹道进行修正。此时亦可采用本文方法进行仿真。计算模型和网格划分如图8所示，网格划分技术参数与4.1节相同，网格数约260万。仿真时，修正组件滚转角固定在图8(a)位置即滚转角为0°；初始高低攻角为0°；后体初始转速1 500 rad/s；其余条件与4.2节相同。图8中：F1和F3表示一对差动舵受到的空气动力；F2和F4表示一对控制舵受到的空气动力。

图8 二维弹道修正弹计算模型和网格Fig.8 Calculation model and meshes of two-dimensional trajectory correction projectile

图9给出了步长0.5 ms下紧耦合法计算的攻角曲线及特征点P1～P6，攻角变化规律为双圆运动，并与解析方法[29]计算的结果对比，两种方法在频率和幅值方面大致相似，但不完全相同。从攻角收敛趋势来看，控制力向上舵偏角产生的攻角向下，与鸭舵布局的尾翼稳定弹相反，这一结论与文献[29]一致。

图9 攻角仿真结果Fig.9 Simulated result of angle of attack

图10 弹丸表面压力分布Fig.10 Pressure distribution on the surface of projectile

为了分析弹丸飞行过程中表面压力变化情况，观察特征点处的流场，如图10所示，视图方向垂直于弹道坐标系中的Otxtyt平面，特征点参数如表2所示。由起始阶段图10(a)可以看出，与超音速来流直接接触的头部和鸭舵稍部压力较高，最大处超过0.17 MPa；前后体连接处压力分布在0.112～0.156 MPa之间，其中鸭舵正后方处在鸭舵激波内部，压力比周围高；图10(b)中弹轴摆动到了左上方，弹头部压力出现不对称，迎风区压力增大；随着弹轴从左上方摆到左下方，弹头迎风区变成右上方，且由于低头使得水平舵有效迎角减小，压力降低，如图10(c)所示；当弹轴从左下方位置摆到图10(d)中下方略偏右位置时，上方差动舵有效迎角变大同时处于迎风区，下方有效迎角减小同时处于背风区，使得前者压力增大后者压力减小；当弹轴从图10(d)位置略微上摆至图10(e)位置时，水平舵有效迎角增大，压力也增大，直至摆动到图10(f)时，弹丸近似完成一个慢圆周期。紧耦合计算中，随着弹轴摆动，弹丸表面迎风区域和背风区域不断变化，形成时变且复杂的压力分布，压力分布进而又影响弹轴摆动，这与弹丸实际飞行过程是类似的。综上所述可知，与普通弹道计算方法相比，紧耦合法计算不仅能向设计人员提供弹道数据，还能提供流场参数，展现更加实际的弹丸飞行过程。

表2 二维弹道修正弹攻角仿真结果Tab.2 Simulated results of angle of attack of two-dimensional trajectory correction projectile

5 结论

耦合计算流体力学和刚体动力学方法是外弹道学一个新的发展方向。本文考虑高速旋转弹丸的特点，利用壁面速度表征旋成体表面的高速旋转，结合高速旋转弹丸弹道模型，将弹轴坐标系平动速度和转动角速度附加于控制体表面，建立了基于ALE形式新的流动模型。在此基础上提出一种适合于高速旋转弹丸的基于4阶Runge-Kutta法流场和弹道耦合数值计算方法。最后通过对比传统方法的仿真结果和解析结果验证了紧耦合计算方法的有效性。得到主要结论如下：

1)在进行耦合计算时，气动力和力矩耦合模型对计算结果影响较大，采用线性插值和提供实时瞬态值即紧耦合时计算结果收敛性好，精度高。

2)时间步长的选取很重要，在合理时间步长下进行紧耦合计算时，可以有效模拟弹丸飞行过程。

3)紧耦合方法模拟具有复杂外形的旋转稳定弹飞行时，能够反映出弹丸的运动特性，流场变化规律符合客观实际。

本文提出的耦合计算方法可为新型旋转稳定弹，如加装了精确制导组件的二维弹道修正弹设计和计算提供一种新的理论和方法。今后可进一步考虑控制部件的运动，为新型旋转稳定弹控制组件的设计提供基础。