fdtd讲座1125

1、FDTD时域有限差分法一.时域有限差分法概述二.时域有限差分法的特点三.非分裂完美吸收边界FDTD方法 的研究一一.时域有限差分法概述时域有限差分法概述1.FDTD法的发展背景法的发展背景 FDTD方法最初由K.S.Yee提出，它是一种电磁场数值计算，通过将时间和空间离散化来直接求解麦克斯韦方程组的方法。对麦克斯韦方程的偏微分形式进行中心差分。 尽管Yee早在1966年发表了他的论文，但是这种方法没有立刻得到认可，这不仅因为当时计算机性能（CPU速度和内存等）的限制，而且缺少一种有效的高精度的吸收边界条件（ABC）来实现自由空间的无反射截断。每一个磁场分量由四个电场分量环绕，同样每一个电场分量

2、由四个磁场分量环绕。电场和磁场在时间顺序上交替抽样，抽样时间间隔彼此相差半个时间步不需要进行矩阵求逆运算，由给定相应电磁问题的初始值，就可以逐步推进地求得以后各个时刻空间电磁场的分布 这种电磁场分量的空间取样方式不仅符合法拉弟感应定律和安培环路定律的自然结构，而且这种电磁场各分量的空间相对位置也适合于麦克斯韦方程的差分计算，能够恰当地描述电磁场的传播特性。 作为一种电磁场的数值计算方法，时域有限差分法具有一些非常突出的特点，也是它的优点。 正是由于这些，使得越来越多的人对它产生了浓厚的兴趣，并得到越来越广泛的应用。 二二.时域有限差分法的特点时域有限差分法的特点 (1)直接时域计算直接时域计算

3、 时域有限差分法直接把含时间变量的 Maxwell旋度方程在Yee氏网格空间中转换为差分方程。 在这种差分格式中每个网格点上的电场(或磁场)分量仅与它相邻的磁场(或电场)分量及上一时间步该点的场值有关。 在每一时间步计算网格空间各点的电场和磁场分量，随着时间步的推进，即能直接模拟电磁波的传播及其与物体的相互作用过程。 这一特点使它能直接给出非常丰富的电磁场问题的时域信息，给复杂的物理过程描绘出清晰的物理图像如果需要频域信息，则只需对时域信息进行傅立叶变换，为获得宽频带的信息，只需在宽频谱的脉冲（高斯脉冲）激励下进行一次计算 (2)广泛的适用性广泛的适用性 由于时域有限差分法的直接出发点是概括电

4、磁场普通规律的Maxwell方程，这就预示着这一方法应具有最广泛的适用性。 从具体的算法看，在时域有限差分法的差分格式中被模拟空间电磁性质的参量是按空间网格给出的，因此，只需设定相应空间点以适当的参数，就可模拟各种复杂的电磁结构。 媒质的非均匀性、各向异性、色散特性和非线性等均能很容易地进行精确模拟。 由于在网格空间中电场和磁场分量是被交叉放置的，而且计算中用差分代替了微商，使得介质交界面上的边界条件能自然得到满足，这就为模拟复杂的结构提供了极大的方便。 (3)节约存储空间和计算时间节约存储空间和计算时间 在时域有限差分法中每个网格电场和磁场的六个分量及其上一时间步的值是必须存储的。此外还有描

5、述各网格电磁性质的参数以及吸收边界条件和连接条件的有关参量，它们一般是空间网格总数N的数倍。 所以，时域有限差分法所需要的存储空间直接由所需的网格空间决定，与网格总数N成正比。 在计算时，每个网格的电磁场都按同样的差分格式计算，所以，就所需的主要计算时间而言，也是与网格总数N成正比 相比之下，若离散单元也是N，则矩量法所需的存储空间与3N的2次方成正比，而所需的CPU时间则与(3N)的2次方至(3N)的3次方成正比 (4)适合并行计算适合并行计算 很多复杂的电磁场问题不能计算往往不是没有可选用的方法，而是由于计算条件的限制。 当代电子计算机的发展方向是运用并行处理技术，以进一步提高计算速度。

