雙色散模型的辛?xí)r域有限差分算法?

王輝1)3) 黃志祥1)? 吳先良1)2) 任信鋼3) 吳博1)

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

2)(合肥師范學(xué)院電子信息工程學(xué)院,合肥 230061)

3)(香港大學(xué)電機(jī)電子工程學(xué)院,香港,薄扶林道)

1 引 言

色散媒質(zhì)在電磁場(chǎng)領(lǐng)域有著廣泛的應(yīng)用[1?3],比如等離子體,雙負(fù)媒質(zhì),吸波材料等都屬于色散媒質(zhì).如果媒質(zhì)的介電常數(shù)或者磁導(dǎo)率是頻率的函數(shù),那么這種媒質(zhì)被稱為色散媒質(zhì).現(xiàn)實(shí)中的很多物質(zhì)都是色散的,比如海水,雪等.由于色散媒質(zhì)的復(fù)雜特性,使得解析法在求解色散媒質(zhì)的電磁場(chǎng)問題上變的十分困難.因此對(duì)它們的研究只能借助于數(shù)值方法以及科學(xué)實(shí)驗(yàn).在諸多數(shù)值方法中,時(shí)域有限差分方法(FDTD)[4,5]是研究色散媒質(zhì)電磁特性的有力工具.因?yàn)镕DTD方法具有寬帶特性,一次運(yùn)行可以得到寬頻結(jié)果.電通量密度D和電場(chǎng)E,磁通量密度B和電場(chǎng)H之間的關(guān)系反應(yīng)了材料的復(fù)雜性,很多時(shí)候并不簡(jiǎn)單的假設(shè)本構(gòu)關(guān)系式D=εE,B=μH中的ε,μ是一個(gè)恒定的標(biāo)量,ε,μ可以是一個(gè)張量,從而可以描述媒質(zhì)的各項(xiàng)異性.當(dāng)研究非線性媒質(zhì)時(shí)也可以把ε,μ寫成非線性函數(shù),當(dāng)然媒質(zhì)參數(shù)也可以是隨時(shí)間變化的函數(shù).介電常數(shù)也可以是位置的函數(shù),即表征空間媒質(zhì)的非均勻性.本文研究的重點(diǎn)是ε,μ是頻率的復(fù)雜函數(shù)的情形.色散媒質(zhì)能更好的描述現(xiàn)實(shí)世界,因此我們有必要對(duì)色散媒質(zhì)的FDTD方法展開研究.為了仿真色散媒質(zhì),研究者們提出了多種色散模型,比如Debye,Drude,Lorentz以及Drude-Lorentz模型[6?8].

辛?xí)r域有限差分算法(SFDTD)[9?11]在穩(wěn)定性,長時(shí)間迭代以及色散誤差等方面比普通FDTD算法有很大的優(yōu)勢(shì).構(gòu)造SFDTD算法需要將包含媒質(zhì)參數(shù)的Maxwell方程寫成矩陣形式,在計(jì)算時(shí)間演化矩陣exp(?tL)時(shí),要求L分裂成的兩個(gè)矩陣U,V滿足矩陣指數(shù)exp(cl?tU),exp(dl?tV)收斂. 但是由于色散媒質(zhì)的電磁參數(shù)相對(duì)復(fù)雜,使得矩陣U,V的構(gòu)造相當(dāng)困難.文獻(xiàn)[12]對(duì)無耗散項(xiàng)的Drude模型構(gòu)造了SFDTD算法.針對(duì)更具有普遍意義的有耗散項(xiàng)Drude-Lorentz模型,本文建立了可以處理雙色散模型的SFDTD算法.

2 理論和數(shù)值推導(dǎo)

2.1 Drude-Lorentz色散模型的辛算法框架

Drude模型是比較常用的色散模型,它可以很好地?cái)M合一些材料的電磁參數(shù),有時(shí)候?yàn)榱诉M(jìn)一步提高對(duì)材料電磁參數(shù)的擬合,需要引入Lorentz修正項(xiàng),即通過Drude和Lorentz模型的疊加來提高擬合的精度.此時(shí)復(fù)的相對(duì)介電常數(shù)ε及復(fù)的相對(duì)磁導(dǎo)率μ可表示為

其中,ε∞及μ∞分別為無限頻率的相對(duì)介電常數(shù)和相對(duì)磁導(dǎo)率,ωeD和ωhD分別為電等離子體角頻率和磁等離子體角頻率,γeD及γhD分別為電碰撞頻率和磁碰撞頻率,?ε及?μ為引入Lorentz修正項(xiàng)后介電常數(shù)及磁導(dǎo)率的相對(duì)改變量,?eL及?hL為非損耗諧振頻率,ΓeL及ΓhL是損耗系數(shù).為簡(jiǎn)單起見,這里只考慮一項(xiàng)修正,對(duì)于多項(xiàng)修正可類推.

無耗,各項(xiàng)同性,均勻介質(zhì)空間中的Maxwell方程組的頻域形式為

考慮到本構(gòu)關(guān)系,引入如下輔助變量:

因此,(3),(4)式可改寫為

