何光峰　遲潔茹

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

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

中圖分類(lèi)號(hào)： O241.8；TB5文獻(xiàn)標(biāo)識(shí)碼： A

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

1CNADFDTD方法

1.1CNADFDTD理論

二維TE波在x，y和z三個(gè)方向上存在Ex，Ey和Hz分量，寫(xiě)成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為時(shí)間步長(zhǎng)，Δx，Δy分別為空間步長(zhǎng)；ε為介質(zhì)介電常數(shù)，μ為磁導(dǎo)系數(shù)。其迭代過(guò)程如下：

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

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

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。

最終迭代過(guò)程為：對(duì)式（3）進(jìn)行迭代；對(duì)式（1b）迭代求解En+1y；對(duì)式（3）和式（1b）更新；對(duì)式（1c）的磁場(chǎng)值迭代更新。

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

根據(jù)V.Neumann穩(wěn)定性分析，二維情況下的電磁場(chǎng)分量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ù)；ε是增長(zhǎng)因子；kx和ky是波常數(shù)；Δx，Δy分別代表沿x方向和沿y方向網(wǎng)格的大小。

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

ε-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）中的增長(zhǎng)因子滿(mǎn)足ε=1，因此時(shí)間步數(shù)無(wú)論取何值，CNADFDTD都是穩(wěn)定的。

1.3CNADFDTD在完全匹配層的公式

以有限空間模擬無(wú)限空間時(shí)，以完全匹配層[1315]（perfectly matched layer，PML）作為吸收邊界，則二維TE波在PML吸收介質(zhì)中，按照CNADFDTD的方法進(jìn)行差分[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取值與左端場(chǎng)分量節(jié)點(diǎn)的空間位置相同；且。

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

將磁場(chǎng)分量代入迭代的點(diǎn)場(chǎng)分量求值中，利用CNADFDTD，先求出En+1x和En+1y，再對(duì)磁場(chǎng)值進(jìn)行迭代。

2數(shù)值分析

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

仿真條件：仿真空間為20 mm×20 mm均勻網(wǎng)格，磁場(chǎng)激勵(lì)源為正弦波，位于仿真空間中心位置，頻率f=1×1012 Hz，c=30×108 m/s，波長(zhǎng)λ=3 mm，空間分辨率Δx=Δy=01 mm，時(shí)間步數(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。

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

3結(jié)束語(yǔ)

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

參考文獻(xiàn)：

[1]葛德彪， 閆玉波. 電磁波時(shí)域有限差分方法[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]王秉中. 計(jì)算電磁學(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.