應(yīng)用FDTD方法解決電磁輻射問(wèn)題

1、 應(yīng)用FDTD方法解決電磁輻射問(wèn)題 自電磁場(chǎng)基本方程以來(lái)，電磁場(chǎng)理論和應(yīng)用的發(fā)展已經(jīng)有一百多年的歷史。目前，電磁波的研究已深入到各個(gè)領(lǐng)域，應(yīng)用十分廣泛，例如無(wú)線電波傳波，光纖通信和移動(dòng)通信，雷達(dá)技術(shù)，微波，天線，電磁成像，地下電磁探測(cè)，電磁兼容等等。在各類復(fù)雜系統(tǒng)中的電磁問(wèn)題，主要依靠各種電磁場(chǎng)數(shù)值計(jì)算方法加以解決。隨著電子計(jì)算機(jī)處理能力和存儲(chǔ)容量的巨大發(fā)展，更促進(jìn)了這些計(jì)算方法在實(shí)際問(wèn)題中的應(yīng)用。目前在電磁場(chǎng)領(lǐng)域應(yīng)用的數(shù)值算法也是種類繁多，各有其優(yōu)缺點(diǎn)，常用的電磁場(chǎng)計(jì)算方法大致有：FDTD Finite difference time domain（時(shí)域有限差分法）TLM Tran

2、smission line method （傳輸線法）FEM Finite element method （有限元法）BEM Boundary element method （邊界元法）MoM Method of moments （矩量法）其中時(shí)域有限差分法(FDTD)理論經(jīng)過(guò)30多年的發(fā)展和完善，已經(jīng)成為時(shí)域電磁場(chǎng)數(shù)值計(jì)算的主要方法之一，并廣泛應(yīng)用各類實(shí)際工程電磁場(chǎng)中。一、 FDTD法簡(jiǎn)介時(shí)域有限差分法以差分原理為基礎(chǔ)，直接從概括電磁場(chǎng)普遍規(guī)律的麥克斯韋旋度方程出發(fā)，將其轉(zhuǎn)換為差分方程組，在一定體積內(nèi)和一段時(shí)間上對(duì)連續(xù)電磁場(chǎng)的數(shù)據(jù)采樣。因此，它是以電磁場(chǎng)問(wèn)題的最原始、最本質(zhì)、最完備的數(shù)值模擬

3、。以它為基礎(chǔ)制作的計(jì)算程序，對(duì)廣泛的電磁場(chǎng)問(wèn)題具有通用性，因此得到了廣泛的應(yīng)用。1. Yee差分算法基本原理考慮空間一個(gè)無(wú)源區(qū)域，其煤質(zhì)參數(shù)不隨時(shí)間變化且各向同性，由Maxwell方程組中的兩個(gè)旋度方程在直角坐標(biāo)系中可導(dǎo)出六個(gè)耦合公式： 其中為介電常數(shù)（F/m）；為磁導(dǎo)率（H/m）；為電導(dǎo)率（S/m）；為磁阻率（）。 按照Yee差分算法，首先在空間建立矩形差分網(wǎng)格，網(wǎng)格節(jié)點(diǎn)的空間坐標(biāo)與一組相應(yīng)的整數(shù)標(biāo)號(hào)一一對(duì)應(yīng)： 而該點(diǎn)的任意函數(shù)F(x,y,z,t)在時(shí)刻的值可以表示為 分別為矩形網(wǎng)格沿x，y，z方向的空間步長(zhǎng)是時(shí)間步長(zhǎng)。時(shí)域差分法的實(shí)質(zhì)就是將計(jì)算問(wèn)題在時(shí)間和空間進(jìn)行離散: 用中心差分取二階精

