周东强，严 艺，杨 易，高振良，刘奕鑫

（北京空间飞行器总体设计部，北京100094）

0 引言

电磁发射是一种全新概念的发射方式。电磁发射技术将电磁能变换为发射物体需要的瞬态动能，可在短距离内将发射物体加速至高速，能有效突破传统发射的速度和能量极限，在军事和民用领域具有巨大的优势和广阔的应用前景[1-2]。电磁发射过程中对发射装置加载脉冲大电流，伴随着磁场扩散、热传导、结构变形等物理过程，是具有瞬态特性的多物理场耦合问题[3]。

电磁发射的数值仿真是发射装置设计、试验的重要支撑手段。国内外对电磁发射的数值仿真与计算都非常重视，开展了大量研究工作，使用多种数值计算方法模拟电磁发射的多物理场相互作用。Powell 和Zielinski[4]提出了计算电磁场、电流密度、导轨和电枢温度等随空间、时间分布的数学模型。Ghassemi 和Pasandeh[5]考虑导轨‒电枢之间的大电流产生的热能对导轨‒电枢组件的电学、热学、力学性能影响，在麦克斯韦方程组中引入能量方程，同时考虑导轨和电枢中的电传导、欧姆热和摩擦热，得到非线性控制微分方程组，计算导轨和电枢中的热和电磁感应分布。随着有关数值仿真应用需求的增长，研究人员逐渐开发出用途各异的软件，如美国的EMAP3D，法国和德国的EMAS，英国的MEGA 等[6-11]，一些大型的商业软件也可以进行电磁轨道发射装置仿真计算。然而，国外部分自用仿真软件不对外公开，且其他大型商业软件针对轨道电磁发射在电枢高速运动多物理场耦合计算方面存在局限性，在处理高速滑动电接触问题上能力不足，无法计算特有的速度趋肤效应等[3,11-15]，因此可以说，对固体电枢电磁发射过程的建模和仿真尚处在不断探索与研究之中[16-20]。

本文基于有限元法的三维电磁、热、力瞬态耦合算法，根据时间步进顺序耦合方法构建电磁、热、力多物理场耦合集成框架，用以指导相关计算程序的编制，实现对电磁发射中电枢高速运动状态下电磁、热、力瞬态耦合计算以及电枢高速滑动等多物理场瞬态状态的定量分析。

1 瞬态电磁分析与有限元模型

轨道电磁发射是瞬态电磁相互作用的过程，伴随着导轨与电枢间的电流传导、电磁感应以及趋肤效应、涡流等现象。麦克斯韦方程是表征所有宏观电磁场现象的基本方程[3,5,18]，电磁分析的实质是求解给定边界条件下的麦克斯韦方程组问题。

电磁发射的三维瞬态电磁计算基于麦克斯韦理论方程和体现发射装置材料特性的本构关系，采用伽辽金方法[21]对计算模型进行离散，使用三维高斯积分方法进行数值积分，通过计算机编程实现瞬态电磁、电位的计算。电磁和电位瞬态分布的仿真计算是电磁发射中多物理场耦合分析的重要环节。

1.1 麦克斯韦方程微分形式与电磁分析本构关系

三维瞬态电磁场分析是以麦克斯韦方程组为基础，根据设定的边界条件求解场域内电场、磁场、电流密度、涡流等状态量随时间和空间变化的过程[6]。麦克斯韦方程组由4个基本定律组成，分别为安培环路定律、法拉第电磁感应定律、高斯电通定律和高斯磁通定律。为了方便使用有限元数值算法，选用其微分形式建立三维瞬态电磁分析的数学模型：

式(1)～式(4)中：H为磁场强度；J为电流密度；D为电位移；E为电场强度；B为磁感应强度；ρ为电荷密度。

不同材质的发射装置对发射中多物理量的作用结果不同，为了反映电磁发射装置的材质构成对电场、磁场、感应电流的相互作用的影响，在电磁发射的数值计算中使用反映材质宏观性质的本构关系来计算介电常数、磁导率、电导率对电磁场的影响和作用[6]：