6、(5)计算程序的通用性计算程序的通用性 由于Maxwell方程是时域有限差分法计算任何问题的数学模型，因而它的基本差分方程对广泛的问题是不变的。 因此一个基础的时域有限差分法计算程序，对广泛的电磁场问题具有通用性，对不同的问题或不同的计算对象只需修改有关部分，而大部分是共同的。 频率升高? (6)简单、直观、容易掌握简单、直观、容易掌握 由于时域有限差分法直接从Maxwell方程出发，不需要任何导出方程，这样就避免了使用更多的数学工具，使得它成为所有电磁场的计算方法中最简单的一种。 其次，由于它能直接在时域中模拟电磁波的传播及其与物体作用的物理过程，所以它又是非常直观的一种方法。 由于它既简单

7、又直观，掌握它就不是件很困难的事情，只要有电磁场电磁场的基本理论知识，不需要数学上的很多准备，就可以学习运用这一方法解决很复杂的电磁场问题。这样，这一方法很容易得到推广，并在很广泛的领域发挥作用。1.麦克斯韦方程组和麦克斯韦方程组和Yee算法算法在无源、线性、各向同性、非色散介质中，麦克斯韦方程组表示为： （1.1a） （1.1b）*HEHt EHEt 三三.非分裂完美吸收边界FDTD方法的研究 在三维直角坐标系中，方程1.1a和1.1b可表示成如下微分形式： （1.2a）1yxzHEEtxz（1.2b）1yxzEEHtyx（1.2c）1yxzEHEtzy1yxzxHEHEtyz（1.2d）1

