一種用于非均勻網(wǎng)格劃分的低數(shù)值色散的AA-ID-FDTD方法

劉桂英 謝孟洪 趙珉 陳偉軍

(嶺南師范學(xué)院信息工程學(xué)院,廣東湛江 524048)

0 引言

時域有限差分(FDTD)方法是電磁器件建模的有效方法。然而,由數(shù)值離散方法引起的數(shù)值色散誤差(NDE)影響計算精度,使得仿真結(jié)果偏離真實值。為了降低NDE,精細網(wǎng)格FDTD被提出,這無疑增大了計算量。為了改善數(shù)值色散,高階和ID方法被引進FDTD,但在非均勻網(wǎng)格劃分技術(shù)中仍然存在各向異性數(shù)值色散誤差。為了避免上述問題,人工各向異性(AA)參數(shù)被引進FDTD來抑制各向異性NDE,提高計算精度。因為AA參數(shù)只用來改變自由空間介電常數(shù),所以該方法不會增加計算量。

本文提出了一種基于各向同性色散(ID)和人工各向異性(AA)參數(shù)的時域有限差分方法來對電磁波在材料中的傳播進行建模,該方法取代了傳統(tǒng)的以空間各向同性方式對場進行采樣的中心有限差分方法。用AA參數(shù)來修正自由空間介電常數(shù)。文中用一個數(shù)值實例對本文提出的方法進行了驗證,結(jié)果表明對于非均勻網(wǎng)格,本文的AAID-FDTD方法的準確性要高于ID-FDTD。

1 數(shù)學(xué)公式

以文獻[1]提出的ID-FDTD方法公式為基礎(chǔ),引入電場人工各向異性參數(shù),2-D TE波的更新函數(shù)可以表示為:

其中ε0和μ0分別表示自由空間的介電系數(shù)和磁導(dǎo)系數(shù)。εx和εy分別是電場Ex和Ey的人工各向異性參數(shù),Δx和Δy分別是x和y方向的網(wǎng)格尺寸。α是權(quán)重系數(shù),角標i和j表示采樣點。對于EM波,2-DTEz波的更新函數(shù)為:

如文獻[2],定義單個FDTD網(wǎng)格的形狀Z和空間分辨系數(shù)R為:

假設(shè)45φ=°,A=1,(8)可以寫成:

利用牛頓迭代公式,由(9)可以得到權(quán)重系數(shù)α的結(jié)果。

圖1為不加AA參數(shù)時權(quán)重系數(shù)相對于空間分辨率R的關(guān)系,由圖1可以看出對于不同的Z,α趨于不同的常數(shù)。圖2為加上AA參數(shù)時權(quán)重系數(shù)相對于空間分辨率R的關(guān)系。由圖2可以看出對于不同的Z,α都趨于同一常數(shù)(0.1667)。

圖1 沒有人工各向異性參數(shù)情況下,權(quán)重因子與空間分辨率系數(shù)的關(guān)系圖(S=0.8)Fig.1 Weighting factor α versus spatial resolution coefficient R without AA parameters. S=0.8 for all of the cases

圖2 有人工各向異性參數(shù)情況下,權(quán)重因子與空間分辨率系數(shù)的關(guān)系圖(S=0.8)Fig.2 Weighting factor α versus spatial resolution coefficient R with AA parameters. S=0.8 for all of the cases

圖3表明相對數(shù)值速度A是波傳播角度φ的函數(shù)。由圖可以看出,在沒加人工各向異性參數(shù)時,在空間分辨系數(shù)R為常數(shù)的情況下,隨著Z的增大,波傳播角度φ的函數(shù)相對數(shù)值速度也變大。當加入人工各向異性參數(shù)時,可以到,當Z取值為3,6,12時,數(shù)值色散誤差都明顯降低。

圖3 在相同的空間分辨率系數(shù)情況下,相對數(shù)值速度與波傳輸角度的關(guān)系圖(R=10和S=0.8)Fig.3 Variation of relative numerical velocity with wave propagation angle in a 2-D ID-FDTD grid for same spatial resolution coefficient. R=10 and S=0.8 for all of the cases

2 數(shù)值結(jié)果

為了驗證本文所提出方法的性能,以如圖4所示的2-D諧振腔為例,計算其中的傳輸場。以正弦調(diào)制的高斯脈沖(14)為信源:

圖4 ID-FDTD的計算模型Fig.4 Diagram of computational domain for ID-FDTD

式中,Td=1/(2fc),Tc=3Td,fc=5GHz,持續(xù)時間Tf=4ns,時間步長Δt=2×10-13s,仿真時長取20000Δt。在該計算區(qū)域內(nèi)采用300×100的網(wǎng)格,其尺寸大小為0.1×0.1mm2。x軸y軸方向采用良好導(dǎo)電性能的金屬邊界。分別采用傳統(tǒng)的ID-FDTD和本文方法計算其橫向電場。

圖5顯示的是在觀察點o (240Δx, 150Δy),用傳統(tǒng)的ID-FDTD和本文方法計算的歸一化的橫向電場的時域波形圖。結(jié)果顯示兩種方法及仿真結(jié)果一致。

圖5 采用文中提出的方法和傳統(tǒng)的ID-FDTD方法計算的瞬態(tài)電場值Fig.5 Normalized transient electric fields calculated with the conventional ID-FDTD and proposed method

表1為不同方法的計算能力對比。運行該仿真的電腦配置為Intel ? Core (TM) i5-4210 and 8 GB RAM。如表1所示,當采用300×9的網(wǎng)格時,采用本文方法和傳統(tǒng)ID方法計算出的相對誤差從0.47%降為0.41%,其中本文方法的AA參數(shù)取值為ex=0.9119,ey=0.9993。此外,為了驗證本文方法的計算精度,用離散傅氏變換求得觀測點的諧振頻率,諧振頻率的相對誤差定義為:|fProposedmethodorID-FDTD-freferencesolution|/freferencesolution×100%,其中freferencesolution是采用均勻網(wǎng)格300×300,由現(xiàn)存的ID-FDTD方法求得的諧振頻率。 fproposedmethod和fID-FDTD分別是在如表所示非均勻網(wǎng)格劃分的情況下,由本文方法和現(xiàn)存ID-FDTD方法所求得的諧振頻率。如表1所示,與現(xiàn)存的方法相比,本文方法求得的相對誤差更小。

表1 計算二維諧振腔結(jié)果比較Tab.1 Comparison of the Computational Efforts for the 2-D cavity

3 結(jié)論

本文提出了人工各向異性參數(shù)修正的AA-ID-FDTD用于非均勻網(wǎng)格的計算,該方法在沒有增加計算時間的情況下有效提高了計算精度。以2-D腔體橫向電場計算為例,該方法的數(shù)值仿真結(jié)果明顯優(yōu)于現(xiàn)存的計算精度較高的ID-FDTD方法。