何光峰　遲潔茹

摘要： 針對傳統(tǒng)的時域有限差分法受Courant穩(wěn)定條件的限制，且存在交替方向隱式時域有限差分法數(shù)值色散較大的問題，本文以TE波為例，研究了CrankNicoloson差分方式的近似去耦時域有限差分法基本原理，并對其穩(wěn)定性進行分析，證明該方法是無條件穩(wěn)定。通過數(shù)值仿真，從運行時間和吸收效果方面與傳統(tǒng)的時域有限差分法和交替方向隱式時域有限差分法進行對比。仿真結(jié)果表明，近似去耦時域有限差分法比交替方向時域有限差分法的吸收效果好，但比傳統(tǒng)的時域有限差分法吸收效果差；近似去耦時域有限差分法比交替方向時域有限差分法運行時間長，但比傳統(tǒng)時域有限差分法運行時間短，說明近似去耦時域有限差分法突破了Courant穩(wěn)定條件的限制，且在吸收效果方面比交替方向時域有限差分法好。該研究具有廣闊的應(yīng)用前景。

關(guān)鍵詞： 時域有限差分法； CrankNicoloson差分格式方案； 無條件穩(wěn)定

中圖分類號： O241.8；TB5文獻標識碼： A

1966年，Yee提出傳統(tǒng)時域有限差分（finitedifference timedomain，F(xiàn)DTD）方法，該方法得到廣泛應(yīng)用與發(fā)展[13]。時域有限差分法是用古典顯式差分的方法對麥克斯韋方程進行差分，進而對麥克斯韋方程進行微分求解。運用顯式差分的FDTD方法必須滿足Courant穩(wěn)定條件，即時間步數(shù)的選取受到空間步長（離散網(wǎng)格大?。┑南拗?。為降低數(shù)值色散誤差，在空間步數(shù)選取時必須滿足遠遠小于其采樣頻率的波長的條件，這就使其在解決實際問題的過程中，為得到更加真實的電磁波特性，不得不采用非常小的空間步數(shù)，進而導(dǎo)致計算量過大而無法實現(xiàn)。為克服傳統(tǒng)FDTD方法的缺陷，T.Namiki等人[47]提出了交替方向隱式（alternative direction implicit，ADI）FDTD方法。由于ADIFDTD方法對麥克斯韋方程組的所有場分量進行交替差分，從而實現(xiàn)了無條件穩(wěn)定，克服了傳統(tǒng)FDTD方法必須滿足Courant穩(wěn)定條件的限制，但ADIFDTD方法在一個時間步數(shù)上需要進行兩次迭代，增加了每步計算的時間量和存儲空間，且由于每步差分方程時間不同步，對計算誤差產(chǎn)生了非常大的影響。2003年，Sun G L等人[811]利用CrankNicoloson半隱式差分方式，對麥克斯韋方程組進行差分，使用近似去耦（approximately decoupling，AD）法而提出的CrankNicoloson差分方式的近似去耦（cranknicoloson approximatelydecoupling，CNAD）FDTD方法。CNADFDTD方法在使用CrankNicoloson半隱式差分方式對二維麥克斯韋方程求解過程中，需要求解一個大型稀疏矩陣方程組，從而占用計算機大量內(nèi)存和運行時間，因此采用近似方法求解過稀疏矩陣，將E在n時刻的值代替n+1時刻的值，從而把對稀疏矩陣求解變成對三對角矩陣求解，使計算過程簡單化，并減小所需的計算機內(nèi)存和運行時間?；诖?，本文以橫電（transverseelectric，TE）波為例，介紹了CNADFDTD方法，并對FDTD、ADIFDTD和CNADFDTD三種方法分別進行仿真實驗。實驗結(jié)果表明，在相同條件下，CNADFDTD比ADIFDTD的吸收效果要好，但與FDTD有一定的差距，CNADFDTD的運行時間比ADIFDTD的運行時間長。

1CNADFDTD方法

1.1CNADFDTD理論

二維TE波在x，y和z三個方向上存在Ex，Ey和Hz分量，寫成CN差分方程[911]分別為

Ex|n+1i+1/2，j=Ex|ni+1/2，j+a1ΓyHz|i+1/2，j+1/2/Δy （1a）

Ey|n+1i，j+1/2=Ey|ni，j+1/2-a1ΓxHz|i+1/2，j+1/2/Δx（1b）

Hz|n+1i+1/2，j+1/2=Hz|n+1i+1/2，j+1/2+a2ΓyEx|i+1/2，j+1/Δy-a2ΓxEy|i+1，j+1/2/Δx（1c）

其中

ΓyHz|i+1/2，j+1/2=Hz|n+1i+1/2，j+1/2-Hz|n+1i+1/2，j-1/2+Hz|ni+1/2，j+1/2-Hz|ni+1/2，j-1/2（2）

式中，a1=Δt/2ε，a2=Δt/2μ；Δt為時間步長，Δx，Δy分別為空間步長；ε為介質(zhì)介電常數(shù)，μ為磁導(dǎo)系數(shù)。其迭代過程如下：

