第33卷第2期 2 0 18年5月
青岛大学学报（工程技术版）
JO U R N A L OF Q IN G D AO U N IV E R S IT Y (E&T)
Vol. 33 No. 2
May 2 0 18
文章编号！1006 -9798(2018)02 -0123 -05; D O I:10.13306).1006-9798.2018.02.025
何光峰，迟洁茹，范昊博，孙桂琪
(青岛大学电子信息学院，山东青岛266071)
摘要：针对传统的时域有限差分法受C o u ra n t稳定条件的限制，且存在交替方向隐式时域有限差
分法数值散较大的问题，本文以T E波为例，研究了C ra n k-N ico lo so n差分方式的近似去耦时域
有限差分法基本原理，并对其稳定性进行分析，证明该方法是无条件稳定。通过数值仿真，从运行 时间和吸收效果方面与传统的时域有限差分法和交替方向隐式时域有限差分法进行对比。仿真结 果表明，近似去耦时域有限差分法比交替方向时域有限差分法的吸收效果好，但比传统的时域有限
差分法吸收效果差；近似去耦时域有限差分法比交替方向时域有限差分法运行时间长，但比传统时
域有限差分法运行时间短，说明近似去耦时域有限差分法突破了C o u ra n t稳定条件的限制，且在吸
收效果方面比交替方向时域有限差分法好。该研究具有广阔的应用前景。
关键词：时域有限差分法；C ra n k-N ico lo so n差分格式方案；无条件稳定
中图分类号！〇241.8;T B5 文献标识码：A
1966年，Y e e提出传统时域有限差分（finite-difference tim e-d o m a in,F D T D)方法，该方法得到广泛应用与发
展13]。时域有限差分法是用古典显式差分的方法对麦克斯韦方程进行差分，进而对麦克斯韦方程进行微分求解。运用显式差分的F D T D方法必须满足C o u m n t稳定条件，即时间步数的选取受到空间步
长（离散网格大小）的限制。为降低数值散误差，在空间步数选取时必须满足远远小于其采样频率的波长的条件，这就使其在解决
实际问题的过程中，为得到更加真实的电磁波特性，不得不采用非常小的空间步数，进而导致计算量过大而无法
实现。为克服传统F D T D方法的缺陷，T.N a m ik i等人[4—7]提出了交替方向隐式（alternative direction im p lic it，A D I)F D T D方法。由于A D J F D T D方法对麦克斯韦方程组的所有场分量进行交替差分，从而实现了无条件稳
定，克服了传统F D T D方法必须满足C o u ra n t稳定条件的限制，但A D S F D T D方法在一个时间步数上需要进行
两次迭代，增加了每步计算的时间量和存储空间，且由于每步差分方程时间不同步，对计算误差产生了非常大的
影响。2003年，S u n G L等人811]利用C ra n k-N ico lo so n半隐式差分方式，对麦克斯韦方程组进行差分，使用近
似去稱（（pproxim ate ly decoupling,A D)法而提出的C rank-N icoloson差分方式的近似去稱（crank-
nicoloson ap-p ro xim a tely-decoupling，C N A D)F D T D方法。C N A D-F D T D方法在使用C rank-N icoloson半隐式差分方式对二
维麦克斯韦方程求解过程中，需要求解一个大型稀疏矩阵方程组，从而占用计算机大量内存和运行时间，因此采
用近似方法求解过稀疏矩阵，将£在^时刻的值代替w+1时刻的值，从而把对稀疏矩阵求解变成对三对角矩阵
求解，使计算过程简单化，并减小所需的计算机内存和运行时间。基于此，本文以横电（transverse-e le ctric，T E)波为例，介绍了C N A D-F D T D方法，并对F D T D、A D I-F D T D和C N A D-F D T D三种方法分别进行仿真实验。实
验结果表明，在相同条件下，C N A D-F D T D比A D I-F D T D的吸收效果要好，但与F D T D有一定的差距，C N A D-
F D T D的运行时间比A D--F D T D的运行时间长。
1 CNAD-FDTD 方法
1.1 CNAD FDTD理论
二维T E波在和^：三个方向上存在E z，&和分量，写成C N差分方程911]分别为
收稿日期：2017 -08- 28；修回日期：2018-01-10
基金项目：山东省自然科学基金资助项目（Z R2016F M11)
作者筒介：何光峰（1993 -)，男，山东人，硕士研究生，主要研究方向为智能信息处理％
通讯作者：迟洁茹（1970 -)，女，副教授，主要研究方向为智能信息处理，医学成像等％ E m a il:q d u c h iie n i%163.c o m
124青岛大学学报（工程技术版）第33卷
*E x |"i t l/2,j—E x |#+1/2,E."\V y ^H z |++1/2，e+ 1/2]/&D
)BD |#.i l/2 —B y |#，E+ 1/2—"l2i L K Z |++i/2，E+ l/2$/
+KT |#.1/2 !.1/2 — K t |#.1/2! + 1/2."!2^ [B i |++1/2，E+1 ]/&y— "2 21 [ By |+.1 ,+ 1/2]/其中(1") (l b) (l c)
式中，"r
2[K|+ 1/2，E+ 1/2]=K|#.1/2,e+1/2—K|#+1/2，E+ 1/2z |#. 1/2，—1/2⑵=&/2s，2=&/2";〜为时间步长，&i，&y分别为空间步长;s为介质介电常数，"为磁导系数。其迭代
过程如下:
1)将式（C)代人式（1")和式（1$)，消去#.1时刻磁场值K"+1，可得#.1时刻电场值B"+1，实现对电场值
。
2)将求得的电场值B#+1代人式（1c)，更新磁场值K#+1。其中，在对电场值进行迭代中，将B y代替B y+1，可将对大型非有限带宽稀疏矩阵的求解简化为对三对角矩阵的求解，从而将式（1a)化简为
(1.2b1 )B i |.2'—b2 (B: |#.11/2，—1+B: |.心.：）=
(1 —2b2)B i |#. 1/2, .$1 (B^ |#. 1/2，—1.B i |#+1/2，+ 1)+2"1 (K」#+1，+ 1/2—K T|#+1,E—1/2)—
2b i b y (B y |#+1,E+ 1/2—B y|#,E+ 1/2—B y |#+1,E—1/2.B y |#! —1/2)D
其中，b y=c&/2&y。
最终迭代过程为：对式（3)进行迭代；对式（l b)迭代求解B y+1；对式（3)和式（l b)更新；对式（1c)的磁场值迭。
1.2 C N A D-F D T D的稳定性分析
根据V.N e u m a n n稳定性分析，二维情况下的电磁场分量B i B y和K/12]分别为
B x=^>a£°exp['+^y A y)]，B y=tpB£n P x p^j(k x&x.k y&y)]，H T=4c£°exp[E(k x&x.k y&y"(4)
式中，办，办，4c分别为各电磁分量初始系数s是增长因子；i和7是波常数；A i，A y分别代表沿i方向和沿y 方向网格的大小。
将式（4)代人式（1")〜（0中，消去系数得关于增长因子的方程为
(―1)2 =—"(s+1)2
\f2
其中，" =^^sin2 (kx A x/2) .^^sin2(k yA y/2)。将"(">0)代人式（5)，可求得s 为s f i A x s f i A y
1 —"士2
1."
由于式（6)中的增长因子满足|s|=1，因此时间步数无论取何值，C N A D-F D T D都是稳定的。
1.3 C N A D-F D T D在完全匹酉E层的公式(5)
(6)
以有限空间模拟无限空间时，以完全匹配层[315](p e r f e c t ly m a tc h e d la y e r，P M L)作为吸收边界，则二维T E 波在P M L吸收介质中，按照C N A D-F D T D的方法进行差分[617]，得
'Q x(m)B x |#+1/2,j=D e z(m)B x |#+l/2,E.r y I H t |,+ 1/2,.1/2 ]
Q; (m)By |#j i l/2(D y(m)By |#，.1/2 +H [H|.1/2，j+ 1/2]
) (7)
H x |#.1/2,j+l/2 — Qz(m) H x |n.1/21/2―D^ (,m)2x L B y |++l,j+l/2
y |#.1/2 ,j.1/2 Q h y(m) H y |#+i1/2,j+l/2.Dy(m)T y [_B x |++1/2,j.l]
其中，m取值与左端场分量节点的空间位置相同；且。
Q; (m) ((|〇 .$y(m) AZ)，D h y(m) =2^(m)A^A^(')将磁场分量代人迭代的点场分量求值中，利用C N A D-F D T D，先求出Bx+1和B f1，再对磁场值进行迭代。2数值分析
为验证C N A D-F D T D方法的正确性，以二维平面上的一个点源的辐射场为例，分别对传统F D T D、A D I-F D T D及C N A D-F D T D方法使用C语言编写程序，在V isual Studio 2015中进行仿真，对得到的数据利用Tec-p lo t 360软件进行成像处理。
仿真条件：仿真空间为20 m m X20 m m均匀网格，磁场激励源为正弦波，位于仿真空间中心位置，频率/=1
第2期
何光峰，等：基于二维T E 波常用时域有限差分算法的分析
125
X 10
12
H z ，c = 3. 0X 108
m /s ，波长A = 3m m ，空间分辨率&r = &；y = 0. 1m m ，时间步数〜 = 0. 155p s ，吸收边界厚
度为1 m m 。P M L 中电导率的分布为
$m /x
| T T 0
$
T &
) =
----(=------
(9)
式中，（是P M L 的厚度+0
为P M L 层靠近F D T D 区的界面位置=
为整数。为了得到更好的吸收效果，$
m /x
的最
佳取值[18—20药绘图
$为
，應— 2)
(10)
其中，*()是理论反射系数值，*() = 10
—6。
通过仿真实验，得到距离中心点源为5个网格处的磁场分量，3种方法的数值比较结果如图1所示。由图1 可以看出，3种方法的结果基本相同，验证了 C N A D -F D T D 方法的正确性；通过仿真计算，得到坐标为（189,189) 处的反射误差如图2所示。由图2可以看出，在相同条件下，传统F D T D 的吸
收效果最好，C N A D -F D T D 的吸收 效果比A D S F D T D 的吸收效果要好，主要由于A D S F D T D 差分等式两边的时间不对称所导致，而C N A D -F D T D 由
一步，且时间上对称，故其
比A D S F D T D
。
图1 3种方法的数值比较结果图2坐标为（189,189)处的反射误差
由于传统
FDTD 受到
Courant
稳定条件的限制，时间步数必须充分小，即必须满足
C FL = &/
。当
CFL=1时，
FDTD 、ADI-FDTD 和
CNAD-FDTD
在相同的时间步数和空间分辨率的条件下，3
种方法的运行仿
真时间由表
1所示
。
表1
相同条件下3种方法的运行仿真时间
间步数F D T D 间/s A D I -F D T D 时间/s C N A D ~F D T D 时间/s 间步数F D T D 间/s A D I -F D T D 时间/s C N A M D T D 时间/s
1 000 2. 211 5.6418. 18650 00073.573204.966317.7645 0008. 07324.25238.474100 000
143.657
378.122
642.723
10 000
14. 806桃花岛奇遇
49 239
75.313
表2
不同C F L 条件下的运行时间
大学生微博
由表1
可以看出，当三者的运行条件相同时，
FDTD
的运行时间最短，而
CNAD-FDTD
的运行时间最长。由
于
CNAD-FDTD 和
ADIFDTD
是无
定，因此可以通过改变
CFL
取值，相应的减少运行时间步数。当其
它条件不变，
CFL= 1，
FDTD 运行步数为
50 000
时，通过改变
CFL
取值，比较ADIFDTD 和
CNAD-FDTD
的仿
真时间。不同CFL 条件下的运行时间如表
2所示。由
表
2可以看出，当
CFL 的取值增大时，
ADIFDTD
和
CNAD-FDTD
都可明显的缩短运行时间。但对于仿真 复杂结构的
，需要进行
划
，
FDTD 可能由于
计算量较大无法实现，
可利用
AD-FDTD
和
C NA
D -FD TD
来实现。
C F L 取
A D I -F D T D 时间/s C N A D ~F D T D 时间/s
2116.906178.175460.33674.20010
25.500
39.230
3结束语
$
本文主要讨论无条件稳定的
CNAD-FDTD
，
由
于其不受稳定条件的限制，所以可通过改变时间步的大小，缩
短
CPU
的运行时间。通过与传统的
FDTD 和
ADI-FDTD
两种方法的对比可以得出，
CNAD-FDTD
的吸收效
果
126青岛大学学报（工程技术版）第33卷比ADI-FDTD的吸收要好，但运行时间比ADI-FDTD的多，且吸收效果相比于FDTD还有一段差距。因此，CNAD-FDTD是一种值得深人研究和推广应用的时域算法。
参考文献：
葛德彪，闫玉波.电磁波时域有限差分方法[M].西安：西安电子科技大学出版社，2011.
[2] Yee K S. Numerical Solution of Initial Boundary Value Problems Involving MaxwelCS Equation in Isotropic Media[J]. IEEE
Transactions on Antennas and Propagation，1966，14(3) ：302 - 307.
#]Schneider J B. A Selective Survey of the Finite-Difference Time-Domain Literature[J]. IEEE Antennas and propagation，1995，37(4)：39 - 57.
#]Namiki T. A New FDTD Algorithm Based on Alternating-Direction Implicit Method [J]. IEEE Transaction on Microwave TheoryandTechniques，1999，47(10)：2003 - 2007.
[5] Namiki T. 3-D ADI-FDTD Method-Unconditionally Stable Time-Domain Algorithm for Solving Full Vector M axw ells Equa
tions [J]IEEE Transaction on Microwave Theory and Techniques，2000，48(10) ：1743 - 1748.
[6] Zheng F, Chen Z，Zhang J. Toward the Development of a Three-Dimensional Unconditionally Stable Finite-Difference Time-
Domain Method[J]. IEEE Transaction on Microwave Theory and Techniques，2000，48(9)：1550 - 1558.
[7] ZhengHongxing，LeungKW. A Nonorthogonal ADI-FDTD Algorithm for Solving Two Dimensional Scattering Problems[J].
IEEE Transaction on Antennas and Propagation，2009，57(12)：3891 - 3902.
[8] Sun G L，Trueman C W. Approximate Crank-Nicolson Schemes for the Finite-Difference Time-Domain Method for TEz Waves
[J]. IEEE Transaction on Antennas and Propagation，2004, 52(11) ：2963 - 2972.
[9] Sun G L，Trueman C W. Unconditionally Stable Crank-Nicolson Scheme for Solving the Two-Dimensional M axw ells Equa-
tions[J]. Electronics Letters，2003，39(7)：595 - 597.
[10] X ie X，Pan G，Hall S. A Crank-Nicholson-Base Unconditionally Stable Time-Domain Algorithm for 2D and 3D problem[J].
Microwave and Optical Technology Letters，2007，49(2)：261 - 265.
[11] Li J X，Jiang H L，Zhao X M，et al. Effective CNAD-and ADE-Based CFS-PML Formulations for Truncating the Dispersive
FDTD Domains[J]. IEEE Antennas and Wireless PropagationLetters，2015，14：1267 - 1270.
[12]王秉中.计算电磁学[M].北京：科学出版社，2002.
[13] Berenger J P. A Perfectly Matched Layer for the Absorption of Electromagnetic Waves[J]. Journal of Computational Phys
2013小企业会计准则ics，1994, 114(9) ：185 - 200.
[14] Berenger J P. Perfectly Matched Layer for the FDTD Solution of Wave-Structure Interaction Problem[J]. IEEE Transaction
高尔基资料
on Antennas andPropagation，1996，44(1) ：110-117.
[15] Berenger J P. Three-Dimensional Perfectly Matched Layer for the Absorption of Electromagnetic Waves[J]. Journal of Com
putational Physics，1996，127(2)：363 - 379.
[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) ：2572 - 2583.
[17] K u safM，Oztoprak A Y. An Unco Nditionally Stable Split-Step FDTD Method for Low Anisotropy [J]. IEEE Microwave
and Wireless Components Letters，2008，18(4)：224 - 226.
[18] Lu Yijun，ShenC Y. A Domain Decomposition Finite-Difference Method for Parallel Numerical Implementation of Time-De
pendent M axw ells Equations[J]. IEEE Transactions on Antennas and Propagation，1997，45(3)：556 - 562.
[19] Park J H，Strikwerda J C. The Domain Decomposition Method for M axw ells Equations in Time Domain Simulations with
Dispersive Metallic Media[J]. Society for Industrial and Applied Mathematics，2014, 32(2) ：684 - 702.
[20] Taflove A，ME Brodwin. Computation of the Electromagnetic Fields and Induced Temperatures Within a Model of the Micro
wave-Irradiated Human Eye[J]. IEEE Transactions on Microwave Theory C Techniques，1975，23(11) ：888 - 896.
着重号用法补充规则
j 不应使用文字下加直线或波浪线等形式表示着重。文字下加直线为专名号形式（4. 16)文字下加
t 波浪线是特殊署名号(A. 13. 6)。着重号的形式统一为相应项目下加小圆点。
I 示例：下面对本文的理解，不正确的一项是(误）
下面对本文的理解，不正确的一项是（正)
第2期何光峰，等：基于二维T E波常用时域有限差分算法的分析127 Analysis of Time-Domain Finite Difference Algorithm Based
on Two-Dimensional TE Wave
HE Guangfeng，CHI Jieru，FAN Haobo，SUN Guiqi
(School of Electronic Inform ation，Q ingdao University，Q ingdao 266071，C hina)
A b stra c t:Because the tra d itio n a l tim e-dom ain fin ite difference m ethod is restricted by Courant s the num erical dispersion of the fin ite difference m ethod is very im p o rta n t.In this paper，the s proxim ate decoupling tim e domain fin ite difference m ethod (C N A D-F D T D)is analyzed，and it i the m ethod is u nconditional and stable.T h e n，by num erical s im u la tio n，fro m tw o aspects of running tim e and absorption effect are compared w ith those of the conventional fin ite difference tim e domain m ethod (F D T D)and
based on alternating im p lic it fin ite difference tim e domain m ethod(A D--F D T D)，it is k n o w n t h a t t h e C N A D-
F D T D has a better absorption effect than A D--F D T D，b ut has a gap co m p a re d w ith the tra d itio n a l F D T D，the
C N A D-F
D T D runs longer than the A D--F D T D，b ut shorter than the tra d itio n a l F D T D.I t ca n b e concluded
th a t the approxim ate decoupling tim e domain fin ite d ifferencem ethod (C N A D-F D T D)breaks the restrictio n of
Courant s ta b le c o n d itio n a n d is b e tte r than th e A D I-F D T D in th e a b s o rp tio n e ffe c t.
K e y w o rd s：tim e domain fin ite difference m e th od；C rank-N icoloson d ifferen tia l scheme；unconditional s ta b ility
(上接第107页）
Study of the Road Roughness Based on Measured Road Elevation Data
XU Zhongyuan，LI Xuefeng，SUN Rongqiang，ZHU Xiaoqing，LIU Dawei (School of Electrom echanic E ngin eering，Qingdao U n iv e rs ity，Qingdao 266071，China)
A b stra c t：In order to learn about the re g u la rity of the self-paved experim ental pavem ent，the leveling in s tru
m ent is adopted to measure the differences o f three-dim ensional elevation in differen t locations of the pavem ent，
and then obtain the three-dim ensional elevation curved surface of w hole pavement by statistical a sides，the em pirical mode decom position (E M D)m ethod is adopted to deal w ith the trend ter re g u la ritie s，a n d o b ta in th e d a ta s e rie so fp a ve m e n t irre g u la ritie s w h ic h a c c o rd w ith the relative g ra d e o fp a v e-
m e n t，and then reconstruct th e m a p o fp a v e m e n t irre g u la ritie s，a n d solvea n da na lyze th e p o w e r s p e c tru m o f the
内蒙古师范大学学报experim entalpavem ent by establishing the A R (autoregressive，A u t o Regressive)m odel.T he re
s u lts o f the
p o w e rs p e c tru m d e n s ity a n a ly s is indicate th a t th e w h o le g r a d e o f th ep ave m en t belongs to class G and canbe
used t o s t udy the vehicle dynamic response.The analysis results show th at in the low f experim ental road displacement p ower spectrum density is larger and in the high frequency section，the experi
m ental road displacement power spectrum density is sm all.It can be seen th at the characteristic the experim ental pavement are sim ilar to the G level pavem ent，and the w hole road surface belongs t el pavement.The research p rovides more accurate info rm a tion about the road roughness fo r the dynamic re
sponse of the vehicle.
K e y w ord s：road roughness；E M D;A R m od e l；power spectral density
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议