电磁发射是一种利用电能替代传统化学能完成发射的全新发射方式^[1]，电磁轨道发射技术在军用和民用领域都有着巨大的潜在优势和广阔的应用前景。在信息化、电能化快速发展的大趋势下，电磁轨道炮有着常规火炮无法比拟的优势，全电化更有利于发射过程的精准控制，超高的发射初始速度令一众常规兵器无法望其项背，更高的能量转化效率以及更大的毁伤动能使得电磁轨道炮成为了未来战场最具潜力的新概念武器之一。早在1901年，挪威物理学家伯克兰就正式提出了电磁发射概念。美国在电磁轨道发射技术方面走在世界的前列。1920年，美国费城的电炮公司研制了用于火炮的电磁加速器^[2]。随着对于电磁炮特性研究的不断深入，电磁炮从特性分析，到样机实验再到实践应用^[3-7]，美国最新公布的实验数据显示，其电磁炮的性能指标已经达到较高水平^[8]。

对于电磁炮物理特性的仿真，温艳玲^[9]基于电流和磁场的多物理场耦合有限元分析得到磁场和电流的分布，并利用电磁场仿真结果实现电流趋肤效应下轨道电阻梯度计算。郭仁荃^[10]采用COMSOL仿真软件建立轨道炮分析发射过程中电枢前端磁场分布的基础上，对不同金属屏蔽方式及屏蔽层厚度比例对屏蔽性能的影响进行了仿真。金亮^[11]给出了基于Kerrisk轨道电感梯度的解析计算方法、电流扩散和速度趋肤效应影响下的电枢控制方程。关晓存^[12]运用通用有限元分析软件ANSYS的耦合计算流程，对单级感应线圈炮中电枢电磁场和温度场进行了仿真。王昊^[13]通过有限元分析软件COMSOL分别建立了单次发射条件下电磁轨道炮三维模型及连续发射模式下电磁轨道炮二维模型。楼宇涛^[14]结合场路模型计算了电磁轨道炮系统的发射效率和涡流能耗，讨论了不同管身结构和材料下的发射效率，进一步分析了管身对电磁轨道炮力学特性的影响。马硕^[15]用COMSOL三维模型计算输入电流大小、电流波形、轨道、电枢几何尺寸和物性参数对轨道中电流密度和温度的影响。关永超^[16]建立了包括脉冲电源模块和电磁轨道炮的全电路模型，实现了从电源放电至电磁轨道炮发射过程的全电路模拟。葛一凡^[17]基于COMSOL有限元分析软件建立了轨道炮二维瞬态有限元模型，考虑电枢高速运动产生的速度趋肤效应的影响，分析了不同驱动电流上升沿时间作用下的电磁热紧耦合轨道的电流与温度分布规律。

综上所述，相关学者基于COMSOL对于电磁轨道炮的物理特性进行了大量的研究，但这些研究在进行网格划分时采用“静网格”。COMSOL通过手动划分网格功能，可以对重点研究的区域的网格进行自主划分。在进行网格划分时，可以对电枢部分进行更细致的网格划分以提高计算准确性。但在求解枢轨电磁特性时，尤其是在时变电流的作用下，不同时刻电枢位所处的位置与枢轨接触面的物理环境会明显不同，如果按照原来的网格进行计算，会对结果造成很大误差，但如果在每一个时间步上对电枢以及轨道都进行动网格的重新划分，会大大增加计算复杂度。基于此，本文提出一种新的网格划分方式，对电枢及枢轨接触部分，进行细网格划分，并采用滑移网格，建立滑移网格移动域，使整个细网格区域进行移动，电枢区域的滑移网格不涉及网格的变形与重构，而轨道部分采用动网格功能，并进行粗网格的划分。数值结果表明，采用动网格划分功能求解瞬时、快速的模型的电磁特性时，其计算性能明显优于传统的静网格划分。

1.   理论分析

电磁轨道炮最基本的工作原理是利用脉冲大电流产生电磁力驱动弹丸，将弹丸加速至超高速发射。图1为电磁轨道炮工作原理示意图，两根导体轨道一根接入电流，一根接地，与接触的电枢形成回路。回路电流流过金属轨道在其周围形成强大的磁场，方向垂直于两轨道所在平面向上。搭载弹丸的载流电枢在轨道强磁场的作用下受到洛伦兹力的推动不断加速，从而使弹丸在极短的时间内获得较高发射速度。

图  1  轨道炮工作原理

Figure  1.  Working principle of railgun

下载: 全尺寸图片 幻灯片