1）將式（1c）代入式（1a）和式（1b），消去n+1時刻磁場值Hn+1，可得n+1時刻電場值En+1，實現(xiàn)對電場值迭代。

2）將求得的電場值En+1代入式（1c），更新磁場值Hn+1。其中，在對電場值進行迭代中，將Eny代替En+1y，可將對大型非有限帶寬稀疏矩陣的求解簡化為對三對角矩陣的求解，從而將式（1a）化簡為

1+2b2xEx|n+1i+1/2，j-b2xEx|n+1i+1/2，j-1+Ex|n+1i+1/2，j+1=1-2b2xEx|ni+1/2，j+b2x（Ex|ni+1/2，j-1+Ex|ni+1/2，j+1）+2a1Hz|ni+1，j+1/2-Hz|ni+1，j-1/2-2bxbyEy|ni+1，j+1/2-Ey|ni，j+1/2-Ey|ni+1，j-1/2+Ey|ni，j-1/2（3）

其中，by=cΔt/2Δx，bx=cΔt/2Δy。

最終迭代過程為：對式（3）進行迭代；對式（1b）迭代求解En+1y；對式（3）和式（1b）更新；對式（1c）的磁場值迭代更新。

1.2CNADFDTD的穩(wěn)定性分析

根據(jù)V.Neumann穩(wěn)定性分析，二維情況下的電磁場分量Ex，Ey和Hz[12]分別為

Ex=ψAεnexp[j（kxΔx+kyΔy）]， Ey=ψBεnexp[j（kxΔx+kyΔy）]， Hz=ψCεnexp[j（kxΔx+kyΔy）]（4）

式中，ψA，ψB，ψC分別為各電磁分量初始系數(shù)；ε是增長因子；kx和ky是波常數(shù)；Δx，Δy分別代表沿x方向和沿y方向網(wǎng)格的大小。

將式（4）代入式（1a）～（1c）中，消去系數(shù)得關(guān)于增長因子的方程為

ε-12=-a（ε+1）2（5）

其中，a=Δt2εμΔx2sin2kxΔx/2+Δt2εμΔy2sin2（kyΔy/2）。將a（a>0）代入式（5），可求得ε為

ε=1-a±2aj1+a（6）

由于式（6）中的增長因子滿足ε=1，因此時間步數(shù)無論取何值，CNADFDTD都是穩(wěn)定的。

1.3CNADFDTD在完全匹配層的公式

以有限空間模擬無限空間時，以完全匹配層[1315]（perfectly matched layer，PML）作為吸收邊界，則二維TE波在PML吸收介質(zhì)中，按照CNADFDTD的方法進行差分[1617]，得

Qex（m）Ex|n+1i+1/2，j=Dex（m）Ex|ni+1/2，j+ΓyHz|i+1/2，j+1/2Qey（m）Ey|n+1i，j+1/2=Dey（m）Ey|ni，j+1/2+ΓxHz|i+1/2，j+1/2Hzx|n+1i+1/2，j+1/2=Qhx（m）Hzx|n+1i+1/2，j+1/2-Dhx（m）ΓxEy|i+1，j+1/2Hzy|n+1i+1/2，j+1/2=Qhy（m）Hzy|n+1i+1/2，j+1/2+Dhy（m）ΓyEx|i+1/2，j+1（7）

其中，m取值與左端場分量節(jié)點的空間位置相同；且。

Qex（m）=（2ε0+σy（m）Δt）， Dhy（m）=12μ0+σmy（m）ΔtΔtΔy（8）

將磁場分量代入迭代的點場分量求值中，利用CNADFDTD，先求出En+1x和En+1y，再對磁場值進行迭代。

2數(shù)值分析

為驗證CNADFDTD方法的正確性，以二維平面上的一個點源的輻射場為例，分別對傳統(tǒng)FDTD、ADIFDTD及CNADFDTD方法使用C語言編寫程序，在Visual Studio 2015中進行仿真，對得到的數(shù)據(jù)利用Tecplot 360軟件進行成像處理。

仿真條件：仿真空間為20 mm×20 mm均勻網(wǎng)格，磁場激勵源為正弦波，位于仿真空間中心位置，頻率f=1×1012 Hz，c=30×108 m/s，波長λ=3 mm，空間分辨率Δx=Δy=01 mm，時間步數(shù)Δt=0155 ps，吸收邊界厚度為1 mm。PML中電導(dǎo)率的分布為

σz（z）=σmaxz-z0mdm（9）

式中，d是PML的厚度；z0為PML層靠近FDTD區(qū)的界面位置；m為整數(shù)。為了得到更好的吸收效果，σmax的最佳取值[1820]為

σmax=-（m+1）ln[R（θ）]2ηd（10）

其中，R（θ）是理論反射系數(shù)值，R（θ）=10-6。