4、度：對(duì)空間離散： 對(duì)時(shí)間離散：為了獲得（1.10）式的精度，并滿足（1.3）至（1.8），可將空間任意矩形網(wǎng)格上的E和H的六個(gè)分量如圖1-1所示放置。 圖1-1 Yee網(wǎng)格單元及電磁場(chǎng)空間離散點(diǎn)關(guān)系每個(gè)磁場(chǎng)分量由四個(gè)電場(chǎng)分量環(huán)繞著；每個(gè)電場(chǎng)分量由四個(gè)磁場(chǎng)分量所環(huán)繞。為獲得（1.11）的精度，可將E和H在時(shí)間上相差半個(gè)步長(zhǎng)交替計(jì)算。按照這些原則，可將式（1.3）（1.8）化為差分方程組。以（1.3）式為例，差分后得到（1.12） 其余的也可以寫出，每個(gè)網(wǎng)格點(diǎn)上的個(gè)場(chǎng)分量的新值依賴于該點(diǎn)在前一時(shí)間步長(zhǎng)時(shí)刻的值即該點(diǎn)周圍臨近點(diǎn)上另一場(chǎng)量在早半個(gè)時(shí)間步長(zhǎng)時(shí)的值。因此任意時(shí)刻可一次算出一個(gè)點(diǎn)，并行算法可

5、計(jì)算出多個(gè)點(diǎn)。通過(guò)這些運(yùn)算可以交替算出電場(chǎng)磁場(chǎng)在各個(gè)時(shí)間步的值。具體參考12，這里不多敘述。2. 數(shù)值穩(wěn)定性條件與網(wǎng)格劃分?jǐn)?shù)值計(jì)算中一個(gè)比較重要的問(wèn)題就是數(shù)值解的穩(wěn)定性問(wèn)題，這在FDTD方程中取決于網(wǎng)格的空間步長(zhǎng)和時(shí)間不長(zhǎng)的關(guān)系。 Ø 數(shù)值穩(wěn)定性條件以上所建立的麥克斯韋旋度方程的有限差分算式是顯式的，時(shí)間步長(zhǎng) ，空間步長(zhǎng) 必須滿足一定關(guān)系，否則會(huì)使得數(shù)值表現(xiàn)不穩(wěn)定。這種不穩(wěn)定表現(xiàn)為在解顯示差分方程時(shí)，隨著時(shí)間步數(shù)的增加，計(jì)算結(jié)果也無(wú)限制的增加。Taflaove等對(duì)Yee差分格式穩(wěn)定性進(jìn)行討論，導(dǎo)出了對(duì)時(shí)間步長(zhǎng)的限制性條件： 當(dāng) 時(shí)，數(shù)值穩(wěn)定條件（12）可簡(jiǎn)化為 一維條件下，這一穩(wěn)定性

6、條件就簡(jiǎn)化為： 上式可以看出，即要求時(shí)間步長(zhǎng)不能大于電磁波傳播一個(gè)空間步長(zhǎng)所需的時(shí)間，否則，就破壞了電磁波傳播的因果關(guān)系12。如果計(jì)算空間中的媒質(zhì)是不均勻的，那么穩(wěn)定條件對(duì)不同的煤質(zhì)區(qū)域是不同的，我們只需取最大的V滿足的條件，在其他區(qū)域中也就得到滿足了。此時(shí)有：其中n為空間維度。Ø 網(wǎng)格劃分由數(shù)值穩(wěn)定性條件可知，網(wǎng)格劃分需要滿足一定的條件，在有足夠條件的情況下，將網(wǎng)格劃分的越細(xì)，計(jì)算的結(jié)果就會(huì)越精確。目前有一款交互式軟件網(wǎng)格劃分軟件GEOM，根據(jù)設(shè)計(jì)需要，可自動(dòng)的用FDTD的方法分析幾何結(jié)構(gòu)，產(chǎn)生離散模型。但是計(jì)算機(jī)計(jì)算能力畢竟是有限的，有時(shí)候提高一定比例的網(wǎng)格大小，并不一定能夠獲得

