《万方数据-数字化期刊群》全文上网期刊
CNKI《中国学术期刊(网络版)》全文收录期刊
《中文科技期刊数据库》(维普网)全文收录期刊
超星期刊域出版平台、博看网全文收录期刊
日本科学技术振兴机构数据库收录

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

面向电磁发射的多场耦合仿真分析方法

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

周东强, 严艺, 杨易, 等. 面向电磁发射的多场耦合仿真分析方法[J]. 航天器环境工程, 2020, 37(6): 561-569 doi:  10.12126/see.2020.06.005
引用本文: 周东强, 严艺, 杨易, 等. 面向电磁发射的多场耦合仿真分析方法[J]. 航天器环境工程, 2020, 37(6): 561-569 doi:  10.12126/see.2020.06.005
ZHOU D Q, YAN Y, YANG Y, et al. Simulation analysis of multi-field coupling field during electromagnetic emission[J]. Spacecraft Environment Engineering, 2020, 37(6): 561-569 doi:  10.12126/see.2020.06.005
Citation: ZHOU D Q, YAN Y, YANG Y, et al. Simulation analysis of multi-field coupling field during electromagnetic emission[J]. Spacecraft Environment Engineering, 2020, 37(6): 561-569 doi:  10.12126/see.2020.06.005

面向电磁发射的多场耦合仿真分析方法

doi: 10.12126/see.2020.06.005
详细信息
    作者简介:

    周东强(1985—),男,硕士学位,高级工程师,从事航天器总体设计。E-mail: zhoudq2004@163.com

  • 中图分类号: V417+.7

Simulation analysis of multi-field coupling field during electromagnetic emission

图(14) / 表ll (1)
计量
  • PDF下载量:  23
  • 文章访问数:  267
  • HTML全文浏览量:  32
文章相关
  • 中图分类号:  V417+.7
  • 收稿日期:  2020-04-09
  • 修回日期:  2020-11-17
  • 网络出版日期:  2020-11-23
  • 刊出日期:  2020-12-25

面向电磁发射的多场耦合仿真分析方法

doi: 10.12126/see.2020.06.005
    作者简介:

    周东强(1985—),男,硕士学位,高级工程师,从事航天器总体设计。E-mail: zhoudq2004@163.com

  • 中图分类号: V417+.7

摘要: 针对电磁发射中高速运动的电枢与导轨间多物理场耦合计算及时间状态量连续传递问题,采用时间步进耦合算法构建多物理场耦合仿真框架,基于瞬态电磁分析、热传导与有限元方法开发出适用于轨道电磁发射的电磁、热、力多物理场瞬态耦合数值计算程序。并以典型轨道电磁炮为例,建立电磁、热、力三维有限元数值模型,采用单脉冲大电流为激励加载源,动态仿真分析轨道电磁发射中电枢与导轨电流密度、磁感应强度、焦耳热、电磁力等时间状态量的分布。计算结果直观反映出速度趋肤效应等物理现象,同时可方便获取电磁发射中多物理场随空间域和时间域变化的瞬态数据,为电磁发射装置的设计、优化及评估提供仿真计算支撑。

English Abstract

