高 林，謝國(guó)大，黃志祥，許 杰，吳先良

(安徽大學(xué) 計(jì)算智能與信號(hào)處理教育部重點(diǎn)實(shí)驗(yàn)室，安徽 合肥 230039)

隨著新材料的出現(xiàn)和應(yīng)用，與頻率有關(guān)的特性材料在寬頻帶上的精確建模引起了現(xiàn)代計(jì)算電磁學(xué)領(lǐng)域研究人員的廣泛關(guān)注[1-2]，研究人員開(kāi)發(fā)了很多數(shù)值方法，以精確描述色散介質(zhì)的電磁特性.基于Yee網(wǎng)格的時(shí)域有限差分 (finite-difference time-domain，簡(jiǎn)稱 FDTD) 算法，將Maxwell方程在時(shí)間和空間上直接差分，具有簡(jiǎn)單、直觀、高精度等優(yōu)點(diǎn)，F(xiàn)DTD是目前使用較多的電磁計(jì)算數(shù)值算法[3].傳統(tǒng)FDTD算法采用的是顯式差分，時(shí)間步長(zhǎng)必須滿足CFL(courant-friedrich-levy)穩(wěn)定條件，這個(gè)條件要求電磁波在數(shù)值網(wǎng)格中的計(jì)算傳播速度不小于物理傳播速度，時(shí)間步長(zhǎng)的最小值由空間網(wǎng)格尺寸決定.在等離子體結(jié)構(gòu)仿真[4]中，須采用高分辨率，然而高分辨率須網(wǎng)格尺寸小、時(shí)間步長(zhǎng)小，這樣就導(dǎo)致仿真時(shí)間長(zhǎng).為了消除這些限制，T.Namiki提出了無(wú)條件穩(wěn)定的交替方向隱式時(shí)域有限差分 (alternating direction implicit finite-difference time-domain，簡(jiǎn)稱ADI-FDTD)算法[5]，此算法把傳統(tǒng)FDTD算法中一個(gè)時(shí)間步分解為2個(gè)子時(shí)間步,然后交替采用顯式和隱式差分得到電場(chǎng)和磁場(chǎng)分量的時(shí)域步進(jìn)方程.ADI-FDTD算法因具有2階精度和無(wú)條件穩(wěn)定特性，在計(jì)算電磁學(xué)中得到了廣泛應(yīng)用[6-8].

筆者提出一種用于精確模擬波在色散金屬結(jié)構(gòu)中傳播的3維交替方向隱式時(shí)域有限差分算法，用廣義關(guān)鍵點(diǎn)(critical points，簡(jiǎn)稱CPs)模型模擬等離子體結(jié)構(gòu)、研究金屬材料的電磁特性.采用輔助微分方程(auxiliary differential equation，簡(jiǎn)稱ADE)[9]描述電極化與電場(chǎng)間的本構(gòu)關(guān)系，降低算法的復(fù)雜度和應(yīng)用難度.最后，通過(guò)數(shù)值仿真驗(yàn)證算法的有效性.

1 在色散介質(zhì)中的ADI-FDTD算法

1.1 色散介質(zhì)中ADI-FDTD算法的Maxwell方程

差分形式的Maxwell方程組用矩陣表示為

(1)

其中：D為電通量密度；E為電場(chǎng)強(qiáng)度；H為磁感應(yīng)強(qiáng)度；μ0為真空中的磁導(dǎo)系數(shù)；t為時(shí)間；矩陣A，B分別為

(2)

采用文獻(xiàn)[5]中的ADI-FDTD算法，建立電通量密度D的更新方程.

在ADI-FDTD算法中，需要將n→n+1的時(shí)間步分解為n→n+1/2和n+1/2→n+1兩個(gè)子時(shí)間分步.在n→n+1/2子時(shí)間步，對(duì)式(1)的E和H在時(shí)間和空間上進(jìn)行中心差分離散，分別得到

(3)

(4)

將式(2)中的矩陣A，B代入式(3)，可得

(5)

將式(2)中的矩陣A，B代入式(4)，可得

(6)

在n+1/2→n+1子時(shí)間步，對(duì)式(1)中E和H在時(shí)間和空間上進(jìn)行中心差分離散，分別得到

(7)

(8)

將式(2)的矩陣A，B代入式(7)，可得

(9)

將式(2)中的矩陣A，B代入式(8)，可得

(10)

(11)

(12)

(13)

(14)

(15)

(16)

由文獻(xiàn)[10]可知，在色散介質(zhì)中D和E間的關(guān)系為

(17)

Pm(ω)=ε0χ(ω)E(ω)，

(18)

(19)

(20)

1.2 CPs模型的輔助微分方程

近年來(lái)，多個(gè)色散模型[11-15]被提出.筆者利用CPs模型研究金屬材料在光學(xué)頻率下的電磁特性.文獻(xiàn)[16]給出的磁化率表達(dá)式為

(21)

其中：φm為相位；ω為等離子體頻率；Cm為振幅；Γm為增寬；Ωm為間隙能；ejφm為時(shí)諧因子，其表達(dá)式為

ejφm=cosφm+jsinφm，e-jφm=cosφm-jsinφm.

(22)

將式(22)代入式(21)，可得

(23)

將式(23)代入式(18)，得到的電極化強(qiáng)度P的表達(dá)式為

(24)

其中