7、更多的精度，同時(shí)計(jì)算的效率卻下降了很多。而且對(duì)于一些復(fù)雜的模型結(jié)構(gòu)，由于介質(zhì)層的物理尺寸很小，較大的網(wǎng)格無(wú)法覆蓋這樣的結(jié)構(gòu)，就導(dǎo)致了劃分后介質(zhì)層的缺損。因此需要在滿足穩(wěn)定性條件的前提下對(duì)劃分的大小進(jìn)行優(yōu)化。最好的辦法就是采用非均勻網(wǎng)格劃分的辦法。合理的利用了計(jì)算區(qū)域內(nèi)不同的計(jì)算需要，對(duì)于一些關(guān)鍵位置的地方，如不同介質(zhì)的交接處，或者是結(jié)構(gòu)細(xì)微處，或者是饋源的地方，對(duì)網(wǎng)格細(xì)化，而在介質(zhì)均一的位置，則放大網(wǎng)格的大小，使得總的網(wǎng)格數(shù)量得到降低。對(duì)于非均勻劃分方案的確定，也并不是盲目的過(guò)程?？梢允紫壤幂^大的網(wǎng)格大小進(jìn)行一下預(yù)模擬，對(duì)模擬的結(jié)果進(jìn)行分析，判斷電場(chǎng)(或磁場(chǎng)等其他場(chǎng)分布)分布的規(guī)律及形態(tài)，對(duì)

8、于關(guān)鍵計(jì)算的部位，再采用較小的網(wǎng)格劃分。其次，可以利用遺傳算法等優(yōu)化算法，對(duì)網(wǎng)格劃分進(jìn)行局部最優(yōu)化的搜索，提高網(wǎng)格劃分的質(zhì)量。下圖1.2是對(duì)手機(jī)天線幾種網(wǎng)格劃分示意圖。 圖1.2a 均勻網(wǎng)格劃分 圖1.2b 非均勻網(wǎng)格劃分 圖1.2c含有亞網(wǎng)格的網(wǎng)格劃分除了非均勻網(wǎng)格劃分的方法，另外一種在這個(gè)基礎(chǔ)上的改進(jìn)網(wǎng)格劃分方法就是亞網(wǎng)格劃分的方法14。亞網(wǎng)格劃分在原有網(wǎng)格劃分的基礎(chǔ)上，在計(jì)算的關(guān)鍵點(diǎn)，或者是重點(diǎn)關(guān)注的部分，或者是介質(zhì)層物理尺寸很薄的地方，將網(wǎng)格進(jìn)一步劃小，而其它部分的網(wǎng)格則保持不變。這樣做的好處在于可以在不影響其它部分的情況下，提高重點(diǎn)計(jì)算區(qū)域的計(jì)算精度。同時(shí)網(wǎng)格的數(shù)量不會(huì)因此激增。3.

9、 吸收邊界條件當(dāng)我們考慮像輻射和散射這種開(kāi)放向題時(shí)，原則上應(yīng)該在無(wú)限空間中計(jì)算電磁場(chǎng)，但是，由于計(jì)算機(jī)存儲(chǔ)量的限制不可能考慮很大的空間,，因而，一個(gè)基本問(wèn)題是如何處理網(wǎng)格的截?cái)嘞蝾}。一種有效的方法是在截?cái)嗵幵O(shè)置一種吸收邊界條件，如果它能吸收全部向外傳播的波，就等效于考慮了一個(gè)無(wú)限大的空間。目前構(gòu)造吸收邊界條件的思路主要有兩種：一種是從外形波方程出發(fā)構(gòu)造的透射邊界條件，如Mur邊界條件，這種類型的透射邊界條件具有構(gòu)造簡(jiǎn)單，內(nèi)存需求小等特點(diǎn)。另一種是在邊界上引入吸收材料，電磁波在無(wú)反射地進(jìn)入吸收材料后被衰減掉，如PML。這種方法構(gòu)造復(fù)雜，內(nèi)存需求量較大，但在很大的入射角度的吸收條件較好。Mur給出

