一种求解细线瞬态散射的稳定算法

韦宇祥, 孙海涛, 马跃华, 王舒申, 胡彦宇

(上海机电工程研究所，上海 201109)

0 引言

电磁脉冲武器能够产生爆发性瞬态电磁脉冲，在金属导体上产生瞬态峰值很高的感应电流，对精密的光电子系统产生干扰、损伤和破坏。美国和俄罗斯等军事强国在电磁脉冲武器领域发展很快[1]，电磁脉冲弹药在近年的局部战争中都有实战应用。因此，针对电磁脉冲的作用机理及防护研究具有重大现实意义。

求解电磁脉冲在目标上产生的宽带时域响应是分析电磁脉冲作用机理的重要手段。常见的时域方法有时域有限差分法(FDTD)、时域有限元法(TDFE)和时域积分方程方法(TDIE)等。其中TDIE方法求解问题的边界条件自动满足，没有吸收边界问题，计算量相对较小，所以TDIE方法成为计算电磁学的重要研究方向。然而，求解TDIE普遍采用的时间步进法(MOT)在计算后期会出现严重的振荡，后期不稳定性问题影响了其应用。Tijhuisd等认为[2]，时间和空间离散会带来误差，在面元上进行内积运算时所作的各种近似也会产生误差，这些误差在时间步进的过程中会逐步积累，最终导致发散。Rynne等认为[3]，除了离散误差积累以外，还可能有以下2种因素：1)目标本身存在某些内谐振频率，MOT算法系统实际上是一个频率低通系统，这些谐振频率如果在MOT系统的通频带以内，就可能会引起谐振，从而发散；2)空间网格和时间离散的不精确，会“诱导”某些本来不会产生谐振的高频发生谐振。Rynne等提出，面元网格的最小尺度和时间步长应满足Courant准则，并采用三步或四步电流系数平均法抑制后期振荡[4-5]。电流系数平均法、或新基函数方法[6-9]和时间外推方法[10]，在一定程度上改善了后期不稳定性，但并不能完全解决问题。

本文采用显式MOT方法计算细导线在外加瞬变电场的辐射下的电流响应，分析其后期振荡问题。提出一种新的电流系数平均法，不仅可以改善前面时间步的电流值，而且能对当前时间步的电流值进行判断，并采取措施预防振荡，有效消除迭代计算的累积误差，改善后期不稳定性问题。

1 MOT方法

长度为2l，半径为a的细导线在入射脉冲电场的激发下将产生感应电流分布I(z,t)，并激发散射场，利用导体表面的切向场为零的条件，可以得到时域电场积分方程：

(1)

(2)

式中，ε0和μ0分别为真空中的介电常数和磁导率，c=1/(ε0μ0)1/2为光速。

假设细线上电流可以展开为：

(3)

这里有：

(4)

式中，Im(t)表示第m子域的电流系数，为待求量，是时间t的函数。

将式(2)离散：

Am,n=Az(zm,tn)

(5)

式(5)中的积分项为几何参数，称为阻抗系数。

(6)

将式(5)分离出k=m的自身单元项后可写为：

Am,n=Im(tn-|zm-zm|)Km,m+

(7)

(8)

将式(1)离散，用中心差分近似可得：

(Am+1,n-2Am,n+Am-1,n/(Δz)2-

(Am,n+1-2Am,n+Am,n-1)/(cΔt)2

项目位于黑龙江省哈尔滨市，建设于2014年。供热面积为20万m2，主要为居民住宅楼冬季采暖供热。供热系统采用原生污水源热泵系统，市政污水管网中的污水经过污水换热器换热后，进入热泵机组，通过消耗少量电能，获得热能。污水量约为2000t/h，温度为11℃。

(9)

Am,n=2Am,n-1-Am,n-2+(Δt)2Fm,n+

(10)

将式(10)代入式(7)可得：

(11)

式(11)表明，当前时间步电流系数Im ,n可以由前面时间步的结果进行求解，可交替应用式(7)、式(8)和式(11)式随时间逐步推进(MOT)进行计算，称为显式解(explicit solution)。MOT方法优点是无需进行矩阵求逆计算，缺点是后期出现振荡。

(12)

三步平均法在一定程度上抑制了振荡，延缓了振荡出现的时间。文献[11～12]提出了先五步平均后三步平均的方法，进一步提高了算法的稳定性，可以获得比三步平均法更好的结果：

(13)

(14)

2 数值计算结果

假设电磁脉冲武器发射的电磁脉冲可以用双指数脉冲模型描述[13]：

(15)

式中，E0=50 kV/m，k0=1.05为修正系数，α、β是表征脉冲前后沿的参数，α=4.0×107，β=4.76×108。

电磁脉冲直接照射到电路的细导线，细导线的半长度为l=0.2 m，半径为a=1.0 mm，分为26段，利用MOT方法求解。细导线上的电流响应如图1所示。

图1 电流响应

从图1可以看出，未经平均(MOT)的电流响应在0.1μs的时候就出现了明显的振荡，经过先五步平均后三步平均(Five steps MOT)的电流响应在0.4 μs后出现明显的振荡；在频域上，如图2所示，尽管经过平均的电流响应频谱在7 GHz附近的高频分量幅值远远小于原值，然而并不能完全消除振荡，所以在时域上反映为延缓了振荡出现的时间。

图2 电流响应的频谱

3 新的电流系数修正方法

下面分析电流振荡的曲线，研究如何采取措施修正电流系数，达到改善甚至消除后期振荡的目的。

在没有振荡的理想的平滑电流响应曲线中，如图3(a)所示，相邻的3个时间步的电流值Im,n-1、Im,n和Im,n+1可以近似认为在同一直线上；当出现振荡的时候，如图3(b)所示，这几个值将不在同一直线上。电流值Im,n和Im,n+1的中心为p，近似认为电流曲线经过Im,n-1和p，以此重新计算Im,n和Im,n+1。中心p的坐标为：

t0=(tn+tn+1)/2

I0=(Im,n+Im,n+1)/2

(16)

图3 电流响应曲线

电流I关于时间t直线方程为：

I=(I0-Im,n-1)/(t0-tn-1)t+Im,n-1-

(I0-Im,n-1)/(t0-tn-1)tn-1

(17)

特别地，需要对振荡情况进行判断：

(18)

应用新方法对细线的电流响应进行计算，结果如图4所示。很明显，新方法得到的电流响应曲线与先五步后三步平均法得到的曲线几乎完全重合，但是，如前所述，在0.4 μs后期，先五步后三步平均法的结果开始出现明显的振荡，而新方法得到结果趋于稳定；在频域上，两者的低频段几乎完全重合，在时域上反映为2种方法获得的初始响应几乎完全一致；在10 GHz附近的高频段，与先五步后三步平均法不同，新方法的结果几乎完全消除了高频分量，也就是说，获得了稳定的后期响应曲线。

4 结束语

本文分析了求解TDIE的普通MOT方法和电流系数先五步后三步平均法的后期稳定性。先五步后三步平均法相对而言改善了稳定性，但只是延缓了后期振荡出现的时间。通过对电流曲线的振荡进行分析，提出的新方法不仅改善了前面时间步的电流值，而且对当前时间步的电流值进行判断，并采取措施预防振荡，计算结果优于先五步后两步平均法的结果，消除高频分量，有效改善后期稳定性。■

图4 电流响应