基于RKADE-SS-FDTD 方法的無(wú)條件穩(wěn)定高階CFS-PML 算法

李建雄，魏之棟

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

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

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

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

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

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

截?cái)喾谴呕入x子體介質(zhì)的HO-CFS-PML 在二維TMz 波下的麥克斯韋微分方程組可以寫為：

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

為了得到Jz分量的值，式（4）采用二階龍格庫(kù)塔格式進(jìn)行離散，以第1 個(gè)子時(shí)間步為例，可以得到：

為了保證局部截?cái)嗾`差盡可能小，各參數(shù)需要滿足如下條件：

將式（12）代入式（6）中，可以得到第1 個(gè)子時(shí)間步中的電磁場(chǎng)分量公式為：

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

使用克蘭克-尼科爾森（Crank-Nicolson）格式對(duì)式（16）進(jìn)行離散，則其表達(dá)式為：

2 數(shù)值算例分析

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

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

相對(duì)反射誤差定義為：

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

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

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

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

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

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

圖3 CFLN=8 時(shí)電場(chǎng)和磁場(chǎng)的波形Fig.3 Waveforms of electric field and magnetic field when CFLN=8

3 結(jié)論

本文提出了一種無(wú)條件穩(wěn)定HO-RKADE-SSCFS-PML 算法，用于截?cái)喾谴呕入x子體。數(shù)值算例結(jié)果表明：

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

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

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