1 引言
电磁轨道炮是一种利用电磁力(能)推进物体到高速或超高速的发射技术[1-3]，目前最高速度可达到5.9 km/s。而传统火药发射技术受火药燃气滞止声速的限制，其炮口初速一般难以突破1.5 km/s。在防空反导系统、反装甲系统等领域[4-5]，作为新型高速发射技术和革新传统火炮的主要技术候选者，电磁轨道炮是当今各国的研究热点。
利用电磁轨道炮发射的弹丸大多有杆式弹丸、制导弹丸、脱壳穿甲弹等类型，其中脱壳穿甲弹能够应用于防空反导以及反装甲领域。脱壳穿甲弹通常由弹体、弹托两部分组成，当弹丸出炮口后，弹托在空气动力与重力的作用下产生一系列气体扰动，流场中产生强烈的激波与膨胀波[6-7]，使弹体的起始扰动增大，进而影响弹体的射击密集度。因此有必要对弹托分离过程开展流体力学和动力学分析，研究弹托能否顺利分离以及对弹体初始运动的影响[8]。
国外的Hayden等[9]通过电磁发射技术将弹体在2～4 km/s超高速的发射条件下成功飞行，并拍摄到试验中弹托分离的过程。Erengil等[10]建立了新的IAT弹托分离弹道模型，预测出弹托表面的压力分布以及弹托脱壳过程中的运动轨迹。Heavey等[11]运用 CFD++、Fluent等软件对脱壳穿甲弹固定位置上的二维流场及三维流场进行仿真，研究了出炮口瞬间弹丸表面的压力分布。国内的李湘平等[12]采用动网格技术耦合流体与动力学方程对一体化弹丸弹托分离进行模拟，获得了不同初速对弹托分离的影响。
目前国内对于电磁发射超高速穿甲弹的脱壳过程的相关研究及文献较少。本文采用动网格技术建立了一种流体力学与6DOF刚体动力学耦合的超高速穿甲弹脱壳数学模型，对电磁发射穿甲弹在不同马赫数下的脱壳过程进行数值模拟，详细分析了超高速条件下脱壳过程中流场的作用特点，获得了不同马赫数下弹托相对弹体的运动轨迹以及弹体和弹托的气动参量随脱壳过程的变化规律，为研究电磁发射超高速穿甲弹弹托分离的动力学问题提供参考。
2 超高速穿甲弹脱壳过程数学模型
电磁轨道炮发射的弹丸出膛口瞬间，若忽略弹托分离所受到的电磁力作用，其主要受到空气动力与重力作用。迎面来的高速空气流会在弹托前端形成正激波，并在前腔形成高压区，弹体会产生斜激波，弹体弹托的激波会相互发生碰撞。弹托在高压和激波的作用下发生分离，流场变化剧烈。弹体会以较高的初速继续飞行，弹托会逐渐偏离原弹道，然后各自飞行。整个脱壳过程涉及到复杂的流场变化以及弹体弹托的刚体运动，基于此，采用动网格技术建立了一种流体力学与6DOF刚体动力学耦合的超高速穿甲弹脱壳数学模型。
2.1 流体力学基本方程
针对所研究的电磁发射超高速穿甲弹脱壳过程有关动网格的流体力学数值计算，采用N-S[13-14]粘性非定常流体控制方程，选择k-ε-Realizable作为湍流模型，能够表征超音速下流体的附着和分离，具体数学模型如下：
质量方程：
(1)
动量方程：
(2)
能量方程：
▽(κ▽T)+▽(τu)
(3)
湍流方程：
Gk-ρε-YM
(4)
(5)
(6)
质量方程中： ρ、u分别为流体的密度与速度矢量；
动量方程中：p、τ分别为流体的压强与平均应力张量；
能量方程中：e、T分别为单位流体体积的内能与温度，其中e=p/(ρ(γ-1))，p=ρRT，γ为理想气体绝热指数，κ为热导率；
湍流方程中：k、ε分别为单位质量湍流脉动能及其耗散率；μ为流体粘度，μt=ρCμ(k2/ε)为湍流粘度；Gk=μtS2为湍动能生成率，其中
为系数；YM=2ρεM2，M为湍流马赫数；σk、σε、C1、C2为湍流模型参数。
2.2 6DOF刚体动力学方程
运用上述流体控制方程求得每一时刻的方程解，并在弹托表面对压力进行积分后，结合6DOF刚体动力学方程求解弹托的分离[15]过程。初始条件需要设置弹托质心的位置，由动力学方程在流场解的基础上求解弹托质心运动的线加速度和角加速度，再求得质心的线速度和角速度，进而得到弹托质心各个时刻的位置和转动的角度。
对于刚体质心的平移运动，在惯性坐标系下求解其方程：
式中：vG表示惯性坐标系下刚体的速度， fG表示惯性坐标系下刚体受到的合外力。
角运动的动力学方程为：
(8)
MB=RMG
(9)
(10)
式中：L为刚体的惯量张量，ωB为体坐标系下刚体的角速度，MB为体坐标系下刚体所受的力矩，MG为惯性坐标系下刚体所受的力矩，R为坐标转换矩阵，ξ、θ、ψ表示绕x、y、z轴转动的角度。
3 基于动网格技术的脱壳过程数值建模
通过Fluent软件进行仿真，模型采用基于密度的耦合求解方式，湍流模型采用k-ε-Realizable模型。对于离散格式，时间项采用基于隐式向后差分格式，空间上采用有限体积法，对流项采用Roe-FDS格式。
3.1 基于有限体积法网格模型建立
电磁发射超高速穿甲弹分为弹体和弹托，弹托由四瓣组成，其形状较为复杂。为了网格生成以及简化计算，省略弹体与弹托之间的螺纹部分，在弹体头部建立坐标系，利用三维软件建成的弹体与弹托的模型如图1所示。建立空气域为矩形，长度为30d，宽度和厚度为15d。
图1 弹丸(a)、弹体与弹托三维模型(b)示意图
Fig.1 Three-dimensional model of projectile body and sabot
首先对弹体及弹托划分表面网格，其次建立空气域的模型后划分空气域与弹体及弹托间的体网格，由于动网格技术的特殊性，故整体用非结构网格进行划分，同时对弹丸附近区域的网格进行加密，总单元数约为70万。弹丸与空气域网格如图2所示，弹体表面网格划分如图3所示，弹托表面网格细节部分如图4所示。
图2 弹丸与空气域网格示意图
Fig.2 Mesh of projectile and air domain
图3 弹体表面网格示意图
Fig.3 Surface mesh of projectile body
图4 弹托表面网格示意图
Fig.4 Surface mesh of sabot
3.2 初始边界条件
物面边界采用无滑移绝热壁面，计算域外边界为压力远场，压力为一个标准大气压，密度为1.177 kg/m3，温度为300 K，重力加速度为9.8 kg/m2。分别计算4Ma、5.3Ma、6.2Ma下的脱壳过程，且攻角为0°，即在对称来流的条件下4个弹托受气动力和重力作用相对于弹体发生6DOF运动。初始边界条件如表1所示。
表1 初始边界条件
Table 1 Initial boundary condition
压力/Pa密度/(kg·m-3)温度/K重力加速度/(kg·m-2)攻角/(°)101 3251.1773009.80
3.3 基于动网格的流体力学与6DOF刚体动力学耦合模型数值解算
动网格技术[16]结合各种网格变形、重构或重叠技术，可实现对流体边界运动或者变形问题的模拟。针对弹托分离过程中流体网格会发生大的变形，选用动网格中的弹簧光顺法和局部网格重构法[17-23]。弹簧光顺法将网格节点间用弹簧相互连接，当节点发生位移后，通过弹簧传递力，使得弹簧连接在一起的节点在新的位置达到平衡，从而模拟流体单元的变形。当位移过大时，网格会产生较大的倾斜变形，网格质量已恶化并有可能影响计算精度，此时再采用局部网格重构法对网格进行重新划分。
在Fluent中设置Dynamic Mesh，动网格区域选择弹托部分，即在整个过程中使弹托的网格进行运动；其次通过UDF编写弹托的质量与转动惯量，使弹托做6DOF运动。
首先求解脱壳过程中每一时刻流体控制方程，得到流体的压力分布，再在弹托表面上对压力进行积分求得气动参数，将其代入6DOF刚体动力学方程，求解出弹托的动力学参量，再由动网格技术中的弹簧光顺法和局部网格重构法计算弹托运动后的位置，从而模拟弹托与弹体的分离过程，求解过程如图5所示。
在动网格求解过程中，对于边界移动的任意控制体积V，φ的通用守恒方程[24-25]的积分形式写为：
(11)
式中：u为流体速度矢量，ug为动网格移动的速度， Γ为耗散系数，Sφ为φ标量的源项，∂V表示控制体积V的边界，A为控制体积面单元向量。
对弹托的6DOF运动可利用4阶多点Adams-Moulton(见式(12))方法[26]对上述式(7)—式(10)离散求得弹托运动的速度与加速度：
(12)
图5 数值求解过程框图
Fig.5 Numerical solution process
4 仿真结果与分析
4.1 压力与密度云图分布
以5.3Ma时的脱壳状态为例，弹丸飞行速度达到1 800 m/s以上，将弹托分离过程分为初始分离阶段、中间分离阶段以及弱耦合与自由飞行阶段。图6至图8采用是x-y剖面的密度等势分布来表示弹托分离过程，本文以y轴正向与z轴正向的弹托为主讨论脱壳过程。图6表示初始分离过程，图7表示中间分离过程，图8表示弱耦合与自由飞行过程。
图6 初始分离过程压力云图
Fig.6 Initial separation process
在初始分离阶段，弹头部分有斜激波，同时弹托前腔的高压使弹托前端形成了脱体激波。弹托在气动力与重力作用下开始分离，弹托前端正激波开始变形为弓形激波(如图6(b)所示)，并且沿着弹体表面由中间部分向尾部延伸，导致弹体表面的压力增加(如图9所示)。弹托前腔的高压区以及弹托其他凸起部分的斜激波能够为弹托提供翻转力矩。
图7 中间分离过程压力云图
Fig.7 Intermediate separation process
图8 弱耦合与自由飞行过程压力云图
Fig.8 Weak coupling with free flight processes
图9 t=0.5 ms时的压力云图
Fig.9 Pressure cloud graph at t=0.5 ms
在中间分离阶段，弹体与弹托之间的斜激波与弹体发生碰撞，且碰撞程度随着弹托与弹体之间间隙的增大而减小。从图7(a)可看出弹体与弹托之间的斜激波与弹体中间部分发生碰撞形成反射激波，而后又与弹托内表面中间部分发生碰撞，此过程中弹托与弹体之间有多道反射激波，使得弹体与弹托内表面压力均增大(如图10所示)，在一定程度上为弹托分离提供径向力，此时弹体气动参数变化剧烈。
图10 t=0.78 ms时的压力云图
Fig.10 Pressure cloud graph at t=0.78 ms
在弱耦合与自由飞行的阶段，弹托内表面逐渐成为主要受力的迎风面，同时受力面的压力逐渐增大。弹托的斜激波与弹体尾部的激波交汇，使弹体周围气流的密度与压强增大，如图8(a)、图11所示。在t=1.98 ms时弹托分离到与弹体近似垂直的姿态，此后弹体的气动参数变化受弹托影响较小，两者之间的干扰消失，弹体与弹托各自做自由飞行。
图11 t=1.78 ms时的压力云图
Fig.11 Pressure cloud graph at t=1.78 ms
4.2 压力矢量分布
图12～图15分别为6.2Ma条件下不同脱壳位置的压力矢量图。图中的截面为xoy面与xoz面的角平分面,来流速度为2 148.23 m/s，压力为0.101 MPa。由图12、图13可得，在脱壳初始阶段主要在弹托迎风槽、弹托尾部与弹体尾部形成涡流，即A、E、G处，压力为：6.17 MPa，0.054 MPa，0.082 MPa。因此涡流处形成较大压差，从而为弹托分离提供翻转力矩。H处的压力达到3.15 MPa，结合图6、图9云图，弹体前半段表面承受较大的压力。
由图14、图15可得，在脱壳中间阶段主要在B、D、G处形成涡流，压力分别为：0.12 MPa，0.081 MPa，0.32 MPa。如图15所示，在E处与G处之间气流压力达到1.23 MPa，结合图7、图10云图，得出激波反射叠加导致弹体尾翼处压力较大，进而影响其气动参数的变化，对弹体飞行产生干扰。
图12 Ma=6.2、t=0.1 ms的压力矢量图
Fig.12 Pressure vector diagram with Ma=6.2 and t=0.1 ms
图13 A处涡流压力矢量图
Fig.13 Local vortex diagram
图14 Ma=6.2、t=0.78 ms压力矢量图
Fig.14 Pressure vector diagram with Ma=6.2 and t=0.78 ms
图15 E、G间的压力矢量图
Fig.15 Vector diagram of the pressure between E and G
4.3 不同马赫数下气动参量随脱壳过程的变化规律
图16、图17分别表示不同马赫数下分离过程中弹托阻力系数与升力系数随时间变化的趋势。由图可知，2个曲线变化趋势基本相同。在初始及中间分离阶段时，由于弹托俯仰角的增加，弹托内表面逐渐成为迎风面，阻力系数逐渐增大，升力系数在一定范围内波动，最大值为3.5～4。在弱耦合阶段时，由于弹托继续向后翻转使得受力面积减小，弹托阻力系数在弹托内表面与来流垂直时增加至峰值，此后开始下降；而弹托翻转导致垂直于弹托方向上的投影越来越小，升力系数逐渐减小，直到弹托内表面与来流垂直时减小为0，此后升力改变方向并逐渐增大。同时，马赫数越大，弹托能更快翻转到与来流垂直的姿态，即弹托阻力系数与升力系数达到峰值的时间越短，且峰值越小。
图16 弹托阻力系数曲线
Fig.16 Sabot drag coefficient vs.time
图17 弹托升力系数曲线
Fig.17 Sabot lift coefficient vs.time
图18表示不同马赫数下弹托俯仰力矩系数随时间变化的趋势。由图可知，变化趋势基本相同，且变化规律的分析与阻力系数和升力系数类似，可得出，马赫数越高，俯仰力矩系数达到峰值的时间越短，且峰值越小。
图18 弹托俯仰力矩系数曲线
Fig.18 Sabot pitch moment coefficient vs.time
图19、图20表明弹托分离对弹体气动参数在初始及中间分离阶段是有波动的，但随着弹托的分离，对弹体气动参数的变化影响减小直至稳定。随着弹托的分离，弹体的阻力系数与升力系数在近似相同时刻出现极值，例如在6.2Ma的条件下，当t=0.7 ms左右时，阻力系数与升力系数均达到极值，分别为0.4、1.8左右。结合图7(a)、图14以及图15，此时弹体与弹托之间存在激波反射，使得弹体后半部分至尾翼处表面压力增大，达到1.23 MPa，导致阻力系数与升力系数增大；之后t=1.7 ms左右，弹体与弹托之间干扰减小，使得弹体尾翼处表面压力下降至0.12 MPa，从而使得气动系数减小至稳定。
图19 弹体阻力系数曲线
Fig.19 Projectile body drag coefficient vs.time
图20 弹体升力系数曲线
Fig.20 Projectile body lift coefficient vs.time
4.4 不同马赫数对弹托分离姿态和轨迹的影响
图21、图22分别表示不同马赫数下弹托重心相对于分离前的轴向距离和径向距离随时间变化的趋势。由图可知，随着时间的增加，弹托加速做轴向运动远离弹体，而近似匀速做径向运动远离弹体。马赫数越高，弹托的轴向速度和径向速度越大，同一时刻的轴向距离和径向距离越大，越容易分离。
不同马赫数下弹托分离过程中俯仰角随时间变化的趋势如图23所示。由图可知，分离时弹托初速越大，在分离过程中同一时刻弹托的俯仰角越大。
图21 弹托重心的轴向距离曲线
Fig.21 The axial distance of sabot’s center of gravity varies with time
图22 弹托重心的径向距离曲线
Fig.22 The radial distance of the sabot’s center of gravity varies with time
图23 弹托俯仰角曲线
Fig.23 Pitch angle of sabot varies with time
图24、图25分别表示不同马赫数下弹托沿X轴方向的速度变化曲线以及绕Y轴方向的角速度变化曲线。由图可知，沿X轴方向的速度在某一时刻前加速增大，之后缓慢增大，绕Y轴方向的角速度先逐渐增大后逐渐减小。
例如在6.2Ma的条件下，由于t=1.7ms左右时弹托的内表面与来流垂直，在此之后升力的方向相反，使合力减小，合加速度减小，因而沿X轴向的速度开始缓慢增加，绕Y轴方向的角速度开始逐渐减小。同时马赫数越高，同一时刻弹托沿X轴方向的速度越大，绕Y轴方向转动的角速度越大，曲线峰值越高。
图24 弹托沿X方向的速度曲线
Fig.24 The velocity of sabot along the X-axis varies with time
图25 弹托绕Y轴方向的角速度曲线
Fig.25 The angular velocity of sabot around the Y-axis varies with time
5 结论
1) 对于弹托分离过程，常规静态网格无法模拟弹托运动引起的流场变化，对弹托建立了动态网格，采用弹簧光顺法传递节点力，通过6DOF动力学与流体力学耦合求解。当位移过大时，通过局部网格重构法对网格重新划分，计算新的弹托6DOF位置及流场特性，使节点在不同时刻不同位置达到力平衡，解决弹托动态运动及流场变化的耦合问题。
2) 初始阶段时弹托迎风槽处的涡流存在压差，压力为弹托分离提供翻转力矩，在压力达到一定阈值时，弹托翻转，气动参数波动；中间阶段时弹体与弹托之间的激波经过多次反射使得弹托表面压力增加，径向力使弹托分离，气动参数增大并逐渐达到峰值。同时弹体表面受到压力，引起其气动参数剧烈变化，对弹体飞行产生干扰。
3) 分离后期阶段时弹托轴向距离与径向距离增加，弹体与弹托之间的激波相互影响减小，弹托表面压力没有太大波动，两者耦合作用降低并做自由飞行。
[1] MA W M,LU J Y.Thinking and study of electromagnetic launch technology[J].IEEE Transactions on Plasma Science,2017,45(7):1071-1077.
[2] LI X,QIU Q X,LI P F.Research on the dynamic properties of electromagnetic railgun launch system[J].Journal of Physics:Conference Series,2019,1345(3):1-8.
[3] ZHU Y W,WANG Y,YAN Z M,et al.Multipole field electromagnetic launcher[J].IEEE Transactions on Magnetics,2010,46(7):2622-2627.
[4] FAIR H D.Advances in electromagnetic launch science and technology and its applications[J].IEEE Transactions on Magnetics,2009,45(1):225-230.
[5] CUI M Y,WANG X Z,BAI J,et al.Comparison and simulation analysis of augmented electromagnetic rail launchers[C],Proceedings of 2019 International Conference on Applied Mathematics,Modeling,Simulation and Optimization.Guilin,China:DEStech Publications,2019.
[6] LI X P,LU J Y,FENG J H.Simulation of sabot discard for electromagnetic launch integrated projectile[J].IEEE Transactions on Plasma Science,2018,46(7):2636-2641.
[7] 王小召,肖林楠,吉继昌,等.一种导弹电磁线圈弹射器方案设计与仿真分析[J].兵工自动化,2020,39(12):72-75，85.
Wang X Z,Xiao L N,Ji J C,et al.Design and simulation analysis of a missile electromagnetic coil ejector [J].Ordnance Industry Automation,2020,39(12): 72-75，85.
[8] 李湘平，鲁军勇，冯军红，等.电磁发射弹丸飞行弹道仿真[J].国防科技大学学报,2019,41(04):25-32.
LI X P,LU J Y,FENG J H,et al.Simulation of flight ballistic of electromagnetic launch projectile[J].Journal of National University of Defense Technology,2019,41(04):25-32.
[9] HAYDEN T E,ELDER D J.Hypervelocity projectile development for electromagnetic guns[J].IEEE Transactions on Magnetics,1991,27(1):452-457.
[10] ERENGIL M E,ZIELINSKI A E.Effect of electromagnetic launch on armature/sabot discard[J].IEEE Transactions on Magnetics,2001,37(1):67-72.
[11] HEAVEY K R,JAMES D,JUBARAJ S.Computational fluid dynamics flow field solutions for a kinetic energy(KE) projectile with sabot[J].Nuclear Instruments & Methods in Physics Research B,2003,268(7/8):1084-1087.
[12] 李湘平，鲁军勇，冯军红，等.采用动网格技术的弹托分离仿真模型[J].国防科技大学学报,2018,40(05):9-13.
LI X P,LU J Y,FENG J H,et al.Simulation model for sabot discard using dynamic mesh technique[J].Journal of National University of Defense Technology,2018,40(05):9-13.
[13] 罗倩，胡常莉.基于局部时均Navier-Stokes模型的分域计算方法[J].兵工学报,2020,1-9.
LUO Q,HU C L.A domain computing method based on partially-averaged Navier-Stokes turbulence model[J].Acta Armamentarii,2020,1-9.
[14] 蒋胜矩，陈小龙，张鹏飞.一种基于动态嵌套网格的膛内弹托分离力预测方法[J].弹箭与制导学报,2019,39(06):98-101.
JIANG S J,CHEN X L,CHEN P F.A prediction method of sabot separating force in barrel based on moving chimera grid[J].Journal of Projectiles Rockets Missiles and Guidance,2019,39(06):98-101.
[15] 罗驭川，黄振贵，杜宏宝，等.空心弹结构设计及弹托分离过程仿真[J].兵器装备工程学报,2020,41(05):5-8.
LUO Y C,HUANG Z G,DU H B,et al.Design of hollow projectile structure and simulation of separation process of sabot[J].Journal of Ordnance Equipment Engineering,2020,41(05):5-8.
[16] LI M J.Interaction between free surface flow and moving bodies with a dynamic mesh and interface geometric reconstruction approach[J].Computers & Mathematics with Applications,2021(81):649-663.
[17] ZENG D H,ETHIER C R.A semi-torsional spring analogy model for updating unstructured meshes in 3D moving domains[J].Finite Elements in Analysis and Design,2005,41(11/12):1118-1139.
[18] DEGAND C,CHARBEL F.A three-dimensional torsional spring analogy method for unstructured dynamic meshes[J].Computers & Structures,2002,80(3):305-316.
[19] MATSUSHIMA K,MURAYAMA M,NAKAHASHI K.Unstructured dynamic mesh for large movement and deformation[C].40th AIAA Aerospace Sciences Meeting & Exhibit.Reno,NV,U.S.A.:American Institute of Aeronautics and Astronautics,2002.
[20] ZHELEZNYAKOVA A L,SURZHIKOV S T.Calculation of a hypersonic flow over bodies of complex configuration on unstructured tetrahedral meshes using the AUSM scheme[J].High Temperature,2014,52(2):271-281.
[21] 姜毅，郝继光，傅德彬.导弹发射过程三维非定常数值模拟[J].兵工学报,2008,29(08):911-915.
JIANG Y,HAO J G,FU D B.3D unsteady numerical simulation of missile launching[J].Acta Armamentarii,2008,29(08):911-915.
[22] LIJEWSKI L E.Comparison of transonic store separation trajectory predictions using the pegasus/DXEAGLE and beggar codes[C].15th Applied Aerodynamics Conference.Atlanta,GA,U.S.A.:American Institute of Aeronautics and Astronautics,1997.
[23] 古刚，李宣.一体化弹丸弹体和弹托分离特性研究[J].舰船科学技术,2020,42(15):32-37.
GU G,LI X.Research on separation characteristic of intergration launch projectile’s body[J].Ship Science and Technology,2020,42(15):32-37.
[24] SNYDER D,EVANGELOS K,JOHN A.Transonic store separation using unstructured CFD with dynamic meshing[C]//33rd AIAA Fluid Dynamics Conference and Exhibit.Orlando,Florida:American Institute of Aeronautics and Astronautics,2003.
[25] PANAGIOTOPOULOS E E,KYPARISSIS S D.CFD transonic store separation trajectory predictions with comparison to wind tunnel investigations[J].International Journal of Engineering,2010,3(6):16.
[26] POSTSDAM M,GURUSWAMY G.A parallel multiblock mesh movement scheme for complex aeroelastic applications[C].39th Aerospace Sciences Meeting and Exhibit.Reno,NV,U.S.A.:American Institute of Aeronautics and Astronautics,2001.