一種提高內存使用效率的時域有限差分算法

張 品 陳亦望 傅 強

（解放軍理工大學工程兵工程學院，江蘇 南京 210007）

1.引 言

時域有限差分方法（FDTD）已經廣泛應用于電磁散射計算領域［1－3］。然而，傳統FDTD方法所需計算資源量非常大，并且受到CFL穩定條件的限制，使得計算效率不高，計算時間較長，這一切都給使用FDTD方法進行電磁仿真計算帶來困難。為了突破CFL穩定條件的限制，Takefumi Namiki提出了交替隱式時域有限差分法（ADI－FDTD）［4－5］，提高了計算速度，但是計算公式較為復雜，占用內存較大。近來，J.Shibayama等人提出的局部一維時域有限差分法（LOD－FDTD）也被應用到電磁計算中［6－7］。由于 LOD－FDTD 方法比 ADI－FDTD 方 法的計算公式簡單，占用內存較少，計算效率得到提升［6］。

雖然LOD－FDTD方法擺脫了CFL條件限制，具有更高的計算效率，但是由于其求解過程中仍然需要解方程組，所以LOD－FDTD方法計算所需內存量依然比傳統FDTD方法要大，相對有限的計算機資源給LOD方法的應用帶來了制約。

減縮時域有限差分法（R－FDTD）是一種有效減少內存使用量的方法［8］。R－FDTD方法依據麥克斯韋方程組中三個電場方程與三個磁場方程之間通過散度方程建立的相互聯系，避免了部分場量值的存儲，使得原本需要儲存的麥克斯韋方程組6個方程量減少為4個方程量。循環使用一維數組儲存場量代替二維數組儲存（二維情況），或二維數組代替三維數組儲存（三維情況），這樣就達到了平均節約33%內存的目的［8－9］。

在文獻［9］中，R－FDTD已經用于改進二維ADI－FDTD方法，但該文假設在無源區電場或磁場仍為零散度關系，這將導致大時間步長計算時計算的發散。不同于利用傳統FDTD計算時電、磁場分量散度關系為零，經過研究發現，在無源場區環境下利用局部一維時域有限差分法（LOD－FDTD）計算時，電、磁場分量的散度非零，并推導了該散度關系的具體表達式。基于此散度關系，提出一種新的時域有限差分方法——局部一維減縮時域有限差分法（LOD－R－FDTD），并給出了具體計算格式。該方法使電、磁場的計算不僅能夠擺脫CFL條件限制、減少計算時間，而且節省平均達1／3的內存使用量。通過數值仿真實驗證明了新方法的準確性和有效性。

2.LOD－FDTD計算時電、磁場分量的空間散度

僅考慮二維TE波模型，但是同樣可以推廣到三維情況［10－11］。

2.1 二維TE波LOD－FDTD公式推導

對于二維TE波情況，有Ez＝Hx＝Hy＝0。Maxwell方程組在各向同性、無損耗介質空間的形式為

在直角坐標系中可表示為

式中U＝ ［Ex，Ey，Hz］T；

對式（1）應用近似：（?U）／?t≈ （Un＋1－Un）／Δt，Un＋（1／2）≈ （Un＋1＋Un）／2，得到Crank－Nicolson形式

式中I為單位矩陣。式（2）近似表示為

式（3）可通過兩個時間步驟進行計算：

子時間步1：n→n＋1／2

子時間步2：n＋1／2→n＋1

將各場分量代入式（4）和（5）得到

子時間步1：n→n＋1／2

子時間步2：n＋1／2→n＋1

文獻［12］已證明，以上方法無條件穩定。

2.2 LOD－FDTD計算時的電、磁場分量空間非零散度關系

基于經典電磁場計算理論，在無空間電荷的區域電、磁場散度為零［8］，即

將傳統FDTD Yee's各場量計算式代入式（12）、（13），可證明當電、磁場的各場量采用傳統FDTD格式計算時，式（12）、（13）的空間離散形式成立。但是，當采用 LOD－FDTD 法計算時，式（12）、（13）的空間離散形式不再成立。下面以TE波為例，推導LOD－FDTD算法下電、磁場分量的空間散度關系。

在時間步n＋1／2時，

對于電場散度有

將式（6）和（7）代入式（14），有

假設在n＜0的時間步，各場分量均為零，得到電場空間散度關系為

在空間網格（i，j）處，n＋1／2時間步的Ey由式（16）計算，具體形式為

在時間步n＋1時，