結(jié)合逆傅里葉變換,(5)–(12)式對(duì)應(yīng)的時(shí)域形式為

方程(13)–(20)可進(jìn)一步改寫成矩陣形式

其中,G=(H E KLJLP Q KDJD)T,

（1）科教部門負(fù)責(zé)科研經(jīng)費(fèi)的管理，負(fù)責(zé)審核和控制各項(xiàng)科研經(jīng)費(fèi)的支出。各專項(xiàng)科研經(jīng)費(fèi)劃撥我院后，科研科根據(jù)醫(yī)院文件規(guī)定，分級(jí)別按比例給予匹配一定的經(jīng)費(fèi)，并按項(xiàng)目負(fù)責(zé)人立戶，一題一本（科研經(jīng)費(fèi)使用記錄本），標(biāo)明項(xiàng)目名稱、項(xiàng)目負(fù)責(zé)人、允許使用經(jīng)費(fèi)的額度。

方程(21)可表示為

其中,

為簡(jiǎn)化記號(hào),引入

引入辛算法近似,考慮時(shí)間步進(jìn)從t=0到t=?t,結(jié)合(24)式有

考慮到U及V不可交換性,利用基本辛矩陣的乘積,矩陣算符可近似為

其中,對(duì)應(yīng)參數(shù)的定義與普通辛算法一致.盡管Uα≠0及Vα0(α ≥ 2),指數(shù)算符exp(dl?tV)及exp(cl?tU)可由表達(dá)式(31),(32)定義并計(jì)算:

至于空間一階偏導(dǎo)數(shù)的離散可借助于四階交錯(cuò)差分,即

因而,結(jié)合時(shí)間及空間離散,最終建立了空間四階時(shí)間任意階的SFDTD算法以處理雙色散媒質(zhì),而且本文構(gòu)建的方法可以直接退化為單色散媒質(zhì)及普通媒質(zhì)的SFDTD算法.

3 數(shù)值算例

3.1 雙色散平板模型的傳輸問題

為了驗(yàn)證本文方法的有效性,首先考慮平板的傳輸問題. 以結(jié)構(gòu)為50 nm厚且為雙色散材質(zhì)的平板作為研究對(duì)象,參考貴金屬金的相對(duì)介電常數(shù),設(shè)置該平板的相對(duì)介電常數(shù)為一個(gè)Drude項(xiàng)加一個(gè)Lorentz修正項(xiàng),其對(duì)應(yīng)的參數(shù)分別為:ε∞=1.17152,?ε=2.233994,γeD=3.21×1013rad/s,ωeD=1.37×1016rad/s,ΓeL=2π×6.22×1013rad/s及ωeL=2π×1.31×1015rad/s.另外設(shè)置該平板的相對(duì)磁導(dǎo)率與相應(yīng)的相對(duì)介電常數(shù)相同,即 μ∞=1.17152,?μ =2.233994,γhD=3.21× 1013rad/s,ωhD=1.37×1016rad/s,ΓhL=2π×6.22×1013rad/s,ωhL=2π×1.31×1015rad/s,為了對(duì)SFDTD和普通FDTD算法進(jìn)行比較,計(jì)算區(qū)域設(shè)為300 nm,平板置于計(jì)算區(qū)域的中心,最外層為20個(gè)網(wǎng)格的PML,入射的平面波源為調(diào)制的高斯激勵(lì)

其中 ω =4π×1014rad/s,τ=5.0×10?15s,t0= τ,其他相關(guān)計(jì)算參數(shù)見表1.圖1給出了該平板在125–275 THz頻段透射系數(shù)的精確解,SFDTD解以及普通FDTD解.由圖1可以看出SFDTD(紅色圓點(diǎn)),FDTD(黑色三角形)和精確解(藍(lán)色實(shí)線)的誤差相當(dāng).表1給出了兩種算法結(jié)果對(duì)應(yīng)的內(nèi)存以及計(jì)算時(shí)間的比較.由表1可以看出,在相同精度下,SFDTD可以設(shè)置較大的穩(wěn)定性系數(shù)以及網(wǎng)格步長,從而可以減少計(jì)算時(shí)間.結(jié)果證明本文建立的Drude-Lorentz雙色散模型的辛算法框架是有效的.

圖1 (網(wǎng)刊彩色)透射系數(shù)比較

表1 算法性能比較

3.2 三維U形開口諧振環(huán)(SRRs)的特性分析

