張 莉, 張 健,2*

(1. 四川師范大學(xué) 數(shù)學(xué)與軟件科學(xué)學(xué)院, 四川 成都 610066； 2. 電子科技大學(xué) 數(shù)學(xué)科學(xué)學(xué)院， 四川 成都 611731)

1 預(yù)備知識

Yee[1]構(gòu)造了一種特殊的網(wǎng)格剖分(后來被稱為交錯網(wǎng)格剖分),對電場E和磁場H的分量在時間和空間上進行離散,將含有時間變量的Maxwell旋度方程轉(zhuǎn)化為一組差分方程,最后按時間層逐步推進求解空間電場和磁場,即時域有限差分方法(FDTD),它是一種高效的算法,后來被應(yīng)用在更加廣泛的問題上[2-4].但是,這種方法是條件穩(wěn)定的,即要滿足Courant-Friedrichs-Lewy(CFL)[2]穩(wěn)定性條件的限制.在二維情況下,

其中,c是光速,Δt是時間步長,Δx和Δy是空間步長,因此,當(dāng)要解決的問題要求空間步長必須足夠小時,按CFL條件需要時間步長相應(yīng)取很小,這使得計算量大幅度增加,計算時間延長,有時候甚至不能實現(xiàn).

以上方法是針對下面的電導(dǎo)率和磁阻率都等于零的特殊Maxwell方程組進行計算模擬的.本文對磁導(dǎo)率不為零的Maxwell方程組構(gòu)造其分裂時域差分格式,構(gòu)造的差分格式是兩步的顯格式.兩步顯格式的差分更有利于推導(dǎo)得到解電場的線性方程組.本文給出2種格式詳細(xì)的誤差分析和數(shù)值穩(wěn)定性分析.

2 Maxwell方程的分裂格式

Maxwell方程

(1)

(2)

(3)

E=(Ex(x,y,t),Ey(x,y,t))表示電場,Hz=Hz(x,y,t)表示磁場，其中ε為介電常數(shù),單位為F/m;μ為磁導(dǎo)系數(shù),單位為H/m;σm為磁導(dǎo)率,σm即為介質(zhì)的磁損耗,單位為Ω/m.(x,y)∈Ω=[0,a]×[0,b],t∈(0,T].?Ω為Ω的邊界,n為?Ω的外法向量,并且滿足理想導(dǎo)體邊界條件

E×n=0, (0,T]×?Ω

(4)

和初始條件

E(x,y,0)=E0(x,y)=(Ex0(x,y),Ey0(x,y)),
Hz(x,y,0)=Hz0(x,y).

(5)

由文獻[12]可知,問題(1)～(3)存在唯一解.為了理論分析,假設(shè)(1)～(3)式的解具有如下正則性：

(6)

(7)

為了簡單起見,只考慮ε、μ、σm為常數(shù)的情況,文中的方法可以推廣到變系數(shù)(ε、μ、σm是關(guān)于x和y的函數(shù))的情形.

對空間區(qū)域Ω和時間區(qū)域[0,T]做如下的網(wǎng)格剖分,這種網(wǎng)格稱為交錯網(wǎng)格剖分:

為書寫簡單,對函數(shù)F(x,y,t),令

(8)

和

(9)

對(8)和(9)式應(yīng)用特殊的差分近似,可得如下差分格式,稱為兩步分裂時域有限差分格式TS-FDTDI(Two-Step Finite-Difference Time-Domain)格式:

第一步

(10)

第二步

(11)

可以看出,每一步包含2個方程,這比ADI-FDTD方法更加簡單,更利于編寫程序,但此格式關(guān)于時間是一階的.為了提高精度,在第一步中加入擾動項,得到修正后的差分格式(S-FDTDII格式):

第一步

(12)

第二步

(13)

由邊界條件可以得到(10)～(13)式的邊界條件為

(14)

其中m=n或m=n+1代表時間層.最后,由初值條件(5)可以得到差分格式的初值為

(15)

因為2個格式的處理方式類似,以S-FDTDII格式為例,推導(dǎo)差分的計算過程.首先，由(12)式的第二個方程可以得到

(16)

然后將(16)式帶入到(12)式的第一個方程,則可以得到

(17)

類似的分析可得

(18)

和

(19)

3 差分格式的誤差分析

(20)

(21)

(22)

(23)

(24)

(25)

由(6)和(7)式可得

(26)

其中

(27)

4 數(shù)值穩(wěn)定性分析

利用Fourier分析方法研究2種差分格式的數(shù)值穩(wěn)定性.2種格式分析方法類似,以修正格式TS-FDTDII格式為例進行說明.

定理4.1TS-FDTDII格式是無條件穩(wěn)定的.

證明假定TS-FDTDII格式的差分解為:

(28)

其中

因為(Ex0,Ey0,Hz0)是非零向量,因此方程組(29)的系數(shù)矩陣的行列式等于零,則

(30)

計算(30)式,化簡整理后可得

(31)

因為

kx=kcos(φ),ky=ksin(φ),
Nλ=λ/h,ω=ck.

(32)

令S=cΔt/h表示CFL數(shù)[10],其中λ是波長,Δx=Δy=h表示空間步長,則Nλ是一個波長內(nèi)的節(jié)點個數(shù),k和φ是K的圓柱坐標(biāo).因為ξ≠0,因此增長因子ξ的方程可以改寫為

(ξ-1)(c3ξ3+c2ξ2+c1ξ+c0)=0,

(33)

其中

(34)

當(dāng)σm=0,方程(33)等價于文獻[10]的方程(4.1),容易得方程(33)的一個實根為ξ=1.

