第十一章时域有限差分方法.docx
- 文档编号:29734352
- 上传时间:2023-07-26
- 格式:DOCX
- 页数:59
- 大小:78.85KB
第十一章时域有限差分方法.docx
《第十一章时域有限差分方法.docx》由会员分享,可在线阅读,更多相关《第十一章时域有限差分方法.docx(59页珍藏版)》请在冰豆网上搜索。
第十一章时域有限差分方法
第十一章-时域有限差分方法
第十一章时域有限差分方法
自从1966年K.S.Yee创建时域有限差分法(FiniteDifferenceTimeDomain,简称FDTD)
[1]以来,已经发展成为一种理论完整、应用广泛的数值方法,并且与矩量法和有限元法一起奠定了计算电磁学的基础。
本章将介绍时域有限差分的基本理论,数值模拟技术,若干相关的专题以及工程实例。
11-1差分的基本概念
时域有限差分法是对微分形式的Maxwell方程进行差分求解的技术。
在详述其之前,首先简单回顾差分的基本概念。
已知分段连续函数在位置处的增量可表示为fxx,,
(11-1-1),,,,,fxfxxfx,,,,,,
其差商为
,,,fxfxxfx,,,,,,(11-1-2),,,xx
x当,0时,fx的导数定义为差商的极限,即,,
,,,fxfxxfx,,,,,,'limlim(11-1-3)fx,,,,,,,,xx00,,xx
x当足够小时,的导数可以近似为fx,,
dff,,(11-1-4)dxx,
根据导数取值位置的不同,差分格式分为前向差分、后向差分和中心差分。
前向差分定义为
fxxfx,,,,,,,,f(11-1-5),,,xxx
后向差分定义为
fxfxx,,,,,,,,f(11-1-6),,,xxx
中心差分定义为
fxxfxx,,,,,22,,,,,f(11-1-7),,,xxx
fxx,,将在点x处展开为Taylor级数,得,,
23dddfxfxfx,,,,,,1123(11-1-8)fxxfxxxx,,,,,,,,,,,,,23d2!
d3!
dxxx
371
23dddfxfxfx,,,,,,1123(11-1-9)fxxfxxxx,,,,,,,,,,,,,23d2!
d3!
dxxx
将方程(11-1-8)和(11-1-9)代入(11-1-5)~(11-1-7)后可以发现,前向和后向差分具有一阶精度,中心差分具有二阶精度。
11-2FDTD概述
时域有限差分法直接将两个旋度Maxwell方程,在空间和时间上用中心差分格式进行离散,从而获得一组递推方程。
中心差分格式能够保证时域有限差分的解具有二阶精度,并且
[2]在满足Courant条件时其结果是稳定的。
在Yee的差分格式中,计算区域在x,y,z方向上分别用直角坐标网格进行离散。
划分电场的网格称为电网格,而划分磁场的网格称为磁网格,所谓时域有限差分网格通常是指电网格。
按Yee的定义,在电网格单元中,电场采样与电网格单元的棱边重合,磁场采样则位于电网格面的中心且与电网格面垂直,如图11-2-1(a)所示。
同样,在磁网格单元中,磁场采样与磁网格单元的棱边重合,电场采样则位于磁网格面的中心且与磁网格面垂直。
电网格单元与磁网格单元的相互位置关系如图11-2-1(b)所示。
磁场单元磁场单元HHEExxzzHHzzHHyy
EEyyHHzz
HHyyEEzzHHyyHHxxEEzzEExxEEHHxxxxzzEEyyEEyyyy电场单元电场单元xx(a)电网格单元(b)电网格和磁网格的位置关系
图11-2-1Yee差分格式中电场和磁场的位置
时域有限差分中的基函数是以采样点为中心的脉冲函数,也就是说我们假设电场在电网格单元的棱边上均匀分布,同时也在磁单元网格面上均匀分布。
磁场基函数的定义与此类似。
nt,,tnt,,12在时域中,电场采样时刻为,为时间采样间隔,且假定电场在时间间隔到,,
nt,nt,,12nt,,1均匀分布。
磁场则在时间间隔到均匀分布。
由第一章获悉,当介质,,,,
中存在电损耗及磁损耗时,两个Maxwell旋度方程可以表示为
Er(,)t(11-2-1a),,,,,,HrEr(,)(,)tt,t
Hr(,)tm,,,,,ErHr(,)(,)tt,,(11-2-1b),t
m,式中,为介电常数,,为电导率,,为磁导率,代表磁损耗,可以称为磁耗率。
在直角坐标系中,将上式写成分量形式,获得的6个相互耦合的偏微分方程为
E,,,H,E1ymxz,,,,H(11-2-2a),,xx,,,tzy,x,,
372
H,E1,E,,ymxz(11-2-2b),,,,Hyy,,,,,,txz,,y
E,,,E,H1ymxz(11-2-2c),,,,H,,zz,,,tyx,,,z
H,,,E,H1yxz,,,,E(11-2-2d),,xx,,,tyz,,,x
E,H1,H,,yxz(11-2-2e),,,,Eyy,,,,,,tzx,,y
H,,,H,E1yxz(11-2-2f),,,,E,,zz,,,txy,,,z
这里将介质参数也写成分量形式,以便可以直接模拟一些简单的各向异性介质,同时也有利于处理计算区域内部的介质界面问题。
方程(11-2-2a)~(11-2-2f)建立了时域有限差分数值方法模拟时变电磁场与三维物体相互作用的基础,它与边界条件结合,即可解决几乎所有的电磁问题。
若以电网格为参考系,离散化的电磁场可以表示为
n(11-2-3a)EijkEixjykznt,,,,,,,12,,12,,,,,,,,,xx
n(11-2-3b)EijkEixjykznt,12,,12,,,,,,,,,,,,,,,yy
nEijkEixjykznt,,12,,12,,,,,,,,(11-2-3c),,,,,,zz
n,12HijkHixjykznt,12,12,12,12,12,,,,,,,,,,(11-2-3d),,,,,,,,,,xx
n,12HijkHixjykznt,,,,,,,,,,12,,1212,,12,12(11-2-3e),,,,,,,,,,yy
n,12HijkHixjykznt,,,,,,,,,,12,12,12,12,,12(11-2-3f),,,,,,,,,,zz
现实世界中的电场和磁场充满整个空间,同时在时间上也是连续的。
但是,在时域有限
nt,,t差分的Yee格式中,电场总是在整时间步上采样,而磁场则在半时间步(n+1/2)上采样。
因此,空间某点在某时刻的电磁场分量需要通过时间和空间上进行插值才能获得。
同样为了计算频域的电磁场,在进行Fourier变换时也应考虑到电场和磁场在时间上具有半个步长的时移。
如果忽略电磁场采样在时间或空间上的错位,将会导致模拟结果在高频段产生误差。
利用时间和空间中心差分公式以及式(11-2-3a)~(11-2-3f),可将Maxwell方程(11-2-2a)~
(11-2-2f)改写成下面递推形式:
m,,,,0.5tnn,,1212xxHijkHijk,12,12,12,12,,,,,,,,,xxm,,0.5t,,xx
nn,EijkEijk,12,1,12,,,,,,,,,,tyy,,m,,,,,0.5tz,xx,
373
nn,EijkEijk,1,12,,12,,,,,,,,zz(11-2-4a),,,y,
m,,,,0.5tyynn,,1212HijkHijk,,,,,12,,1212,,12,,,,yym,,0.5t,,yy
nn,EijkEijk,,,,1,,12,,12,,,,,tzz,,m,,,,,0.5txyy,
nn,EijkEijk,,,,12,,112,,,,,,xx(11-2-4b),,,z,
m,,,,0.5tnn,,1212zzHijkHijk,,,,,12,12,12,12,,,,,zzm,,0.5t,,zz
nn,EijkEijk,,,,12,1,12,,,,,,,txx,,m,,,,,0.5tyzz,
nn,EijkEijk,,,,1,12,,12,,,,,yy(11-2-4c),,,x,,
,,,0.5tnn,1xxEijkEijk,,,12,,12,,,,,,xx,,0.5t,,xx
nn,,1212,HijkHijk,,,,,12,12,12,12,,,,,,tzz,,,,,,,0.5tyxx,
nn,,1212,HijkHijk,,,,,12,,1212,,12,,,,yy(11-2-4d),,,z,,
,,,0.5tyynn,1EijkEijk,12,,12,,,,,,,,yy,,0.5t,,yy
nn,,1212,HijkHijk,,,,,12,12,12,12,,,,,,tzz,,,,0.5tx,,,yy,
nn,,1212,HijkHijk,12,12,12,12,,,,,,,,,xx(11-2-4e),,z,,
,,,0.5tnn,1zzEijkEijk,,12,,12,,,,,,,zz,,0.5t,,zz
nn,,1212,HijkHijk,,,,,12,,1212,,12,,,,,tyy,,,,,,,0.5tx,zz,
374
nn,,1212,HijkHijk,12,12,12,12,,,,,,,,,xx(11-2-4f),,,y,式中介质参数的空间坐标序号和对应的场分量序号相同。
为了求解方程组(11-2-4),需用适当的边界条件截断时域有限差分的计算区域。
时域有限差分法涉及的边界条件通常包括理想电导体(PerfectElectricConductor,简称PEC),理想磁导体(PerfectMagneticConductor,简称PMC),吸收边界件(AbsorbingBoundaryCondition,
[4-12][13-15]简称ABC)和周期边界件(PeriodicBoundaryCondition,简称PBC)。
除了截断计算区域外,还必须对计算区域内部的非均匀介质边界进行处理。
如果非均匀介质边界与网格重合,根据电场和磁场的位置定义,磁场采样位于两个网格中心的连线上,与其相应的有效磁介质参数应该是两个对应网格磁参数的加权平均,而电场位于电网格棱边上。
根据法拉第电磁感应定律,磁场的积分环路最多可以跨越电场周围的4个网格,因此与电场相联系的有效电介质参数应该是相关网格参数的加权平均。
弯曲的金属或介质边界需要用共
[16-24]形技术处理。
由上递推方程可见,电场和磁场在时间和空间上都是交错的。
因此,我们只需要存储前一时刻的电磁场值,而当前时刻的电(磁)场值将由前一时刻的电(磁)场值和前半时刻磁(电)场值求得。
此外,三维电场和磁场数组的大小与空间格点数量相同,但是在时间上只需要保留最后一个时刻的电磁场值即可。
这些特点对于实际编程是很有价值的。
定义电磁场的位置序号的方法很多,通常采用的方法如11-2-2图所示。
y(n-1,n)EE(0,n)xx
,,(0,n-1)Ey,,E(n,n-1)y(n-1,n-1)HH(0,n-1)zz
E(1,1)yE(0,1)yE(0,1)E(1,1)xx(n-1,0)HzH(0,0)z,,(0,0)EE(1,0)E(n,0)y,yy,x0E(n-1,0)E(0,0)xxE(1,0)x
图11-2-2在Yee差分各式中电磁场的相对位置
11-3网格数值色散
已知在非色散媒质中,不同频率的电磁波具有相同的传播速度,此时电磁波角频率和波数的关系为
2,,,222kkk,,,(11-3-1),,xyzc,,
k式中k、和k分别是沿x、y和z方向的波数。
然而,对于时域有限差分来说,即使在非色yxz
散媒质中,电磁波在不同方向上仍然可能具有不同的传播速度,这种特性称为网格数值色散。
网格数值色散与网格大小、形状以及差分格式有关。
下面讨论在直角坐标系中时域有限差分的网格色散误差。
375
设一个平面波的波函数可以表示为
j(),tkxkykz,,,xyz(11-3-2)xyzt,,,e,,,,,0
x,t令、、和分别为空间和时间的采样间隔,那么,上式的离散形式可以表示为,y,z
j(),ntkIxkJykKz,,,,,,,nxyz(11-3-3)IJK,,e,,,,,0
式中,I,J,K是时间和空间采样点的序号。
n
已知自由空间中电磁场各分量满足下述波动方程
2222,,,,,,1(11-3-4),,,,,0,,22222,,,,xyzct,,
式中c为真空光速。
利用时间和空间的中心差分公式离散上式,求得
nnn,,,ijkijkijk,,,,1,,2,,1,,,,,,,,
2,x
nnnijkijkijk,1,2,,,1,,,,,,,,,,,,,,,2,y(11-3-5)nnnijkijkijk,,12,,,,1,,,,,,,,,,,,,,2,z
nnn11,,ijkijkijk,,2,,,,,,,,,,,,,,,,2,t
将(11-3-3)代入上式,得
2222,,ky,,,,,,,,kx,,1111,kz,,t,,,,,,yxz(11-3-6)sinsinsinsin,,,,,,,,,,,,,,,,,,,2222ctxyz,,,,,,,,,,,,,,,,,,,,
,x0,,z0当,,时,上式即变为式(11-3-2),可见网格数值色散随网格的,,y0
尺寸减小而减弱。
此外,网格数值色散误差还与传播方向有关。
为了说明简单起见,令
,,xy。
那么,一根z方向的无限长线源,产生的电场分量在x-y平面的分布特性如图E,,z
11-3-1所示。
图中的白色圆环代表柱面波的波面。
由此可见,电磁波在=45:
,135:
,225:
和315:
的对角线方向上,网格数值色散误差最小。
而且,网格尺寸越小,网格数值色散误差
也越小。
(a)(b)(c),,,,xy,5,,,,xy,10,,,,xy,20
图11-3-1时域有限差分的网格色散误差随网格尺寸的变化
376
11-4稳定性分析
对于以时间步进为基础的数值方法,解的稳定性是衡量时域方法性能的主要判据之一。
时域方法的稳定性与其物理模型、差分格式以及网格结构有关。
直接求解式(11-3-6),得
,ky,,,kx,kz2111,,,,,y222xz,,(11-4-1),,,,,arcsinsinsinsinct,,,,,,222,,,,,,txyz222,,,,,,,,
若上式中为虚数,由式(11-3-2)可见,电磁场的解将随时间快速衰减或发散。
为了保证为,,
实数,式(11-4-1)中括号内的部分应该满足下述条件
ky,,,kx,kz,111,,,,y222xz(11-4-2)ct,,,,sinsinsin1,,,,,,222222,,,xyz,,,,,,
因为上式根号下的正弦平方项最大值为1,所以时间步长,t的取值应该满足下式
1,,t(11-4-3)111c,,222,,,xyz
该式为获得时域有限差分稳定解的必要条件,即Courant稳定条件,又称为Courant-Friedrichs
[25]-Lewy稳定判据。
式(11-4-3)表明,时域有限差分的时间步长取决于3个方向上的网格尺寸和传播速度。
为了进一步理解时域有限差分稳定条件(11-4-3)的意义,将此式在一维情况下简化为ctx,,,。
那么,电磁波在自由空间中沿x方向从第n个网格传播到第n+1个网格所需的时间为。
所以在时域有限差分模拟中,若时间步长,那么当电磁波传播到,,,txc,,,txc
第n+1个网格之前即获得场值。
显然,这将违背因果关系,导致不稳定解。
11-5截断时域有限差分网格的边界条件
在时域有限差分法的求解中,吸收边界条件极为重要。
对于一个开域问题,由于计算机内存和速度的限制,时域有限差分法的计算区域不可能延伸到无穷远。
时域有限差分方法早在1966年就已创建,但是真正得到快速发展还是在1981年荷兰科学家G.Mur建立吸收边界
[5]条件以后。
利用Mur的吸收边界条件,时域有限差分法才有可能解决实际工程问题。
此后,梅冠香(K.K.Mei)和方家园(J.Y.Fan)等学者在Mur的吸收边界条件基础上提出了所谓的超[3,6]吸收边界条件。
美国Illinois大学周永祖(W.C.Chew)教授发展了中国科学家廖振鹏等提出
[7]的吸收边界条件,使其成为著名的廖氏吸收边界条件。
与Mur吸收边界条件相比,廖氏吸收边界条件具有更好的吸收特性。
但是,其主要缺点是高阶吸收边界条件不稳定。
正当人们试图改进高阶廖氏吸收边界条件的稳定特性的时候,1994年法国学者J.P.Berenger提出完全
[811]匹配层(PerfectlyMatchedLayer,简称PML)边界条件,人们的注意力开始转向PML边~
界条件。
PML边界条件的创建引起了吸收边界条件领域y内一场革命,这种边界条件在理论上可以吸收来自不同
方向、不同频率的电磁波。
本节将首先简单介绍理想导,,电体(PEC)、理想导磁体(PMC)和Mur吸收边界条件,
然后重点讨论PML边界条件。
虽然PML被称为边界条()H不可利用zHz,,,件,实际上应该被看作一种各向异性的吸波材料。
Ey
在介绍各种边界条件之前,首先说明为什么时域有xEx限差分法仿真需要边界条件。
差分方程式(11-2-4)仅给出()H不可利用z,,
图11-5-1计算区域边界的电磁场分布377
在空间某点某一时刻电场和磁场之间的相互关系,并不包含与计算区域相关的边界信息,因此需要与边界条件结合才能完整解决电磁问题。
边界条件包括截断计算区域的边界条件和不同介质界面上的边界条件,本节主要介绍截断计算区域的边界条件。
按照式(11-2-4)给出的时域有限差分递推关系,图11-5-1中计算区域内部的电场和磁场分别可用周围的磁场和电场,以及同一位置前一时刻的场值求得。
因为计算区域外部的磁场是未知的,所以计算区域边界上的切向电场值并不能由式(11-2-4)的递推关系导出。
边界条件的任务是寻求一种无须任何外部信息即能求解计算区域边界上切向电场的方法。
11-5-1PEC和PMC边界条件
理想导电体(PEC)是一种理想边界条件,它能够完全反射电磁波。
应用PEC边界条件截断时域有限差分网格时,只需简单地强迫理想导电体表面上总电场的切向分量为零即可。
PEC边界条件一般用于波导、谐振腔、微带电路和微带天线底板的近似。
为了节省内存和缩短计算时间,PEC边界条件可以在对称面上截断具有对称性问题的计算区域,此外,还可用于截断平面波正入射时的对称周期结构。
理想导磁体也是一种理想边界条件,也能够完全反射电磁波。
PMC边界条件一般用于截断对称结构问题以减小计算区域,也可用作平面波正入射时截断对称周期结构的边界条件。
PMC与PEC边界条件结合可以模拟平面波正入射时
二维对称周期结构问题。
应用PMC边界条件截断时EE,HijkH(i+1/2,j,k)yz2z2,,,,,,域有限差分网格时,可以强迫理想导磁体上的总磁,H(i-1/2,j,k)z1
场切向分量为零。
,,,,,,由于时域有限差分中的切向磁场位于半个网格
上,PMC边界条件无法与物体表面重合,实际应用计算区域计算区域
很不方便。
一种比较好的实现方式是使PMC边界条PMCPMC边界边界件不与切向磁场重合,而像PEC边界条件一样与计图11-5-2PMC边界条件的应用
算区域边界上的切向电场重合,然后用镜像原理实
现PMC边界条件,如图11-5-2所示。
这样设置的PMC边界条件,其切向电场仍然处在PMC边界上。
PMC边界上切向电场的计算需要边界两边的切向磁场和计算区域边界上法向磁场的信息。
根据镜像原理,PMC边界两边的切向磁场分量和数值相等但是方向相反,即HHz1z2
nn,,1212HijkHijk,,,,,,12,12,12,12,,,,,zz
EPMC边界上的切向电场HH由磁场和利用时域有限差分递推方程求得,即yz1z2
,,,0.5t,tyynn,1EijkEijk,12,,12,,,,,,,,,yy,,,,0.50.5tt,,,,yyyy
,12nn,12,,212,12,Hijk,,,,,Hzx(11-5-1),,,,,,0.510.51,,,,,,,,zkzkxixi,,,,,,,,,,,,,,
E由式(11-5-1)可见,计算切向电场不仅需要计算区域边界外的切向磁场,而且还需要计y
算区域边界上的法向磁场。
也就是说,计算区域边界上的法向磁场也需要在每一时间步递推。
11-5-2Mur吸收边界条件
在PML边界条件出现之前,Mur吸收边界条件对于时域有限差分的发展曾发挥过重要作用。
至今,对于要求不高的求解,如波导、微带天线及微波电路等问题,使用Mur吸收边界
378
条件仍然能够获得满意的结果。
与PML边界条件相比,Mur吸收边界条件节省内存,且计算速度较快。
考虑到仅在一些特殊情况下使用Mur吸收边界条件,本节只介绍一阶Mur吸收边界条件。
二阶Mur吸收边界条件的推导过程留给读者自己完成。
[4,5]Mur边界条件是对Engquist-Majda吸收边界条件的一种差分近似,因此在介绍Mur边界条件之前,需要回顾一下Engquist-Majda吸收边界条件。
通常,沿方向传播的平面波,x可定义为,可见当时间t增加时x应该减小,以保持相位恒定的波面。
该平面波满足,xct,,,
的波动方程为
,1,,(11-5-2),,,xt,0,,,,,,xct,,
满足上式的平面波解只能沿方向传播。
图11-5-3,x,,xctxct,,,,,,表示一个沿着-x方向传播的平面波,除了计算区域边界点边界点
左边界上的场值未知外,内部的场值均可通过时域,,,,,,,,xx有限差分递推关系获得。
计算区域内部每一点的场xx=0=0xx=1=1值都满足方程(11-5-2),下面讨论边界上场值的计算
图11-5-3吸收边界条件的图形解释方法。
x,0由于在计算区域外没有任何已知的场信息,在边界点处的场不能够通过时域有限差分递推关系得到。
Engquist-Majda吸收边界条件使用计算区域内部的已知信息表示边界点x,0上的场值。
式(11-5-2)对于垂直入射平面波严格成立,如果对式(11-5-2)在时间和空间上
x,0作前向差分近似,那么在边界点处,时刻的场值可表示为nt,,1,,
ctct,,,,nnn,1(11-5-3),,,1,,,xxx,,,001,,,,xx,,
由上式可见,计算区域边界上的场值可用前一时刻的内部场信息表示。
因为式(11-5-2)为一维波方程,式(11-5-3)仅对垂直入射的情况成立。
Mur对式(11-5-2)进行时间和空间上的中心差分近似,得
11,,11nn,,nn,1,,22,,,,,,,(11-5-4),,xx,,1011,,xx,,,,xct,,22,,
上式具有二阶精度。
然而,式(11-5-4)没有包含计算区域边界上在nt,,1时刻的场值。
Mur,,
nn,12进一步对该式中和取时间和空间的平均,即,,x,1/2x,0
1nnn,,121,,,,,(11-5-5a),,xxx,,,0002
1nnn,,,,,(11-5-5b),,110,,xx,x22
将上述两式代入方程(11-5-4),化简后求得
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第十一 时域 有限 方法