中國電子科技集團(tuán)公司第二十研究所 李宏亮 黨 杰

時間域電磁波有限差分方法(FDTD)是一種常用的電磁波數(shù)值模擬方法，通過將電磁波進(jìn)行網(wǎng)格化處理，使得電場和磁場均勻的分布在網(wǎng)格上，從而求解電磁波方程。本文基于麥克斯韋方程組，開展了基于頻率域電磁波有限差分的相關(guān)研究。為了驗證頻率域有限差分的可行性以及提高數(shù)值模擬精度，主要采用了兩種方法：第一種方法是在麥克斯韋方程組的基礎(chǔ)上直接推導(dǎo)出頻率域電磁波方程，避免了時間域遞推造成的累計誤差；第二種方法是采用優(yōu)化25點(diǎn)差分方法構(gòu)建頻率域電磁波方程算子，不僅保證了模擬的精度，而且能有效的抑制數(shù)值頻散。此外，為了驗證算法的可行性，在PML邊界條件下進(jìn)行了頻率域電磁波的數(shù)值模擬。

隨著電磁波理論研究與實際應(yīng)用的推廣，麥克斯韋方程組的數(shù)值解法也被廣泛的關(guān)注和研究。已經(jīng)提出的方法有頻率域的矩量法、有限元法、有限差分法以及時間域的傳輸線矩陣法、有限差分等方法。各種方法都有自己的優(yōu)勢以及應(yīng)用上的局限性，要么計算過于復(fù)雜，要么占用內(nèi)存大、計算速度慢，在實際應(yīng)用中往往需要將幾種方法混合使用。

有限差分方法的中心思想是用差分代替微分，通過將研究區(qū)域網(wǎng)格化處理，使電場和磁場有規(guī)律的分布在網(wǎng)格點(diǎn)處，從而得到網(wǎng)格區(qū)域的電場和磁場值。本文通過推導(dǎo)麥克斯韋方程組得到頻率域二維電磁波方程，然后整理成矩陣方程形式，采用優(yōu)化25點(diǎn)有限差分方法構(gòu)建頻率域差分算子，對算子進(jìn)行稀疏化存儲，使得頻率域數(shù)值模擬更加簡便快捷，同時也證明了該方法的有效性。

1 頻率域電磁波優(yōu)化25點(diǎn)有限差分

1.1 頻率域電磁波方程推導(dǎo)

麥克斯韋方程組的微分形式是這樣的，如式(1)：

式(1)中，E表示電場強(qiáng)度，B表示磁場強(qiáng)度；ρ表示電荷密度，J表示電流密度；ε0和μ0分別表示真空中的介電常數(shù)和磁導(dǎo)率；是矢量微分算子，和分別表示求取散度和旋度。

為簡化后續(xù)推導(dǎo)，考慮真空、均勻介質(zhì)(即ρ=0，J=0)情況下，麥克斯韋方程組變?yōu)槭?2)：

式(2)方程兩側(cè)同時取旋度，并作進(jìn)一步簡化，得到式(3)：

考慮二維情況，可以將式(3)改寫為式(4)：

式(4)兩側(cè)對時間做傅氏變換，即：

式(5)中，ω表示角頻率；和分別是電場強(qiáng)度和磁場強(qiáng)度的頻率域表示。

1.2 基于優(yōu)化25點(diǎn)有限差分方法的算子離散化

圖1所示為優(yōu)化25點(diǎn)有限差分方法中各個權(quán)系數(shù)分布示意圖。在求解時，需要用到該點(diǎn)周圍的25個點(diǎn)，而這25個點(diǎn)對于要求解的二階差分算子“貢獻(xiàn)程度”是不同的。如圖1(a)，每一行需要用到5個點(diǎn)來構(gòu)造2個差分算子，根據(jù)每一行中各列元素距離待計算點(diǎn)的距離遠(yuǎn)近，設(shè)置加權(quán)系數(shù)c、d，將他們結(jié)合起來作為這一行的差分算子，因此可以得到5個差分算子。然后，再將這5個差分算子通過加權(quán)系數(shù)b1、b2、b3進(jìn)行加權(quán)平均即可得到最終的離散形式。同理，也通過類似的方法可以獲得，如圖1(b)所示。圖1(c)是方程左端加速項的差分格式，其權(quán)系數(shù)以差分點(diǎn)為圓心，分布在一組同心圓上。

圖1 優(yōu)化25點(diǎn)差分示意圖

為了使網(wǎng)格色散和數(shù)值各向異性盡可能最小，必須確定一組最優(yōu)的加權(quán)系數(shù)使相速度和群速度盡可能的接近。求取這些系數(shù)可以采用梯度法或牛頓法，通過給定一組初始值，不斷地迭代搜索以使該系數(shù)能夠同時使相速度和群速度接近。

本文采用Min利用高斯牛頓法計算出來的最優(yōu)加權(quán)系數(shù)：

1.3 吸收邊界條件及離散方程

電磁波在空間中傳播時相當(dāng)于是在一個半無窮的介質(zhì)中傳播，沒有邊界的影響。但在數(shù)值模擬過程中，由于計算區(qū)域有限，而且人為的引入了反射界面，這就會導(dǎo)致出現(xiàn)一些在實際中不會出現(xiàn)的反射波，影響最終的模擬結(jié)果。因此，必須構(gòu)造一個吸收邊界，以期望盡可能的抑制這種邊界反射，從而更加真實的模擬電磁波的傳播過程，研究電磁波的傳播特性。

