復雜媒質目標電磁散射問題的有限元分析

彭文峰 宛 汀 郁美艷

（南京理工大學，江蘇 南京210094）

引 言

隨著微波和光學工程的迅速變革，以及人工合成媒質技術的不斷發展，復雜電磁媒質的研究受到越來越多的重視.選擇和運用恰當的電磁場問題的數值計算方法對復雜電磁媒質進行研究有著重要的理論意義和工程應用價值.在一些傳統的數值計算方法中，使用積分方程法如矩量法的算法，往往導致復雜的格林函數，有時候格林函數甚至是不可得的.使用時域有限差分法（FDTD）時，因為需要計算額外的時間導數而變得不穩定.用于分析平面分層結構的譜域法與直線法也需要重新計算格林函數.雖然利用場分解特性［1］，可以將手征材料等效成兩個沒有耦合的簡單材料進行計算，但是這種方法卻不能應用到更具有通用形式的一般雙各向異性媒質中.眾所周知，有限元方法［2］是一種通用性很強的數值算法，它能靈活地適應物體幾何結構和材料的變化.與FDTD相比，有限元法更適用于分析任意不均勻復雜結構.而與矩量法相比，使用有限元算法計算雙各向異性媒質的優點在于推導雙各向異性的泛函公式相對比較簡單，有較強的復雜媒質的處理能力.而且，有限元方法生成的線性方程組的系數陣為稀疏矩陣，節省了內存.因此，本文研究將有限元方法發展應用到含有復雜媒質結構的電磁散射問題中，包括各向異性媒質、雙各向同/異性媒質，并驗證其準確性和通用性.

本文建立了具有一般通用形式的線性媒質——雙各向異性媒質［3］的有限元泛函公式模型，采用完全匹配層［4］來截斷散射問題開放區域，并推導了有限元矩陣方程的具體表達式.然后，運用該方法對四種典型的復雜媒質結構散射問題進行了數值計算，并將計算結果和相關文獻進行比較.結果證明，有限元方法能很好地解決含有各向異性、雙各向同性以及雙各向異性的各種復雜媒質目標電磁散射的數值仿真問題.本文的研究工作具有較強的通用性，能夠對多種復雜媒質目標進行分析，在研究雷達目標隱身和反隱身技術、復雜天線系統設計、現代電子系統的電磁兼容性分析等領域均能有效地發揮作用.

1 基本理論

對于復雜媒質，其本構關系為

圖1給出了一個典型的采用有限元方法分析電磁散射問題的示意圖.我們采用完全匹配層（PML）作為區域截斷邊界條件.在PML區域內，有效磁導率和介電常數為對角矩陣的形式.對于z方向的PML，強加如下約束條件

圖1 三維目標的散射

假定波的傳播方向為z，當bc＝1，a＝b時，可以得到0反射.為了使入射波充分的衰減，a、b、c應取作復數.令a＝b＝1/c＝α-jβ，α決定了此介質中的波長，β決定了波的衰減程度.則PML內的相對磁導率與電導率可以表示為

PML參量的取值依賴于波的傳播方向和PML層本身的位置（PML層可以設在面、邊或者角落）.針對二維和三維問題時，我們選擇不同的ε＝r和＝μr應用于不同的方向.棱邊和體角落含有特殊的張量，例如放于棱邊，則對于z方向上的一條棱邊，應選擇r和r的值為r＝r＝［μr］x×［μr］y；放在體角上，則為三個方向上的各張量的乘積.

由復雜媒質的本構關系以及麥克斯韋方程組推導出雙各向異性媒質的電場波動Helmholtz方程為

對于復雜媒質，有限元泛函可以寫為［2］

由于式（6）的工作變量是總場，在散射問題中有總場E＝Esc＋Einc，其中Esc為散射場，Einc為入射場，因此式（6）也可以寫為

又因Einc為已知場，在變分公式求偏導時等于零，因此式（7）中只含Einc的項可以去掉，應用第一矢量格林定理，可得

