局部特征值解的無條件穩定FDTD高效實施方案

趙斯晗，魏 兵，2，何欣波，2

(1.西安電子科技大學 物理與光電工程學院，陜西 西安 710071；2.西安電子科技大學 信息感知技術協同創新中心，陜西 西安 710071)

傳統時域有限差分方法(Finite Difference Time Domain，FDTD)由于其簡單、高效等特點在計算電磁學中占據重要地位，然而傳統有限差分方法需滿足CFL條件以保證解的穩定性[1-3]。當分析含有精細結構目標的電磁問題時，需要用更小的網格對其進行精確描述，這使得時間步長隨之變小，完成一次計算需要大量的時間步數，導致計算效率較低[4]。

為了解決這一問題，JIAO D團隊提出了基于濾除空間不穩定模的顯式無條件穩定時域有限差分方法(explicit and Unconditionally Stable Finite Difference Time Domain，US-FDTD)[5-6]。由于數值系統的不穩定性是由較大特征值對應的特征模引起的[7]，而這些不穩定模式對所關心頻段內的場值貢獻很小，因此將其消除不會影響計算結果的精確性。然而，當未知量個數或不穩定模式很多時，求解全域矩陣特征值問題及場值迭代的計算成本仍然很高。

由于系統中的不穩定模式是由空間中的細網格和與細網格緊鄰的粗網格造成的[8]，因此筆者給出了一種基于局部特征值求解的顯式無條件穩定時域有限差分方法(Unconditionally Stable Finite Difference Time Domain based on Local eigenvalue solution，USL-FDTD)的高效方案用以解決上述計算成本高的問題。筆者將計算域分為兩部分，細網格和與之緊密相鄰的粗網格為區域I，其余粗網格為區域Ⅱ。于是原始系統矩陣相應地被分為四個局域矩陣塊，僅需找到區域I所對應局域矩陣中的不穩定模式即可獲得全域不穩定模式。這大大降低了特征值求解及矩陣向量相乘的計算復雜度，降低了內存需求和計算時間。

1 US-FDTD算法基礎

為了便于分析，此處考慮無耗問題。麥克斯韋微分方程組差分離散后可寫為如下形式[9]：

(1)

(2)

其中，{h}表示計算區域中磁場未知量構成的向量；{e}表示計算區域內電場未知量構成的向量；{j}為激勵源矢量；Δt表示時間步長；n，n+1，n±1/2表示時刻；Se表示對電場的旋度算符；Sh表示對磁場的旋度算符；ε表示介電常數；μ表示磁導率；D1/ε表示元素為1/ε的對角矩陣；D1/μ表示元素為1/μ的對角矩陣[9]。

將式(1)代入式(2)，可得

(3)

將式(2)中，n+1→n，n→n-1，n+1/2→n-1/2，得

(4)

聯立式(3)和式(4)，可得到電場{e}關于時間的二階方程：

(5)

其中，{f}n=-(jn+1/2-jn-1/2)/Δt，為-?j/?t在n時刻的中心差分[9]。

式(5)實際上為式(6)的中心差分離散：

(6)

系統中的不穩定模式可通過求解全域矩陣S的特征值問題獲得：

SΦi=λiΦi，i=1，2，3，…，n，

(7)

其中，λi和Φi分別表示S中第i個特征值和特征向量，n代表S矩陣的維度。若時間步Δt按照采樣定理確定，則大于4/Δt2的特征值所對應的特征向量Φi稱為不穩定模式[11]。

設Φh為所有不穩定模式構成的矩陣，則式(6)修改為

(8)

對式(8)采用中心差分離散后可得US-FDTD方法的迭代公式[9]：

(9)

2 USL-FDTD方案

2.1 USL-FDTD理論

在計算含有精細結構目標的電磁問題時，由于實際物理尺寸的限制，其空間離散尺度遠遠小于計算上限頻率準確計算所需的空間離散尺度[12]。為滿足CFL條件，用于仿真計算的時間步長就會很小。若時間步長不滿足穩定性條件，則會導致計算的不穩定。

圖1 分區域編號示意圖

不穩定模式是由空間中的細網格和與細網格緊鄰的粗網格造成的[8]，若對計算域未知量進行分組編號，如圖1，子域I包含所有細網格和與細網格直接相鄰的粗網格對應的未知量；子域Ⅱ包含剩余粗網格對應的未知量。

分別對區域Ⅰ、Ⅱ進行編號。Ⅰ區域粗細網格交界處棱邊的電場，如圖2所示，需做插值處理。圖2中h1和h2分別代表1號和2號面片的磁場，e1～e7代表棱邊的電場。“×”代表h2關于e2對稱的位置，h×代表未知磁場。l1和l2分別代表h1和h2到“×”處的距離。

圖2 粗細網格交界處示意圖

h×可用已知磁場表示為

(10)

e2可表示為

(11)

其中，

(12)

將式(12)代入式(11)中，得

(13)

對邊界處棱邊所對應的Sh矩陣元素按照式(13)進行更新。

由于Sh經過了插值處理，這時S矩陣不再是一個對稱矩陣，需要對不穩定模式Φh進行正交化處理，以滿足條件[13]：

(14)

這時S矩陣可自然地被分為4個小矩陣塊，分塊后的S矩陣和電場未知量{e}向量相乘如圖3所示。

圖3 系統矩陣S和電場矢量相乘

圖3中，S11矩陣對應子域Ⅰ，包含所有不穩定模式；S22矩陣對應子域Ⅱ，無不穩定模式；S12矩陣和S21矩陣為耦合矩陣，分別代表區域Ⅰ對區域Ⅱ的場值貢獻和區域Ⅱ對區域Ⅰ的場值貢獻，無不穩定模式。

