李建雄，魏之棟

（1.天津工業(yè)大學(xué)電子與信息工程學(xué)院，天津 300387；2.天津工業(yè)大學(xué) 天津市光電檢測技術(shù)與系統(tǒng)重點實驗室，天津 300387）

隨著計算機技術(shù)與電磁基礎(chǔ)理論的結(jié)合日益緊密，時域有限差分法（finite-difference time-domain，F(xiàn)DTD）在電磁學(xué)的各個領(lǐng)域得到了廣泛的應(yīng)用。然而，傳統(tǒng)FDTD 方法的計算速度受到柯朗-弗里德里希斯-列維（Courant-Friedrichs-Lewy，CFL）條件的限制。因此，學(xué)者們提出了一些無條件穩(wěn)定的FDTD 方法，如交替方向隱式FDTD（alternatingdirectionimplicitFDTD，ADI-FDTD）[1-3]。但是，當(dāng)時間步長較大時，ADI-FDTD方法存在較大的離散誤差和色散誤差[4]。分裂步FDTD（split-step FDTD，SS-FDTD）方法[5-7]也是一種無條件穩(wěn)定的算法。SS-FDTD 方法根據(jù)“Strang”分裂格式將1 個時間步分裂為3 個子步。與ADI-FDTD 方法相比，該方法具有更好的精度[8]。

等離子體作為一種色散介質(zhì)，在天線、等離子體光子晶體器件和隱身技術(shù)中發(fā)揮著重要作用[9-12]。而二維場景下的等離子體算法也在很多領(lǐng)域中有廣泛的應(yīng)用，如確定反射計系統(tǒng)的各種參數(shù)以及應(yīng)對航天器返回時的通信中斷問題[13-14]等。到目前為止，學(xué)者們已經(jīng)提出了幾種解決色散介質(zhì)的方法。在這些方法中，輔助微分方程（auxiliary differential equation，ADE）方法[15-17]避免了復(fù)雜的轉(zhuǎn)換和卷積，是一種簡單有效的色散介質(zhì)仿真方法。近年來，Song 等[18-19]學(xué)者提出了一種龍格-庫塔輔助微分方程（Runge-Kutta ADE，RKADE）的算法。該算法采用龍格-庫塔技術(shù)代替線性差分技術(shù)對非磁化等離子體的微分方程進行離散。與傳統(tǒng)的ADE 方法相比，RKADE 具有更高的精度，并且不會增加存儲和計算負擔(dān)。但是，Song 等僅將該算法和ADI 算法進行了結(jié)合，依然存在ADI 算法在大時間步下計算精度較低的問題。

在模擬開放或半開放的電磁輻射和散射問題時，為了用有限的網(wǎng)格空間模擬無限空間，必須在計算區(qū)域的截斷邊界處給出吸收邊界條件（absorbing boundary condition，ABC）。在實際計算中，完美匹配層（perfectly matched layer，PML）是目前最常用也是吸收效果最好的一種ABC 算法。其中，復(fù)頻率偏移PML（complex frequency-shifted PML，CFS-PML）具有實現(xiàn)簡單、對倏逝波吸收效果好等優(yōu)點，深受學(xué)者們的歡迎。此外，由于高階PML（higher-order PML，HO-PML）的吸收效果優(yōu)于一階PML（first-order PML，F(xiàn)O-PML），近年來CFS-PML，特別是HO-CFS-PML 受到越來越多的關(guān)注以及廣泛的研究[20-21]。但是目前缺少對結(jié)合RKADE算法的PML 算法研究，這限制了RKADE 算法在開放空間下的應(yīng)用。

本文提出了一種截斷等離子體介質(zhì)的二維無條件穩(wěn)定HO-RKADE-SS-CFS-PML 算法，將RKADE算法與SS-FDTD 算法相結(jié)合，消除了CFL 穩(wěn)定性條件的限制；并且通過采用ADE 方法在計算域邊界處實現(xiàn)HO-CFS-PML。通過數(shù)值算例驗證該算法的有效性，并與FO-RKADE-SS-CFS-PML 算法和HORKADE-ADI-CFS-PML 算法相比較。

1 HO-RKADE-SS-CFS-PML 算法公式推導(dǎo)

截斷非磁化等離子體介質(zhì)的HO-CFS-PML 在二維TMz 波下的麥克斯韋微分方程組可以寫為：

式中：Ez為電場分量；Hx、Hy為磁場分量；ε0為真空下的介電常數(shù)；μ0為真空下的磁導(dǎo)系數(shù)；Jz為極化電流密度分量；v 為電子碰撞頻率；ωp為等離子體角頻率；Sη（η=x，y）為HO-CFS-PML 的拉伸坐標(biāo)變量，定義為：

為了得到Jz分量的值，式（4）采用二階龍格庫塔格式進行離散，以第1 個子時間步為例，可以得到：