在轨道以及电枢区域，满足电流守恒条件

式中：${\boldsymbol{J}}$为电流密度，${Q_{j,v}}$为总电荷量，$\sigma$为导体的电导率，${\boldsymbol{E}}$为电场强度，${\boldsymbol{D}}$为电位移矢量，V为两根轨道的电势差。得到电流密度方程后，根据Maxwell方程求解电流产生的感应磁场

通过安培定律可得到电枢的受力方程

在电磁力F的作用下，电枢获得一定的加速度a，其动力学方程满足牛顿第二定律$F = ma$。电枢在轨道运动的过程中，受到轨道摩擦力的作用，同时，强大的电流会产生较高的热量从而对轨道的几何形状、枢轨之间的压力以及电阻都会产生很大的影响。

2.   基于COMSOL的电磁轨道炮模型的建立

2.1   模型参数的选取

在基于COMSOL进行轨道炮模型的构建时，首先要明确电枢和轨道涉及到的方程以及边界条件。建立的轨道炮模型参数如表1所示。

表  1  电磁轨道炮模型参数

Table  1.  Model parameters of electromagnetic railgun

track length/mm	track width/mm	track thickness/mm	center to center spacing/mm	orbital conductivity/(S·m−1)	armature length/mm	armature width/mm	armature conductivity/(S·m−1)
900	40	20	50	$5.998 \times {10^7}$	50	30	$3.774 \times {10^7}$

下载: 导出CSV 

| 显示表格

构建的电磁轨道炮模型总共划分为56个区域，在COMSOL中进行物理场设置时，需要对物理场的求解环境与求解区域进行设置，并且对求解区域所满足的物理场方程以及边界条件、初始值进行设定。在电枢与轨道接触面之间，由于电枢和轨道的不完全接触会导致电流在传输过程中分布状态产生很大变化，进而对求解的瞬态磁场值产生很大影响。所以需要设置接触面以及材料的一些信息。表2给出了在进行轨道炮特性仿真时用到的一些参数。

表  2  电枢模型参数取值

Table  2.  Parameter values of armature model

average height of the rough surface ${\sigma _{{\mathrm{asp}}}}$/μm	average slope of the rough surface ${m_{{\mathrm{asp}}}}$/μm	microhardness of solids ${H_{\mathrm{c}}}$/Pa	coefficient of friction $\;{ \mu _{\mathrm{f}}}$	viscous coefficient of friction ${C_{\mathrm{f}}}$	rectangular pulse current intensity ${I_0}$/MA
1	0.4	$3 \times {10^9}$	0.11	0.03	0.7

下载: 导出CSV 

| 显示表格

2.2   全局常微分方程以及动网格设置

针对电磁轨道发射系统的理论方程，调用软件中的全局常微分方程，对动网格进行编辑处理。当电枢位置发生改变后，需要重新嵌套描述系统特性的数学模型，然后再次对系统进行计算。因此，在计算电磁发射的过程中，调用全局常微分方程，即ODE方程，对电枢及其与轨道接触的区域设置移动网格。使用全局常微分和微分代数方程接口和逻辑表达式来指定用于描述速度和位置的微分方程。根据牛顿第二定律，可得电磁力

式中：a为电枢运动的加速度；m为电枢的质量。COMSOL全局方程需要定义电枢的速度v，根据能量守恒原理，电磁发射系统的电磁能增量为

电磁力F的表达式为

式中：x为电枢位移量，$L(x)$为电感量；$I(t)$为电流，$L'$为电感梯度。式（6）表明，弹丸受力与电流的平方和电感梯度有关，不随弹丸位置的变化而变化。因此，在COMSOL中进行全局方程定义

式中：t为脉冲电流的作用时间。

图2为滑移网格和动网格划分方式图示。COMSOL的移动网格接口可用于所有或部分几何相对于绝对参考坐标系移动的模型，中间电枢及其接触的轨道部分设置为滑移网格，其网格大小不会改变，不需要考虑电枢的自身变形，在计算过程中，只需要将网格部分整体平移，不需要进行网格的重新构建。在轨道与电枢不接触的部分，进行变形网格的划分，COMSOL的变形网格接口可用于几何形状发生变形的模型，使材料区域跟随几何变形刷新。在实际的计算中，电枢及其与轨道接触部分的网格随着计算的进行向右移动，而左侧的网格被拉伸，右侧的网格被压缩。COMSOL中，参考坐标系称为空间坐标系，相对材料静止的坐标系称为材料坐标系，默认情况下，变形网格接口使空间坐标系中的网格按照坐标变换所定义的进行变形

