FDTD原理及例子

1、FDTD数值分析法数值分析法 目录目录1.麦克斯韦方程的基础知识2.一维和三维Maxwell方程的Yee算法3.数值稳定性分析4.吸收边界条件5.波源的设置6.编程思路基础知识基础知识麦克斯韦方程微分形式：0DHJtBEtBD BPM方式基础知识基础知识麦克斯韦方程微分形式：0DHJtBEtBD FDTD方式将时间进行差分，并且磁场与电场交替迭代更新,eemmmeDHJJEtBEJJsHtBD 时谐场形式：00()( , , , )( , , , )( , , , )( , , , )j tj tHjjEEjHE x y z tE x y z t eH x y z tHx y z t e 对于

2、有耗媒质：一维Maxwell方程的Yee算法一维Maxwell方程，介质参数和场量均与x，y无关(无损，电导率和磁导率为0)利用一阶导数的二阶中心差分近似，上面的方程变为zHtEyrx01zEtHxry01zkHkHktkEkEnynyrnxnx)()()(1)()(212/1212/101zkEkEktkHkHnxnxrnyny)() 1()(1)()(210212/1212/1一维Maxwell方程的Yee算法采用归一化磁场 使得电场与归一化磁场有相同的数量级，于是可以得到FDTD迭代公式为式中， 为自由空间中的光速。HH00)() 1()()21()21(212/12/1kEkEkztc

3、kHkHnxnxrnyny)21()21()()()(2/12/11kHkHkztckEkEnynyrnxnx001c一维Maxwell方程的Yee算法 用计算机语言表示的FDTD公式式中，时间变量已隐含在迭代公式中，以及 只要给定了所有空间点上电/磁场的初值，就可以一步一步地求出任意时刻所有空间点上的电/磁场值。 1*kExkExkcakHykHy 1*kHykHykcbkExkEx)(;)(;);21(21kztckcbkztckcakEkExkHkHyrrxy一维Maxwell方程的Yee算法xExExEyHyH0123k0n2n2/3n1n2/1n)() 1()()21()21(212

4、/12/1kEkEkztckHkHnxnxrnyny)21()21()()()(2/12/11kHkHkztckEkEnynyrnxnx三维Maxwell方程的Yee算法Ey( i, j, k )xyzEzEyEyEzEzHxExExExHzHyFDTD离散中的Yee元胞1.每一个磁场分量由四个电场分量环绕；每一个电场分量由四个磁场分量环绕2.这种空间取样方式符合法拉第感应定律和安培环路定律3.电场和磁场在时间上交替抽样，抽样时间彼此相差半个时间步三维Maxwell方程的Yee算法Maxwell旋度方程为tHHE 1xzyxHyEzEtHyxzyHzExEtH1 1zyxzHxEyEtHxyz

5、xEzHyHtE1yzxyExHzHtE1zxyzEyHxHtE1tEEH三维Maxwell方程的Yee算法11221111,2222211,22nnnxxi jki jkxi jkHHHOttt2,21,1,21,21,21,zOzEEzEnkjiynkjiynkjiy221,21, 1,21,21,yOyEEyEnkjiznkjiznkjiz采用时间平均近似22121,21,2121,21,21,21,2tOHHHnkjixnkjixnkjixynkjiEnkjiEznkjiEnkjiEkjiDnkjiHkjiDnkjiHzzyybHxxaHxx, 1, 1,1,最后得：可以看到有个2阶小

6、量，这个要忽略掉，因此是近似的三维Maxwell方程的Yee算法同理，可以得到其他2个磁场分量的FDTD方程znkjiEnkjiExnkjiEnkjiEkjiDnkjiHkjiDnkjiHxxzzbHyyaHyy, 1, 1,1,xnkjiEnkjiEynkjiEnkjiEkjiDnkjiHkjiDnkjiHyyxxbHzzaHzz, 1, 1,1,所在空间位置wHaHwttD22所在空间位置wHbHwttD22三维Maxwell方程的Yee算法利用对偶原理： ，并注意到E与H在时间上差半个步长，可以直接从磁场FDTD公式得到电场的FDTD公式。如：，HEEH,bbaaCDCD,ynkjiHn

7、kjiHznkjiHnkjiHkjiCnkjiEkjiCnkjiEzzyybExxaExx1,1, 1,1,1, 1,1,所在空间位置wEaEwttC22所在空间位置wEbEwttC22三维Maxwell方程的Yee算法媒质参数赋值在所有空间点给电磁场分量赋初值求所有空间离散点上n+1时间步的磁场求所有空间离散点上n+1时间步的电场n=n+1nnmax结 束NoYes三维Maxwell方程的Yee算法 介绍了求解矢量Maxwell方程的FDTD Yee算法，归纳起来，Yee算法的主要特点有： 1）Yee算法采用耦合的Maxwell旋度方程，同时在时间和空间求解电场和磁场，而不是采用波动方程只求

8、解电场或磁场。 2） Yee网格在三维空间这样安排E和H分量，使得每一个E或H分量由四个H或E循环的分量所环绕。 3） Yee算法以蛙跳算法在时间上安排E和H分量。在某一时刻，使用前一时刻的E数据计算所有H分量。然后，再使用刚计算的H数据计算所有的E分量。如此循环，直至完成时间步进过程。数值稳定性问题（1）FDTD计算中每一步都是有误差的，随着时间步进，误差会不断积累。如果误差的积累不会造成总误差的增加，就成FDTD法是稳定的，否则成为不稳定的。数值不稳定性会造成计算结果随时间步进无限增加。（2）FDTD法是有条件稳定的，即：时间步必须必须小于一定值以避免数值不稳定性。（3）数值稳定性分析方法

9、是建立在Courant等人几十年前提出的经典方法基础上。这种方法首先把有限差分算法分解为相互分离的时间和空间本征值问题。 数值稳定性问题22111yxct二维的Yee算法数值Courant稳定性条件三维的Yee算法数值Courant稳定性条件2221111zyxct便于理解，当Yee元胞为立方体时，此两式表明时间间隔必须等于或小于波以光速通过yee元胞对角线长度1/3或1/2所需时间差分近似所带来的数值色散一维标量波动方程为例22222xuctu设在离散空间点 ,离散行波解为nitx ,xiktnjninietxuu,将上式代入差分方程得：22221112122tOxOctuuuuuxtcun

10、inininininitjxk jxk jtjeeextce222最后得色散关系1cos1cos121ttcxxk我们想要k与没加元胞之前的k一样，但是差分近似后k变的不一样差分近似所带来的数值色散将平面波带入差分方程所出现的稳定性、色散、以及各向异性，这些特性并非由介质的物理特性所引起，而是数值计算中的差分近似所致，在FDTD数值计算中，稳定性、色散、各向异性将影响计算精度。吸收边界条件 理论上说，求解空间是无限大的，但是由于计算的数据容量问题，需要在有限空间的周围做特殊处理，使得向边界面行进的波在边界处保持“外向行进”，无明显的反射现象，并且不会使内部空间的场产生畸变。 1）Mur吸收边界

11、条件 2）Berenger完美匹配层（PML） 3）各向异性匹配层吸收边界条件Berenger完美匹配层（PML） 只考虑二维TE情况，对于二维TM和三维情况可采用类似方法进行分析。 分析范围在PML层内yHEtEzxx0 xHEtEzyy0 xEyEHtHyxzz*0*分别表示自由空间中的电导率和磁阻率 分别表示自由空间中的电导率和磁阻率 Berenger为了引入规定损耗和阻抗匹配的新自由度，将 Hz分裂为两个分量 和 ，即 同时引入了新的电导率 和磁损耗 ，并规定TE情形的四个场分量（而不是通常的3个）由下列方程耦合在一起：zxHzyzxzHHHzyHHzHzzHyx,*,yxBereng

12、er完美匹配层（PML）yHHEtEzyzxxyx0 xHHEtEzyzxyxy0 xEHtHyzxxzx*0yEHtHxzyyzy*0换句话说，Berenger构造了一般新的非物理媒质（称为PML媒质），在该媒质中场满足左端的方程（并不一定是Maxwell方程） 经过一系列公式推导：最后得出无反射匹配条件0*00*0 yyxxGHEZ000022sincosyxWWG00Z表明PML媒质中的波阻抗与入射角无关，与真空中波阻抗相同，这意味着从真空中的任意角度入射到PML媒质交界面时将会无反射地进入PML媒质中，并在PML媒质中衰减地传播 Berenger完美匹配层（PML）1.若分界面垂直于x

13、轴（y轴），要求二者具有相同的横向电导率和磁导率，且横向和纵向电导率、磁导率均满足阻抗匹配条件。Berenger完美匹配层（PML）如何计算：任何普通介质均可以视为特殊的PML介质，但麦克斯韦方程只需计算三个分量，而不是PML的四个分量。另外由于在介质中电磁波衰减的很快，常规FDTD中Yee的差分格式已不再适用，需要指数差分0011122,1,1/2,1/2,11xxnnnnttyyzzi jijijijxEeEeHHx0011122,1 2,1 2,1,11yynnnnttxxzzi ji ji ji jyEeEeHHy波源的设置F 为了用FDTD法模拟电磁场工程问题，必须在FDTD网格中引

14、入电磁波激励源。F 常用的波源种类： 平面波源用于电磁散射问题； 导波源用于微波网络参数计算； 电流源或电压源用于微波电路或天线的激励。F 常用的设置方式： 初始条件； 硬源； 电流源； 总场/散射场公式硬波源F 硬波源可以通过规定FDTD空间网格中电场或磁场分量满足所希望的时间函数进行简单地设置。如：在网格点（is,js,ks)建立Ez的硬源为:F 硬波源相当于点源，模拟从源点向外辐射的具有与源函数性质相同的数值波。如果材料结构放在离波源一定距离的地方，辐射的数值波最终传播到这个结构，一部分传输过去，一部分反射。在原理上，时间步进可以连续进行，直到所有的瞬态波消失。0n 0,tnfEEnkj