通過仿真實驗，得到距離中心點源為5個網(wǎng)格處的磁場分量，3種方法的數(shù)值比較結(jié)果如圖1所示。由圖1可以看出，3種方法的結(jié)果基本相同，驗證了CNADFDTD方法的正確性；通過仿真計算，得到坐標為（189，189）處的反射誤差如圖2所示。由圖2可以看出，在相同條件下，傳統(tǒng)FDTD的吸收效果最好，CNADFDTD的吸收效果比ADIFDTD的吸收效果要好，主要由于ADIFDTD差分等式兩邊的時間不對稱所導(dǎo)致，而CNADFDTD由于只需要一步，且時間上對稱，故其效果比ADIFDTD要好。

3結(jié)束語

本文主要討論無條件穩(wěn)定的CNADFDTD，由于其不受穩(wěn)定條件的限制，所以可通過改變時間步的大小，縮短CPU的運行時間。通過與傳統(tǒng)的FDTD和ADIFDTD兩種方法的對比可以得出，CNADFDTD的吸收效果比ADIFDTD的吸收要好，但運行時間比ADIFDTD的多，且吸收效果相比于FDTD還有一段差距。因此，CNADFDTD是一種值得深入研究和推廣應(yīng)用的時域算法。

參考文獻：

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

[2]Yee K S. Numerical Solution of Initial Boundary Value Problems Involving MaxwellS Equation in Isotropic Media[J]. IEEE Transactions on Antennas and Propagation， 1966， 14（3）： 302307.

[3]Schneider J B. A Selective Survey of the FiniteDifference TimeDomain Literature[J]. IEEE Antennas and propagation， 1995， 37（4）： 3957.

[4]Namiki T. A New FDTD Algorithm Based on AlternatingDirection Implicit Method [J]. IEEE Transaction on Microwave Theory and Techniques， 1999， 47（10）： 20032007.

[5]Namiki T. 3D ADIFDTD MethodUnconditionally Stable TimeDomain Algorithm for Solving Full Vector Maxwells Equations [J]IEEE Transaction on Microwave Theory and Techniques， 2000， 48（10）： 17431748.

[6]Zheng F， Chen Z， Zhang J. Toward the Development of a ThreeDimensional Unconditionally Stable FiniteDifference TimeDomain Method[J]. IEEE Transaction on Microwave Theory and Techniques， 2000， 48（9）： 15501558.

[7]Zheng Hongxing， Leung K W. A Nonorthogonal ADIFDTD Algorithm for Solving Two Dimensional Scattering Problems[J]. IEEE Transaction on Antennas and Propagation， 2009， 57（12）： 38913902.

[8]Sun G L， Trueman C W. Approximate CrankNicolson Schemes for the FiniteDifference TimeDomain Method for TEz Waves[J]. IEEE Transaction on Antennas and Propagation， 2004， 52（11）： 29632972.

[9]Sun G L， Trueman C W. Unconditionally Stable CrankNicolson Scheme for Solving the TwoDimensional Maxwells Equations[J]. Electronics Letters， 2003， 39（7）： 595597.

[10]Xie X， Pan G， Hall S. A CrankNicholsonBase Unconditionally Stable TimeDomain Algorithm for 2D and 3D problem[J]. Microwave and Optical Technology Letters， 2007， 49（2）： 261265.

[11]Li J X， Jiang H L， Zhao X M， et al. Effective CNADand ADEBased CFSPML Formulations for Truncating the Dispersive FDTD Domains[J]. IEEE Antennas and Wireless Propagation Letters， 2015， 14： 12671270.

[12]王秉中. 計算電磁學(xué)[M]. 北京： 科學(xué)出版社， 2002.

[13]Berenger J P. A Perfectly Matched Layer for the Absorption of Electromagnetic Waves[J]. Journal of Computational Physics， 1994， 114（9）： 185200.

[14]Berenger J P. Perfectly Matched Layer for the FDTD Solution of WaveStructure Interaction Problem[J]. IEEE Transaction on Antennas and Propagation， 1996， 44（1）： 110117.

[15]Berenger J P. ThreeDimensional Perfectly Matched Layer for the Absorption of Electromagnetic Waves[J]. Journal of Computational Physics， 1996， 127（2）： 363379.

[16]Stanislav O， Pan G W. An Updated Review of General Dispersion Relation for Conditionally and Unconditionally Stable FDTD Algorithms[J]. IEEE Transaction on Antennas and Propagation， 2008， 56（8）： 25722583.

[17]Kusaf M， Oztoprak A Y. An Unco Nditionally Stable SplitStep FDTD Method for Low Anisotropy[J]. IEEE Microwave and Wireless Components Letters， 2008， 18（4）： 224226.

[18]Lu Yijun， Shen C Y. A Domain Decomposition FiniteDifference Method for Parallel Numerical Implementation of TimeDependent Maxwells Equations[J]. IEEE Transactions on Antennas and Propagation， 1997， 45（3）： 556562.

[19]Park J H， Strikwerda J C. The Domain Decomposition Method for Maxwells Equations in Time Domain Simulations with Dispersive Metallic Media[J]. Society for Industrial and Applied Mathematics， 2014， 32（2）： 684702.

[20]Taflove A， ME Brodwin. Computation of the Electromagnetic Fields and Induced Temperatures Within a Model of the MicrowaveIrradiated Human Eye[J]. IEEE Transactions on Microwave Theory & Techniques， 1975， 23（11）： 888896.