图  2  滑移网格+动网格图示

Figure  2.  Diagram of sliding mesh and dynamic mesh

下载: 全尺寸图片 幻灯片

式中：变量$x,y,{\textit{z}}$是空间坐标系中网格点的坐标。$X,Y,Z$是材料坐标系网格点的坐标，$f,g,h$的物理意义是空间坐标系中某点与材料坐标系中某点形成映射，在计算时，空间中该点的属性即为材料坐标中对应点的属性，从而实现材料在空间中的“移动或变形”。动网格第n+1步的体积为

即下一步的运动情况是由当前时间步的计算结果确定的，各个时间步的体网格更新是基于边界条件新的位置由COMSOL自动来完成的。而滑移网格模型的网格是刚性的，第$n + 1$步的体积${V_{n + 1}} = {V_n}$。在数据传递方面，滑移网格模型是通过滑移面（交界面）传递区域计算电枢速度以及电流密度结果数据，而动网格模型直接通过变形网格传递数据。由于该计算模型所有的尺寸比较常规，且不存在极细结构，所以本次计算将滑移网格区域设置为三角形与四边形网格结合，并设置为超细化网格。变形网格部分设置为长方体粗网格。相对于传统“静网格”而言，一般仅在有运动物体参与的仿真模型中使用。引入动网格的概念之后，仿真流程就不再是单独的一条流水线，而变成了循环迭代模型。动网格能够对模型进行实时的网格划分，以满足对运动的模型的物理环境的更新。

图3-图5分别给出了三种不同的静态网格划分示意图。表3为不同网格划分方式的具体数据，从表3可以看出，采用静网格方式对电磁炮模型进行剖分，采用较粗的网格，计算的域单元个数较少，但这种网格划分方式的区域中包含比指定的最小单元更小的区域，即某些计算区域小于最小网格尺寸，会导致其计算结果不准确。这些比网格最小尺寸还要小的计算区域主要来自于电枢。若想要得到超细化划分状态下的精确解，则需要计算的网格数量很大，其计算区域的单元个数约为粗网格的16.3倍。采用轨道为粗网格划分，电枢为细网格的划分方式进行网格划分，其计算区域的个数缩减为全部细网格划分的1/3。这种划分方式解决了电枢处计算区域小于最小网格划分区域的问题，但对于高速运动的电枢而言，如果在枢轨接触区域仍采用粗网格划分，在接触面的电流传递的计算中会带来较大误差。而滑移网格+动网格的划分方式，将电枢及其接触到的轨道部分看成一个整体，在这个区域内的网格采用极细化的划分方式，而其余轨道部分的网格划分由自由空间的坐标到材料空间的坐标映射来完成对网格的重新划分，这种划分方式不但能大大减少对于细化网格划分的区域，而且能保证枢轨接触区域附近的网格较小，使得计算结果更精确。

表  3  三种不同静网格划分方式所划分的网格信息

Table  3.  Mesh information divided by three different static meshing methods

mesh information	domain units number	boundary elements number	edge elements number	maximum mesh size/mm	minimum mesh size/mm
mesh①	57500	17172	1973	90	16.2
mesh②	935392	79766	4279	18	0.18
mesh③	326996	55474	5206	90	0.18
mesh④	372512	53554	4254	90	0.18

下载: 导出CSV 

| 显示表格

图  3  全部为细网格划分

Figure  3.  All divided into fine meshs

下载: 全尺寸图片 幻灯片

图  4  全部为粗网格划分

Figure  4.  All are coarse mesh divisions

下载: 全尺寸图片 幻灯片

图  5  轨道为粗规网格划分，电枢为细网格划分

Figure  5.  Coarse track mesh and fine armature mesh

下载: 全尺寸图片 幻灯片

2.3   激励电流的设置

由公式（6）可知，激励电流的选取对电枢电磁力的求解有很大关联，端子处外加的驱动电流的波形、峰值以及上升时间对电枢运动也会产生很大的影响。COMSOL软件可以进行电流函数的自定义，通过添加函数类型可对外加激励电流进行设置，通过插值，可以设定激励电流大小随时间的变化。通过设置矩形脉冲电流和高斯脉冲电流两种不同的激励电流的形式，求解在不同电流驱动下电枢以及轨道的电流密度分布与电枢运动速度。比较不同激励电流下不同网格划分方式的计算结果。

在实际的电磁轨道炮应用中，激励电流能达到kA～MA量级，设置矩形脉冲电流的峰值大小为0.7 MA，脉宽为2 ms，通过插值的方法设置脉冲电流形式如图6所示，脉宽为2 ms，峰值为0.7 MA。

