摘要：针对单元烧蚀算法精度较低的问题,基于导轨式电磁发射装置的二维控制方程,推导出采用求解量移动法处理运动问题的伽辽金有限元离散方程,建立了电磁-温度-运动耦合场有限元计算模型。以单元节点为烧蚀判断对象,采用节点烧蚀算法处理电枢材料烧蚀问题,建立了电枢熔化波模型,对运动电枢的烧蚀速度进行了数值计算。模型的计算结果表明:在激励电磁感应强度为40 T、电枢运动速度为150 m/s的情况下,采用节点烧蚀法得到的熔化波烧蚀速度为Barber理论模型计算值的94%,而相关文献采用单元烧蚀法的计算值为Barber理论模型计算值的73%。因此,与采用单元烧蚀法相比,采用节点烧蚀法的计算值与Parks和Barber两个经典理论模型的计算值更加相近,验证了该方法的正确性。
关键词：导轨式电磁发射装置;运动边界;电枢;熔化波;烧蚀速度;耦合场;相变
在脉冲电流的驱动下,导轨式电磁发射装置的电枢在两轨道间做加速运动,电枢与轨道存在着滑动电接触,电枢的高速运动对轨道及电枢内部的磁场、电流分布具有显著的影响[1-2],这种现象称为速度趋肤效应。文献[3-4]采用有限差分法建立了二维电磁场、温度场数学离散模型,文献[5-6]则采用有限元法建立计算模型,得到了电枢运动下的电流分布特点。
由于电流集中于电枢尾翼,所以在焦耳热的作用下,电枢材料最先在尾翼处发生熔化,熔化部分在电磁力的作用下迅速从电枢本体剥离,随后熔化区域在新的电枢边界处重新形成,这种熔化形式称为熔化波。对于电枢熔化烧蚀问题,国外学者Parks和Barber分别从理论上推导出了二维情况下的电枢熔化烧蚀速度[7-8]。此后,国内外学者分别从一维、二维模型出发,采用有限差分法、有限元法建立了基于单元烧蚀算法的电枢烧蚀模型,得出了不同电枢速度下的熔化波烧蚀速度[9-12]。从公开文献中的分析结果可知:采用单元烧蚀算法时,若电枢运动速度较大,则其计算结果与理论值偏差过大。本文认为:选取单元作为烧蚀判断对象时,实际上是取单元内节点的温度算术平均值作为单元的温度,在单元内温度梯度较大的情况下,势必会产生较大的误差,甚至可导致相变过程的数值计算方法发生失效。
为此,本文首先建立了考虑电枢运动及材料相变的电磁-温度耦合场的二维有限元模型,并针对单元烧蚀算法应用在电枢熔化烧蚀计算中时精度较差的问题,采用节点烧蚀法建立电枢熔化波模型,对运动电枢的烧蚀速度进行数值计算,得到了与经典理论计算值更加接近的电枢熔化烧蚀速度。
1 导轨式电磁发射装置控制方程
假设电枢与轨道间的接触为理想接触,不考虑电枢与轨道间的接触电阻,当电枢材料发生熔化时,电枢材料立即在电磁力的作用下脱离电枢,电枢的实际物理边界发生相应变化,且忽略电枢非熔化区域的形变。
1.1 电磁场控制方程
求解域如图1所示,采用静止坐标系,二维下电枢与轨道内的磁扩散方程[2-3]为
(1)
求解域边界条件为
(2)
式中:σ、μ、vx分别为材料的电导率、磁导率及电枢的运动速度,均是与区域相关的参数,假设轨道沿x轴的负方向运动,则vx为轨道运动速度,而如果电枢静止不动,则vx为0;Bz为电磁感应强度,对于二维问题,仅存在z方向的分量;J为施加的电流密度,其方向平行于求解域平面,对于边界S2、S3,J仅存在x方向的分量,而对于边界S1、S7,J仅存在y方向分量。模型中边界S1、S2处施加的电磁感应强度用于等效引入电流密度激励。
图1 导轨式电磁发射装置二维计算模型
1.2 温度场控制方程
如图1所示,在二维情况下,电枢与轨道内的热传导方程满足
(3)
式中:κ、c、ρ、EL分别为材料的热导率、比热容、质量密度和熔化潜热。熔化潜热即材料发生熔化相变过程中单位质量所吸收的热量。
在温度场中,求解域的边界条件为
(4)
1.3 网格移动法处理运动问题
对于运动电磁场问题,为消除运动项系数矩阵对称性的影响,通常采用独立坐标系。求解域的运动部分与静止部分分别采用不同的坐标系,将控制方程中的运动项消去,通过网格的实际移动来考虑运动项的影响[13]。
1.4 求解量移动法处理运动问题
假设电枢静止,轨道沿x负方向运动,通过移动求解量考虑运动项的影响,其等效性说明如下。
矢量B是时间和空间上的函数,在二维情况下
(5)
由式(1)、式(5)可得
(6)
对于运动部分,运动坐标系下的时间离散表达式为
(7)
因此,在t时刻的节点求解量Bz计算完成后,运动部分运动的距离为ds。在进行t+Δt时刻的计算时,需用到节点在t时刻时所处位置的求解量。采取的处理方法为:网格不发生移动,t时刻计算所得的节点求解量沿运动方向平移ds后作为t+Δt时刻计算时所用的求解量数值,即通过求解量的移动来模拟电枢的运动。
1.5 有限元离散方程
根据求解区域的电磁场和温度控制方程,采用移动求解量法处理运动问题,最终得到导轨式电磁发射装置二维耦合场计算模型的伽辽金有限元离散方程为
(8)
(9)
式中:ke、me为电磁场中的系数矩阵
为温度场中的系数矩阵;Fe为温度场中的右端向量。
采用四节点四边形单元,则系数矩阵和右端向量在局部坐标系下的形式为
(10)
(11)
(12)
(13)
(14)
式中
(15)
下角标i、j为单元节点编号;ξ、η为局部坐标系下的坐标轴;J为雅克比矩阵,表示整体坐标系与局部坐标系之间的转换关系。
对于形如式(10)～式(14)中的数值积分,采用高斯积分法进行计算。对于瞬态场问题,采用Crank-Nicholson法对时间进行离散。
2 电枢材料烧蚀的处理
电枢在加速运动过程中,由于电流在电枢尾翼处的集中,可导致此处的材料发生熔化,造成电枢的烧蚀,因此在温度场计算过程中,需考虑材料的相变过程。本文采用热源法,为避免单元温度梯度较大导致的计算精度降低、甚至相变过程的数值计算失败,以节点为判断对象:当单元内的节点达到材料熔点时,节点温度不再发生变化;当其吸收的热量达到熔化潜热时,节点处发生相变,认为节点及其附近发生了烧蚀。当节点发生熔化时,其对应位置的电枢材料在电磁力的作用下迅速脱离电枢本体,磁场迅速扩散至电枢的非熔化区域边界,即电磁感应强度的边界条件发生了改变。由式(16)可知:磁场在材料中的渗透深度δ与材料电阻率ρ的1/2次方成正比。
(16)
在程序中增大熔化位置的电阻率,可增加磁场的渗透深度,使磁场较快地渗透至电枢的非熔化边界,相当于更新了Bz的边界条件。在程序编写过程中发现:将熔化位置的电阻率增大4个数量级时即可得到较为准确的结果。因此,对于存在烧蚀节点的单元,采用非均匀材质单元的处理方法,在高斯积分的过程中,将熔化节点的电阻率增大4个数量级,这种方法称为节点烧蚀法。
具体的处理过程描述如下:图2为正处于相变过程的单元,其中节点1、4为已发生烧蚀的节点,节点2、3为非烧蚀节点,在电磁场计算过程中,其单元系数矩阵ke的元素为
(17)
在温度场计算过程中,其单元的右端向量Fe中的元素为
(18)
式中:n为每一维的高斯积分点数,本文取n=2;σe(ξp,ηq)为单元内高斯积分点(p,q)处的电导率,其具体数值与单元内节点的相态有关,对于烧蚀节点1、4,其对应的高斯积分点g1、g4的电导率减小4个数量级。至此,完成了电枢烧蚀的数值模拟计算。
图2 局部坐标系下的单元高斯积分点示意图
考虑材料相变的二维电枢熔化波有限元计算流程如图3所示。
图3 电磁-温度-运动耦合场有限元计算流程图
3 有限元计算结果分析与验证
为方便对比分析,本文模型所用的材料参数及求解域均与文献[12]相同。求解域如图1所示,所用材料的参数见表1。本文模型与文献[12]中Stefani模型的求解区域单元剖分密度相同(每个单元边长为0.2 mm)。
在激励磁感应强度Bz为40 T(对应的区域为S1、S2边界)、电枢运动速度vx为100 m/s的情况下,电枢的熔化波烧蚀演变如图4所示,对应的电流密度云图如图5所示。当电枢材料发生熔化烧蚀时,电流集中点从电枢尾翼处移至熔化波前端,形成了烧蚀间隙在电枢轨道接触面由电枢尾翼至电枢前端的扩展过程,烧蚀间隙的扩展速度即为熔化波烧蚀速度。
表1 轨道电枢材料的参数
在Bz分别为32和40 T、vx分别为100和150 m/s的情况下,通过本文模型计算得出的熔化波烧蚀深度(hm)随时间的变化关系如图6所示。为了验证本文模型的计算结果,将根据本文模型计算出的电枢熔化波烧蚀速度与经典Parks、Barber理论模型的计算值及Stefani模型的计算值进行了对比,结果见表2。Parks和Barber根据其理论模型建立的二维固体电枢熔化波烧蚀速度公式分别为
。可以看出,电枢熔化波烧蚀的理论速度是与时间无关的参数。
图4 电枢熔化波烧蚀演变图
图5 熔化波烧蚀过程的电流密度
表2 不同模型计算的电枢烧蚀速度对比
由图6a可知:对于特定的某一工况,本文模型计算得出的熔化波电枢烧蚀深度与时间近似成线性关系,即不同时刻的电枢熔化波烧蚀速度基本不变,这与Parks和Barber的理论公式相吻合。由图6b可知:用Stefani模型计算出的烧蚀深度与时间成明显的非线性关系,且在计算前期,电枢熔化波烧蚀速度计算值明显偏小。
由表2的对比可知:当Bz=40 T、vx=100 m/s时,本文模型计算出的熔化波烧蚀速度略大于Stefani模型的熔化波烧蚀速度计算值,但小于Parks、Barber理论公式的计算值;当Bz=40 T、vx=150 m/s时,本文模型计算出的熔化波烧蚀速度为Barber理论公式计算值的94%,而Stefani模型的熔化波烧蚀速度计算值为Barber理论公式计算值的73%。
(a)本文模型计算结果
(b)Stefani模型计算结果[12]
图6 熔化波烧蚀深度与时间的关系
显然,在烧蚀计算的整个过程中,Stefani模型的计算结果均偏小,主要原因在于:Stefani模型选取单元作为烧蚀判断对象,实际上是取单元内节点的温度算术平均值作为单元的温度,当单元内节点温度差异较大时,会导致其确定的单元烧蚀时刻滞后,最终造成计算前期材料烧蚀时刻延迟,且电枢熔化波烧蚀速度计算值偏小。
图7 不同电枢速度下的熔化波烧蚀速度
取Bz=40 T,并增加求解域的单元剖分密度(每个单元边长为0.05 mm),计算不同电枢速度下的熔化波烧蚀速度,并与Parks和Barber理论公式的计算值进行对比,结果如图7所示。由图7可知:在不同电枢速度下,本文模型计算出的熔化波烧蚀速度值介于Parks与Barber理论公式计算值之间。
由以上对比分析可知:与单元烧蚀法相比,在激励磁感应强度场较大、电枢运动速度较高的情况下,采用本文提出的节点烧蚀法可得到稳定的、且与经典理论公式计算值较吻合的烧蚀速度。
4 结 语
本文通过求解量移动的方法处理运动边界,采用节点烧蚀法处理电枢的烧蚀问题,建立了考虑材料相变的导轨式电磁发射装置二维电磁-温度-运动耦合场计算模型,并进行了电枢熔化波烧蚀速度的数值计算。计算结果表明:与单元烧蚀法相比,采用本文模型得到的熔化波烧蚀速度与Parks和Barber两个经典模型得到的理论计算结果更接近,间接验证了本文方法的正确性。
