雙色散模型的辛時域有限差分算法?

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

1)(安徽大學計算智能與信號處理教育部重點實驗室,合肥 230039)

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

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

1 引 言

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

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

2 理論和數值推導

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

Drude模型是比較常用的色散模型,它可以很好地擬合一些材料的電磁參數,有時候為了進一步提高對材料電磁參數的擬合,需要引入Lorentz修正項,即通過Drude和Lorentz模型的疊加來提高擬合的精度.此時復的相對介電常數ε及復的相對磁導率μ可表示為

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

無耗,各項同性,均勻介質空間中的Maxwell方程組的頻域形式為

考慮到本構關系,引入如下輔助變量:

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

結合逆傅里葉變換,(5)–(12)式對應的時域形式為

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

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

（1）科教部門負責科研經費的管理，負責審核和控制各項科研經費的支出。各專項科研經費劃撥我院后，科研科根據醫院文件規定，分級別按比例給予匹配一定的經費，并按項目負責人立戶，一題一本（科研經費使用記錄本），標明項目名稱、項目負責人、允許使用經費的額度。

方程(21)可表示為

其中,

為簡化記號,引入

引入辛算法近似,考慮時間步進從t=0到t=?t,結合(24)式有

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

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

至于空間一階偏導數的離散可借助于四階交錯差分,即

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

3 數值算例

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

為了驗證本文方法的有效性,首先考慮平板的傳輸問題. 以結構為50 nm厚且為雙色散材質的平板作為研究對象,參考貴金屬金的相對介電常數,設置該平板的相對介電常數為一個Drude項加一個Lorentz修正項,其對應的參數分別為:ε∞=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.另外設置該平板的相對磁導率與相應的相對介電常數相同,即 μ∞=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,為了對SFDTD和普通FDTD算法進行比較,計算區域設為300 nm,平板置于計算區域的中心,最外層為20個網格的PML,入射的平面波源為調制的高斯激勵

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

圖1 (網刊彩色)透射系數比較

表1 算法性能比較

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

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

圖2 (網刊彩色)SRRs結構示意圖

圖3 (網刊彩色)垂直入射時SRRs結構的透射系數,反射系數和吸收系數

圖3中,根據透射系數的變化可以看出該SRRs結構的諧振頻率對應的波長位于1850 nm附近,吸收峰位于1740 nm附近.

表2 算法性能比較

4 結 論

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

附錄A 矩陣指數展開

以(32)式為例,下面可以證明其展開后的計算結果是收斂的,展開后矩陣指數的收斂這是本文構建的算法的基礎.矩陣U的表達式為(28)式,根據矩陣的乘法運算,可以得到任意整數個U矩陣乘積的結果:

令

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

簡單記α=cl?t,于是

同樣可以得到

[1]Ai F,Bai Y,Xu F,Qiao L J,Zhou J 2008 Acta Phys.Sin.57 4189(in Chinese)[艾芬,白洋,徐芳,喬利杰,周濟2008物理學報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物理學報61 154101]

[3]Wang W S,Zhang L W,Ran J,Zhang Y W 2013 Acta Phys.Sin.62 184203(in Chinese)[王五松,張利偉,冉佳,張冶文2013物理學報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