一種用于非均勻網格劃分的低數值色散的AA-ID-FDTD方法

劉桂英 謝孟洪 趙珉 陳偉軍

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

0 引言

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

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

1 數學公式

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

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

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

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

利用牛頓迭代公式,由(9)可以得到權重系數α的結果。

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

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

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

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

圖3 在相同的空間分辨率系數情況下,相對數值速度與波傳輸角度的關系圖(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 數值結果

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

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

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

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

圖5 采用文中提出的方法和傳統的ID-FDTD方法計算的瞬態電場值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的網格時,采用本文方法和傳統ID方法計算出的相對誤差從0.47%降為0.41%,其中本文方法的AA參數取值為ex=0.9119,ey=0.9993。此外,為了驗證本文方法的計算精度,用離散傅氏變換求得觀測點的諧振頻率,諧振頻率的相對誤差定義為:|fProposedmethodorID-FDTD-freferencesolution|/freferencesolution×100%,其中freferencesolution是采用均勻網格300×300,由現存的ID-FDTD方法求得的諧振頻率。 fproposedmethod和fID-FDTD分別是在如表所示非均勻網格劃分的情況下,由本文方法和現存ID-FDTD方法所求得的諧振頻率。如表1所示,與現存的方法相比,本文方法求得的相對誤差更小。

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

3 結論

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