式(5)～式(7)中：ε为介电常数、μ为磁导率、σ为电导率，各向同性材料用标量表示，各向异性材料用张量表示；v为电枢发射速度。

为了求解麦克斯韦方程，通过定义标量电势Φ与矢量磁势A，将包含电场与磁场2个场量的一阶微分方程转化为2个独立的包含1个场量的二阶微分方程[5]。将电位标量势E= -∇Φ与本构关系D=εE代入高斯电通定律（式(3)）中推导可得电位标量势方程

磁位矢量势满足高斯磁通定律，由本构关系B=μH知μH= ∇×A，代入安培环路定律（式(1)）中推导可得磁位矢量势方程

1.2 电位标量与磁位矢量有限元列式

为了使用计算机编程实现电磁发射中电磁、电位的数值积分计算，本文使用伽辽金方法建立电位标量势方程的三维有限元格式，构造8节点六面体等参单元，其形函数N为三维空间参考坐标(ξ,η,ζ)的函数，

2 热传导分析与有限元模型

电磁发射过程中产生热的途径主要有电枢与导轨摩擦生热、脉冲激励电流生热及感生涡流生热等；热扩散途径主要为热传导和热对流。本文主要考虑瞬态热传导问题的基本方程与有限元方法在热传导中的应用。

2.1 热传导与热对流方程

2.2 瞬态热传导有限元列式

与电磁计算相似，将计算域离散成有限个8节点六面体单元，单元的温度场可以用节点温度Ti插值得到，即

在瞬态热传导计算中此温度场Ti是与时间t有关的函数，Ni为8节点六面体形函数。使用伽辽金积分法建立三维瞬态热传导的等效积分形式[5]

将节点温度插值代入热传导方程（式(15)），得到有限元方程

由此，热传导方程可以转换为以时间t为独立变量的线性常微分方程组。式(19)中：C为热容矩阵；K为热传导矩阵；P为温度载荷列阵；T˙为节点温度对时间偏倒列阵。

3 瞬态多场耦合计算程序

计算程序对时间单元进行离散，以电磁、热、力收敛为准则进行编制与开发，构建瞬态多场耦合计算程序框架如图1所示。计算数据信息主要包括控制信息与有限元前处理信息：控制信息主要为物理计算时间与计算时间步长以及输出控制；有限元前处理信息主要包含节点信息、单元拓扑关系、材料信息和边界条件等。

图1 瞬态多场耦合计算程序框架Fig.1 Framework of the coupled transient multi-field calculation program

计算程序中的各子模块功能分别为：

1）电场分析子模块：电流曲线插值，组装电场刚度矩阵；

2）磁场分析子模块：组装磁场刚度矩阵、磁场远场边界；

3）热传导分析子模块：组装热传导刚度矩阵、主动冷却与被动冷却边界条件；

4）动态响应子模块：结构动态响应计算，电枢速度与位移计算；

5）电接触子模块：电枢与导轨接触面准匹配网格搜索，瞬变状态量插值。

轨道电磁炮发射激励源是由外在电流加载到2根导轨端部形成回路，电枢与导轨接触是电场分布计算的先决条件；当电枢发射脱离导轨后回路中断，外在激励源不再维持回路电磁感应，计算时仅需考虑热场传导和结构动态效应，不再考虑电磁感应，因此电枢出膛后不再调用电场和磁场分析子模块。

使用FORTRAN90搭建程序主体框架，各子模块与关键部分均采用子程序模式，主体框架与子程序间为调用与被调用关系，主体框架代码如下：

4 典型轨道电磁发射计算与分析

基于已完成的多物理场耦合有限元计算程序，以典型轨道电磁炮发射过程为例，分析计算电流密度、电磁场、热场以及发射出口速度等多物理场的状态量分布情况。发射装置主要由导轨、电枢、聚碳酸酯绝缘体、约束槽等组成，为了同时保证计算精度、速度与效率，对计算模型进行合理简化。

4.1 轨道电磁发射模型

