基于TDIE?FDTD混合算法的電磁環境效應分析

摘 要： 電磁環境效應是在現代航空、航海、航天以及地面系統設計中必須考慮的因素。針對帶有線纜結構的復雜系統的電磁環境效應問題，發展了一種時域積分方程（TDIE）和時域有限差分（FDTD）的混合方法，采用時域平面波（PWTD）加速的TDIE方法求解場問題，高階FDTD（2，4）方法求解線纜響應。數值算例表明了方法的有效性。

關鍵詞： 電磁環境效應； 線纜； 時域積分方程； 時域有限差分； 時域平面波

中圖分類號： TN01?34 文獻標識碼： A 文章編號： 1004?373X（2013）21?0157?06

0 引 言

隨著數字設備和集成電路的廣泛應用，現代電子信息系統對雷電、高功率脈沖武器等產生的強電磁環境極為敏感。這種瞬態電磁場會通過天線、孔縫或外部電纜，耦合到信息設備機箱內部線纜網絡從而對終端電路產生電磁干擾（EMI），甚至造成損傷。

在系統設計階段就對其進行電磁環境效應分析無疑是一種有效手段。但對于典型的系統級分析會面臨多尺度問題，電子信息系統的載體一般是電大或超電大尺寸，系統的線纜網絡在一維上（傳輸方向）電大、另一維（橫截面）電小。

為了解決極其復雜的電磁環境效應分析問題，20世紀70年代Baum提出了電磁拓撲理論[1]（Electromagnetic Topology，EMT），其基本思想是將復雜系統分割為相互聯系的空間，各空間單獨求解，再將結果進行綜合，并在此基礎上發展了BLT（Baum?Liu?Tesche）方程方法[2]，但EMT理論人為地對系統進行分割，分析誤差較大。傳統的數值計算方法，如矩量法[3]（Method of Moment，MoM）、有限元法（Finite Element Method，FEM） [4]、時域有限差分方法（Finite Difference Time Domain，FDTD）[5]、傳輸線矩陣法（Transmission Line Matrix，TLM）[6]、時域積分方程法[7]等等，以及其加速、改進算法等[8?12]，解決多尺度問題、線纜加載問題仍有較大難度。一種可行的解決方案是將多種方法進行混合，已有的代表性的成果包括：Ba?c?采用FFT加速的TDIE方法求解帶線纜結構飛機的電磁環境效應問題[13]；Xie采用FDTD方法和SPICE等效電路模型求解線纜網絡的響應[14?15]； Liu采用FDTD方法和MNA方法混合求解帶線纜和非線性負載結構的系統響應[16?18]。另外，一些電磁場商業軟件（如CST、HFSS等）也紛紛推出復雜系統的電磁兼容和電磁環境效應仿真解決方案。

本文發展了一種基于TDIE?FDTD的混合方法，采用PWTD加速的TDIE方法求解場問題，高階FDTD（2，4）方法求解線纜響應，結合兩種方法的優勢對系統的電磁環境效應進行有效預測，數值算例表明了方法的有效性。

1 PWTD加速的TDIE方法求解場問題

基于Maxwell方程組，金屬問題的時域電場積分方程表達式為：

[n×n×Eir，t=n×n×μ04πs1R?Jr，τ?tds-14πε0?sds-∞τ??Jr，tRdt r，r∈s]

式中：[Eir，t]表示入射場；[Jr，t]為表面電流；[n×n]表示取矢量的表面負切向分量；[?]作用于場[r]上；[?]作用于源[r]上；[s]表示電流源分布表面；[τ=t-Rc]表示時間延遲。

傳統的求解方法是時間遞推方法：

式中：[Ij]為第[j]時刻電流系數向量；[Vj]為第[j]時刻激勵向量；[Zi]為阻抗矩陣。