8、yxzyEHHEtzx（1.2e）1yxzzHHEEtxy（1.2d）1.1 1.1 一维一维MaxwellMaxwell方程的方程的FDTDFDTD算法（算法（1 1）一维自由空间Maxwell方程利用一阶导数的二阶中心差分近似，上面的方程变为01 (1.1)yxHEtz 01 (1.2)yxHEtz (1.3) )()(1)()(212102/12/1zkHkHtkEkEnynynxnx )() 1(1)()(2/12/1021211zkEkEtkHkHnxnxnyny(1.4) n 时间步长k 空间步长1.1 1.1 一维一维MaxwellMaxwell方程的方程的FDTDFDTD算法（

9、算法（2 2）电场与磁场分量的空间-时间分布图2/1nxE2k1kk1k2k2/1nxE2k1kk1k2k2/3k2/1k2/1k2/5k3/ 2k nxH0 xtn c 0023xxtcc 稳定条件 时间步长t与空间步长x，y和z之间必须满足一定条件，否则将出现数值不稳定性1.1 1.1 一维一维MaxwellMaxwell方程的方程的FDTDFDTD算法（算法（3 3）001/c 1/21/211220()()( )( )1nnnnyyxxHkHkEkEktx 11/21/211220()()(1)( )1 nnnnyyxxHkHkEkEktx 1/21/211220( )( )()()n

10、nnnxxyytEkEkHkHkx11/21/211220()()(1)( ) nnnnyyxxtHkHkEkEkxEE00由于介电常数和磁导率不在一个数量级1/21/2112200( )( )()()nnnnxxyytEkEkHkHkx 11/21/2112200()()(1)( ) nnnnyyxxtHkHkEkEkx 1.1 1.1 一维一维MaxwellMaxwell方程的方程的FDTDFDTD算法（算法（3 3） 用计算机语言表示的FDTD公式式中，时间变量已隐含在迭代公式中, 只要给定了所有空间点上电和磁场的初值，就可以一步一步地求出任意时刻所有空间点上的电场和磁场值。 1*5 .

11、 0kExkExkHykHy 0.5* 1Ex kEx kHy kHy k(1.7) (1.8) 1/21/21122( )( )0.5 ()()nnnnxxyyEkEkHkHk11/21/21122()()0.5 (1)( ) nnnnyyxxHkHkEkEk一维边界条件 102nxnxEE0023xxtcc 20 xtc距离稳定性条件电磁波传播距离设置的边界条件一维电磁波在介质中传播1/21/2112200( )( )()()nnnnxxyyrtEkEkHkHkx 11/21/2112200()()(1)( ) nnnnyyxxtHkHkEkEkx 1*5 . 0kExkExkHykHy

12、* 1Ex kEx kCb kHy kHy k 0.5/Cb kepsilon (1.20a)EHJt 01 (1.20b)HEt JE在有耗媒质中001rrEHEt 00( )( )1( )yxxrrHtE tE ttz 000( )( )1( )yxxrrHtE tE ttz 00( )( )1yxHtE ttz 1/21/21/21/2000(1/ 2)(1/ 2)( )( )( )( )12nnnnnnyyxxxxrrHkHkEkEkEkEktx 我们对方程进行时域和空间的近似一维电磁波在有耗介质中传播一维电磁波在有耗介质中传播00112tx 1/21/2001/ 2( )1( )1(

13、1/ 2)(1/ 2)22nnnnxxyyrrrttEkEkHkHk 1/21/21/21/2000(1/ 2)(1/ 2)( )( )( )( )12nnnnnnyyxxxxrrHkHkEkEkEkEktx 1/21/21/21/20( )( )1( )( )(1/ 2)(1/ 2)22nnnnnnxxxxyyrrEkEktEkEkHkHk 1/21/2000(1)21/ 2( )( )(1/ 2)(1/ 2)(1)(1)22nnnnrxxyyrrrtEkEkHkHktt 因此我们得到计算机的计算公式 *(1 )ex kca kex kcb khy khy k 0.5*( 1)hy khy

14、kex kex k一维电磁波在有耗介质中传播一维电磁波在有耗介质中传播1/21/2000(1)21/ 2( )( )(1/ 2)(1/ 2)(1)(1)22nnnnrxxyyrrrtEkEkHkHktt 因此我们得到计算机的计算公式 *(1 )ex kca kex kcb khy khy k 0.5*( 1)hy khy kex kex k*/(2*)eafdt sigmaepsz epsilon (1.)/(1.)ca keafeaf 0.5/(*(1.)cb kepsiloneaf一维电磁波在有耗介质中传播一维电磁波在有耗介质中传播对于任意媒质，麦克斯韦方程组的微分形式为 )()()(*0

15、EDr1.2 1.2 一维有耗媒质的一维有耗媒质的FDTDFDTD算法算法(2.1a) (2.1b) (2.1c) （2.1b）写成了频域形式DHt 01HEtEE00DD001HtD001)()()(*EDrEtH100令这样和我们上一节提到的公式是一样的 ，只是D变成了E。1.2 1.2 一维有耗媒质的一维有耗媒质的FDTDFDTD算法算法(2.2a) (2.2b) (2.3a) (2.3b) (2.3c) *0()rrj 对于有耗媒质 ,我们可表示为下式0()()()rDEEj00( )( )( )trD tE tE tdt从傅立叶变换理论，我们知道频域1/j是时域的积分，由于我们是对频

16、域的采样，上式变为00nnnirtDEE1.2 1.2 一维有耗媒质的一维有耗媒质的FDTDFDTD算法算法(2.4) (2.5) 2.4代入2.3b(2.6) 0crj，我们分离电场,得到公式1000nnnnirttDEEE(2.7) 1000nninrtDEEt(2.8a) E n是当前的电压值和当前的D值和以前E的和我们令00nnitIE(2.8b) 公式2.8可简化为0nnnrDIEt10nnntIIE(2.8c) 1.2 1.2 一维有耗媒质的一维有耗媒质的FDTDFDTD算法算法现在,FDTD公式可以表示为:dxk=dxk+0.5*(hyk-1-hyk)exk=gaxk *(dxk

17、-ixk)ixk=ixk+gbxk*exkhyk=hyk+0.5*(exk-exk+1)其中,gaxk=1/(epsilon+(sigma*dt/epsz)gbxk=sigma*dt/epsz所有记录媒质信息都包含在公式2.9b,2.9c中.对于自由空间,gax=1 和gbx=0；对于有耗媒质，gax和gbx由公式2.10计算。公式2.9a和2.9d，其中包含空间差分，并不随媒质而改变。1.2 1.2 一维有耗媒质的一维有耗媒质的FDTDFDTD算法算法(2.9a) (2.9b) (2.9c) (2.9d) (2.10a) (2.10b) 0tIDErnnn10nnntIIE0( )()( )

18、rDEj我们可以避免在时间域上处理棘手的卷积积分, 我们可以利用Z变换,10( )( )( )( )rtD zE zz I zE z01/( )( )( )1rtD zE zE zz1010/( )( )( )( )1ttI zE zz I zE zz10( )( )( )rD zz I zE zt10nnnrDIEt10nnntIIE1.3 1.3 利用利用Z Z变换公式变换公式HtD001)()()(*EDrEtH100(2.3a) (2.3b) (2.3c) 1.3 1.3 二维二维MaxwellMaxwell方程的方程的FDTDFDTD算法算法 (1)(1)对于二维仿真，可以用TE和T

19、M波表示。对于TM波对于二维仿真，可以用TE和TM波表示。对于TM波，公式2.3可以简化为下式)(100yHxHtDxyz)()()(*zrzEDyEtHzx001xEtHzy001，其中02 cxt1.3 1.3 二维二维MaxwellMaxwell方程的方程的FDTDFDTD算法算法 (1)(1)(3.1a) (3.1b) (3.1c) (3.1d) 我们可以得到公式： dz(i,j)=dz(i,j)+0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1);ez(i,j)=gaz(i,j)*(dz(i,j)-iz(i,j); iz (i,j)=iz(i,j)+gb

20、z(i,j)*ez(I,j)hx(i,j)= hx(i,j)+0.5*( ez(i,j)-ez(i,j+1);hy(i,j)= hy(i,j)+0.5*( ez(i+1,j)-ez(i,j);xH和zD与一维的简单的有损介质表达方式是一样的 1.3 1.3 二维二维MaxwellMaxwell方程的方程的FDTDFDTD算法算法 yH(3.2a) (3.2a) (3.2c) (3.2d) (3.2e) 基本理论：如果一束在媒质A中的波，接触的媒质B，反射系数的大小通过两种媒质的阻抗描述通常我们假设磁导率为常数，所以，当波照射从 到 时将会产生反射。但当 随介电常数改变时， 将保持为常数， 为零

21、，没有反射发生。但这并不能解决我们的问题，进入媒质的波会继续传播。这就需要一种有耗媒质，波在到达边界之前，在这种媒质中很快衰减。我们可以通过设置磁导率和介电常数为复数，因为虚数部分代表着衰减。 BABA1.4 PML1.4 PML吸收边界吸收边界 （1 1）121.4 PML1.4 PML吸收边界吸收边界 （2 2）)(0yHxHcDjxyz)()()(*zrzEDyEcHjzx0 xEcHjzy0我们把公式3.1转化为复数形式的麦克斯韦方程(3.3a) (3.3b) (3.3c) (3.3d) 对于3.3 式，我们引入有耗媒质，其介电常数和磁导率定义为 *Fz*Fx*Fy)()()(0*yH

22、xHcyxDjxyFzFzz)()()(*zrzEDyEcyxHjzFyFyx0*)()(xEcyxHjzFyFyy0*)()(（1） 只与电位移D有关，与E无关；（2）我们在PML中增加的有耗媒质与真实的介电 常数没有关系。 1.4 PML1.4 PML吸收边界吸收边界 （3 3）(3.4) (3.5a) (3.5b) (3.5c) (3.5d) F(引入的有耗媒质只对吸收边界起作用） 1.4.1 实现PML的两个基本条件 (1)Sack提出的构造PML的两个基本条件1、对于PML来说，媒质1的阻抗应该是一常数，1*0FxFx阻抗为1是归一化 2、垂直边界的方向，介电常数和相对磁导率的值必须

23、是相反的 *1FyFx*1FyFx0*mDFmFm0*mHFmFm我们假设介电常数和相对磁导率复数形式为：下面的公式满足3.7式，其中1FmFm000DmHmD(3.6) (3.7a) (3.7b) (3.8a) (3.8b) (3.9a) (3.9b) 引入的参数1.4.1 实现PML的两个基本条件 (2)我们把3.9代入3.61/ )(1/ )(100*0jxjxFxFx满足第一个条件。如果 值随着进入PML区而逐渐增大，那么电场和磁场将逐渐衰减。如果我们只在X方向应用PML，公式可以简化只有X方向的)()(0*yHxHcxDjxyFzzxEcxHjzFyy0*)(*Fx*FxyEcxHj

24、zFyx0*)(3.10a) (3.10b) (3.10c) 1.4.1 实现PML的两个基本条件 (3)公式3.9代入3.10)()(1 (00yHxHcDjxjxyzDyEcHjxjzxD010)(1 (xEcHjxjzyD00)(1 (我们注意到公式3.11b和3.11cHx和公式Hy的磁导率相反，满足第二条件。下面我们把公式3.11a左边转化为差分形式：2),(),()(),(),()(2/12/102/12/10jiDjiDitjiDjiDDitDnznzDnznzzDz2)(1 1),(2)(1 1),(02/102/1titjiDtitjiDDnzDnz(3.11a) (3.11

25、b) (3.12) 求Dz与Dz-1的关系我们令02/)(11)(2tiigiD002/)(12/)(1)(3titiigiDD同样，我们可以处理3.11c，给出下式其中),2/1(5 .0)(2),()(3),(2/12/1jiHigijiDigijiDnynznz)2/1,()2/1,(),2/1(jiHjiHjiHnxnxny),(), 1(5 . 0) 2/ 1( 2),()( 3),(2/12/11jiEjiEifijiHifijiHnznznynz， 02/)2/1(11)2/1(2tiifiD002/)2/1(12/)2/1(1)2/1(3titiifiDD其中 i+1/2，是磁

26、场和电场相差1/2网格1.42 PML公式的推导 (1)(3.12a) (3.12b) (3.13) (3.14) (3.15a) (3.15b) 公式3.11b转化为下式00( )1zDzxExEjHcyjy ，其中xecurlxjiEjiEyEnznzz_),()1,(2/12/1_)(_)2/1,()2/1,(0001TnDnxnxxecurltxxecurlctjiHjiH) 2/ 1,(2)(_) 2/ 1,() 2/ 1,(2/ 1001jiItxecurlxtcjiHjiHnHxDnxnx2/10 xtc02)()()( 1tiixnifiD在时域代表积分，j/1j在时域代表差分

27、，上式可化为) 2/ 1,()( 1_5 . 0) 2/ 1,() 2/ 1,(2/ 11jiIifiecurljiHjiHnHxnxnx1.42 PML公式的推导 (2)(3.16) (3.17) )()(1 (00yHxHcDjxjxyzDyEcHjxjzxD010)(1 (xEcHjxjzyD00)(1 (1/2( ,)(1 / 2,)( ,1 / 2)( ,1 / 2)nnnnzyxxDijHijHijHij),(), 1(5 . 0) 2/ 1( 2),()( 3),(2/12/11jiEjiEifijiHifijiHnznznynz) 2/ 1,()( 1_5 . 0) 2/ 1,

28、() 2/ 1,(2/ 11jiIifiecurljiHjiHnHxnxnx最后引入匹配参数后的二维差分方程为了计算f 和g ，我们没有必要去计算实际的电导率，引入辅助参数，3)_(*333. 0pmllengthixn pmllengthi_,.,2 , 1)()(1ixnifi)(11()(2ixnigi)(1)(1()(3ixnixnigi)(1 ifi)(2 igi)(3 igi这些辅助参数有一定的变化范围0 0.333 1 0.75 1 0.5对于主要问题空间， =0， =1， =1。因此从主要问题空间到 PML空间是无缝转换。)( 1 ifi)(2 igi)( 3igi1.42 P

29、ML公式的推导 (3)(3.18) (3.19a) (3.19b) (3.19c) PML区折叠部分12 ,3 , 2 , 3fjfj fj g j g j减小，增大 PML 边界参数关系 以上我们应用PML在方向。显然，在把y方向加到其中。)()(1)()(1 (000yHxHcDjyjxjxyzDDyEcHjyjxjzzDD0010)(1 ()(1 (xEcHjyjxjzyDD0100)(1)()(1 (和前面的处理过程相同，使用下面公式代替3.13), 2/ 1(5 . 0)( 2)( 2),()( 3)( 3),(2/12/1jiHjgjigijiDjgjigijiDnynznz)2/

30、1,()2/1,(), 2/1(jiHjiHjiHnxnxny1.42 PML公式的推导 (4)(3.20a) (3.20b) (3.20c) (3.21) 1.42 PML公式的推导 (5)， 在Y方向的可给出如下式), 2/1()( 1_5 . 0)2/1(2), 2/1()2/1(3), 2/1(2/11jiIjfjecurlifijiHifijiHnHynynyecurljiIjiInHxnHx_)2/1,()2/1,(2/12/1),() 1,(_2/12/1jiEjiEecurlnznz在X方向的可给出如下式 ecurljiIjiInHxnHx_)2/1,()2/1,(2/12/1

31、)2/1,()( 3)2/1,(1jiHjfjjiHnxnx(3.22a) (3.22b) (3.22c) )2/1,()( 1_5 . 0)(22/1jiIifiecurljfjnHx(3.23a) (3.23b) (3.23c) 1.5 1.5 三维三维MaxwellMaxwell方程的方程的FDTDFDTD算法（算法（1 1） Yee 首先将空间按立方体分割，电磁场的六个分量在空间的取样点分别放在立方体的边沿和表面中心点上，电场与磁场分量在任何方向始终相差半个网格步长。Ey( i, j, k )xyzEzEyEyEzEzHxExExExHzHy下面是三维麦克斯韦方程 、 )(100zHy

32、HtDyzx)(100 xHzHtDzxy)(100yHxHtDxyz)(100yEzEtHzyx)(100zExEtHxzy)(100 xEyEtHyxz对于三维，和二维处理PML边界完全相同，只是我们处理三个方向公式3.11变成下式)()(1)()(1)()(1 (01000yHxHcDjzjyjxjxyzzyx)()(1 ()(1)()(1 (0000yHxHjzcDjyjxjxyzzyx1.5 1.5 三维三维MaxwellMaxwell方程的方程的FDTDFDTD算法（算法（1 1）(5.1) (5.2) 比二维麦克斯韦方程 ，多出一项我们定义1.5 1.5 三维三维MaxwellM

33、axwell方程的方程的FDTDFDTD算法（算法（1 1）hcurljzchcurlcz_1)(_000)()(1 ()(1)()(1 (0000yHxHjzcDjyjxjxyzzyxhcurljIDz_1)(_()(1)()(1 (0000DzzzyxIzhcurlcDjyjxj公式右边比二维多出一项(5.5) (5.4) (5.3) 1.6总场散射场方法 总场/散射场方法是基于Maxwell方程的线性特性和下列电磁场的分解 totincscattotincscatEEEHHH式中下标inc表示不存在任何材料时的入射波场。下标scat表示散射波场，最初它们是未知的，它们是入射波与材料相互作用产生的场。下标tot表示总场。根据Maxwell方程的线性特性，无论是入射场、散射场或总场都满足Maxwell方程，所以FDTD法可以独立地应用于入射场、散射场和总场。如图4-3所示，将计算域分成两个区：区域1和区域2。区域1中包含了所有散射体。在区域1中用FDTD法模拟总场，称为总场区。区域2中为自由空间的一部分