為了保證局部截斷誤差盡可能小，各參數(shù)需要滿足如下條件：

將式（12）代入式（6）中，可以得到第1 個子時間步中的電磁場分量公式為：

式中：Fzxn和Gyxn（n=1，2）均為輔助變量。下面以Fzxn為例對輔助變量的更新公式進行推導(dǎo)：

使用克蘭克-尼科爾森（Crank-Nicolson）格式對式（16）進行離散，則其表達式為：

2 數(shù)值算例分析

為了展示HO-RKADE-SS-CFS-PML 算法的有效性，采用了一個二維空間下的數(shù)值算例。在該算例中，激勵源設(shè)為調(diào)制高斯脈沖，最大頻率為70 GHz，中心頻率為35 GHz。每個波長下的采樣點為100。計算域為100Δx×100Δy 的二維網(wǎng)格空間，空間步長取Δx=Δy=Δ=4.3×10-5m。計算域邊界設(shè)置10 層的PML 以截斷非磁化等離子體介質(zhì)。傳統(tǒng)FDTD 算法的最大時間步長為=10-13s，而無條件穩(wěn)定的算法時間步長為Δt=CFLN·，CFLN 為正數(shù)。等離子體參數(shù)為ωp=2π×28.7 Grad/s，v=20 GHz。觀察點設(shè)置在計算域左下角，距離兩邊PML 僅1 個空間單位，算例的模型圖如圖1 所示。

圖1 數(shù)值算例草圖Fig.1 Sketch of numerical example

相對反射誤差定義為：

圖2 各CFS-PML 算法的相對反射誤差Fig.2 Relative reflection errors of different CFS-PML algorithms

由圖2（a）可知，在CFLN=1 時，HO-RKADE-SSCFS-PML 算法擁有最小的相對反射誤差，為-116 dB，但相比其他算法優(yōu)勢不夠明顯。由圖2（b）可知，隨著CFLN 的增大，RKADE-ADI-CFS-PML 算法的相對反射誤差會明顯增大，而HO-RKADE-SS-CFS-PML 算法則基本保持不變，其BRRE 值在不同CFLN 的時候均在-116 dB 附近。這說明HO-RKADE-SS-CFS-PML算法不僅擁有更好的電磁波吸收效果，而且不會隨著CFLN 的增大而減弱，因此，該算法在CFLN 較大時會有更明顯的優(yōu)勢。

表1 給出了各個PML 算法所需的CPU 時間、相對傳統(tǒng)CFS-PML 算法的時間減少率以及BRRE 值。

表1 各CFS-PML 算法的CPU 時間、時間減少率以及BRRETab.1 CPU time，time reduction and BRRE of different CFS-PML algorithms

由表1 可知，RKADE-SS-CFS-PML 和RKADEADI-CFS-PML 算法都可以通過增大時間步長來減少計算時間，當(dāng)CFLN=8 時，各算法消耗的時間相差不大，均能減少80%左右的時間。但是此時HO-RKADESS-CFS-PML 算法的BRRE 為-115 dB，遠遠小于其他的算法。由此說明，HO-RKADE-SS-CFS-PML 算法在大時間步長下?lián)碛懈玫碾姶挪ㄎ招省?/p>

算法的穩(wěn)定性可以通過波形的長時間表現(xiàn)[21]來證明。圖3 所示為在CFLN=8 時，HO-RKADE-SSCFS-PML 算法的Ez、Hx和Hy分量。該仿真的時間擴大為上述算例的1 000 倍。由圖3 可知，各分量的波形均長期收斂于零點，證明該算法具有無條件穩(wěn)定性，并且PML 具有有效的吸收效果。

圖3 CFLN=8 時電場和磁場的波形Fig.3 Waveforms of electric field and magnetic field when CFLN=8

3 結(jié)論

本文提出了一種無條件穩(wěn)定HO-RKADE-SSCFS-PML 算法，用于截斷非磁化等離子體。數(shù)值算例結(jié)果表明：

（1）該算法克服了CFL 穩(wěn)定條件的限制，在CFLN=8 時可以節(jié)省79.2%的計算時間，與傳統(tǒng)RKADECFS-PML 算法相比，提高了計算效率。

（2）在CFLN 相同時，本算法的BRRE 值小于FO-RKADE-SS-CFS-PML 算法和2 種RKADE-ADICFS-PML 算法，證明本文所提出的算法擁有更強的電磁波吸收能力。

（3）與2 種RKADE-ADI-CFS-PML 方法不同的是，HO-RKADE-SS-CFS-PML 的吸收性能不會隨著CFLN 的增加而產(chǎn)生明顯的下降。在CFLN=8 時，本算法的BRRE 值可以達到-115 dB，遠遠大于HORKADE-ADI-CFS-PML 算法的-88 dB，證明其在大時間步長下?lián)碛懈蟮膬?yōu)勢。