轨道电磁发射仿真模型包含有限电磁传播域网格、导轨网格、电枢网格和电流激励片。如图2所示，电磁传播域网格单元、导轨网格单元、电枢网格单元为8 节点六面体体单元，电流激励片网格单元为四边形面单元。金属导轨尺寸为15 mm×20 mm×500 mm，磁场域尺寸为200 mm×200 mm×1000 mm。

图2 轨道电磁发射计算模型Fig.2 Model for calculating orbit electromagnetic emission

根据电磁发射装置设计通常使用的材质情况，在数值计算中对各计算区域进行材料参数赋值。如表1所示：导轨采用铜制材料；电枢采用铝制材料；导轨与电枢以外电磁域媒质选用真空媒质材料参数。

表1 计算所用的相关材料参数Table 1 Material parameters for the electromagnetic calculation

4.2 载荷与边界条件

根据轨道电磁炮电磁、热、结构计算的载荷与边界情况，定义激励载荷1个（脉冲电流荷载），边界条件3个（无限电磁域、热对流边界、导轨位移约束）。

电流加载方式采用瞬态脉冲式电流，在模拟电磁发射过程时选择历时3 ms的单脉冲式电流加载。选择导轨与电枢形成的导体回路的2个端头分别为电流输入端和电流输出端（参见图3）；在输入端与输出端建立加载单元，通过导轨输入端头与输出端头加载单元将激励电流分配到导轨加载单元节点上（参见图4），使电流均匀加载到激励单元，按照激励单元各节点电流贡献值分配电流承载值。

图3 轨道电流加载方向示意Fig.3 Schematic diagram of current loading direction

图4 电流加载单元与加载节点示意Fig.4 The current loading unit and theengaging nodes

4.3 多场耦合计算结果与分析

4.3.1 脉冲电流对电枢速度的影响

为了比对脉冲电流峰值对电枢发射速度、加速度的影响规律，分别对电流脉冲峰值为3.0 kA、3.2 kA、3.4 kA、3.6 kA 共4种激励工况进行计算，脉冲加载历程皆为3 ms，峰值状态均持续0.5 ms。

图5～图7分别为电枢的发射速度曲线、发射位移曲线和发射加速度曲线，各图中曲线C1～C4分别代表电流脉冲峰值为3.0 kA、3.2 kA、3.4 kA、3.6 kA 的4种激励工况。可以看出：从开始加载电流到0.6 ms 时刻期间，电枢发射速度明显提升，0.6 ms之后速度增长趋缓；曲线C4与曲线C1相比，出口速度提高68.61 m/s，电枢脱离导轨时间缩短0.6 ms；电枢加速度变化与加载电流变化高度相关，电枢脱离导轨时电枢与导轨间接触有效电阻发生变化，加速度会出现微小波动。

图5 计算出的电枢发射速度曲线Fig.5 Curve of calculated launching speed of the armature

图6 计算出的电枢位移曲线Fig.6 Curve of calculated launching distance of the armature

图7 电枢加速度计算曲线Fig.7 Curve of calculated acceleration of the armature

4.3.2 导轨与电枢中电流密度分布

分析电磁发射全过程导轨和电枢中的电流密度分布情况：电枢在导轨中滑行时，导轨的最高电流密度约为平均电流密度的1.39～1.42倍，电枢的最高电流密度约为平均电流密度的1.53～1.57倍；当电枢即将脱离导轨时，电枢与导轨的接触面较小，电流集中明显增强。

图8显示了0.4 ms时刻电流密度集中区域，其中电流密度最大区域在电枢尾部与导轨接触位置，体现了电枢高速滑动尾部趋肤效应。为了观察电流集中区内电流密度的时间历程曲线，并且能够与非电流集中区内电流密度的时间历程曲线对比，选取7个检测点作为电流密度时间历程取值点（见图9）：node1为导轨网格单元附属节点，处在导轨外载电流加载端头内侧；node2、node3、node5、node6位于电枢与导轨接触面；node4与node7位于电枢尾部与头部对称位置。

图8 计算出的0.4 ms时刻电流密度分布云图Fig.8 Cloud map of calculated current density after 0.4 ms of current loading

