復(fù)雜媒質(zhì)目標(biāo)電磁散射問(wèn)題的有限元分析

彭文峰 宛 汀 郁美艷

（南京理工大學(xué)，江蘇 南京210094）

引 言

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

本文建立了具有一般通用形式的線性媒質(zhì)——雙各向異性媒質(zhì)［3］的有限元泛函公式模型，采用完全匹配層［4］來(lái)截?cái)嗌⑸鋯?wèn)題開(kāi)放區(qū)域，并推導(dǎo)了有限元矩陣方程的具體表達(dá)式.然后，運(yùn)用該方法對(duì)四種典型的復(fù)雜媒質(zhì)結(jié)構(gòu)散射問(wèn)題進(jìn)行了數(shù)值計(jì)算，并將計(jì)算結(jié)果和相關(guān)文獻(xiàn)進(jìn)行比較.結(jié)果證明，有限元方法能很好地解決含有各向異性、雙各向同性以及雙各向異性的各種復(fù)雜媒質(zhì)目標(biāo)電磁散射的數(shù)值仿真問(wèn)題.本文的研究工作具有較強(qiáng)的通用性，能夠?qū)Χ喾N復(fù)雜媒質(zhì)目標(biāo)進(jìn)行分析，在研究雷達(dá)目標(biāo)隱身和反隱身技術(shù)、復(fù)雜天線系統(tǒng)設(shè)計(jì)、現(xiàn)代電子系統(tǒng)的電磁兼容性分析等領(lǐng)域均能有效地發(fā)揮作用.

1 基本理論

對(duì)于復(fù)雜媒質(zhì)，其本構(gòu)關(guān)系為

圖1給出了一個(gè)典型的采用有限元方法分析電磁散射問(wèn)題的示意圖.我們采用完全匹配層（PML）作為區(qū)域截?cái)噙吔鐥l件.在PML區(qū)域內(nèi)，有效磁導(dǎo)率和介電常數(shù)為對(duì)角矩陣的形式.對(duì)于z方向的PML，強(qiáng)加如下約束條件

圖1 三維目標(biāo)的散射

假定波的傳播方向?yàn)閦，當(dāng)bc＝1，a＝b時(shí)，可以得到0反射.為了使入射波充分的衰減，a、b、c應(yīng)取作復(fù)數(shù).令a＝b＝1/c＝α-jβ，α決定了此介質(zhì)中的波長(zhǎng)，β決定了波的衰減程度.則PML內(nèi)的相對(duì)磁導(dǎo)率與電導(dǎo)率可以表示為

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

由復(fù)雜媒質(zhì)的本構(gòu)關(guān)系以及麥克斯韋方程組推導(dǎo)出雙各向異性媒質(zhì)的電場(chǎng)波動(dòng)Helmholtz方程為

對(duì)于復(fù)雜媒質(zhì)，有限元泛函可以寫為［2］

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

又因Einc為已知場(chǎng)，在變分公式求偏導(dǎo)時(shí)等于零，因此式（7）中只含Einc的項(xiàng)可以去掉，應(yīng)用第一矢量格林定理，可得

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

在獲得整個(gè)分析區(qū)域的泛函之后，接下來(lái)要進(jìn)行的是區(qū)域離散的工作.選用靈活的四面體單元對(duì)整個(gè)區(qū)域進(jìn)行離散，生成網(wǎng)格之后，根據(jù)選定的切向連續(xù)的矢量基函數(shù)（9）來(lái)建立有限元方程.

式中：i＝1，2，…，6，i表示第i條邊；i1i2表示第i條邊兩個(gè)端點(diǎn)的編號(hào)；Li為體積坐標(biāo).

選擇相同的基函數(shù)和加權(quán)函數(shù)Ni，可以得到復(fù)雜媒質(zhì)的有限元散射公式為

根據(jù)伽遼金方法，對(duì)于其單元每一棱邊元的殘數(shù)加權(quán)積分為0，結(jié)合式（10），得到整個(gè)e單元的矩陣形式為

式中：

其中，

式中，

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

式中：S是包圍散射體的任意閉合曲面；Er為等效的無(wú)窮遠(yuǎn)處入射的平面波，其表達(dá)式為

