時域有限元法在三維電磁輻射和散射問題中的應(yīng)用

時域有限元法在三維電磁輻射和散射問題中的應(yīng)用

時域有限元法近年來，通信系統(tǒng)和雷達系統(tǒng)的工作頻率范圍有所增加，這要求天線和相關(guān)微波設(shè)備具有頻帶特性。超寬帶移動通信設(shè)備基本工作在微波波段,內(nèi)部幾何結(jié)構(gòu)和材料分布十分復(fù)雜。很多先進設(shè)計還需要外加電子材料,如人工電磁帶隙材料、左手材料等,以期進一步提高其性能。如此一來,對這一類具有復(fù)雜結(jié)構(gòu)和非均勻介質(zhì)的超寬帶天線、超寬帶雷達的瞬態(tài)輻射、目標的瞬態(tài)散射問題進行建模與仿真提出了新要求。十多年來,有大量的文獻采用不同的時域分析計算方法,例如時域有限差分法(FDTD)、時域積分法(TDIE)和時域有限元法(TDFEM)等,對各種不同的結(jié)構(gòu)進行了理論分析。FDTD法因其理論簡單、方法通用而成為流行的數(shù)值模擬工具。然而,它不太適用于復(fù)雜幾何結(jié)構(gòu)和介電材料分布。TDIE法在分析無界均勻介質(zhì)中的金屬或均勻介質(zhì)目標的散射和輻射時具有優(yōu)勢,原因有兩方面:一是它把未知量限制在目標表面積上;二是,它通過格林函數(shù)自然地滿足輻射條件。然而,當數(shù)值模擬的區(qū)域包含復(fù)雜非均勻介質(zhì)時,TDIE法就會面臨非均勻介質(zhì)體內(nèi)帶來的大型矩陣皆為滿陣而非稀疏矩陣,從而需要消耗大量的計算資源的困難。相反地,在眾多的方法中,TDFEM因其繼承了FEM特別適用于復(fù)雜的幾何結(jié)構(gòu)和介電特性分布的優(yōu)點及在非均勻介質(zhì)體內(nèi)產(chǎn)生的也是稀疏矩陣的長處,因而在解決寬帶特性的電磁場問題上彌補了FDTD和TDIE法的某些不足而被人們漸漸關(guān)注。時域有限元法起源于上世紀80年代,最初是由A.C.Cangellaris等人提出的微分形式的點匹配時域有限元法。該方法結(jié)合了FDTD顯式積分法的簡明性和傳統(tǒng)有限元的靈活性,但是由于采用的是節(jié)點元,在強加邊界條件時存在困難;又因同時采用電場和磁場,變量增加了一倍;此外還會出現(xiàn)源于電場和磁場點在網(wǎng)格上的交替出現(xiàn)而產(chǎn)生的類似FDTD法的越級(leap-frog)問題。為了避免或解決這些問題,1990年代提出無條件穩(wěn)定的隱式時域有限元法[11,12,13,14,15,16],該方法基于電場(或磁場)的矢量有限元空間基函數(shù)和時間基函數(shù)展開的二階矢量波方程,其優(yōu)點是簡單而穩(wěn)定,缺點是每個時間步必須求解一個矩陣方程;對于開域空間問題,采用傳統(tǒng)的吸收邊界條件,有限元計算空間十分龐大,這個問題更加嚴重,以致無法進行。由于這些原因,這類時域有限元法發(fā)展緩慢。直到2000年代初,由于完全匹配層(PML)和正交基函數(shù)概念的提出和應(yīng)用,時域有限元法才得到較大的發(fā)展。目前時域有限元法在國外,特別是以J.M.Jin的科研組在PML的應(yīng)用、正交基函數(shù)的引進及其在二維和三維電磁輻射問題的應(yīng)用方面,作了若干開創(chuàng)性的工作;在國內(nèi)除外,僅見二維時域有限元法關(guān)于穩(wěn)定性分析的報道。以下就時域有限元法應(yīng)用于三維電磁輻射和散射問題的原理、發(fā)展中遇到的挑戰(zhàn)和解決方案、尚待解決的問題和發(fā)展趨勢作簡要介紹,以期推動它在我國的應(yīng)用和發(fā)展。全文共分5個部分,本節(jié)是引言;第二部分給出基本原理和公式;第三部分論述時域有限元法在應(yīng)用于三維電磁問題時遇到的主要問題以及相應(yīng)的解決方法;第四部分討論時域有限元法在該領(lǐng)域的應(yīng)用和發(fā)展的前景;最后一節(jié)為結(jié)論。1理論發(fā)展和公式1.1單元狀態(tài)方程的匹配和點匹配時域有限元法起源于上世紀80年代,由A.C.Cangellaris等人首次提出,它將FDTD法作有限元插值擴展得到了節(jié)點時域有限元算法,稱為點匹配FEM。很明顯,這種方法的提出和發(fā)展的源動力是希望結(jié)合FDTD顯式積分法的簡明性和傳統(tǒng)有限元空間離散過程的靈活性,其原理和公式如下:與FDTD方法相似,點匹配時域有限元法計算區(qū)域被2個互補的網(wǎng)格(分別是電場和磁場)離散化,網(wǎng)格的構(gòu)造規(guī)則是使每個磁場網(wǎng)格的一個單元包含另一個單元的一個節(jié)點,反之亦然。按照有限元法的過程,電場E和磁場H在計算區(qū)域上可以被近似表達為:E(r,t)=Μ∑i=1?i(r)Ei(t),Η(r,t)=Ν∑j=1?j(r)Ηj(t)?(1)這里M和N分別是電場和磁場節(jié)點的數(shù)目,Ei和Hj分別是電場和磁場在節(jié)點上的值,?i,?j是已知的有限元插值函數(shù)。利用(1)式可以表達在單元內(nèi)部的任意一點的場變量值。然后,將Maxwell旋度方程化為下面的狀態(tài)方程組:μdΗjdt=-L∑l=1?φl(rj)×El(t)j=1,2,?Ν?σEi+εdEidt=L∑l=1?φl(ri)×Ηl(t)i=1,2,?Μ?(2)其中φl(r)是有關(guān)單元節(jié)點l的單元形狀函數(shù),L是單元的節(jié)點數(shù)目。?φl(rj)表示在r=rj處?φl(r)的值。將(2)式中的時間微分用差分代替,就得到微分方程形式的點匹配TDFEM計算公式。值得注意的是此方法主要缺陷在于:每個場的三個分量都位于相同的節(jié)點上,這就導(dǎo)致了在媒質(zhì)的發(fā)生電磁特性突變的邊界表面上強加合適的邊條件變得十分復(fù)雜。1.2whit充換電子束考慮到點匹配有限元法的缺陷,1990年之后有人提出另一途徑得到的TDFEM,此類TDFEM是離散化二階矢量波方程(也即旋度-旋度方程),從Maxwell’s方程組中消去一個場變量而得到的[11,12,13,14,15,16],其中典型的是StephenD.Gedney等人提出的一種無條件穩(wěn)定的隱式FETD方法,其原理和公式簡述如下:時變電場必須滿足依賴時間的非齊次波動方程:?×μ-1r?×E+μ0σ?E?t+εrc20?2E?t2=-μ0?J?t(3)式中,c0是自由空間光速;σ是電導(dǎo)率;εr和μr是相對介電常數(shù)和磁導(dǎo)率。J是有限域Ω的表面?Ω上的電流密度。所求問題可以簡化成把?Ω界定為是一個理想導(dǎo)電面(PEC)或理想導(dǎo)磁面(PMC)的電磁場邊值問題。作(4)式與測試函數(shù)T的內(nèi)積,應(yīng)用Green第一恒等式可以導(dǎo)出(3)式的弱形式?Ω[μ-1r(?×Τ)(?×E)+μ0σΤ?E?t+εrc20Τ?2E?t2]dΩ=-?Ωμ0Τ?J?tdΩ(4)將矢量場用Whitney邊棱元基函數(shù)作權(quán)重函數(shù)Wi展開,展開系數(shù)ej是時間連續(xù)函數(shù),同時,測試函數(shù)也取為Whitney基。代回(4)式,并由此導(dǎo)出二階微分方程[Τε]1c20d2edt2+[Τσ]η0c0dedt+[S]e=-f(5)這里[Tε]、[Tσ]和[S]是不依賴時間的矩陣:[Τε,σ]i,j=?Ω{εrσ}Wi?WjdΩ(6)[S]i,j=?Ω1μr?×Wi??×WjdΩ(7)fi=?Ωμ0Wi??J?tdΩ(8)將(5)式中的時間微分用Newmark-β差分公式展開,可以得到積分形式的隱式TDFEM迭代公式,據(jù)此,由第n-1和第n個時間步的場值,可計算出第n+1個時間步的場值。2積分形式的隱式時域有限元法上一節(jié)中闡述的第一類直接求解Maxwell方程組的顯式時域有限元法,優(yōu)點是不需要在每個時間步上求解一個矩陣方程;但是它一般同時在電場和磁場上工作,變量增加了一倍,而且需要不同種類的基函數(shù)來展開它們;另外,由于電場和磁場點在網(wǎng)格上的交替出現(xiàn),會導(dǎo)致類似FDTD法的越級(leap-frog)問題。由于這些原因,這一方法未得到推廣應(yīng)用,人們轉(zhuǎn)而尋求另外的解決方案——積分形式的隱式時域有限元法。積分形式的隱式時域有限元法的優(yōu)點是簡單而且穩(wěn)定,其缺點是只適用于理想導(dǎo)電面或理想導(dǎo)磁面;對于無界空間的開域問題,采用傳統(tǒng)的吸收邊界條件,會使計算區(qū)域變得十分龐大;同時,它在每個時間步上需要求解一個空間矩陣方程。由于上述原因,加之當時計算機速度和內(nèi)存的限制,時域有限元法長期不能廣泛流行。這種情況直到完全匹配層(PML)的出現(xiàn)和應(yīng)用才得到改變。下面分兩節(jié)闡述時域有限元法是如何解決上述問題的,首先闡述無界自由空間的截斷問題——完全匹配層法,然后闡述龐大計算量問題的兩個解決方案——正交基函數(shù)和子域分解并行計算。2.1時域有限元法微波輻射和散射問題的重要特點是開域,因此在用時域有限元法來分析和計算電磁場時,所面臨的一個實際問題就是如何確定網(wǎng)格劃分的邊界和如何處理場在這些邊界上的值。不像FDTD法在吸收邊界條件(ABC)上有很好的發(fā)展,很長一段時間里,很少有對時域有限元法的精確ABCs進行研究的,也僅僅實現(xiàn)了一階和二階ABCs。直到把完全匹配層(perfectlymatchedlayer:PML)的概念引入到時域有限元法,這一問題才得到很好的解決?？梢哉f,Berenger在1994年提出的完全匹配層把這項研究推向了新階段。PML擺脫了以往吸收邊界條件研究中拘泥于實際可實現(xiàn)的物理吸收邊界的框框,提出了一種基于各向異性介質(zhì)的數(shù)學(xué)邊界。這種吸收邊界的性能比以往提出的邊界優(yōu)越,而且為進一步研究能與不同介質(zhì)相匹配的完全匹配層提供了新的思路。1995年,ZacharySacks等提出了一種易于實現(xiàn)和調(diào)控的PML。這種各向異性介質(zhì)的性能參數(shù)[ε]和[μ]為復(fù)對角張量,為了使PML層的內(nèi)表面成為一個完全的無反射界面,對角張量應(yīng)具有以下形式:[μ]=μ0[Λ][ε]=ε0[Λ](9)Λ(r,ω)=diag{γyγzγx,γzγxγy,γxγyγz},γξ=1+σξ(ξ)jωε0ξ=x,y,z(10)這樣電場在PML中滿足?×μ-1(Λ-1??×E)-ω2εΛ?E=0(11)(11)式看似簡單,但是由于Λ是張量,而且是頻率ω的函數(shù),在進行時域計算時,必須利用Laplace或Fourier變換將(11)式轉(zhuǎn)化為時域下的方程,數(shù)學(xué)實現(xiàn)十分繁瑣。因此,雖然PML是在十多年前提出的,并已成為FDTD通常用的截斷方法,并已應(yīng)用于頻域有限元法,但其實現(xiàn)于時域有限元法在2002年之后才有報道。為了得到(11)式對應(yīng)的時域離散矩陣方程,我們采用Laplace變換及其反變換。同時,將電場根據(jù)空間基函數(shù)Nj展開,并利用測試基函數(shù)Ni得到時域版本的半離散的矩陣方程:[Τ]d2dt2{e}+[B]ddt{e}+[S]{e}-{h}+{g}=0(12)Τij=?VεΝi?ΝjdV(13)Bij=?Vεε0Νi?J?ΝjdV(14)Sij=?V1μ(?×Νi)?(?×Νj)dV+?Vεε20Νi?Κ1?ΝjdV(15)hi=?Vεε30Νi?∑j(Κ2?uj?Νj)dV(16)gi=?V1με0(?×Νi)?∑j[(L1?u+j+L2?u++j)?(?×Νj)]dV(17)其中J,K1,K2,L1,L2都是對角矩陣,矩陣元是σξ(ξ=x,y,z)的函數(shù);并且uj,u+j,u++j是含有卷積的對角矩陣diage-σξε0t*ej(t),e-σξε0t*ej(t),e-σξε0t*ej(t)(18)其中σ的下標ξ代表{x,y,z}的輪換對稱表達,在uj中分別對應(yīng)于{x,y,z};uj+中分別對應(yīng)于{y,z,x};而在uj++中分別對應(yīng)于{z,x,y}。在利用PML和TDFEM求解電磁輻射和散射問題時,首先在輻射體或散射體的外側(cè)設(shè)置封閉的PML層,在PML層內(nèi)部,求解由(12)式所代表的TDFEM矩陣方程,在PML層內(nèi)表面所包圍包括輻射體和散射體內(nèi)的空間,則求解由(5)所代表的TDFEM矩陣方程。由于PML層很薄(約n分之一波長)而且可以離輻射體或散射體很近,因此與ABCs相比極大的壓縮了計算空間。2.2計算域和大型天線陣列盡管TDFEM應(yīng)用了PML算法后大大縮小了計算區(qū)域,然而,對于模擬電大尺寸天線和大型天線陣列,計算域和未知場變量個數(shù)十分龐大的問題依然存在;另外在每個時間步上都要求解一個階數(shù)極高的矩陣方程,這大大增加了計算的復(fù)雜度。為解決計算內(nèi)存和時間問題人們進一步采取了各種辦法,比較典型的有引入正交基函數(shù)和子域分解技術(shù):2.2.1第i號邊的面基函數(shù)有學(xué)者提出構(gòu)造正交基函數(shù)的辦法來解決在每個時間步上都要求解一個矩陣方程問題,通過正交基函數(shù)可以使基于Newmark-β的時域有限元法的隱式計算公式中的系數(shù)矩陣成為對角陣或高度稀疏的矩陣,從而迅速縮小每個時間步的計算時間。下面簡要介紹它的構(gòu)造原理:考慮一個矢量波動方程,用Whitney基函數(shù)對電場和磁場進行逼近。對于二維問題,原始的線性Whitney邊基函數(shù)定義為Wi=ξj?ξk-ξk?ξj(19)這里的ξk是與節(jié)點k有關(guān)的線性節(jié)點基函數(shù),且第i號邊連接節(jié)點j和k,這些函數(shù)自然滿足單元邊界的切向場的連續(xù)性,對法向分量的不連續(xù)性沒有影響。對于三維問題,面基函數(shù)一般定義為:Ui=2(ξj?ξk×?ξm+ξk?ξm×?ξj+ξm?ξj×?ξk)(20)其中,編號為i的面由節(jié)點j,k和m連接而成,ξj,ξk,ξm表示體積坐標。Ui的法向分量在單元邊界是連續(xù)的,并且它的通量僅在穿過編號為i的面時為1,其他為0。顯然這兩組邊基函數(shù)都不是正交的。然后對已有的基函數(shù)進行擴充,可構(gòu)造出新的一組正交基函數(shù)。首先我們要定義新的內(nèi)積形式:?u,v?=∑i=13αiu(mi)?v(mi)(21)這里mi代表第i號邊的中點,ml代表編號為i的面的中心點。對于二維的正交基函數(shù),要再定義3個輔助的基函數(shù):Bi=ξjξkn^ii=1,2,3(22)ξi是第i號邊的單位法向矢量。之所以構(gòu)造這3個輔助基函數(shù)是因為它們具有以下性質(zhì):1)沿剖分的棱邊的切向分量為零;2)Bi僅在第i號邊上有非零法向分量;3)根據(jù)(21)式的內(nèi)積定義,這些基函數(shù)是兩兩正交的。那么我們可以定義新的一組基函數(shù)Z:Ζi=Wi-∑j=13?Wi,Bj??Bi,Bj?Bj.(23)由它的定義可知,新的邊基函數(shù)Zi是兩兩正交的。同樣的思路,文獻中構(gòu)造了一組三維正交完備的基函數(shù)。2.2.2子域和局部時間步另一種有效解決TDFEM龐大計算域的方法是子域分解技術(shù),即將原始的計算域分成幾個較小的子域并采用并行計算。每個子域問題包含了縮小規(guī)模的計算域,可以用稀疏矩陣求解器直接求解,并且與原來的單一區(qū)域問題相比,它減少了整體的計算復(fù)雜性和計算時間。文獻中提出了一種新型的子域分解方案。以一個一般的計算域V為例,以金屬表面S1和阻抗表面S2為邊界。通過任意的人工邊界SD將V分裂成V1和V2兩個子域。子域上的電場和磁場值用Ei和Hi表示,i=1,2。在子域上二階矢量波方程可以寫成如下的形式:?×(μr-1??×Ei)+1c02εr?2Ei?t2+μ0σ?Ei?t=-μ0?Ji?t(24)?×(εr-1??×Ei)+1c02μr?2Ηi?t2+μεrσ?Ηi?t=?×(εr-1Ji)(25)S1上的邊界條件是n^×E=0,n^×(?×Η)=0(26)S2上的邊界條件是n^×(μ-1?×E)+γe??t(n^×n^×E)=0n^×(ε?×Η)+γh??t(n^×n^×Η)=0(27)如果取γe=ε/μ和γh=μ/ε就是ABCs,對(24)(25)式以矢量基函數(shù)Ni測試得到弱形式的波方程,強加邊條件并離散化,就可以得到矩陣形式的方程并求解。該方法是基于每個子域的雙場二階矢量波方程,對同一空間網(wǎng)格卻不同的時間域上求解電場和磁場。每個子域問題可以小到足以靠一個稀疏矩陣求解器直接求解。這樣,通過對每個子域的預(yù)先分解,并在時間步進之前將信息存儲在內(nèi)存。在每個時間步,子域的問題靠高效利用局部的預(yù)先分解矩陣解決。它的獨特之處在于并不需要解全局的矩陣方程,缺陷是由于該方法將原來計算域分解成了幾個單獨的子域,而各個子域問題由于離散誤差,將不再與原有的問題精確匹配,導(dǎo)致數(shù)值色散誤差的增加。3正演中間層存在的問題如上節(jié)所述,時域有限元法在PML的應(yīng)用、正交基函數(shù)的引進、子域分解技術(shù)等方面取得了一些進展,但尚有若干問題需要解決,還存在很大的發(fā)展空間,例如:3.1混合方法3.1.1求解各區(qū)域的邊界積分文獻提出了一種新型混合時域有限元-邊界積分方法(TDFE-BI)分析三維電磁散射問題。它使用了一種人工邊界將無限的解區(qū)域劃分成內(nèi)部和外部區(qū)域,分別用時域有限元(TDFEs)和邊界積分(BIs)求解,兩者通過在人工邊界上的場連續(xù)性連接起來。它的優(yōu)點是可以使用一種快速算法,即多層平面波時域(PWTD)方法來求解邊界積分。在考慮一個電大尺寸目標時,引用這個方法大大降低了計算費用;缺點是每個時間步上需要求解一個矩陣方程。有限元-邊界積分混合方法(FEM-BI)在頻域有著較好的應(yīng)用,在時域則剛開始起步,尚待發(fā)展。3.1.2tdfem-floquet定理混合算法周期性結(jié)構(gòu)在各種電磁散射和輻射問題的應(yīng)用有著悠久的歷史。特別是新型的、人工構(gòu)造的周期結(jié)構(gòu)——電磁帶隙材料和左手材料(其介電常數(shù)和磁導(dǎo)率均為負值)的提出,掀起了一股研究周期性結(jié)構(gòu)的熱潮。為了分析這類周期結(jié)構(gòu)的電磁散射和輻射特性,文獻、提出了在二維和三維問題上實現(xiàn)TDFEM與FloquetABC混合的方法:當一個周期性結(jié)構(gòu)包含許多個周期單元晶格時,一般認為對該結(jié)構(gòu)在周期重復(fù)的方向上無限延展的近似分析是合理的。TDFEM-Floquet定理混合算法通過引入一個轉(zhuǎn)換場變量使得在單元晶格的四個面上的周期性邊條件很容易實現(xiàn),能從對一個周期單元晶格的分析得到無限周期結(jié)構(gòu)的特性,該方法的優(yōu)點是將有限元計算區(qū)域限制在一個周期的范圍,大大縮小了計算量。缺點是在每個時間步上也需要解一個矩陣方程,而且轉(zhuǎn)換場變量的引入對時域穩(wěn)定性的影響還沒有數(shù)學(xué)上嚴格的證明。3.1.3新型混合方法PSTD法借助平面波譜展開理論,利用快速傅里葉變換(FFT)及其逆變換代替麥克斯韋方程中的空間差商,該法最初應(yīng)用于FDTD方法,在FDTD方法中PSTD法在處理復(fù)雜的曲線邊界的時候采用階梯近似,有較大的精度誤差。時域有限元法(TDFEM)使用于處理復(fù)雜的曲線邊界,但是,它需要更多的內(nèi)存。為了收益這兩種方法的好處,文獻提出了一種新型時域擬譜法(PSTD)和時域有限元法(TDFEM)混合方法。相比TDFEM-FDTD方法,TDFEM-PSTD大大緩和了空間離散上的條件限制。3.2層級單元的求解除了混合算法的研究,我們也可以通過自適應(yīng)FEM網(wǎng)格剖分技術(shù)改進TDFEM的精度。通常的改善有限元分析精度的方法是利用高階基函數(shù),這些基函數(shù)可分為兩類:插值基函數(shù)和層級基函數(shù)。層級單元可以通過兩種方式來提高精度:一種方式是在誤差估計大的區(qū)域中提高基函數(shù)的階(p-細化);另一種方式是保持基函數(shù)的階數(shù)不變,來細分單元(h-細化)。hp自適應(yīng)有限元法(hp-FEM),顧名思義,是允許h和p都能在網(wǎng)格剖分上獨立地改變。這就使得用更少的未知量獲得想要的精度成為可能。然而,一個主要的難題是如何決定h和p的分布,以及如何確定“細化的方向”。另外網(wǎng)格的局部加密從原理上來說是很簡單的事情,但是在三維的網(wǎng)格