對于電場散度有

將式（9）和（10）代入式（18），有

將式（16）代入式（19），得到電場空間散度關系為

在空間網格（i，j）處，n＋1時間步的Ex由式（20）計算，具體形式為

同理，對于二維TM波可得到類似的電、磁場散度的空間表達形式。

2.3 LOD－FDTD與R－FDTD結合

在二維情況下以TE波為例，設計算區域為Nx×Ny的無電荷矩形空間區域。采用LOD－FDTD計算時，所計算場量為Ex、Ey和Hz三個分量。在n＋1／2和n＋1時間步共需6個二維數組存儲整個計算區域的分量值。將LOD－FDTD與R－FDTD結合，應用非零散度關系式，在n＋1／2時間步，在整個計算區域存儲Ex、Hz分量，不在整個計算區域存儲Ey分量，而僅用一維數組（i）循環存儲j＝c（c＝1，…，Ny）的網格上的Ey分量；在n＋1時間步，在整個計算區域存儲Ey、Hz分量，不在整個計算區域存儲Ex分量，而僅用一維數組（j）循環存儲i＝c（c＝1，…，Nx）的網格上的Ex分量。這樣一來，僅需要4個二維數組存儲整個計算區域的分量值，LOD－R－FDTD方法所需內存僅為 LOD－FDTD 所需內存的2／3。

理論上可以任意選定一個電場量不使用二維數組存儲其在整個計算空間的數值，而是利用非零散度關系式（式（16）、（20））得到某個計算列上的數值，并存儲在一維數組中，在隨后的計算過程中沿著與此計算列垂直的方向循環更新。但是，如果在計算區域中選定空間距離最長的方向作為一維數組更新方向，則可以最大程度地節約內存。

仍以二維 TE波情況為例。使用 LOD－RFDTD算法，選定（i）、（j）兩個分量不需要使用二維數組儲存所有計算空間的場值，而僅需要一維數組儲存分別沿x軸、y軸方向上不斷更新的場值和.為了更清晰地描述該算法的計算過程，具體算法的偽代碼見圖1。

3.試 驗

為了評估LOD－R－FDTD方法的有效性和準確性，本文采用提出的LOD－R－FDTD方法分別計算TE波的二維自由空間和三維自由空間傳播情況，并將計算結果與傳統FDTD以及LOD－FDTD方法的計算結果相比較。

3.1 二維TE波波導模型

平行波導模型如圖2（a）所示。波導的兩平行表面均為理想電導體（PEC邊界條件）。在波導兩端使用8層完美匹配層（PML吸收層）。計算空間共劃分為75×20的網格，在自由空間網格大小為1.0mm×1.0mm.波導內有兩個帶有狹縫的薄膜層。電介質薄膜厚度為0.5mm，εr＝2.0，μr＝2.0，σ＝15.0S／m.狹縫寬度為0.3mm.薄膜處的網格大小為0.1mm×0.1mm.

使用傳統FDTD方法時，由于CFL條件的限制，時間步長設為0.23ps.當使用LOD－FDTD和LOD－R－FDTD方法時，時間步長設為3.33ps.設左端電場垂直分量為激勵源，右端電場垂直分量為輸出信號觀察點。激勵脈沖波方程為

觀察點處的振幅歸一化輸出信號如圖2（b）所示。觀 察 可 知 LOD－FDTD、LOD－R－FDTD 與 傳 統FDTD方法的計算結果符合較好。LOD－R－FDTD與LOD－FDTD方法計算結果一致，說明LOD－RFDTD方法繼承了LOD－FDTD的準確性。表1提供了不同方法仿真計算時使用計算資源的比較。使用LOD－R－FDTD 和 LOD－FDTD 方法的時間步長均為傳統FDTD方法的5倍，dtLOD－R＝5dt.比較后可發現，使用LOD－R－FDTD方法的CPU計算時間降低為傳統FDTD方法的72.7%。LOD－R－FDTD方法使用的內存為LOD－FDTD方法的78.9%。比較 結 果 表 明：LOD－R－FDTD 方 法 繼 承 了 LODFDTD方法大時間步長的優點，減少了計算時間；同時也繼承了R－FDTD方法的優點，大大減少了計算所需內存。

圖1 LOD－R－FDTD算法偽代碼

圖2 二維TE波模型及計算比較

表1 不同算法間比較

3.2 三維計算模型