10、了適合于在時(shí)域有限差分應(yīng)用的吸收邊界條件的一階和二階的近似形式，具體推導(dǎo)過(guò)程可見(jiàn)1617。 應(yīng)該注意的是Mur的局限性：（1） 邊界面上任意點(diǎn)的插值都是在其領(lǐng)域的三維空間上進(jìn)行的，不允許靠近邊界的區(qū)域出現(xiàn)介質(zhì)的不連續(xù)；（2） 在垂直入射的情況下的計(jì)算精度僅相當(dāng)于一階吸收邊界條件。完全匹配層PML是通過(guò)在截?cái)噙吔缤夥胖锰摂M的吸收介質(zhì)來(lái)系數(shù)外形波的，PML邊界由硬邊界構(gòu)成，在計(jì)算PML邊界時(shí)，就可以不涉及邊界外的未知量，大幅度提高了計(jì)算精度，詳細(xì)過(guò)程見(jiàn)18與Mur邊界條件相比，PML邊界具有更高的計(jì)算精度而且任意方向入射時(shí)均無(wú)反射。在二維自由空間，PML邊界吸收的反射系數(shù)可低于-70dB，比其他各

11、種邊界條件改善約40dB。這種人工涉及的完全匹配層由有耗導(dǎo)電、導(dǎo)磁煤質(zhì)組成，可吸收人以入射角、任何頻率、任意偏振態(tài)的入射電磁波，實(shí)用性更強(qiáng)，它內(nèi)存大，盡管有的文獻(xiàn)19提出的各向異性媒質(zhì)的PML吸收邊界條件所需計(jì)算機(jī)內(nèi)存量減少了一半，但是其實(shí)際計(jì)算量并沒(méi)有減少多少，在不連續(xù)介質(zhì)接觸邊界面和三維區(qū)域的6個(gè)面不全為吸收邊界時(shí)也難以應(yīng)付，具體算法可見(jiàn)19。目前各種吸收邊界條件都存在這一些缺點(diǎn)，吸收邊界條件反射過(guò)大，所適應(yīng)的入射角的范圍不夠?qū)拸V，提高計(jì)算精度等都應(yīng)該是我們努力的方向。二、 電磁輻射問(wèn)題電磁輻射是由空間共同傳送的電能量和磁能量所組成，而該能量是由電荷移動(dòng)所產(chǎn)生的。時(shí)變的電荷和電流是激發(fā)電磁

12、波的源。為了有效使電磁波能量按所要求的方向輻射出去，時(shí)變的電荷和電流必須按某種特殊的方式分布，天線就是設(shè)計(jì)成按規(guī)定的方式有效的輻射電磁波能量的裝置。其實(shí)質(zhì)就是一種導(dǎo)行波與自由空間的轉(zhuǎn)換器件或者換能器，用于發(fā)送和接受自由空間的電磁波。舉例說(shuō)，正在發(fā)射訊號(hào)的射頻天線所發(fā)出的移動(dòng)電荷，便會(huì)產(chǎn)生電磁能量，電磁包括形形色色的電磁輻射，從低頻電磁輻射到高頻電磁輻射，兩者之間還有無(wú)線電波、微波、紅外線和紫外光等。任何形式和形狀的天線都可以看成是有無(wú)限多個(gè)載有交流信號(hào)的基本小線元組成，這些基本元通常被稱為電偶極子，電偶極子產(chǎn)生的電磁場(chǎng)分析計(jì)算是線性天線工程的基礎(chǔ)?？筛鶕?jù)天線上的電流分布來(lái)計(jì)算由其產(chǎn)生的電磁場(chǎng)，

13、由給定的求出A，再根據(jù)求得B，最后由 求得。舉例如下：在自由空間取一個(gè)觀察點(diǎn)P。如圖23所示。用直角坐標(biāo)系表示。其中假設(shè)介質(zhì)的電磁常數(shù)為 、，只要知道表面電流就可以求出自由空間的電磁場(chǎng)分布。 圖2.1 帶有表面電流的理想導(dǎo)體再由電場(chǎng)的矢量波動(dòng)方程可以由Maxwell方程組的第一個(gè)方程解出 在FDTD差分格式中，由于網(wǎng)格中的電流密度取的是平均值，因此其中，由可得： 的解為：在電磁學(xué)中，將自由空間按照輻射距離r劃分為三個(gè)區(qū)域近場(chǎng)反應(yīng)區(qū)，近區(qū)場(chǎng)，遠(yuǎn)區(qū)場(chǎng)。當(dāng)即的區(qū)域稱為近區(qū)，即的區(qū)域稱為遠(yuǎn)區(qū)。近區(qū)場(chǎng)求解過(guò)程見(jiàn)2021。這里以遠(yuǎn)場(chǎng)區(qū)為例，代入（2.4），對(duì)于遠(yuǎn)區(qū)場(chǎng)，和相較于來(lái)說(shuō)可以忽略。 在遠(yuǎn)場(chǎng)區(qū)域，電

