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

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

傳統(tǒng)時域有限差分方法(Finite Difference Time Domain，F(xiàn)DTD)由于其簡單、高效等特點在計算電磁學(xué)中占據(jù)重要地位，然而傳統(tǒng)有限差分方法需滿足CFL條件以保證解的穩(wěn)定性[1-3]。當(dāng)分析含有精細(xì)結(jié)構(gòu)目標(biāo)的電磁問題時，需要用更小的網(wǎng)格對其進(jìn)行精確描述，這使得時間步長隨之變小，完成一次計算需要大量的時間步數(shù)，導(dǎo)致計算效率較低[4]。

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

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

1 US-FDTD算法基礎(chǔ)

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

(1)

(2)

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

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

(3)

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

(4)

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

(5)

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

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

(6)

系統(tǒng)中的不穩(wěn)定模式可通過求解全域矩陣S的特征值問題獲得：

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

(7)

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

設(shè)Φh為所有不穩(wěn)定模式構(gòu)成的矩陣，則式(6)修改為

(8)

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

(9)

2 USL-FDTD方案

2.1 USL-FDTD理論

在計算含有精細(xì)結(jié)構(gòu)目標(biāo)的電磁問題時，由于實際物理尺寸的限制，其空間離散尺度遠(yuǎn)遠(yuǎn)小于計算上限頻率準(zhǔn)確計算所需的空間離散尺度[12]。為滿足CFL條件，用于仿真計算的時間步長就會很小。若時間步長不滿足穩(wěn)定性條件，則會導(dǎo)致計算的不穩(wěn)定。

圖1 分區(qū)域編號示意圖

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

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

圖2 粗細(xì)網(wǎng)格交界處示意圖

h×可用已知磁場表示為

(10)

e2可表示為

(11)

其中，

(12)

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

(13)

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

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

(14)

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

圖3 系統(tǒng)矩陣S和電場矢量相乘

圖3中，S11矩陣對應(yīng)子域Ⅰ，包含所有不穩(wěn)定模式；S22矩陣對應(yīng)子域Ⅱ，無不穩(wěn)定模式；S12矩陣和S21矩陣為耦合矩陣，分別代表區(qū)域Ⅰ對區(qū)域Ⅱ的場值貢獻(xiàn)和區(qū)域Ⅱ?qū)^(qū)域Ⅰ的場值貢獻(xiàn)，無不穩(wěn)定模式。

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

(15)

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

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

(16)

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

設(shè)Vh代表所有特征值大于4/Δt2對應(yīng)的特征向量構(gòu)成的矩陣。由式(6)，通過對S11矩陣的修改，即可消除整個系統(tǒng)的不穩(wěn)定模式：

(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 復(fù)雜度分析

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

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

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

3 數(shù)值結(jié)果

3.1 有擋板的PEC諧振腔

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

圖4 帶有擋板的PEC諧振腔

圖5 場值計算結(jié)果

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

表1 性能對比

3.2 介質(zhì)板的透射系數(shù)

介質(zhì)板相對介電常數(shù)為15，厚度為0.6 cm，如圖6所示。

圖6 介質(zhì)板計算模型

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

圖7 透射系數(shù)

由圖7可看出4種方法的計算結(jié)果一致。在配置為Intel Core 64×4核處理器，運行頻率3.2 GHz，內(nèi)存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。傳統(tǒng)FDTD方法的時間步長為1.67×10-12s，迭代耗時11.16 s。以上表明所給方案在處理含精細(xì)結(jié)構(gòu)目標(biāo)的電磁問題時，不僅保證了計算精度，且其計算性能優(yōu)于傳統(tǒng)FDTD和US-FDTD。

4 結(jié)束語

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