計算一個2cm×1cm×2cm的矩形真空腔體，網格設為Δx＝Δy＝Δz＝1mm，將計算空間劃分為20×10×20個網格。時間步dtLOD－R＝5.8dtCFL.仿真時間長度為8000時間步。在4個腔體壁處采用8層PML吸收邊界。圖3為原始FDTD、LODFDTD和LOD－R－FDTD的時域電場波形。

由圖3（a）可以看到：LOD－R－FDTD 方法是無條件穩定的。由圖3（b）可知，LOD－R－FDTD方法與其他方法的計算結果具有較好的一致性。本例計算中 LOD－R－FDTD 計算所用內存比 LOD－FDTD計算所用內存節約了13.7%。

圖3 電場時域圖像

4.結 論

本文證明了即使在無源區域，局部一維時域有限差分方法（LOD－FDTD）的電、磁場量也不滿足零散度關系，同時推導了該散度關系的具體表達式。基于推導出的非零散度關系，提出了一種新的局部一維縮減時域有限差分方法（LOD－R－FDTD），并給出了具體計算公式。通過二維、三維問題的實驗計算，證明該方法不僅保留了LOD－FDTD的無條件穩定性、增大了時間步長、減少了總計算時間，而且使計算所需內存比LOD－FDTD減少平均可達1／3。該方法與LOD－FDTD的計算結果一致，具有良好的計算穩定性。提出的新方法可用于計算區域較大、計算內存相對有限的電磁計算問題。

［1］SINGH G，TAN E L，CHEN Z N.A split－step FDTD method for 3－D Maxwell’s equations in general anisotropic media［J］.IEEE Transactions on Antennas and Propagation，2010，58（11）：3647－3657.

［2］湯 煒，李清亮，焦培南，等.ADI－FDTD在二維散射問題中的應用［J］.電波科學學報，2003，18（6）：620－624.

TANG Wei，LI Qingliang，JIAO Peinan，et al.Twodimension scattering analysis using ADI－FDTD method［J］.Chinese Journal of Radio Science，2003，18（6）：620－624.（in Chinese）

［3］黃斌科，蔣延生，汪文秉.二維無條件穩定時域有限差分方法［J］.電波科學學報，2003，18（5）：534－539.

HUANG Binke，JIANG Yansheng，WANG Wenbing.Two－dimensional unconditionally stable FDTD method［J］.Chinese Journal of Radio Science，2003，18（5）：534－539.（in Chinese）

［4］NAMIKI T.A new FDTD algorithm based on alternating－direction implicit method［J］.IEEE Transactions on Microwave Theory and Techniques，1999，47（10）：2003－2007.

［5］ROUF H K，COSTEN F，GARCIA D G.Reduction of numerical errors in frequency dpendent ADI－FDTD［J］.IEEE Electronics Letters，2010，46（7）：489－490.

［6］SHIBAYAMA J，MURAKI M，YAMAUCHI J，et al.Efficient implicit FDTD algorithm based on locally one－dimensional scheme［J］. Electronics Letters，2005，41（19）：1046－1047.

［7］TAN E L.Acceleration of LOD－FDTD method using fundamental scheme on graphics processor units［J］.IEEE Microwave and Wireless Components Letters，2010，20（12）：648－650.

［8］KONDYLIS G D，FLAVIIS F D，POTTIE G J，et al.A memory－efficient formulation of the finite－difference time－domain method for the solution of Maxwell equation［J］.IEEE Trans.MTT，2001，（7）：1310－1320.

［9］LIU B，GAO B Q，TAN W，et al.A new FDTD Aalgorithm－ADI／R－FDTD［C］／／Electromagnetic Compatibility，2002 3rd International Symposium on，May 21－24，2002，IEEE 0－7803－7277－8／02：250－253.

［10］AHMED I，CHUA E K，LI E P，et al.Development of the three－dimensional unconditionally stable LODFDTD method［J］.IEEE Transactions on Antennas and Propagation，2008，56（11）：3596－3600.

［11］TANE L.Unconditionally stable LOD－FDTD method for 3－D Maxwell’s equations［J］.IEEE Microw.Wireless Compon.Lett.，2007，17（12）：85–87.

［12］劉國勝，張國基.高階LOD－FDTD方法的數值特性研究［J］.電子與信息學報，2010，32（6）：1384－1387.

LIU Guosheng，ZHANG Guoji.Study for the numerical properties of the higher－order LOD－FDTD methods［J］.Journal of Electronics ＆ Information Technology，2010，32（6）：1384－1387.（in Chinese）