图  6  脉冲激励电流分布示意图

Figure  6.  Schematic diagram of pulse excitation current distribution

下载: 全尺寸图片 幻灯片

3.   结果与分析

图7为矩形脉冲电流激励下不同网格划分时电枢运动速度。为表述方便，记全部采用粗网格方式为网格①，全部采用细网格方式为网格②，轨道采用粗网格、电枢采用细网格划分方式为网格③，电枢与轨道接触部分采用细网格，网格①、②、③全部为“静网格”，其余部分采用粗网格并采用进行粗网格重构的划分方式为网格④。当外加激励电流源为矩形脉冲电流时，电枢受到的合外力近似为恒力，电枢运动的速度与时间的二次方成正比。在采用网格①进行划分时，电枢的速度与其余网格划分方式求解的电枢速度差别较大，网格②和网格④的划分方式求解的电枢速度结果符合得较好。

图  7  不同网格划分方式在矩形脉冲电流作用下电枢速度随时间变化图

Figure  7.  Armature velocity changes with time under the action of rectangular pulse current

下载: 全尺寸图片 幻灯片

图8为高斯脉冲电流激励下不同网格划分时电枢运动速度。当外加激励电流源为时变电流时，轨道以及电枢的电流密度分布发生改变进而导致电枢受力产生变化，网格①计算电枢速度的结果与网格②划分时相差较大，在电枢运动到轨道末端时，二者产生的速度差值为27.48 ${\text{m/s}}$。从计算结果可以看出，网格③在电流为时变电流时计算结果仍会产生较大误差，而网格②与网格④的计算结果在电流变为时变电流时仍然符合得较好。所以网格④的划分方式，不但能大大减少细网格划分求解时求解域的数量，其准确率仍能保持一个较高水平。

图  8  不同网格划分方式在高斯脉冲电流作用下电枢速度随时间变化图

Figure  8.  Armature velocity vs time under Gaussian pulse current

下载: 全尺寸图片 幻灯片

图9为脉冲电流作用下，四种不同网格划分方式求解的电枢在不同位置时中心点处电流密度分布结果。由图可知，中心位置处的电流密度随着电枢的运动呈现出指数型的增加，在运动到${\textit{z}} = 200{\text{ mm}}$时达到最高点，在外加激励电源与电流传输共同的影响下，中心处的电流密度呈现出先上升后下降的趋势，在运动到${\textit{z}} = 600{\text{ mm}}$时，到达轨道末端，电流密度骤减。网格①划分下求解的电流密度与网格②划分时求解的电流密度差距较大，网格④求解的电流密度的结果与网格②差距较小。但计算时间约为网格②求解时间的1/6。本模型利用COMSOL仿真软件，在搭载Intel(R)Xeon(R)Silver 4214R CPU@ 2.40 GHz 2.39 GHz的双核处理器上进行求解，四种不同的网格划分方式与不同的激励电流作用下求解的时间如表4所示。

图  9  不同网格划分方式在电枢中心位置电流密度变化图

Figure  9.  Diagram of current density at the center of the armature with different meshing methods

下载: 全尺寸图片 幻灯片

表  4  不同网格划分方式与激励电流作用下计算时间

Table  4.  Calculation time under different grid partitioning methods and excitation currents

	calculation time/s
	rectangular pulse current	Gaussian pulse current
mesh①	2491	3907
mesh②	59277	78209
mesh③	8602	12611
mesh④	11113	15507

下载: 导出CSV 

| 显示表格

以上分析表明，采用滑移网格+动网格的划分方式可以十分有效的减少计算时间，且相比于轨道粗网格、电枢细网格划分对于电枢速度的求解更精确。随着电枢速度的增加，网格①与网格③求解的结果与网格②求解的结果差距越来越大。这是由于在枢轨接触区域，由于不完全接触所带来的电流损失在每一次求解过程中产生误差的累计。而网格①与网格③的划分方式部分最小网格尺寸仍大于模型几何尺寸，这对求解的结果会产生较大误差。

4.   结　论

本文基于COMSOL动网格功能，提出一种新的网格划分方式—滑移网格+动网格划分方式，对不同激励电流作用下电磁轨道炮的电枢速度以及电磁炮内部电流密度进行了仿真分析，所提出的网格划分方式能在保持较高准确性的情况下，更快速的对的电磁轨道炮的物理特性进行求解。以上分析表明，采用动网格划分功能求解瞬时、快速的模型的电磁特性时，其计算性能明显优于传统的静网格划分。