三維柱坐標系中PML—FDTD算法及其奇點處理方法

歐陽 單志勇

摘 要：將FDTD由直角坐標推廣到三維柱坐標系;為保證其求解的穩定性，給出了更為嚴格的時間步長計算公式，改善了中心軸奇點處理問題;給出了磁場徑向奇點處理的推導公式，邊界條件采用PML完全匹配層。計算機仿真驗證了算法的正確性

關鍵詞：時域有限差分;完全匹配層;電磁輻射;奇點處理

引言

時域有限差分法（FDTD）[1]作為研究電磁場傳播的有效工具，在微波電路的時域分析、天線輻射特性、電磁場散射及生物醫學工程等領域有著廣泛的應用。它是由麥克斯韋方程推得的中心差分方程，對空間和時間步長的選取嚴格，并需要邊界吸收條件來模擬電磁場在無限空間的傳播。尤其關于吸收邊界條件的設置[2]，是時域差分算法的重要一環。FDTD經過多年的發展研究，在處理邊界條件上也越來越完善，其中一種完全匹配層（PML）吸收邊界條件[3]，它可以近似于對各個方向的來波進行完全吸收而無反射.這種邊界處理方法計算效率很高，并且很容易推廣到曲面坐標。由于散射體結構的復雜性，在笛卡兒坐標[4]中采用立方或正方體網格單元來模擬復雜物體表面時會帶來一定的誤差，比如許多生物醫學應用中，柱面幾何是經常遇到的問題，如果FDTD算法可以在柱面坐標系下求解。這就避免了笛卡爾FDTD算法中的數值誤差。因此對 FDTD算法的研究也從直角坐標系推廣到圓柱坐標系[5]。在圓柱坐標系下為保證計算的穩定性，對空間步長和時間步長的要求更加嚴格，沿同一徑向，在角度步長相同時;橫向長度單元是不同的，因此，為了保證數值的收斂性，要取最小的值來計算。以往只研究軸對稱幾何，將三維圓柱問題簡化為二維問題r?z平面[6][7]。然而，對于一般的非軸對稱結構，就需要一個完整的柱面坐標系下三維FDTD格式。但三維圓柱坐標系在中心軸計算會出現奇點問題，這些奇點并不是麥克斯韋方程本身所固有的，而是由微分方程轉變差分方程帶來的，因此處理奇點的各種方法被相繼提出，如基于基數展開法[8]，利用直角坐標系和圓柱坐標系公式轉換法[9]，

本文推導了基于三維擴展柱坐標系的麥克斯韋方程，采用PML作為吸收邊界，同時利用輔助方程重新推到了PML區域內的迭代公式，并給出了處理磁場徑向奇點的方法并推導出相應的公式，電場奇點問題采用文獻[10]。通過編程建模和數值實驗驗證，證明了本文提出處理奇點方法的正確性。

1 柱坐標系下麥克斯韋差分方程及時間步長公式

在擴展坐標下令：

在三維柱面坐標中，與極軸相關的數值奇異性不能用以前開發的級數多項式展開來近似。相反，與極軸相關的場分量都簡單地近似于基于安培定律的表達式，它規范了數值解，并確保在奇異節點附近有一個表現良好的解。我們還注意到，這些公式忽略了外加電流 。如果該區域附近有有源線圈，則應添加電流源。

4三維圓柱下的數值驗證

三維圓柱坐標下，徑向網絡尺寸 ，方位角 ，軸向尺寸 。空間網格數80×80×144，邊界條件采用PML吸收，場計算不涉及遠場計算，當 取不同值時，測試結果如下圖：

圖5，6，7，8，分別表示在自由空間中， 時中心軸上 的值。從圖中可以發現， 誤差較大，隨著 的不斷減小，誤差越來越小，在 時，即 時，誤差可以忽略不計。經過實驗測試表面本文推導的處理磁場奇點 的公式，具有可行性，有效性。

5總結

在三維圓柱坐標系下給出了FDTD有限差分格式并編程實現，并給出了對于圓柱坐標系帶來的奇點問題解決方法，推導出對于解決磁場奇點問題的公式。利用輔助方程法重新推導出了PML中差分公式，通過編寫程序，仿真測試。結果表明了本文處理奇點的可行性，正確性。

參考文獻

[1] 葛德彪. 電磁波時域有限差分方法[M].西安：西安電子科技大學出版社，2011

[2] YEE K S. Numerical solution of initial bounary value problem involving Maxwell’s equations in isotropic media [ J ]. IEEE Trans Antennas and Propagation，1966，14：302 - 307.

[3] BERENGER J P. A perfectly matched layer for the absorption of electro-magnetic vaves [ J ]. Computational Physics，1994，114（2）：185 - 200.

[4] ZHAO L，CANGELLAR IS A C. GT - PML：Generalized theory of perfectly matched layers and its application to the reflectionless truncation of finite-difference time-do-main grids [ J ]. IEEE Trans Microwave Theory and Techniques，1996，44（12）：2555 - 2563.

[5] GHEONJ IAN A，JOBAVA R. Non2uniform conforming mesh generator for FDTD scheme in 3D cylindrical coor-dinate system [A]. IEEE D IPED Proceedings[C]. [ s.l. ]：[ s. n. ]，2000：41 - 44

[6] Y.Chen，R.Mittra，and P.Harms，“Finite-difference time-domain algorithm for solving Maxwell’s equations in rotationally symmetric geometries，” IEEE Trans. Microwave Theory Tech.，vol.44，pp.832-839，June 1996.

[7] D.W.Prather and S.Shi，“Formulation and application of the finite-difference time-domain method for the analysis of axially symmetric diffractive optical elements，” J.Opt. Soc.Amer.A.，vol. 16，pp.1131-1141，1999.

[8] F. Liu，S. Crozier，“An FDTD model for calculation of gradient-induced eddy currents in MRI system，” IEEE Trans Appl Supercond.，vol. 14，pp. 1983-1989，2004.

[9] Zilvinas Kancleris，“Handling of Singularity in Finite-Difference Time-Domain Procedure for Solving Maxwell’s Equations in Cylindrical Coordinate System，” IEEE Trans. Antennas Propagat.，vol. 56，pp. 610-613，2008.

[10] Jieru Chi，Feng Liu，Ling Xia，Tingting Shao，David G. Mason，“An Improved Cylindrical FDTD Method and its Application to Field-tissue Interaction Study in MRI，”IEEE trans on antennas and propagation，vol. 47，no. 2，rebruary 2011