式中：Vsc表示散射體的體積；Ssc表示散射體的表面；Einc為激勵平面波.

在獲得整個分析區域的泛函之后，接下來要進行的是區域離散的工作.選用靈活的四面體單元對整個區域進行離散，生成網格之后，根據選定的切向連續的矢量基函數（9）來建立有限元方程.

式中：i＝1，2，…，6，i表示第i條邊；i1i2表示第i條邊兩個端點的編號；Li為體積坐標.

選擇相同的基函數和加權函數Ni，可以得到復雜媒質的有限元散射公式為

根據伽遼金方法，對于其單元每一棱邊元的殘數加權積分為0，結合式（10），得到整個e單元的矩陣形式為

式中：

其中，

式中，

由方程（11）解得散射場Esc（近場）之后，在遠處r的散射場Esc（r）（遠場）就可以使用等效原理求解，即

式中：S是包圍散射體的任意閉合曲面；Er為等效的無窮遠處入射的平面波，其表達式為

其中α為極化角，當α＝0時入射波為θ極化，α＝π/2時入射波為φ 極化.kinc是傳播矢量，θinc、φinc為激勵平面波的入射角，有

式中：Ni為基函數；ai為有限元方程已求出的基函數的系數.

將式（17）代入式（15）得到

求出了遠處r的散射場，就可以用雷達散射截面（RCS）表示物體的散射特性，雷達散射截面的定義為

由于入射場是平面波，因此有｜Einc（r）｜＝1，則式（19）可以寫為

2 數值算例

下面計算幾個典型的例子來證明本文方法的正確性和有效性.本文均采用PML作為區域截斷邊界條件，為保證PML良好的截斷效果，算例中PML的厚度均取為0.25λ0，距離散射體為0.3λ0，λ0為自由空間波長，參數值取為α＝β＝1.5.

第一個例子分析的是電尺寸為k0a＝0.2π的鐵氧體球.直角坐標系的原點位于球心，設外加偏置磁場在方向得到，即H0＝H0，入射的平面波沿z軸方向，其極化方向沿正x軸方向.各向異性鐵氧體的介質參數為ε0（真空中的介電常數），其相應的導磁率為［5-6］，可表示為

本例中，鐵氧體的參數為μ1＝5μ0，μ2＝jμ0，μ3＝7μ0.圖2為本文計算結果與文獻［7］的比較，可以看出曲線有細小差別，其原因是本程序采用PML作為邊界截斷條件，一般情況下PML并不百分之百吸收，并且-35dB已經是一個很小的值了，所以該誤差在允許的范圍之內.

第二個例子是電尺寸為k0a＝0.5的等離子球，入射的平面波沿z軸正方向入射，其極化方向沿正x軸方向.各向異性等離子體的導磁率為μ0（真空中的導磁率），電容率可表示為

圖3為本文結果與文獻［8］的比較.可以看出，本程序計算的結果和文獻結果在-40dB之上還是非常吻合的.由于本程序采用PML作為邊界截斷條件，一般來講，PML是不能達到百分百吸收的，另外，-40dB是可以忽略的十分小的誤差.

第三個例子是尺寸為k0a＝1.5的手征介質球.為了與文獻對比，手征的本構關系可以寫成如下形式

根據文獻［3］，手征參數設為εDBF＝4ε0、μDBF＝μ0、ε＝εDBF/（1-k2r）、μ＝μDBF/（1-k2r）.入射的平面波沿z軸正方向，其極化方向沿正x軸方向.圖4中的曲線對比再次驗證了本程序的正確性.

第四個例子是雙各向異性媒質圓柱，底面半徑為a＝0.5λ0，高度為h＝0.2λ0，雙各向異性媒質參數具有如下張量形式入射的平面波沿z軸正方向，其極化方向沿正x軸方向.Ω 媒質參數為ε1＝2.0，ε2＝3.0，ε3＝2.0，μ1＝1.2，μ2＝1.2，μ3＝1.0，Ω＝0.0，0.5，1.0.圖5中與文獻［9］曲線的對比驗證了本程序計算雙各向異性媒質的正確性.

