第十一章-时域有限差分方法

371第十一章时域有限差分方法自从1966年K.S.Yee创建时域有限差分法(FiniteDifferenceTimeDomain，简称FDTD)以来[1]，已经发展成为一种理论完整、应用广泛的数值方法，并且与矩量法和有限元法一起奠定了计算电磁学的基础。本章将介绍时域有限差分的基本理论，数值模拟技术，若干相关的专题以及工程实例。11-1差分的基本概念时域有限差分法是对微分形式的Maxwell方程进行差分求解的技术。在详述其之前，首先简单回顾差分的基本概念。已知分段连续函数fx在位置x处的增量可表示为fxfxxfx(11-1-1)其差商为fxfxxfxxx(11-1-2)当x0时，fx的导数定义为差商的极限，即00'limlimxxfxfxxfxfxxx(11-1-3)当x足够小时，fx的导数可以近似为ddffxx(11-1-4)根据导数取值位置的不同，差分格式分为前向差分、后向差分和中心差分。前向差分定义为xfxxfxfxx(11-1-5)后向差分定义为xfxfxxfxx(11-1-6)中心差分定义为22xfxxfxxfxx(11-1-7)将fxx在点x处展开为Taylor级数，得232323ddd11d2!d3!dfxfxfxfxxfxxxxxxx(11-1-8)372232323ddd11d2!d3!dfxfxfxfxxfxxxxxxx(11-1-9)将方程(11-1-8)和(11-1-9)代入(11-1-5)~(11-1-7)后可以发现，前向和后向差分具有一阶精度，中心差分具有二阶精度。11-2FDTD概述时域有限差分法直接将两个旋度Maxwell方程，在空间和时间上用中心差分格式进行离散，从而获得一组递推方程。中心差分格式能够保证时域有限差分的解具有二阶精度，并且在满足Courant条件时其结果是稳定的[2]。在Yee的差分格式中，计算区域在x，y，z方向上分别用直角坐标网格进行离散。划分电场的网格称为电网格，而划分磁场的网格称为磁网格，所谓时域有限差分网格通常是指电网格。按Yee的定义，在电网格单元中，电场采样与电网格单元的棱边重合，磁场采样则位于电网格面的中心且与电网格面垂直，如图11-2-1(a)所示。同样，在磁网格单元中，磁场采样与磁网格单元的棱边重合，电场采样则位于磁网格面的中心且与磁网格面垂直。电网格单元与磁网格单元的相互位置关系如图11-2-1(b)所示。时域有限差分中的基函数是以采样点为中心的脉冲函数，也就是说我们假设电场在电网格单元的棱边上均匀分布，同时也在磁单元网格面上均匀分布。磁场基函数的定义与此类似。在时域中，电场采样时刻为nt，t为时间采样间隔，且假定电场在时间间隔12nt到12nt均匀分布。磁场则在时间间隔nt到1nt均匀分布。由第一章获悉，当介质中存在电损耗及磁损耗时，两个Maxwell旋度方程可以表示为(,)(,)(,)ttttErHrEr(11-2-1a)m(,)(,)(,)ttttHrErHr(11-2-1b)式中为介电常数，为电导率，为磁导率，m代表磁损耗，可以称为磁耗率。在直角坐标系中，将上式写成分量形式，获得的6个相互耦合的偏微分方程为m1yxzxxxEHEHtzy(11-2-2a)图11-2-1Yee差分格式中电场和磁场的位置(a)电网格单元(b)电网格和磁网格的位置关系zxyxExHyEzEyHzH磁场单元电场单元xHxEzEyExHyHyHzEyEzHzxyxExHyEzEyHzH磁场单元电场单元xHxEzEyExHyHyHzEyEzH373m1yxzyyyHEEHtxz(11-2-2b)m1yxzzzzEEHHtyx(11-2-2c)1yxzxxxHEHEtyz(11-2-2d)1yxzyyyEHHEtzx(11-2-2e)1yxzzzzHHEEtxy(11-2-2f)这里将介质参数也写成分量形式，以便可以直接模拟一些简单的各向异性介质，同时也有利于处理计算区域内部的介质界面问题。方程(11-2-2a)~(11-2-2f)建立了时域有限差分数值方法模拟时变电磁场与三维物体相互作用的基础，它与边界条件结合，即可解决几乎所有的电磁问题。若以电网格为参考系，离散化的电磁场可以表示为12,,12,,,nxxEijkEixjykznt(11-2-3a),12,,12,,nyyEijkEixjykznt(11-2-3b),,12,,12,nzzEijkEixjykznt(11-2-3c)12,12,12,12,12,12nxxHijkHixjykznt(11-2-3d)1212,,1212,,12,12nyyHijkHixjykznt(11-2-3e)1212,12,12,12,,12nzzHijkHixjykznt(11-2-3f)现实世界中的电场和磁场充满整个空间，同时在时间上也是连续的。但是，在时域有限差分的Yee格式中，电场总是在整时间步nt上采样，而磁场则在半时间步(n+1/2)t上采样。因此，空间某点在某时刻的电磁场分量需要通过时间和空间上进行插值才能获得。同样为了计算频域的电磁场，在进行Fourier变换时也应考虑到电场和磁场在时间上具有半个步长的时移。如果忽略电磁场采样在时间或空间上的错位，将会导致模拟结果在高频段产生误差。利用时间和空间中心差分公式以及式(11-2-3a)~(11-2-3f)，可将Maxwell方程(11-2-2a)~(11-2-2f)改写成下面递推形式：m1212m0.5,12,12,12,120.5nnxxxxxxtHijkHijktm,12,1,12,0.5nnyyxxEijkEijkttz374,1,12,,12nnzzEijkEijky(11-2-4a)m1212m0.512,,1212,,120.5yynnyyyytHijkHijktm1,,12,,120.5nnzzyyEijkEijkttx12,,112,,nnxxEijkEijkz(11-2-4b)m1212m0.512,12,12,12,0.5nnzzzzzztHijkHijktm12,1,12,,0.5nnxxzzEijkEijktty1,12,,12,nnyyEijkEijkx(11-2-4c)10.512,,12,,0.5nnxxxxxxtEijkEijkt121212,12,12,12,0.5nnzzxxHijkHijktty121212,,1212,,12nnyyHijkHijkz(11-2-4d)10.5,12,,12,0.5yynnyyyytEijkEijkt121212,12,12,12,0.5nnzzyyHijkHijkttx1212,12,12,12,12nnxxHijkHijkz(11-2-4e)10.5,,12,,120.5nnzzzzzztEijkEijkt121212,,1212,,120.5nnyyzzHijkHijkttx3751212,12,12,12,12nnxxHijkHijky(11-2-4f)式中介质参数的空间坐标序号和对应的场分量序号相同。为了求解方程组(11-2-4)，需用适当的边界条件截断时域有限差分的计算区域。时域有限差分法涉及的边界条件通常包括理想电导体（PerfectElectricConductor，简称PEC），理想磁导体（PerfectMagneticConductor，简称PMC），吸收边界件(AbsorbingBoundaryCondition，简称ABC)[4-12]和周期边界件(PeriodicBoundaryCondition，简称PBC)[13-15]。除了截断计算区域外，还必须对计算区域内部的非均匀介质边界进行处理。如果非均匀介质边界与网格重合，根据电场和磁场的位置定义，磁场采样位于两个网格中心的连线上，与其相应的有效磁介质参数应该是两个对应网格磁参数的加权平均，而电场位于电网格棱边上。根据法拉第电磁感应定律，磁场的积分环路最多可以跨越电场周围的4个网格，因此与电场相联系的有效电介质参数应该是相关网格参数的加权平均。弯曲的金属或介质边界需要用共形技术处理[16-24]。由上递推方程可见，电场和磁场在时间和空间上都是交错的。因此，我们只需要存储前一时刻的电磁场值，而当前时刻的电（磁）场值将由前一时刻的电(磁)场值和前半时刻磁(电)场值求得。此外，三维电场和磁场数组的大小与空间格点数量相同，但是在时间上只需要保留最后一个时刻的电磁场值即可。这些特点对于实际编程是很有价值的。定义电磁场的位置序号的方法很多，通常采用的方法如11-2-2图所示。11-3网格数值色散已知在非色散媒质中，不同频率的电磁波具有相同的传播速度，此时电磁波角频率和波数的关系为2222xyzkkkc(11-3-1)式中xk、yk和zk分别是沿x、y和z方向的波数。然而，对于时域有限差分来说，即使在非色散媒质中，电磁波在不同方向上仍然可能具有不同的传播速度，这种特性称为网格数值色散。网格数值色散与网格大小、形状以及差分格式有关。下面讨论在直角坐标系中时域有限差分的网格色散误差。Ex(0,0)Ey(0,0)Hz(0,0)Ey(0,1)Ex(0,1)Ex(1,0)Ex(1,1)Ey(1,0)Ey(1,1)Ey(0,n-1)Ey(n,n-1)Ey(n,0)Hz(0,n-1)Hz(n-1,n-1)Ex(0,n)Ex(n-1,n)Ex(n-1,0)Hz(n-1,0)xy0图11-2-2在Yee差分各式中电磁场的相对位置376设一个平面波的波函数可以表示为j()0,,,exyztkxkykzxyzt(11-3-2)令x、y、z和t分别为空间和时间的采样间隔，那么，上式的离散形式可以表示为j()0,,exyzntkIxkJykKznIJK(11-3-3)式中n，I，J，K是时间和空间采样点的序号。已知自由空间中电磁场各分量满足下述波动方程22222222210xyzct(11-3-4)式中c为真空光速。利用时间和空间的中心差分公式离散上式，求得2221121,,2,,1,,,1,2,,,1,,,12,,,,1,,2,,,,nnnnnnnnnnnnijkijkijkxijkijkijkyijkijkijkzijkijkijk