其中α為極化角，當(dāng)α＝0時(shí)入射波為θ極化，α＝π/2時(shí)入射波為φ 極化.kinc是傳播矢量，θinc、φinc為激勵(lì)平面波的入射角，有

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

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

求出了遠(yuǎn)處r的散射場(chǎng)，就可以用雷達(dá)散射截面（RCS）表示物體的散射特性，雷達(dá)散射截面的定義為

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

2 數(shù)值算例

下面計(jì)算幾個(gè)典型的例子來(lái)證明本文方法的正確性和有效性.本文均采用PML作為區(qū)域截?cái)噙吔鐥l件，為保證PML良好的截?cái)嘈Ч?，算例中PML的厚度均取為0.25λ0，距離散射體為0.3λ0，λ0為自由空間波長(zhǎng)，參數(shù)值取為α＝β＝1.5.

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

本例中，鐵氧體的參數(shù)為μ1＝5μ0，μ2＝j(luò)μ0，μ3＝7μ0.圖2為本文計(jì)算結(jié)果與文獻(xiàn)［7］的比較，可以看出曲線有細(xì)小差別，其原因是本程序采用PML作為邊界截?cái)鄺l件，一般情況下PML并不百分之百吸收，并且-35dB已經(jīng)是一個(gè)很小的值了，所以該誤差在允許的范圍之內(nèi).

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

圖3為本文結(jié)果與文獻(xiàn)［8］的比較.可以看出，本程序計(jì)算的結(jié)果和文獻(xiàn)結(jié)果在-40dB之上還是非常吻合的.由于本程序采用PML作為邊界截?cái)鄺l件，一般來(lái)講，PML是不能達(dá)到百分百吸收的，另外，-40dB是可以忽略的十分小的誤差.

第三個(gè)例子是尺寸為k0a＝1.5的手征介質(zhì)球.為了與文獻(xiàn)對(duì)比，手征的本構(gòu)關(guān)系可以寫成如下形式

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

第四個(gè)例子是雙各向異性媒質(zhì)圓柱，底面半徑為a＝0.5λ0，高度為h＝0.2λ0，雙各向異性媒質(zhì)參數(shù)具有如下張量形式入射的平面波沿z軸正方向，其極化方向沿正x軸方向.Ω 媒質(zhì)參數(shù)為ε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中與文獻(xiàn)［9］曲線的對(duì)比驗(yàn)證了本程序計(jì)算雙各向異性媒質(zhì)的正確性.

3 結(jié) 論

本文介紹了分析復(fù)雜媒質(zhì)散射問(wèn)題的有限元方法.有限元法適用于分析任意非均勻復(fù)雜結(jié)構(gòu)，有較強(qiáng)的復(fù)雜媒質(zhì)的處理能力，生成的線性方程組的系數(shù)陣為稀疏矩陣，節(jié)省了存儲(chǔ)量.算例的數(shù)值分析表明，該方法的計(jì)算結(jié)果與文獻(xiàn)結(jié)果吻合較好，從而證明了其正確性和有效性.此外，該方法適用范圍廣，不僅能夠用來(lái)計(jì)算各向異性媒質(zhì)的散射特性，對(duì)雙各向同性和雙各向異性媒質(zhì)同樣適用，是一種通用性較強(qiáng)的方法.

［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］金建銘，著.電磁場(chǎng)有限元方法［M］.王建國(guó)，葛德彪，譯.西安：西安電子科技大學(xué)出版社，1998.

［3］周 平.脊加載橢圓波導(dǎo)傳輸特性的矢量有限元法分析［J］.電波科學(xué)學(xué)報(bào)，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］孫向陽(yáng)，聶在平，李愛(ài)勇，等.用高階疊層矢量有限元法計(jì)算隨鉆測(cè)井的三維電磁響應(yīng)［J］.電波科學(xué)學(xué)報(bào)，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］耿友林，吳信寶，官伯然.導(dǎo)體球涂覆各向異性鐵氧體介質(zhì)電磁散射的解析解［J］.電子與信息學(xué)報(bào)，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］耿友林.球矢量波函數(shù)在各向異性介質(zhì)散射中的應(yīng)用［D］.西安：西安電子科技大學(xué)，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，F(xiàn)AN 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.