3 結 論

本文介紹了分析復雜媒質散射問題的有限元方法.有限元法適用于分析任意非均勻復雜結構，有較強的復雜媒質的處理能力，生成的線性方程組的系數陣為稀疏矩陣，節省了存儲量.算例的數值分析表明，該方法的計算結果與文獻結果吻合較好，從而證明了其正確性和有效性.此外，該方法適用范圍廣，不僅能夠用來計算各向異性媒質的散射特性，對雙各向同性和雙各向異性媒質同樣適用，是一種通用性較強的方法.

［1］AL-KANHAL M A，ARVAS E.Electromagnetic scattering from a chiral cylinder of arbitrary cross section［J］.IEEE Trans Antennas and Propagation，1996，44（7）：1041-1048.

［2］金建銘，著.電磁場有限元方法［M］.王建國，葛德彪，譯.西安：西安電子科技大學出版社，1998.

［3］周 平.脊加載橢圓波導傳輸特性的矢量有限元法分析［J］.電波科學學報，2009，24（6）：1164-1167.ZHOU Ping.Propagation characteristics analysis for ridged elliptical waveguide by vector finite-element method［J］.Chinese Journal of Radio Science，2009，24（6）：1164-1167.（in Chinese）

［4］孫向陽，聶在平，李愛勇，等.用高階疊層矢量有限元法計算隨鉆測井的三維電磁響應［J］.電波科學學報，2009，24（2）：273-279.SUN Xiangyang，NIE Zaiping，LI Aiyong，et al.The modeling of logging-while-drilling tool’s three-dimensional electromagnetic response using the high order hierarchical vector finite element method［J］.Chinese Journal of Radio Science，2009，24（2）：273-279.（in Chinese）

［5］EDWARD K N，YUNG H Y，DING R S，et al.Finite-difference analysis of H-Plane waveguide Y-junction circulars［J］.Microwave and Optical Technology Letters，1999，20（6）：414-422.

［6］耿友林，吳信寶，官伯然.導體球涂覆各向異性鐵氧體介質電磁散射的解析解［J］.電子與信息學報，2006，28（9）：1740-1743.GENG Youlin，WU Xinbao，GUAN Boran.The analytical solution to the electromagnetic scattering by an anisotropic ferrite-coated conducting sphere［J］.Journal of Electronics ＆Information Technology，2006，28（9）：1740-1743.（in Chinese）

［7］耿友林.球矢量波函數在各向異性介質散射中的應用［D］.西安：西安電子科技大學，2006.GENG Youlin.The application of the sphere vector wave function for the electro-magnetic scattering from the anisotropic media［D］.Xi’an：Xidian University，2006.（in Chinese）

［8］ZHANG Y，WEI X，LI E.Electromagnetic scattering from three-dimensional bi-anisotropic objects using hybrid finite element-boundary integral method［J］.J of Electromagn Waves and App，2004，18（11）：1549-1563

［9］OZGUN O，KUZUOGLU M.A non-iterative domain decomposition method for finite element analysis of 3D electromagnetic scattering problems［C］//IEEE Antennas and Propagation Society International Symposium.San Diego，July，5-11，2008：40-50.

［10］ZHU J，CHEN R S，FAN Z H，et al.Efficient preconditioner for the finite-element boundary integral analysis of electromagnetic scattering in the half space［J］.IET Microwaves Antennas ＆ Propagation，2010，4（3）：374-380.

［11］JACKSON S A，VOUVAKIS M N.Octree-based finite element method for electromagnetic scattering problems［C］//IEEE Antennas and Propagation Society International Symposium.Toronto，July，11-17 2010：1-4.

［12］ZHANG Ruigang，CHEN Huanzhen.An immersed finite element method for anisotropic flow models in porous medium［C］//International Conference on Information Science and Technology.Nanjing，March 26-28，2011：168-175.