一種求解細線瞬態散射的穩定算法

韋宇祥, 孫海濤, 馬躍華, 王舒申, 胡彥宇

(上海機電工程研究所，上海 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 電流響應