周东强, 严艺, 杨易, 等. 面向电磁发射的多场耦合仿真分析方法[J]. 航天器环境工程, 2020, 37(6): 561-569 doi:  10.12126/see.2020.06.005
引用本文: 周东强, 严艺, 杨易, 等. 面向电磁发射的多场耦合仿真分析方法[J]. 航天器环境工程, 2020, 37(6): 561-569 doi:  10.12126/see.2020.06.005
ZHOU D Q, YAN Y, YANG Y, et al. Simulation analysis of multi-field coupling field during electromagnetic emission[J]. Spacecraft Environment Engineering, 2020, 37(6): 561-569 doi:  10.12126/see.2020.06.005
Citation: ZHOU D Q, YAN Y, YANG Y, et al. Simulation analysis of multi-field coupling field during electromagnetic emission[J]. Spacecraft Environment Engineering, 2020, 37(6): 561-569 doi:  10.12126/see.2020.06.005
    • 电磁发射是一种全新概念的发射方式。电磁发射技术将电磁能变换为发射物体需要的瞬态动能,可在短距离内将发射物体加速至高速,能有效突破传统发射的速度和能量极限,在军事和民用领域具有巨大的优势和广阔的应用前景[1-2]。电磁发射过程中对发射装置加载脉冲大电流,伴随着磁场扩散、热传导、结构变形等物理过程,是具有瞬态特性的多物理场耦合问题[3]

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

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

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

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

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

      $$\nabla \times {{H}} = {{J}} + \frac{{\partial {{D}}}}{{\partial t}}\text{;}$$ (1)
      $$\nabla \times {{E}} = \frac{{\partial {{B}}}}{{\partial t}}\text{;}$$ (2)
      $$\nabla \cdot {{D}} = \rho\text{;} $$ (3)
      $$\nabla \cdot {{B}} = 0\text{。}$$ (4)

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

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

      $${{D}} = \varepsilon {{E}}\text{;}$$ (5)
      $${{B}} = \mu {{H}}\text{;}$$ (6)
      $${{J}} = \sigma ({{E}} + {{v}} \times {{B}})\text{。}$$ (7)

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

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

      $${\nabla ^2}\varPhi = {\text{-}} \frac{\rho }{\varepsilon }\text{。}$$ (8)

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

      $${\nabla ^2}{{A}} - \mu \varepsilon \frac{{\partial {{A}}}}{{\partial t}} ={\text{-}}\mu {{J}}\text{。}$$ (9)
    • 为了使用计算机编程实现电磁发射中电磁、电位的数值积分计算,本文使用伽辽金方法建立电位标量势方程的三维有限元格式,构造8节点六面体等参单元,其形函数N为三维空间参考坐标(ξ, η, ζ)的函数,

      $${N_i} = \frac{1}{8}(1 + {{\textit{ξ}} _0})(1 + {\eta _0})(1 + {\zeta _0})\text{。}$$ (10)

      式中:ξ0=ξiξη0=ηiηζ0=ζiζ

      以形函数N为权函数对电位标量势方程在计算域内积分,电位标量势的积分形式为

      $$\int_\varOmega {{{{N}}^{\rm{T}}}\left({\nabla ^2}\varPhi + \frac{\rho }{\varepsilon }\right){\rm{d}}\varOmega } = 0\text{;}$$ (11)

      通过三维高斯数值积分,电位标量势的有限元列式为

      $${{{K}}_{\rm{e}}}{{\varPhi}} = {{{R}}_{\rm{e}}}\text{。}$$ (12)

      同理,以形函数N为权函数对磁位矢量势方程在计算域内积分,磁位矢量势的积分形式为

      $$\begin{split} & \int_\varOmega {{{{N}}^{\rm{T}}}\nabla \cdot \nabla {{A}}(t + \Delta t){\rm{d}}\varOmega } - \int_\varOmega {{{{N}}^{\rm{T}}}\frac{{\mu \varepsilon {{A}}(t + \Delta t)}}{{\Delta t}}} {\rm{d}}\varOmega {\rm{ + }}\\ & \int_\varOmega {{{{N}}^{\rm{T}}}\frac{{\mu \varepsilon {{A}}(t)}}{{\Delta t}}} {\rm{d}}\varOmega {\rm{ + }}\int_\varOmega {{{{N}}^T}{{{J}}_{\rm{s}}}(t + \Delta t)} {\rm{d}}\varOmega = 0\text{;} \end{split}$$ (13)

      通过三维高斯数值积分,磁位矢量势的有限元列式为

      $${{{K}}_A}{{A}}(\Delta t + t) = {{{R}}_J}\text{。}$$ (14)
    • 电磁发射过程中产生热的途径主要有电枢与导轨摩擦生热、脉冲激励电流生热及感生涡流生热等;热扩散途径主要为热传导和热对流。本文主要考虑瞬态热传导问题的基本方程与有限元方法在热传导中的应用。

    • 为了模拟电磁发射过程的热分布与热传导,使用数值方法求解计算模型的热传导和热对流方程。

      热传导基本方程[5]

      $$c\frac{{\partial T}}{{\partial t}} + \nabla \cdot ({\text{-}}{\textit{λ}} \;\nabla T) = q\text{。}$$ (15)

      式中:c为计算域比热容;λ为计算域导热系数;T是温度场;q为热源体积密度。考虑电磁计算域内电流生热,热源体积密度可以写为q=J2/σ,其中,电流密度J包括外载电流、感应电流和涡流。

      热对流方程为

      $${\textit{λ}} \frac{{\partial T}}{{\partial s}}n = {\text{-}}h(T - {T_{\rm{a}}})\text{。}$$ (16)

      式中:h为热对流系数;Ta为外界热场。

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

      $$T = \sum\limits_{i = 1}^8 {{N_i}{T_i}(t)}\text{。} $$ (17)

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

      $$\int_\varOmega {{{{N}}^{\rm{T}}}\left(c\frac{{\partial T}}{{\partial t}} + \nabla \cdot ({\text{-}}{\textit{λ}} \;\nabla T) - q\right){\rm{d}}\varOmega } = 0\text{。}$$ (18)

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

      $${{C\dot T}} + {{KT}} = {{P}}\text{。}$$ (19)

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

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

      图  1  瞬态多场耦合计算程序框架

      Figure 1.  Framework of the coupled transient multi-field calculation program

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

      1)电场分析子模块:电流曲线插值,组装电场刚度矩阵;

      2)磁场分析子模块:组装磁场刚度矩阵、磁场远场边界;

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

      4)动态响应子模块:结构动态响应计算,电枢速度与位移计算;

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

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

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

      program main                !电磁发射过程计算主体程序  ……  call draw_screen               !输入数据界面  call input                  !节点、单元、边界数和电介质材料  mainloop : do mt=1,time/step          !移动单元和参考点      call sort_electrical_element       !导电单元      if (judge_elec_continu==1) then     !判断导体的连续性(电接触模块判断)          call electrical_element     !筛选导电单元          call no_zero_value       !选择非零数          ……          call form_elect        !组装全局单元刚度矩阵          call voltage_condition      !定义电压边界条件(约束电势边界条件                            消除刚度奇异矩阵)          call current_distribution     !求解电流分布          call magnetic_distribution    !求解电磁场分布          call thermal_distribution     !求解热场分布          call motion_velocity      !计算运动速度和位移(动态响应计算)          call output          !输出计算结果      else          call motion_velocity      !运动单元速度和位移          call output          !输出计算结果        ……  end do mainloop  ……

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

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

      图  2  轨道电磁发射计算模型

      Figure 2.  Model for calculating orbit electromagnetic emission

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

      表 1  计算所用的相关材料参数

      Table 1.  Material parameters for the electromagnetic calculation

      组件材料磁导率/
      (×10-6 mm·g·ms-2·A-2)
      电导率/
      (×105 mm-3·g-1·ms3·A2)
      导热系数/
      (×10-3 N·ms-1·K-1)
      导轨1.256 60.54401
      电枢1.256 60.34237
      磁场域真空1.256 600.023
    • 根据轨道电磁炮电磁、热、结构计算的载荷与边界情况,定义激励载荷1个(脉冲电流荷载),边界条件3个(无限电磁域、热对流边界、导轨位移约束)。

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

      图  3  轨道电流加载方向示意

      Figure 3.  Schematic diagram of current loading direction

      图  4  电流加载单元与加载节点示意

      Figure 4.  The current loading unit and the engaging nodes

    • 为了比对脉冲电流峰值对电枢发射速度、加速度的影响规律,分别对电流脉冲峰值为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  计算出的电枢发射速度曲线

      Figure 5.  Curve of calculated launching speed of the armature

      图  6  计算出的电枢位移曲线

      Figure 6.  Curve of calculated launching distance of the armature

      图  7  电枢加速度计算曲线

      Figure 7.  Curve of calculated acceleration of the armature

    • 分析电磁发射全过程导轨和电枢中的电流密度分布情况:电枢在导轨中滑行时,导轨的最高电流密度约为平均电流密度的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时刻电流密度分布云图

      Figure 8.  Cloud map of calculated current density after 0.4 ms of current loading

      图  9  电流密度时间历程曲线取值点

      Figure 9.  The location points for investigating the time profile of current density

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

      图  10  导轨测试点电流密度时间历程曲线

      Figure 10.  Time history of current density at typical positions of the rails

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

      图  11  电枢测试点电流密度时间历程曲线

      Figure 11.  Time history of current density at typical positions of the armature

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

      图  12  导轨与电枢电磁力计算结果矢量分布

      Figure 12.  Distributions of calculated electromagnetic force vectors of the rails and the armature

      图  13  z方向截面磁势分布云图

      Figure 13.  Cloud diagram of magnetic potential on cross section

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

      图  14  不同时刻的热密度分布云图

      Figure 14.  Cloud diagram of heat density at different moments

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

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

参考文献 (21)

目录

    /

    返回文章
    返回