高 林，謝國大，黃志祥，許 杰，吳先良

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

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

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

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

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

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

(1)

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

(2)

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

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

(3)

(4)

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

(5)

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

(6)

在n+1/2→n+1子時間步，對式(1)中E和H在時間和空間上進(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模型的輔助微分方程

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

(21)

其中：φm為相位；ω為等離子體頻率；Cm為振幅；Γm為增寬；Ωm為間隙能；ejφm為時諧因子，其表達(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)

其中

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

(25)

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

(26)

(27)

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

(28)

(29)

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

(30)

(31)

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

(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時間步的方程為

(35)

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

(36)

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

(37)

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

(38)

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

(39)

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

(40)

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

(41)

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

(42)

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

(43)

電場的其他分量同理可得.

2 數(shù)值仿真

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

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

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

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

假定有一充滿色散介質(zhì)的矩形腔體，腔體尺寸為Lx=Ly=Lz=40 mm,網(wǎng)格大小為Δx=Δy=Δz=1 mm，網(wǎng)格涉及的計算區(qū)域為40×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時刻.圖2為FDTD算法和5種不同的CFLN值的ADI-FDTD算法在探測點(diǎn)(3,4,4)處的電場波形圖，圖3為FDTD算法和5種不同的CFLN值的ADI-FDTD算法在探測點(diǎn)(3,3,4)處電場波形圖.由圖2，3可知，不同CFLN值的ADI-FDTD算法和普通FDTD算法，在同一點(diǎn)的電場波形基本吻合.

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

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

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

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

3 結(jié)束語

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