三維柱坐標(biāo)系中PML—FDTD算法及其奇點處理方法

歐陽 單志勇

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

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

引言

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

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

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

在擴(kuò)展坐標(biāo)下令：

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

4三維圓柱下的數(shù)值驗證

三維圓柱坐標(biāo)下，徑向網(wǎng)絡(luò)尺寸 ，方位角 ，軸向尺寸 ?？臻g網(wǎng)格數(shù)80×80×144，邊界條件采用PML吸收，場計算不涉及遠(yuǎn)場計算，當(dāng) 取不同值時，測試結(jié)果如下圖：

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

5總結(jié)

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

參考文獻(xiàn)

[1] 葛德彪. 電磁波時域有限差分方法[M].西安：西安電子科技大學(xué)出版社，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，F(xiàn)eng 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