15、izsss硬波源F 常用源波形有以下三种F 正弦波源为时谐源，提供了频率为f0的正弦波。F FDTD模拟时，一般采用Gaussian脉冲信号激励，因为它可以提供宽频带特性。 Gaussian脉冲信号的频谱也是Gaussian函数。脉冲带通脉冲正弦波Gaussian 2sin Gaussian 2sin00)()(02020tnnfeetnftnfdecaydecaynnnnnn)4/(exp)(222tntnGdecaydecay硬波源 但是，在连续时间步进时，当反射波到达波源位置，由于波源处总场已被规定，并没有考虑网格中可能的反射波（所以称为硬源），硬源对这些反射波就会造成寄生的非物理的再次

16、反射。当规定了一个表面上的总场而没有考虑反射场值时，寄生反射就一定会发生。 解决这一问题的一个简单方法是当反射波到达波源位置时去掉波源。这就要求源离散射体足够远，使散射波到达源位置时源入射波已基本为零，所以这种方法不适合正弦入射波源。编程思路（二维TE+PML边界条件）第一步：基本参数设置，包扩介电系数，磁导系数，频率，波长等编程思路（二维TE+PML边界条件）第二步：网格参数设置编程思路（二维TE+PML边界条件）第三步：材料属性以及激励源设置编程思路（二维TE+PML边界条件）第四步：初始场设置（无激励），矩阵形式编程思路（二维TE+PML边界条件）第五步：中心场区的系数设置编程思路（二维TE+PML边界条件）第六步：PML区的设置，指数差分系数设置以Front层为例，按照相似的办法定义back，left，right，PML边界层编程思路（二维TE+PML边界条件）第七步：movie初始化编程思路（二维TE+PML边界条件）第八步：电磁场随时间迭代运算，之前都是初始化从电场迭代到磁场迭代过程，中间还需要将ex在PML层更新编程思路（二维TE