图9 电流密度时间历程曲线取值点Fig.9 The location pointsfor investigating the timeprofile of current density

图10为导轨各取值点电流密度的时间历程曲线，可以看到，node2＞node3＞node1＞node5＞node6。表明：导轨上电枢尾角与导轨接触位置电流密度最大；在电枢与导轨接触面上沿着电枢发射方向电流密度逐渐变小，电枢头部电流密度小于外载电流加载端电流密度。

图10 导轨测试点电流密度时间历程曲线Fig.10 Time history of current density at typical positions of the rails

图11为电枢各取值点电流密度的时间历程曲线，可以看到，node2＞node3＞node4＞node5＞node7＞node6。表明：电枢尾角处电流密度最大；电枢与导轨接触面电流密度大于电枢体非接触面的电流密度；电枢头部电流密度小于电枢尾部电流密度。

图11 电枢测试点电流密度时间历程曲线Fig.11 Time history of current density at typical positions of the armature

4.3.3 空间计算域电磁场分布

图12为导轨与电枢的电磁力计算结果矢量分布，可以看到，2根导轨主要承受y方向排斥力，并且电磁力主要集中在电流回路，导轨内侧电磁力大于外侧电磁力。截取磁场域z方向截面磁势分布云图（图13）可以看到，两导轨之间磁势大约为导轨两侧磁势的2倍，电枢x方向尾部磁势大约为头部磁势的2倍。磁势最大位置在电枢尾部-x方向，最大磁势区域与电枢所在区域不存在交集，鉴于此可以改善电枢结构，优化最大磁势分布区域。

图12 导轨与电枢电磁力计算结果矢量分布Fig.12 Distributions of calculated electromagnetic force vectorsof the railsand the armature

图13 z 方向截面磁势分布云图Fig.13 Cloud diagram of magnetic potential on cross section

4.3.4 导轨与电枢中热场分布

导轨与电枢是发热最严重的部件，其热源主要为外在电流生热，热源体积密度为q=J2/σ。电磁发射过程中电流生热会在导体中累积，且由于电流传递不均匀，导轨与电枢中的热分布梯度非常明显。图14 为不同时刻导轨与电枢的热密度分布云图。从图中可以看到：0.3 ms时刻热密度分布与电流密度分布情况相近，热密度较高区域与电流密度集中区域相似，在导轨与电枢接触面上热密度较大，电枢尾部热密度比电枢头部热密度大，最大值为1.5×1012J/m3；随后逐步出现热密度累积现象，在0.6 ms时刻电流密度最大的位置为与导轨接触面的电枢尾部，最大热密度为2.4×1012J/m3。图14(c)与图14(d)分别代表电枢离开初始位置一段距离后的热密度分布情况，最大值分别为3.0×1012J/m3和2.8×1012J/m3。可以看出，热密度集中区域仍然为电枢在初始位置尾部的导轨处，2.0 ms时刻导轨电流加载端头整个区域都成为热密度集中区。也就是说，在重复发射情况下，电枢初始位置将会成为热烧蚀集中发生部位。

图14 不同时刻的热密度分布云图Fig.14 Cloud diagram of heat density at different moments

5 结束语

本文基于时间步进耦合算法和瞬态电磁、热传导的有限元算法，研究了适用于电磁发射的多物理场耦合计算程序设计方法，对电磁轨道炮进行三维瞬态模拟计算，得到磁感应密度、电流密度、热密度的时域和空间域分布。计算结果表明：时间步进耦合算法可以较方便地实现电磁发射中多物理场瞬态耦合仿真计算；电磁发射过程中电流密度分布呈动态扩散特点，在电枢与导轨接触面存在明显的速度趋肤效应，最大电流密度分布在电枢尾部与导轨接触位置；导轨中电枢发射的初始位置存在热密度累积现象，该位置也是电磁发射短时间内连续使用较容易发生热烧蚀和热损伤之处。

以上研究有助于对电磁发射中的多物理场耦合演化过程的了解，可为发射装置设计提供辅助支撑。但计算程序中暂未考虑接触摩擦磨损等因素，须对计算模型进一步优化和完善。