直线推进电磁能装备是近年来迅速发展的一种全新概念的加速方式，具有速度快、效率高、启动迅速、可靠性好等技术优势，在军用和民用领域都有着广泛的应用前景[1-3]．直线推进电磁能装备可在数十毫秒内将电枢加至高速，导致其实验环境相对复杂，研究比较依赖于数值模拟计算[4-7]，因此要实现进一步的推广和应用，须提高仿真模型的计算精度实现全真模拟计算．直线推进电磁能装备在实际工作过程中，轨道通入大电流，当大电流通过枢轨接触面时，产生巨大温升，严重降低了推进电磁能装备的效率和性能[8-12]，因此必须在仿真计算中综合考虑接触电阻所产生的影响，才能使其更加接近实际推进工况下的结果．作为表征直线推进电磁能装备性能的两个关键参数，电枢速度和电磁推力均与电感梯度有着直接关系，目前研究主要集中在动态电感梯度的计算上．早期研究认为电感梯度仅与轨道尺寸相关，当轨道尺寸不变时，电感梯度也不变[13-14]，文献[15]分析了电感梯度随频率的分布．文献[16]考虑速度趋肤效应，通过引入速度频率的方式对动态电感梯度进行了计算．文献[17]提出一种综合考虑导轨尺寸和电流扩散影响的电感梯度解析计算方法．在实际工作过程中，接触电阻会使直线推进电磁能装备电感梯度受到影响，上述文献对直线推进电磁能装备进行了充足的理论分析，但均未考虑动态接触电阻所产生的影响．本研究首先通过有限元方法分析了接触电阻、电感梯度、电枢速度等物理量之间的相关性，建立动态接触电阻改进模型；然后采用改进粒子群优化-遗传算法(particle swarm optimization- genetic algorithm，PSO-GA)进行参数辨识，得到考虑动态接触电阻的电感梯度及电枢速度预测函数；最后，通过实验测试数据验证了改进模型的有效性，为直线推进电磁能装备全真模拟研究提供了参考．1 直线推进电磁能装备有限元分析1.1　有限元分析方法电磁能装备直线推进过程可以看作是含有运动导体的三维瞬态涡流场求解问题，直线推进电磁能装备主要包括内外轨道及电枢两个部分．忽略位移电流，采用A-φ法可以推导出电枢与轨道的扩散方程分别如下[18]B=∇×A;E=-∂A∂t+∇φ;∇2A-μaσa(∇φ+∂A∂t-v∇×A)=0 ;   ∇2A-μrσr(∇φ+∂A∂t)  =0   ,                            (1)式中：B为磁感应强度；E为电场强度；∇为哈密顿算子；A为磁矢位；φ为电标位；μa与μr分别为电枢和轨道的磁导率；σa与σr分别表示电枢和轨道的电导率；v为电枢速度；t为时间．考虑到磁感应强度在z方向相互抵消，仅在y方向存在分量，故磁矢位和电流密度仅在x和z方向存在分量．除此之外，电枢速度仅考虑x方向分量，因此忽略∂Ax/∂y，∂Az/∂y，Ay，vy，vz，将方程组(1)写成如下分量形式：∂2Ax∂x2+∂2Ax∂z2-μaσa∂φ∂x+∂Ax∂t=0;∂2Az∂x2+∂2Az∂z2-μaσa∂φ∂z+∂Az∂t-vx∂Ax∂z-∂Az∂x=0;∂2Ax∂x2+∂2Ax∂z2-μrσr∂φ∂x+∂Ax∂t=0;∂2Az∂x2+∂2Az∂z2-μrσr∂φ∂z+∂Az∂t=0.计算域的边界条件，主要包括电枢-轨道接触面条件、轨道-空气域接触面条件以及电枢-轨道接触面条件nij×(Ei-Ej)=0;nij×(Hi-Hj)=0;nij⋅(Ji-Jj)=0;nij⋅(Bi-Bj)=0,式中：nij(i，j=1，2，3)为分界面的法向量；i=1，j=2为电枢-轨道接触面条件；i=1，j=3为电枢-空气域接触面条件；i=2，j=3为轨道-空气域接触面条件；H为磁场强度；J为电流密度．根据有限元网格剖分、单元离散得到有限元矩阵方程[19]KU+M∂U∂t=F，(2)式中：K，M和F分别为总体刚度矩阵、总体质量矩阵和总体载荷向量；U为总体未知函数向量．利用欧拉差分方法对方程(2)进行时间离散，得到如下形式KUt+Δt+MUt+Δt-UtΔt=Ft+Δt，式中Δt为时间步长．电磁推力密度为fe=J×B．1.2　有限元仿真结果实验采用的激励电流如图1所示，拟合得到的激励电流(I)波形也示于图1，表达式为：I=I0 sin(2.5πt)                (0 ms≤t0.2 ms);I0e-0.61(t-0.2)             (0.2 ms≤t≤0.6 ms),      (3)式中I0=318 kA为激励电流峰值．10.13245/j.hust.240842.F001图1实验电流及其拟合的电流波形根据电磁能装备直线推进模型，接触电阻可以表示为Rc=Um/Il-Ra，(4)式中：Um为炮口电压；Il为回路电流；Ra为电枢电阻．依据仿真计算得到的炮口电压和回路电流值，结合式(4)，计算得到电枢接触电阻随时间的变化曲线如图2所示．可以看到：在推进的初始阶段，枢轨接触电阻较大，然后急剧减小，最后趋于稳定．这是由于电枢开始滑动瞬间，微观接触面积急迅速减少，因此接触电阻较大．开始滑动后，接触面积逐渐增加，导致电阻急剧衰减．0.2 ms后，由于枢轨接触面之间产生了熔膜，接触面积保持稳定，进而接触电阻缓慢减小．10.13245/j.hust.240842.F002图2枢轨间接触电阻随时间的变化曲线根据图2中接触电阻随时间变化关系，定义如下拟合函数Rc(t)=∑i=0muiti，式中：Rc(t)为接触电阻随时间(t)变化的函数；m为多项式最高阶数；ui为定义的拟合系数．实际推进过程中，考虑到电枢所受空气阻力及摩擦力影响，列写电枢驱动力方程如下mdvafdt=Fe-fvaf2-μ0m0g，(5)式中：vaf为考虑空气阻力及摩擦力影响的电枢速度；Fe为电枢所受电磁推力；f为风阻系数，取7.9×10-5 kg/m；μ0为滑动摩擦系数，取0.5；m0为电枢质量，取10 g；g为重力加速度，取9.8 m/s2．求解微分方程(5)，得到考虑空气阻力及摩擦力影响下的电枢速度为vaf=Fe-μ0m0gf∙exp[2t((Fe-μ0m0g)f/m0)]-1exp[2t((Fe-μ0m0g)f/m0)]+1．(6)当空气阻力与摩擦力远远小于电磁推力，即fvaf2 /Fe→0，μ0m0g/Fe→0时，式(6)变为vaf=Fe/(m0t)．在推进过程中，采用有限元法计算得到电枢所受到电磁推力最高可达到5×104 N，空气阻力与摩擦力与其相比十分微小，因此可以忽略其影响．为了给后续改进模型的建立提供数据，选取不同的接触电阻值进行计算，得到不同接触电阻下的电枢速度(v)随时间变化曲线如图3所示．10.13245/j.hust.240842.F003图3不同接触电阻下的电枢速度随时间变化曲线电枢推进过程中，受趋肤效应的影响，电流主要集中在电枢中心处与轨道电枢尾翼接触处，而非均匀分布，导致装置电感梯度随着时间产生动态变化．除此之外，接触电阻的存在，导致直线推进电磁能装备所产生的电磁推力有所减小，直线推进电磁能装备整体效率受到影响，电感梯度因此受到影响．为了给后续改进模型的建立提供数据，选取不同的接触电阻值进行计算，得到不同接触电阻下电感梯度(L)随时间变化曲线，如图4所示．10.13245/j.hust.240842.F004图4不同接触电阻下电感梯度随时间变化曲线2 直线推进电磁能装备模型建立进行上述有限元计算时，默认将接触电阻设为常数．在实际推进过程中，接触电阻的动态变化会影响电感梯度以及电枢速度大小，因此采用参数辨识的方法建立考虑动态接触电阻的改进模型．根据总结出的规律建立动态接触电阻改进模型，并通过参数辨识得到改进模型预测函数．2.1　电枢速度预测函数不同时刻下电枢速度随接触电阻变化曲线如图5所示．可以看出：在同一时刻电枢速度随着接触电阻的上升逐渐下降，基本呈现线性变化．10.13245/j.hust.240842.F005图5不同时刻下电枢速度随接触电阻变化曲线综上，可采用近似描述电枢速度随接触电阻的变化关系vi(Rc(t))=rRc(t)+s，式中：vi为受接触电阻影响的电枢速度影响因子；r和s为定义的拟合系数．根据牛顿第二定律公式及电感梯度计算电磁推力公式Fe=m0dv0/dt;Fe=L0I2/2, (7)式中：v0为电枢速度解析值；L0取1.0 μH/m．结合式(3)激励电流公式，采用解析法计算得到电枢速度解析值如下[20]：v0(t)=L0I02m014t-110πsinπt0.2               (0 ms≤t≤0.2 ms);L0I02m01-e-1.22t-0.22.44+0.05      (0.2 mst≤0.6 ms). (8)电枢速度解析值如图6所示．可以看出：电枢速度仿真值与解析值曲线相似，因此可以通过解析值乘以一个系数的方式，将电枢速度仿真值表达出来，定义电枢速度仿真值函数关系v(t,Rc(t))=k0(t,Rc(t))v0(t)，(9)式中：k0为定义的倍乘系数．10.13245/j.hust.240842.F006图6电枢速度解析值不同接触电阻条件下k0随时间变化曲线如图7所示．可以看出：相同接触电阻下，k0随时间动态变化，在0.0～0.2 ms上升，0.22～0.6 ms逐渐趋于平稳．同一时刻，k0值随接触电阻的上升逐渐减小．因此，可以将k0表示为一个仅与时间相关的函数与一个仅与接触电阻相关的函数相乘的形式，定义k0如下k0(t,Rc(t))=k(t)vi(Rc(t))，(10)式中k(t)是一个仅与时间相关的函数，10.13245/j.hust.240842.F007图7不同接触电阻下k0值随时间变化曲线k(t)=∑j=0nwjtj，(11)其中，n为多项式最高阶数，wj为定义的拟合系数．结合式(7)~(11)函数关系，定义改进模型电枢速度预测函数如下：v(t,Rc(t))=k0(t,Rc(t))v0(t);k0(t,Rc(t))=k(t)vi(Rc(t));k(t)=∑j=0nwjtj;vi(Rc(t))=rRc(t)+s. (12)2.2　电感梯度预测函数0.2 ms后，电感梯度逐渐趋于稳定，近似为一条直线，因此可以采用以下公式描述电感梯度随时间变化关系Lt=a1t+b1，(13)式中a1和b1为定义的拟合系数．不同时刻电感梯度随接触电阻变化情况如图8所示．可以看出：电感梯度随着接触电阻上升逐渐下降，二者近似为线性变化，因此可采用如下公式近似描述电感梯度随接触电阻变化关系L(Rc(t))=a2Rc+b2，(14)式中a2和b2为定义的拟合系数．10.13245/j.hust.240842.F008图8不同时刻电感梯度随接触电阻变化曲线结合式(13)和(14)，定义改进模型电感梯度预测函数如下      L(t,  Rc(t))=L(t)+L(Rc(t))=a1t+a2Rc(t)+b0, (15)式中b0=b1+b2．2.3　直线推进电磁能装备模型参数辨识与验证联立式(12)和(15)，即可得到考虑动态接触电阻的改进模型，此模型可以将电感梯度与电枢速度考虑为与接触电阻及时间有关的二元函数．为了使改进模型预测结果有效，所建立的改进模型预测值与仿真值在相同接触电阻下误差应尽可能小，定义ev=1M1N∑i=0M∑j=0Nvi,j-vi,j'vi,j';eL=1M1N∑i=0M∑j=0NLi,j-Li,j'Li,j',式中：ev为电感梯度平均相对误差；eL为电枢速度平均相对误差；M为接触电阻仿真组数；N为时间节点数；vi,j和Li,j分别为对应组的电感梯度及电枢速度预测值；vi,j'和Li,j'分别为对应组的电感梯度及电枢速度仿真值．通过对电枢速度求导的方式得到电枢加速度avt,Rc(t)=dvt,Rc(t)dt，除此之外，还可以电磁推力计算公式计算电枢加速度，即aLt,Rc(t)=I2(t)2m0Lt,Rc(t)，解析计算过程中，二者相等．但在实际计算过程中，由于考虑了接触电阻的影响，因此二者不等．对上述公式进行离散处理，得到：       avk+1,Rck+1={vk+1,Rck+1-vk,Rck}/ΔT;aLk+1,Rck+1=I2k+12m0Lk+1,Rck+1．(16)根据式(16)，定义ea=1T∑k=0Tavk,Rck-aLk,Rckavk,Rck．设定ev，eL和ea为目标函数，对所建立改进模型进行参数辨识，使目标函数为最小．考虑到目标函数中含有一阶微分项，且为非线性，一般辨识方法容易陷入局部最优．因此采用改进的PSO-GA混合算法进行参数辨识[21]，在PSO算法的基础上，通过GA算法的基因序列交叉、变异选出每次迭代的最优个体并保存信息，可有效解决群体因计算趋于稳定而陷入局部最优的问题．PSO-GA混合算法流程如图9所示.10.13245/j.hust.240842.F009图9PSO-GA混合算法流程当m，n=6时，计算得到改进模型参数辨识值如表1所示．采用上述辨识值所得到的Rc(t)与k(t)拟合曲线的R2值分别为0.987 1与0.996 0，说明拟合程度良好．10.13245/j.hust.240842.T001表1改进模型参数辨识值参数辨识值参数辨识值参数辨识值u00.007w00.556r-0.034 6u1-0.029w19.279s1.090 7u2-0.431w2-71.375a1-0.125 0u34.166w3302.410a2-0.033 3u4-13.751w4-707.478b01.185 0u519.910w5848.244u6-10.710w6-405.440优化过程中的迭代误差(E)如图10所示，p为迭代次数．可以看出：当计算开始时，模型初始的迭代误差约为16.64％，随着迭代的进行，最终误差在970代左右逐渐稳定在1.52%．10.13245/j.hust.240842.F010图10参数辨识过程迭代误差对所建立的改进模型进行验证，另取10组有限元计算值与预测值进行比较，计算得到相对误差如表2和3所示．可以看出，改进模型预测值与同组仿真值相对误差均控制在3%以内．因此，可以认为本研究所建立的改进模型能够较为准确地预测电感梯度与电枢速度的大小．10.13245/j.hust.240842.T002表2电枢速度相对误差组号仿真值/(m∙s-1)预测值/(m∙s-1)相对误差/%1849.27835.241.6521 488.101 477.670.7031 720.341 708.570.6841 889.991 861.981.4852 121.962 069.252.4810.13245/j.hust.240842.T003表3电感梯度相对误差组号仿真值/(μH∙m-1)预测值/(μH∙m-1)相对误差/%10.960 20.971 21.1521.073 81.044 72.7131.086 31.093 90.7041.098 71.113 61.3551.111 21.135 52.183 电磁能装备直线推进实验验证为了验证所建立的改进模型的有效性，搭建了电磁能装备直线推进实验平台，如图11所示．脉冲电源通过充电机充电，产生激励电流，作用于电枢产生电磁推力推动电枢发射．10.13245/j.hust.240842.F011图11直线推进电磁能装备实验平台原理图采用光幕靶来测量电枢速度，直线推进电磁能装备实验平台测速原理如图12所示．电枢速度测试方法为：光幕之间距离固定为505 mm，两个光幕把分别与示波器连接，当电枢受到电磁推力飞出并通过光幕时，两个光幕会在示波器中产生信号．实验后，读取两个光幕把产生信号的时间间隔，即为电枢飞出后通过两个光幕的时间间隔，根据时间间隔与光幕之间距离即可算出电枢的出口速度．10.13245/j.hust.240842.F012图12直线推进电磁能装备实验平台测速原理根据所搭建的实验平台，试验得到实际工况下电枢速度，并与考虑动态接触电阻前后的数值进行比较，误差分析结果如表4所示．可以看出：考虑动态接触电阻的改进模型预测值误差更小，更加接近于实际推进情况，能够较为准确地预测实际工况下的输出情况．10.13245/j.hust.240842.T004表4实验值与改进前后计算值相对误差分析值类型电枢速度/(m∙s-1)相对误差/%实验值2 295.45改进前1 986.1213.48改进后2 430.195.874 结论a．采用有限元方法计算得到了枢轨间接触电阻、电感梯度及电枢速度曲线，为电磁能装备直线推进过程中的数值计算问题提供了一定的理论依据．b．研究了动态接触电阻影响下的直线推进电磁能装备电感梯度、电枢速度变化规律．同一时刻下电感梯度与电枢速度值随着接触电阻的上升逐渐下降，基本呈现线性关系．c．根据电感梯度、电枢速度变化规律，建立了考虑动态接触电阻的直线推进电磁能装备改进模型，采用课题组提出的PSO-GA混合算法对所建立模型进行参数辨识，并通过数据验证了改进模型的正确性．d．搭建了直线推进电磁能装备实验平台，对电枢速度进行测试．实验结果表明：建立的改进模型预测值与实验值误差更小，进一步证明了该模型的有效性．改进模型可在一定程度上提高计算精度，为后续电磁能装备直线推进全真模拟研究提供了一定的理论依据．