本文采用PML（如圖2所示）吸收邊界來抑制邊界反射。當(dāng)?shù)卣鸩▊鞑サ竭吔鐣r，隨著進(jìn)入到吸收層內(nèi)傳播距離的增加，電磁波不斷地衰減，這樣當(dāng)波傳播到最外層邊界時，其能量被完全衰減就不會產(chǎn)生邊界反射。當(dāng)電磁波傳播到區(qū)域1時，沿x和z方向都衰減；當(dāng)波傳播到區(qū)域2時，沿x方向不衰減，沿z方向衰減；當(dāng)波傳播到區(qū)域3時，沿x方向衰減，沿z方向不衰減。

圖2 PML吸收邊界示意圖(引自王守東,2003)

具體的，將這種邊界施加到電磁波方程時，需要做一個坐標(biāo)拉，這樣施加PML吸收邊界的電磁波方程就可以表示為式(9)：伸的變換

用d來表示衰減因子，使得PML吸收區(qū)域內(nèi)d≠0，而有效計算區(qū)域內(nèi)d=0。定義，得：

則：

其中：

式中A是在0.5～1.79之間取值的常數(shù)，通常需要通過不斷地試驗確定最佳的參數(shù)，后文中數(shù)值模擬時采用A=1.3；f0是主頻；LPML是PML層的厚度；lf是PML邊界與內(nèi)部有效計算區(qū)域的距離。

將衰減因子代入電磁波方程(5)式，得到：

根據(jù)文獻(xiàn)，為簡化后續(xù)推導(dǎo)，通常省略式(12)中一階偏導(dǎo)項，亦即：

將差分加權(quán)系數(shù)以及式(10)～式(13)帶入到式(9)中，施加PML邊界條件，合并同類項后得到最終的離散方程為：

式中參數(shù)為：

方便起見，將式(14)和式(15)聯(lián)合起來整理成矩陣方程形式：

式(16)中，A是由差分系數(shù)構(gòu)成的一個具有25條對角線的系數(shù)矩陣，該矩陣是一個大型的、非對稱的稀疏矩陣；U是由和構(gòu)成的一個列向量；F是施加的激發(fā)源。若網(wǎng)格區(qū)域為N=m×n，則A是一個2N×2N的矩陣，U是2N×1的向量，F(xiàn)是2N×1的向量。

2 數(shù)值模擬

2.1 數(shù)值模擬

為了驗證各種激發(fā)源加載的有效性，在均勻各向同性介質(zhì)中進(jìn)行數(shù)值模擬實驗。模型的網(wǎng)格大小為160×160，x和z方向的網(wǎng)格間距為dx=dz=0.005m。圖3是激發(fā)源主頻為3GHz，在模型中心位置處激發(fā)得到的單頻切片。

圖3 集中力源單頻切片

2.2 精度對比

為了測試差分點(diǎn)數(shù)對數(shù)值模擬精度的影響，在相同參數(shù)下，采用網(wǎng)格大小201×201，網(wǎng)格間距Δx= Δz=0.005m，激發(fā)源主頻為3GHz，放置在模型中心處。根據(jù)式(16)分別采用9點(diǎn)差分方法和25點(diǎn)差分方法構(gòu)建系數(shù)矩陣，進(jìn)行電磁波正演模擬。圖4是將頻率域各個頻點(diǎn)單頻切片匯總后進(jìn)行反傅里葉變換得到的時間域某時刻波場快照?？梢钥闯觯ㄟ^優(yōu)化25點(diǎn)差分方法構(gòu)建正演算子得到的模擬結(jié)果較之常規(guī)方法而言能夠很好地抑制數(shù)值頻散。

圖4 波場快照

結(jié)論與展望：本文在麥克斯韋方程組的基礎(chǔ)上，推導(dǎo)出了真空中電磁波的頻率域表達(dá)式，并采用優(yōu)化25點(diǎn)有限差分方法對電磁方程進(jìn)行了差分，得到了基于優(yōu)化25點(diǎn)有限差分方法的離散電磁場方程。通過添加PML邊界條件，進(jìn)行了頻率域電磁波數(shù)值模擬，得到了以下認(rèn)識和結(jié)論：

（1）采用25點(diǎn)有限差分方法進(jìn)行電磁波數(shù)值模擬具有更高的精度；

（2）由于頻率域系數(shù)矩陣只與差分系數(shù)有關(guān)，故系數(shù)矩陣A只需要計算一次，這使得頻率域電磁波數(shù)值模擬計算速度快，很適合并行計算，為理論研究和工程實踐提供了一個較好的實現(xiàn)思路，為電磁波成像和反演奠定了一定的基礎(chǔ)。

當(dāng)然，這種方法從目前來看仍有一定的限制，比如從二維模型轉(zhuǎn)向?qū)嶋H三維模型時系數(shù)矩陣需要的巨大內(nèi)存或者說是系數(shù)矩陣稀疏化如何處理，以及實際工程應(yīng)用中大量的頻率域數(shù)據(jù)如何轉(zhuǎn)化等問題都有待進(jìn)一步考慮。