MOT方法的主要計算量在于式右端的矩陣矢量乘積求和運算，它表示先前各時刻表面電流源對當前時刻各場點的貢獻，其運算量為[O（N2s）]。由于上述求和運算需要進行[Nt]個時間步，經典MOT算法的總計算量為[O（NtN2s）]。實際的計算實體（飛機、艦船）基本都是電大尺寸問題，計算量和存儲量都十分巨大，因此必須采用加速方法來進行。

時域平面波（PWTD）方法的機理在于瞬態遠場的平面波展開，通過這種展開可以減少式右端的求和運算。對于空間中某點的散射場，其貢獻來源可分為兩部分：一部分由該場點所在組的近場組[NFP（α）]中的源所產生；另一部分由遠場組[FFP（α）]中的源產生。因此，利用PWTD加速MOT算法求解電磁散射的迭代公式可表達如下：

[Z0Ij=Vj-α∈NFP（α）l=0j-1Zααj-lIαl- α′∈FFP（α）n∈α′（n）fm（r），LeJn（r，t）t=lMtΔt+iΔt] （3）

式中：[Le·]表示電場積分算子；[Mt]為基本子信號離散序列的點數，[i]為整數，且[0≤i

PWTD的計算復雜度為這兩部分計算量之和，采用多層或加窗的PWTD總的計算量可以降至最低[ONtNslogNs]能夠極大地提高計算速度。

2 高階1D?FDTD求解線纜響應

外部電磁場在傳輸線上會激勵其傳輸線模（差模）和天線模（共模）電流，但對于傳輸線終端有影響的只會是傳輸線模電流，另外，在有接地板存在的情況下（實際情況也如此），產生的天線模電流極小，基本可忽略，因此，采用基于傳輸線理論的方法來求解線纜響應。作為一種準靜態方法，傳輸線理論忽略了縱向場分量的影響，在求解波長與線纜橫截面尺寸可比擬的情況下（頻率很高）會有誤差，但在一般情況下還是能夠保證較高精度。求解傳輸線方程的頻域方法的優勢是計算速度快，能夠處理頻變的參數；但其劣勢是難以進行非均勻線和非線性負載的計算，因此，本文采用高階FDTD（2，4）（時間具有二階精度，空間具有四階精度）方法求解線纜響應，由于高階算法的低數值色散性，采用較少的離散點就能獲得滿意的數值精度，提高了計算效率。

復雜系統中，為了減少強電磁脈沖對系統內部電纜線路的輻射耦合，通常采用屏蔽電纜。因為屏蔽的不完整性（能量對材料的穿透、編織屏蔽層存在孔縫、兩端屏蔽層未360°搭接等），芯線也能感應到電流，這一過程一般用屏蔽層的轉移阻抗來表示。其定義是單位長的電纜中有單位電流流過屏蔽層時在電纜芯線與屏蔽層間所形成的開路電壓。實際情況中，電纜的轉移阻抗是與頻率有關的函數，且隨著頻率增加而減小。而為了能夠直接進行時域分析，本文忽略轉移阻抗的頻率相關性，這種假設導致阻抗值要比實際值偏大，可用于電磁兼容的最壞情況估計。

考慮如圖1所示的多芯電纜結構，定義其軸向沿著坐標[x]方向。圖中屏蔽電纜被等效成兩個回路模型，即電纜屏蔽層與參考地之間構成的外部回路以及與內部芯線之間構成的內部回路。忽略線纜內部耦合到屏蔽層外部的能量，當線纜附近的入射場在屏蔽層上激勵起外部電壓源[VEM-lines]與電流源[IEM-lines]時，外部回路可用下面方程進行描述：

[dVsoutdx+ZsoutIsout=VEM-lines] （4）

[dIsoutdx+YsoutVsout=IEM-lines] （5）

式中：[Vsout，][Isout]表示屏蔽層上電壓與電流；[Zsout=][R+jωL]為屏蔽層單位長阻抗；[Ysout=G+jωC]為屏蔽層單位長導納；[R，][L，][G，][C]分別為屏蔽層單位長電阻、電感、電導和電容參數。而內部回路可通過轉移阻抗[Zt]和轉移導納[Yt]與外部回路耦合，因此芯線電壓[Vi]與電流[Ii]滿足如下方程：

式中：[Zt=Rdc+jωLt，Yt=jωCt]表示芯線導體的單位長阻抗矩陣和導納矩陣；[Lt]為轉移電感矩陣；[Ct]為轉移電容矩陣。

根據式（4）～（7）可寫出外、內部回路的時域表達式為：

在實際系統中，電纜大多采用導電細絲組成的編織網構成屏蔽體。對于編織層的屏蔽特性可以用屏蔽體半徑[b]，編織束數[l，]編織束內細線根數[Nl，]細線直徑[d，]細線電導率[σ]以及編織角[?]來描述。該電纜的轉移電阻、轉移電感及轉移電容可用下式計算得到：

式中：[κ]表示投影覆蓋率，定義為[κ=2F-F2]，其中[F=Nldl4πbcos?]為編織網參數填充因子。[K（e），][E（e）]分別表示第一類第二類橢圓積分，[e]為橢圓離心率。[Cout，][Cin]分別為電纜屏蔽層內部電容以及對參考地的外部電容。[εout]和[εin]分別表示電纜外部區域介電常數及內部電介質介電常數。

為了求解方程（8）～（11），通常采用二階的FDTD（2，2）方法，將電壓電流分別在空間與時間上以[Δx2]和[Δt2]為間隔交替離散取點，按中心差分公式將上述方程離散，得：

FDTD（2，4）的穩定性條件為：

這種算法的數值色散關系可以寫為：

從式（27）可以看出，當[Δx=λ2]時，仍然不會造成因為離散而帶來的數值色散。與傳統FDTD方法相比（[Δx≤λ12]），FDTD（2，4）放寬了空間離散間隔所帶來的限制，因而對傳輸線進行粗網格取點能夠有效保證數值精度。

但采用四階中心差分離散得到的迭代方程組不能獲取線纜端點附近的電特性。因此為了保證線纜上電壓電流的連續性，計算[n=1，][N]位置的電流和[n=2，][N]位置處的電壓時，采用式（17）～（20）求解。

3 TDIE?FDTD混合方法求解復雜系統電磁

環境效應

TDIE和FDTD方法混合的關鍵在于電磁場對線纜的耦合。由激勵源的不同，場線耦合理論可分為Taylor模型（磁通量和電通量產生分布電壓源和分布電流源）[20]、Agrawal模型（沿導體的切向入射電場產生分布電壓源）[21]以及Rashidi模型（入射磁場分量產生分布電流源）[22]三種，三種方法在理論上是共通的，只要運用得當都能夠得到同一個解[23]。本文基于電場積分方程求解場，因此采用Agrawal模型求解線纜分布源。

由PWTD方法計算得到系統電磁結構表面電流密度[JEM。]采用如下電場積分方程，獲取線纜位置的場分布：

為了將TDIE場求解器與FDTD線求解器有效結合，需要對線纜位置處的場分布在整個計算時間內進行內插，以獲取FDTD各個不同時間離散點的激勵場。

4 算 例

考慮如圖2所示的穿艙線纜耦合結構，其中線纜為RG?58型同軸電纜，艙內艙外電纜長度均為0.5 m，電纜外屏蔽層對地電阻為100 Ω，芯線與屏蔽層間電阻為50 Ω。矩形理想金屬屏蔽腔尺寸為0.8 m×0.5 m×0.6 m，腔體接地，腔體頂面孔縫中心距左側頂邊0.1 m，左側面三個矩形窄縫長為0.3 m，寬為0.02 m，相鄰孔縫中心距為0.05 m。

為了模擬雷電脈沖對該型電纜的輻射耦合影響，本文采用IEC標準的雙指數脈沖形式的輻射源，具體表達式為：

[E（t）=kE0e-βt-e-αt] （30）

式中：[k=1.3，][E0=50]kV，[α=6.0×108 s-1，][β=][4.0×107 s-1]，波形如圖3（a）所示，其極化角[α=π/4]，入射角[φ=π4，][θ=0。]求解線纜終端[V1]和[V2]端口的電壓響應。

仿真中，先利用兩層PWTD加速的TDIE計算腔體內線纜位置處的場分布。其中艙外的屏蔽線纜將其視為電磁結構與整個腔體一起進行表面剖分，結果如圖3（b）所示。腔體表面被剖分成2 906個三角形，一共生成4 328個未知數。獲取到艙外電纜屏蔽層表面電流后，結合FDTD場線耦合模型，求得內部芯線上的感應電壓、電流。

＼13年21期＼Image＼13t3.tif>

在這個算例中，輻射源耦合到線纜終端負載的路徑有兩條：

（1）入射波直接輻射到腔體外的同軸線，然后通過傳導耦合到內部線路；

（2）入射波通過腔體頂部與側面的孔縫，進入到腔體內部，將能量耦合到線纜負載上。

設定兩種情形進行仿真：

（a）忽略孔縫輻射耦合效應；

（b）綜合考慮傳導與孔縫耦合兩條路徑。

仿真的結果如圖4和圖5所示，對比兩種情形，可以看出內部芯線負載端電壓波形存在振蕩，這是由孔縫存在導致腔體中非均勻場耦合到內部線纜所造成的影響?？偟膩碚f，綜合兩種情形對應的負載端產生的感應電壓波形幾乎重合，說明負載端受外部電磁場直接輻射產生的感應電壓要遠遠大于通過孔縫耦合產生的感應電壓。

5 結 語

本文方法適合分析帶有線纜結構系統的多尺度問題，數值算例說明了方法的有效性。

當然，本文方法有其局限性：

（1）場計算時忽略了線纜感應電流生成的場反過來對已有場的互作用影響，在諧振效應較強的情況下準確性會降低。

（2）線纜計算時采用的是基于傳輸線理論的方法，其本身是一種準靜態方法，忽略了場的縱向分量的影響，在頻率較低時（一般在GHz以下）很準確，但在頻率較高時（線纜束橫截面積與波長可比擬）準確度會降低。

參考文獻

[1] BAUM C E. How to think about EMP interaction [C]// Proceedings of the 1974 Spring FULMEN Meeting. Kirtland AFB， England： [s.n.]， 1974.

[2] BAUM C E， LIU T K， TESCHE F M. On the analysis of general multi?conductor transmission?line networks [J]. Interaction Notes， 1978， 350.

[3] HARRINGTON R F. Filed computation by moment method [M]. Marlabar， FL： Krieger， 1982.

[4] MARTIN H C， CAREY G F. Introduction to finite element analysis： Theory and application [M]. New York： McGraw Hill， 1973.

[5] YEE K. Numerical solution of initial boundary value problems involving Maxwell′s equations in isotropic media [J]. IEEE Transactions on Antennas Propagation， 1966， 14： 302?307.

[6] HOEFER W. The transmission?line?matrix method： theory and applications [J]. IEEE Transactions Microwave Theory Technology， 1985， 33（10）： 882?893.

[7] RAO S M， WILTON D R， GLISSON A W. Electromagnetic scattering by surfaces of arbitrary shape [J]. IEEE Transactions on Antennas propagation， 1982， 30（3）： 409?418

[8] AYGUN K， SHANKER B， ERGIN A， et al. A two?level plane wave time?domain algorithm for fast analysis of EMC/EMI problems [J]. IEEE Transactions on Electromagnetic Compatibility， 2002， 44 （1）： 152?164.

[9] DONEPUDI K C， JIN J M， CHEW W C. A higher order multilevel fast multipole algorithm for scattering from mixed conducting/dielectrics bodies [J]. IEEE Antennas Propagation Magazine， 2003， 51（10）： 2814?2821.

[10] KONDYLIS G D， FLAVIIS F D， POTTIE G J， et al. A memory?efficient formulation of the finite?difference time?domain method for the solution of Maxwell equations [J]. IEEE Transactions on Microwave Theory Technology， 2001， 49（7）： 1310?1320.

[11] JIAO D， ERGIN A， SHANKER B， et al. A fast time?domain higher?order finite?element?boundary?integral method for 3?D electromagnetic scattering analysis [J]. IEEE Transactions on Antennas propagation， 2002， 50（9）： 1192?1202.

[12] TRENKIC V， CHRISTOPOULOS C. Development of a general symmetrical condensed node for the TLM method [J]. IEEE Transactions on Microwave Theory Technology， 1996， 44（12）： 2129?2135.

[13] BA?CI H. Fast and rigorous analysis of EMC/EMI phenomena on electrically large and complex cable?loaded structures [J]. IEEE Transactions on Electromagnetic Compatibility， 2007， 49（2）： 361?381.

[14] XIE H， WANG J， FAN R， et al. Study of loss effect of transmission lines and validity of a SPICE model in electromagnetic topology [J]. Progress in Electromagnetics Research， 2009， 90： 89?103.

[15] XIE H， WANG J， FAN R， et al. A hybrid FDTD?SPICE method for transmission lines excited by a nonuniform incident wave [J]. IEEE Transactions on Electromagnet Compatibility， 2009， 51（3）： 811?817.

[16] LIU Q F， YIN W Y， MAO J F， et al. Accurate characterization of shielding effectiveness of metallic enclosures with thin wires and thin Slots [J]. IEEE Transactions on Electromagnet Compatibility， 2009， 51（2）： 293?300.

[17] LIU Q F， YIN W Y， XUE M F et al. Shielding characterization of metallic enclosures with multiple slots and a thin wire antenna loaded： multiple oblique EMP incidence with arbitrary polarization [J]. IEEE Transactions on Electromagnet Compatibility， 2009， 51（2）： 284?292.

[18] LIU Q F， YIN W Y， LIU P G， et al. Time?domain investigation on cable?induced transient coupling into metallic enclosure [J]. IEEE Transactions on Electromagnet Compatibility， 2009， 51（4）： 953?962.

[19] WANG W， QIN Y J， LIU P G， et al. The solution of transient response for complex transmission line network by hybrid higher order finite difference time domain?modified nodal analysis method [J]. Journal of Electronics and Information Technology， 2012， 34（12）： 2999?3005.

[20] TAYLOR C D. The Response of a Terminated two?wire transmission line excited by a nonuniform electromagnetic Field [J]. IEEE Transactions on Antennas propagation， 1965， 13（6）： 987?989.

[21] AGRAWAL A K. Transient Response ofmulticonductor transmission lines excited by a nonuniform electromagnetic field [J]. IEEE Transactions on Electromagnetic Compatibility， 1980， 22（2）： 119?129.

[22] RASHIDI F. Formulation of field?to?transmission line coupling equations in terms of magnetic excitation field [J]. IEEE Transactions on Electromagnetic Compatibility， 1993， 35（3）： 404?407.

[23] PAUL， C R. Analysis of multiconductor transmission lines [M]. New York ： Wiley， 1994.

[24] WANG W， LIU P G， QIN Y J. An unconditional stable 1D?FDTD method for modeling transmission lines based on precise split?step scheme [J]. Progress in Electromagnetics Research， 2010， 135： 245?260.

[25] WANG W， ZHOU D M， LIU P G， et al. Transient analysis of transmission line with arbitrary loads based on the split?step Crank?Nicolson?FDTD method [J]. Journal of Electronics and Information Technology， 2013， 35（1）： 209?214.

[26] 王為，覃宇建，劉培國，等.基于高階時域有限差分法與改進節點分析法混合求解傳輸線網絡瞬態響應[J].電子與信息學報，2012，34（12）：2999?3005.

[27] 王為，周東明，劉培國，等.端接任意負載傳輸線的分步CN?FDTD分析方法[J].電子與信息學報，2013，35（1）：209?214.