14、場(chǎng)體現(xiàn)出球面波的特性，更加一般的輻射體包括磁感表面 ，可寫出三維坐標(biāo)中電場(chǎng)強(qiáng)度，。過(guò)程中求出的，按式進(jìn)行離散差分。用FDTD方法解決電磁輻射問(wèn)題的一般過(guò)程：1) 用空間上相鄰的節(jié)點(diǎn)，時(shí)間上當(dāng)前以及上一個(gè)時(shí)間表示的當(dāng)前時(shí)間當(dāng)前節(jié)點(diǎn)的電磁場(chǎng)分布表達(dá)式;2) 在足夠的觀測(cè)時(shí)間長(zhǎng)度內(nèi)，計(jì)算在網(wǎng)格邊界處的電磁場(chǎng)，大??；3) 將時(shí)域表達(dá)的，通過(guò)傅立葉變換，轉(zhuǎn)換到頻域表達(dá)的，；4) 通過(guò)計(jì)算得到的場(chǎng)分布，得到網(wǎng)格邊界處的等效表面電流，；5) 運(yùn)用輻射方程式，計(jì)算得到遠(yuǎn)場(chǎng)分布。經(jīng)過(guò)上述幾個(gè)步驟，就可以完成電磁場(chǎng)的時(shí)域差分計(jì)算。三、 FDTD方法的改進(jìn)應(yīng)用傳統(tǒng)的時(shí)域有限差分法分析、計(jì)算偽低頻電磁場(chǎng)問(wèn)題時(shí)，由于空

15、間步長(zhǎng)遠(yuǎn)小于電磁波波源波長(zhǎng)，滿足Courant-Friedrich-Levy穩(wěn)定條件的時(shí)間步長(zhǎng)，遠(yuǎn)小于電磁波波源的周期。這樣將導(dǎo)致很大的時(shí)間步數(shù)。為消除Courant-Friedrich-Levy穩(wěn)定條件限制，在滿足計(jì)算精度的條件下采用較大的時(shí)間步長(zhǎng)以提高算法的計(jì)算效率方面，世界上許多學(xué)者做了大量的探索工作，提出了幾類切實(shí)可行的解決方案，其中具有代表性的主要有兩類。一類是使用改進(jìn)的Lanczoc方法,另一類是采用交替方向隱格式時(shí)域有限差分法(ADI-FDTD)2627。改進(jìn)的Lanczoc方法和ADIFDTD法都是以1966年Yee提出的網(wǎng)格體系為基礎(chǔ)提出的新方法。ADI-FDTD法保留了時(shí)域

16、有限差分法計(jì)算格式簡(jiǎn)單等大部分優(yōu)點(diǎn)，是一種無(wú)條件穩(wěn)定的差分方法。可任意設(shè)定計(jì)算時(shí)間步長(zhǎng)，但隨著時(shí)間步長(zhǎng)的增加計(jì)算精度會(huì)不斷降低。因此三維交替方向隱格式時(shí)域有限差分法的計(jì)算精度尚有待于進(jìn)一步完善。由于應(yīng)用FDTD方法時(shí)受到計(jì)算機(jī)存儲(chǔ)容量的限制，計(jì)算區(qū)域網(wǎng)格空間不能無(wú)限制增加，所以采用并行計(jì)算技術(shù)來(lái)提高時(shí)域有限差分計(jì)算能力也是解決問(wèn)題的有效途徑之一。從串行計(jì)算到并行計(jì)算，可以大幅度提高計(jì)算能力，通過(guò)并行FDTD算法可以解決大內(nèi)存的要求，可以大大縮短時(shí)間。參考文獻(xiàn)1 A. Taflove and K. R. Umashankar, “Review of FD-TD numerical modelin