設{e1}、{e2}分別代表區域Ⅰ和區域Ⅱ的電場未知量，{f}代表激勵源項。分組后的S矩陣、未知量{e}和激勵源項可寫成如下形式：

(15)

式(15)中，只有S11矩陣與細網格及相鄰粗網格相關，即不穩定模式僅存在于S11矩陣中。求解S11的特征值問題即可求得全域不穩定模式：

S11Vi=λiVi，i=1，2，3，…，l，

(16)

其中，λi和Vi分別表示S11中第i個特征值和特征向量，l代表S11矩陣的維度。

設Vh代表所有特征值大于4/Δt2對應的特征向量構成的矩陣。由式(6)，通過對S11矩陣的修改，即可消除整個系統的不穩定模式：

(17)

對式(17)采用中心差分離散，可得

(18)

{e2}n+1=2{e2}n-{e2}n-1-Δt2S22{e2}n-Δt2S21{e1}n-D1/ε{f2}n。

(19)

在每個時間步后，需要對{e1}n+1刪除額外的零空間模式以保證正確性[14]：

(20)

2.2 復雜度分析

下面從特征值求解及場值迭代兩方面分析US-FDTD和USL-FDTD的計算復雜度。

文獻[13]和[15]中，US-FDTD方法計算全域矩陣特征值問題的復雜度為O(k2N)，k為不穩定模個數，N為全域矩陣S的維度；而USL-FDTD中，局部特征值問題的復雜度可降低為O(k2n)，n為矩陣S11的維度(n?N)。因此USL-FDTD提高了特征值求解的計算效率，節省了計算時間。

在式(9)中，全域矩陣和全域電場未知量相乘的復雜度為kO(N)；而在式(18)中，S11矩陣和局域電場未知量相乘的復雜度為kO(n)。因此USL-FDTD降低了矩陣矢量相乘的復雜度，節省了迭代計算的時間。

3 數值結果

3.1 有擋板的PEC諧振腔

在4.00 m×4.64 m諧振腔上下兩側分別加寬度為0.40 m的PEC擋板，構成一個有窄縫的PEC諧振腔，窄縫寬度為0.04 m。由于窄縫很小，因此必須采用細網格模擬。粗網格尺寸為0.10 m，細網格尺寸為0.01 m，網格離散如圖4所示，共離散為4 090條棱邊。區域I為圖中虛線框表示的部分。區域Ⅱ為其余粗網格區域。其中區域I共有526個未知量，區域Ⅱ有3 564個未知量。在(2.00 m，2.32 m)處加微分高斯脈沖線電流源，脈沖寬度為τ=6.666 7×10-9s，t0=5.333 3×10-9s。本例中用傳統FDTD、US-FDTD及USL-FDTD分別計算了(3.5 m，0.5 m)處y方向的電場值，結果如圖5所示。表1給出了USL-FDTD與傳統FDTD和US-FDTD方法的性能對比，使用電腦配置為Intel Core 64×4核處理器，運行頻率3.2 GHz，內存8.0 GB。

圖4 帶有擋板的PEC諧振腔

圖5 場值計算結果

由圖5可看出3個計算結果吻合得很好。表1中，US-FDTD和USL-FDTD方法的時間步長是傳統FDTD方法的10倍，因此US-FDTD和USL-FDTD的迭代時間相比傳統FDTD方法顯著減少。US-FDTD和USL-FDTD需求解矩陣特征值，因此所耗內存比傳統方法大；而USL-FDTD僅需求解局域矩陣特征值，矩陣維度大大降低，所以使用的內存比US-FDTD少，且特征值的求解時間大幅減小，計算性能相比前兩種數值方法優勢明顯。

表1 性能對比

3.2 介質板的透射系數

介質板相對介電常數為15，厚度為0.6 cm，如圖6所示。

圖6 介質板計算模型

粗網格尺寸為1 cm。由于介質板厚度小于一個粗網格的尺寸，為精確模擬，豎直方向需用細網格進行剖分，細網格尺寸為0.1 cm。將計算域離散成640條棱邊，其中區域I共有220條棱邊，區域Ⅱ共有420條棱邊。在介質板上方加一排水平方向的微分高斯脈沖線電流源，脈沖寬度τ=10-9s，t0=8×10-10s。四周采用Mur一階吸收邊界。分別監測無介質板時Q點波形E0和有介質板時Q點透射波形Et，將二者分別做傅里葉變換后的比值作為透射系數。并將USL-FDTD的計算結果與解析法、傳統FDTD方法以及US-FDTD方法進行對比，結果如圖7所示。

圖7 透射系數

由圖7可看出4種方法的計算結果一致。在配置為Intel Core 64×4核處理器，運行頻率3.2 GHz，內存8.0 GB的硬件條件下，USL-FDTD的時間步長為1.67×10-11s，矩陣維度為220，特征值求解花費0.02 s，迭代花費0.055 s；US-FDTD采用的時間步長為1.67×10-11s，矩陣維度為640，特征值求解花費0.46 s，迭代花費0.069 s。傳統FDTD方法的時間步長為1.67×10-12s，迭代耗時11.16 s。以上表明所給方案在處理含精細結構目標的電磁問題時，不僅保證了計算精度，且其計算性能優于傳統FDTD和US-FDTD。

4 結束語

筆者給出了一種基于局部特征值求解的顯式無條件穩定時域有限差分方法的高效實現方案。該方案和US-FDTD均采用大時間步進行迭代，但并不需要求解全域矩陣的特征值，因此相比US-FDTD方法大大降低了矩陣維度，提高了特征值問題的求解效率，降低了內存需求。通過對諧振腔的場值及介質板透射率的計算，說明在處理含有精細結構目標的電磁問題時，該方案的優勢比傳統FDTD和US-FDTD方法更為突出。