下面主要考慮

c3ξ3+c2ξ2+c1ξ+c0=0.

(35)

可以看出增長因子ξ是關(guān)于S、φ、Nλ和σm的函數(shù).下面借助Matlab軟件求出方程(35)的近似解.圖1(a)是表示當(dāng)S=1.4、Nλ=40、σm=2時,方程(35)的根|ξ|隨φ的變化曲線.圖1(b)表示當(dāng)S=1.4、Nλ=40、σm=2時方程(35)的復(fù)數(shù)根|ξ|的曲線圖,可以看出|ξ|=1+O(Δt).同樣可以由圖1看出方程(35)的一個實數(shù)根|ξ|=0+O(Δt).這就表明修正格式TS-FDTDII是無條件穩(wěn)定的.

圖 1 增長因子隨著角度φ的變化情況

圖2表示當(dāng)φ=35°、S=1.5、σm=2時,方程(35)的所有根的模隨Nλ的變化曲線.

圖3的2個圖分別表示當(dāng)Nλ=60、φ=35°、σm=2時,方程(35)的所有根的模隨S的變化曲線.類似圖1和圖2的分析,可得相同的結(jié)果,因為增長因子ξ≠0,所以滿足條件的方程(33)的根的模為|ξ|=1+O(Δt),修正格式TS-FDTDII是無條件穩(wěn)定的.

5 結(jié)論

本文對磁導(dǎo)率不等于零的Maxwell方程進行研究,結(jié)合分裂算子,建立了2個不同精度的差分格式,并給出誤差的詳細(xì)分析和數(shù)值穩(wěn)定性證明.由截斷誤差表達(dá)式可以得到TS-FDTDI關(guān)于時間是1階的,但是關(guān)于空間是2階收斂;修正格式TS-FDTDII關(guān)于時間和空間都是2階收斂的.由數(shù)值穩(wěn)定性分析可以看出,修正格式是無條件穩(wěn)定的.下一步工作是針對一般的Maxwell方程(電導(dǎo)率磁導(dǎo)率都不為零)構(gòu)造一種能量守恒的差分格式,繼續(xù)探討高效的算法來探索電磁場的本質(zhì).

圖 2 增長因子隨著Nλ的變化情況

圖 3 增長因子隨著S的變化情況

[1] YEE K. Numerical solution of inition of initial boundary value problems involving Maxwell’s equations in isotropic media[J]. IEEE Trans Antennas Propagation,1996,14(3):302-307.

[2] TAFLOVE A, BRODWIN M. Numerical solution of steady-state electromagnetic scattering probelems using the time-dependent Maxwell’s equations[J]. IEEE Trans Microwave Theory Tech,1975,23(8):623-630.

[3] TAFLOVE A, HAGNESS S. Computational Electrodynamics:The Finite-Difference Time-Domain Method[M]. 2nd ed . Boston:Artech House,2000.

[4] XIE Z, CHAN C, ZHANG B. An explict four-order staggered finite-difference time-domain method for Maxwell’s equations[J]. Comput Appl Math,2002,147(1):75-98.

[5] NAMIKI T. A new FDTD algorithm based on alternating direction implicit method[J]. IEEE Trans Microwave Theory Tech,1999,47(10):2003-2007.

[6] NAMIKI T, ITO K. Investigation of numerical errors of the two-dimensional ADI-FDTD method[J]. IEEE Trans Microwave Theory Tech,2000,48(11):1950-1956.

[7] ZHENG F, CHEN Z, ZHANG J. Toward the development of a three-dimensional unconditionally stable finite-difference time-domain method[J]. IEEE Trans Microwave Theory Tech,2000,48(9):1550-1558.

[8] ZHENG F, CHEN Z. Numerical dispersion analysis of the unconditionally stable 3-D ADI-FDTD Method[J]. IEEE Trans Microwave Theory Tech,2001,49(5):1006-1009.

[9] ZHAO A. Analysis of the numerical dispersion of the 2-D alternating-direction implicit FDTD method[J]. IEEE Trans Microwave Theory Tech,2002,50(4):1156-1164.

[10] GAO L, ZHANG B, LIANG D. The splitting finite-difference time-domain methods for Maxwell’s equations in two dimensions[J]. Comput Appl Math,2007,205(1):207-230.

[11] CHEN W, LI X, LIANG D. Energy-conserved splitting FDTD methods for Maxwell’s equations[J]. Numer Math,2008,108(3):445-485.

[12] LEIS R. Initial Boundary Value Problems in Mathematical Physics[M]. New York:John-Wiley,1986.

[13] PEACEMAN D, RACHFORD H. The numerical solution of parabolic and elliptic difference equations[J]. SIAM J Soc Indust Appl Math,1955,3(1):28-41.

[14] GEDNEY S, LIU G, RODEN J. Perfectly matched layer media with CFS for an unconditional stable ADI-FDTD method[J]. IEEE Trans Antennas Propagation,2001,49(11):1554-1559.

[15] GAO L P, ZHANG B, LIANG D. Analysis of an ADI finite difference method for the time dependent Maxwell’s equations in 3-D[C]//Advances in Scientific Comput Appl. Hong Kong:Science Press,2004:171-180.

[16] CHEN G, FENG M, XIE C. A new projection-based stabilized method for steady convection-dominated convection-diffusion equations[J]. Appl Math Comput,2014,239(15):89-106.

[17] CHEN G, FENG M. Explicitly uncoupled variational mutiscale for characteristic finite element methods on the unsteady Navier-Stokes equations with high Reynolds number[J]. Appl Math Mode,2015,39(14):4202-4212.