17、g of electromagnetic wave scattering and radar cross section,”IEEE Proc., vol. 77, pp. 682499, May 1989.2 A. Taflove and K. R. Umashankar, “The finite-difference time-domain method for numerical modeling of electromagnetic wave interactions,”Elecrromugn., vol. 10, pp. 105-126, Jan. 1990.3 A. Taflove

18、, “Computational Electrodynamics: The Finite-Difference Time-Domain Method,” Artech House, Norwood, MA, 1995.4 K. R. Demarest, “A finite-difference time-domain technique for modeling narrow apertures in conducting scatterers,” fEEE Trans. Antennas Propagat., vol. 35, pp. 826-831, July 1987.5 P. A. T

19、irkas and C. A. Balanis, “Finite-difference time-domain method for antenna radiation,” IEEE Trans. Antennas Propagat., vol. 40, pp. 334-340, Mar. 1992.6 K. S. Kunz and K. Lee, “A three-dimensional finite-difference solution of the external response of an aircraft to a complex transient EM environmen

20、t: Part IILComparison of predictions and measurements,”IEEE Trans. Electromagn. Compat., vol. 20, pp. 333-341, May 1978.7 Z Li.Y Rahmat-Samii Whip-PIFA combination in wireless handset application:A hybridcircuit model and full wave analysis 20048 葛德彪.閆玉波電磁波時(shí)域有限差分方法 20029 陳軍. 自然使用下手機(jī)對(duì)人體頭部電磁輻射的SAR計(jì)算學(xué)位

21、論文200910 郭佳. 基于FDTD算法的手機(jī)天線對(duì)人體頭部電磁輻射的仿真研究學(xué)位論文200611 盧玲.“時(shí)域有限差分(FDTD)法及其吸收邊界條件的研究，”浙江大學(xué)碩士學(xué)位論文2000 12 KSKunz and RJLuebber“The FiniteDifference TimeDomain Method Electromagnetics”， IEEE TransCRC Press, 199313 K. S. Yee, “Numerical solution for initial boundary value problems involving Maxwells equation

22、s in isotropic media,” IEEE Trans. Antennas Propagat., vol. AP-14, pp. 302307, May 1966.14 MOkoniewski，EOkoniewska，and MAStuehly，“Threedimensional subgriddingalgorithm for FDTD”，IEEE TransAntennas Propagat，1997 Mar v0145，42242915 G. Mur, “Absorbing boundary conditions for the finite-difference appro

23、ximation of the time-domain electromagnetic field equations,” IEEE Trans. Electromagn. Compar., vol. 23, pp. 377-382, Nov. 1981.16 L. Halpern and L. N. Trefethen, “Wide-angle one-way wave equations,”J. Acoustical Soc. America, vol. 84, pp. 1397-1404, Oct. 1988.17 B. Engquist and A. Majda, “Absorbing

24、 boundary conditions for the numerical simulation of waves,” Math. Compuf., IEEE Transvol. 3 1, pp. 629-65 1, July 1977.18 J P Berenger，”A perfectly matched layer for the absorption of electromagnetic waves，”。Comput Physl CS，vol 1 14，PP l 85-200，October 199419 J Fang，and Z Wu，”Generalized perfectly

25、Matched Layer for absolption of propagating and Evanescent waves in lossless and lossy media，”IEEE Trans Microwm，e Theory Tech，vol 44，no 1 2，PP 22 1 6-2222，Dec 1 99620 石新軍. 電偶極子天線電磁場(chǎng)輻射特性，蘭州大學(xué)學(xué)報(bào)，Vol. 45 Supp. Jun. 200921 王長(zhǎng)青. 陳金元，祝西里 Use of the Finite_Difference Time-Domain Method for calculating Near-Field of Electromagetic Radition Systems, 電子學(xué)報(bào), Vol. 18 Nel Jan. 199022 金夢(mèng)筆. 手機(jī)電磁輻射的時(shí)域有限差分法分析及應(yīng)用學(xué)位論文碩士 200623 G. Mur, “Absorbing boundary conditions for the finite-difference approximation of the time-domain elec