U形SRRs環(huán)具有周期性結(jié)構(gòu),它在開關(guān),傳感器等方面有著廣泛的應(yīng)用.作為第二個(gè)例子,利用本文構(gòu)建的色散模型SFDTD算法,計(jì)算分析U形SRRs環(huán)的透射系數(shù),反射系數(shù)以及吸收系數(shù).U形SRRs環(huán)結(jié)構(gòu)如圖2所示,該結(jié)構(gòu)在x,y方向上是周期,且是周期長為p=250 nm的正方形結(jié)構(gòu).其他結(jié)構(gòu)參數(shù)為a=150 nm,w=50 nm,d=75 nm,hg=60 nm,hs=40 nm.下層藍(lán)色區(qū)域?yàn)榻橘|(zhì)基板(GaAs材料介電常數(shù)εr=11.0);上層黃色區(qū)域?yàn)橘F金屬銀,考慮到研究的波長范圍為[1200,2400]nm,用一個(gè)Drude項(xiàng)擬合銀在該波段的介電常數(shù),其相應(yīng)的Drude參數(shù)為ε∞=1.0,ωeD=1.37×1016rad/s,γeD=2.73×1013rad/s.電磁波垂直于SRRs平面入射,其電極化方向平行于y軸,z方向上設(shè)置15層厚度的高階分裂場(chǎng)PML作為吸收邊界.圖3是仿真后得到的透射系數(shù)(黑色實(shí)線),反射系數(shù)(藍(lán)色實(shí)線)和吸收系數(shù)(綠色實(shí)線),而且給出了FDTD算法結(jié)果的對(duì)比.表2是SFDTD和FDTD計(jì)算參量的對(duì)比,可以看出,由于SFDTD算法可以將CFL以及網(wǎng)格步長設(shè)置的較大,從而在計(jì)算時(shí)間和內(nèi)存上比FDTD占優(yōu).

圖2 (網(wǎng)刊彩色)SRRs結(jié)構(gòu)示意圖

圖3 (網(wǎng)刊彩色)垂直入射時(shí)SRRs結(jié)構(gòu)的透射系數(shù),反射系數(shù)和吸收系數(shù)

圖3中,根據(jù)透射系數(shù)的變化可以看出該SRRs結(jié)構(gòu)的諧振頻率對(duì)應(yīng)的波長位于1850 nm附近,吸收峰位于1740 nm附近.

表2 算法性能比較

4 結(jié) 論

SFDTD算法是一種能量守恒,條件穩(wěn)定且能得到寬頻帶特性的有效數(shù)值方法,本文基于矩陣分裂,辛算子技術(shù)以及ADE方法構(gòu)建了Drude-Lorentz電磁雙色散模型的辛算法框架,該算法在時(shí)間上為任意階,可以靈活處理.通過嚴(yán)格的公式推導(dǎo),得到了一個(gè)Drude和一個(gè)Lorentz項(xiàng)疊加的雙色散模型的SFDTD算法迭代公式,對(duì)于多個(gè)Drude或者Lorentz項(xiàng)疊加的色散模型,需要引入更多的中間變量,可以類似推導(dǎo).數(shù)值仿真結(jié)果證明了構(gòu)建的雙色散模型的SFDTD算法的可行性和有效性.該算法為色散材料的研究提供了有力的算法支持.

附錄A 矩陣指數(shù)展開

以(32)式為例,下面可以證明其展開后的計(jì)算結(jié)果是收斂的,展開后矩陣指數(shù)的收斂這是本文構(gòu)建的算法的基礎(chǔ).矩陣U的表達(dá)式為(28)式,根據(jù)矩陣的乘法運(yùn)算,可以得到任意整數(shù)個(gè)U矩陣乘積的結(jié)果:

令

其中I為單位陣.那么B矩陣中的非零元素可計(jì)算得到

簡(jiǎn)單記α=cl?t,于是

同樣可以得到

[1]Ai F,Bai Y,Xu F,Qiao L J,Zhou J 2008 Acta Phys.Sin.57 4189(in Chinese)[艾芬,白洋,徐芳,喬利杰,周濟(jì)2008物理學(xué)報(bào)57 4189]

[2]Liu R,Shi J H,Plum E,Fedotov V,Zheludev N 2012 Acta Phys.Sin.61 154101(in Chinese)[劉冉,史金輝,E.Plum,V.A.Fedotov,N.I.Zheludev 2012物理學(xué)報(bào)61 154101]

[3]Wang W S,Zhang L W,Ran J,Zhang Y W 2013 Acta Phys.Sin.62 184203(in Chinese)[王五松,張利偉,冉佳,張冶文2013物理學(xué)報(bào)62 184203]

[4]Ta flove A,Hagness S C 2005 Computational Electrodynamics:The Finite-Difference Time-Domain Method,third ed.(Boston:Artech House)

[5]Sullivan D M,Hagness S C 2005 Electromagnetic Simulation Using the FDTD Method(New York:IEEE Press)

[6]Prokopidis K P,Tsiboukis T D 2007 Appl.Comput.Eletron.22 287

[7]Deinega A,John S 2012 Opt.Lett.37 112

[8]Vial A,Grimault A-S,Macías D,Barchiesi D,Chapelle M L 2005 Phys.Rev.B 71 085416

[9]Sha W,Huang Z X,Chen M S,Wu X L 2008 IEEE Trans.Antennas Propag.56 493

[10]Hirono T,Lui W,Seki S,Yoshikuni Y 2001 IEEE Trans.Microw.Thory Tech.49 1640

[11]Wang H,Huang Z X,Wu X L,Ren X G 2011 Chin.Phys.B 20 114701

[12]Ren X G,Huang Z X,Wu X L,Lu S L,Wang H,Wu L,Li S 2012 Comput.Phys.Commun.183 1192