利用時(shí)域和頻域間的轉(zhuǎn)換關(guān)系，將式(24)轉(zhuǎn)化為時(shí)域，得到的微分方程為

(25)

為避免ADI-FDTD算法出現(xiàn)不穩(wěn)定情況[17]，筆者增加輔助變量Q, 用2個(gè)1階微分方程替代2階輔助微分方程.替代式(25)的2個(gè)1階微分方程為

(26)

(27)

方程(26)，(27)在n+1/4時(shí)間步上經(jīng)時(shí)域有限差分和平均離散后，分別得到

(28)

(29)

由式(28)，(29)，得到的P和Q在n+1/2時(shí)間步的差分方程為

(30)

(31)

同理，可得P和Q在n+1時(shí)間步的差分方程為

(32)

(33)

其中

fm=b1,m+4b2,m/Δt+Δtb0,m/4,f1,m=2/fm,

f2,m=-[b1,m-(4/Δt)b2,m+(Δt/4)b0,m]/fm,f3,m=[a1,m+(Δt/4)a0,m]/fm,

f4,m=[-a1,m+(Δt/4)a0,m]/fm,h1,m=b1,m-4b2,m/Δt,h2,m=b1,m+4b2,m/Δt.

(34)

式(17)中電通量密度D在n+1時(shí)間步的方程為

(35)

式(17)中的電通量密度D在n+1/2時(shí)間步的方程為

(36)

式(17)中的電通量密度D在n時(shí)間步的方程為

(37)

將式(30)代入式(36)，可得

(38)

將式(37)，(38)代入式(11)，得到的第1個(gè)子時(shí)間步Ex的時(shí)域步進(jìn)方程為

(39)

將式(37)，(38)代入式(12)，得到的第1個(gè)子時(shí)間步Ey的時(shí)域步進(jìn)方程為

(40)

將式(37)，(38)代入式(13)，得到的第1個(gè)子時(shí)間步Ez的時(shí)域步進(jìn)方程為

(41)

將式(32)代入式(35)，可得

(42)

將式(37)，(42)代入式(14)，得到的第2個(gè)子時(shí)間步Ex時(shí)域步進(jìn)方程為

(43)

電場(chǎng)的其他分量同理可得.

2 數(shù)值仿真

為了驗(yàn)證ADI-FDTD算法模擬CPs色散介質(zhì)的優(yōu)越性，對(duì)典型的等離子體模型和金屬諧振腔結(jié)構(gòu)進(jìn)行數(shù)值仿真.

2.1 具有多層金屬結(jié)構(gòu)的電磁波傳播特性

圖1 解析法和3種CFLN值的ADI-FDTD算法的傳輸系數(shù)比較

2.2 基于矩形諧振腔結(jié)構(gòu)驗(yàn)證ADI-FDTD算法的準(zhǔn)確性

假定有一充滿色散介質(zhì)的矩形腔體，腔體尺寸為L(zhǎng)x=Ly=Lz=40 mm,網(wǎng)格大小為Δx=Δy=Δz=1 mm，網(wǎng)格涉及的計(jì)算區(qū)域?yàn)?0×40×40，色散模型的相關(guān)參數(shù)為a0=6.989×1018s-2，a1=0,b1=6.989×108s-1，b2=1,b0=0,ε∞=1.脈沖函數(shù)為E(t)=E-(t-t0)2/(Ts)2， 其中Ts=79.6 ps，脈沖峰值出現(xiàn)在t=t0=4Ts時(shí)刻.圖2為FDTD算法和5種不同的CFLN值的ADI-FDTD算法在探測(cè)點(diǎn)(3,4,4)處的電場(chǎng)波形圖，圖3為FDTD算法和5種不同的CFLN值的ADI-FDTD算法在探測(cè)點(diǎn)(3,3,4)處電場(chǎng)波形圖.由圖2，3可知，不同CFLN值的ADI-FDTD算法和普通FDTD算法，在同一點(diǎn)的電場(chǎng)波形基本吻合.

圖2 FDTD算法和5種不同的CFLN值的ADI-FDTD算法在探測(cè)點(diǎn)(3,4,4)處的電場(chǎng)波形圖

圖3 FDTD算法和5種不同的CFLN值的ADI-FDTD算法在探測(cè)點(diǎn)(3,3,4)處電場(chǎng)波形圖

表1為FDTD和ADI-FDTD算法的時(shí)間步、仿真時(shí)間比較，由表1可以看出，CFLN>8時(shí)ADI-FDTD算法比FDTD快.仿真結(jié)果表明：FDTD，ADI-FDTD算法占用內(nèi)存大小分別為8.608 6，8.767 3 MB，即二者幾乎相同.

表1 FDTD和ADI-FDTD算法的時(shí)間步、仿真時(shí)間比較

3 結(jié)束語(yǔ)

筆者提出了一種用于精確模擬波在色散金屬結(jié)構(gòu)中傳播的3維無(wú)條件穩(wěn)定的ADI-FDTD算法.使用ADI-FDTD算法研究矩形諧振腔中的電磁波，利用基于廣義CPs模型建立的輔助微分方程研究金屬材料在光學(xué)頻率下的電磁特性.ADI-FDTD與FDTD算法的比較結(jié)果表明，二者結(jié)果吻合較好， 但ADI-FDTD算法具有高的效率.因此，ADI-FDTD算法是替代FDTD的有效算法.