1、时域有限差分法时域有限差分法第1讲 一维标量波动方程引言（引言（1）1966年,K.S.Yee（美籍香港人）首先提出了Finite-Difference Time-Domain Method,并用于柱形金属柱电磁散射分析。由于当时计算机技术还比较落后，这一方法并未引起重视。1972年,A.Taflovey应用FDTD研究了UHF和微波对人类眼睛的穿透，以了解“微波白内障”的成因。Taflove成功地应用和发展了Yee的FDTD算法。80年代后期，随着高速大容量计算机的普及，FDTD法得到了迅速发展。如今已应用于涉及波动现象的任何领域。至今，FDTD法的研究与应用仍方兴未艾。引言（引言（2）本课
2、程采用研讨班形式。教师讲授FDTD的基本知识，学生针对某一方向进行较深入的研究。本讲我们考虑描述波动现象的最基本偏微分方程：一维标量波动方程的数值FDTD解，为以后二维、三维Maxwell方程的FDTD分析奠定基础 课程内容取自下列的参考书和近年来相关的一些文献 1A.Taflove,Computational Electrodynamics The Finite-Difference Time-Domain Method,Artech Hourse，1995.2高本庆，时域有限差分法，国防工业出版社，1995.3葛德彪，闫玉波，电磁场时域有限差分法，西电出版社，20021.1 1.1 差分近
3、似（差分近似（1 1）一维标量波动方程 （1-1）上式的解为 （1-2）采用Taylor 展开（1-3）22222xuctu)()(ctxGctxFunninininittxtxtxtxnixuxxuxxuxxuxutxxu,444,333,222,1,!4!3!2),(1.1 1.1 差分近似（差分近似（2 2）于是，有 (1-4）同理，有 (1-5）上式称为二阶偏导数的二阶中心差分格式。将它们代入(1-1)，得（1-6）忽略高次项，便可得到求解的差分迭代公式。2211,222xOxuuuxunininitxni2211,222tOtuuutunininitxni22221112122tOx
4、Octuuuuuxtcunininininini1.1 1.1 差分近似（差分近似（3 3）NoYesn=0在所有空间点给uin,uin-1(i=1:imax)赋初值 n=n+1由(1-6)在所有空间点求uin+1(i=1:imax)结 束nnmax?图1.1 一维波动方程FDTD流程图 1.1 1.1 差分近似（差分近似（4 4）应当注意，在一般情况下(1-6)对时间或空间具有二阶精度。但对 于 的特殊情况，根据解(1-2)，可以证明 于是 所以，（1-6）中的两个剩余项抵消，得到了精确的数值差分公式 （1-7）正因为有这样的奇妙特性，为“魔时间步”(Magic time step).,44
5、4,44inxttucxu2,442,442221212tOtutxuxcxOcinxt1xtcninininiuuuu11111xtc1.2 1.2 数值色散关系数值色散关系(1)(1)色散关系定义为行波的波长随频率的变化关系。为方便起见，色散关系也常表示为行波的波数关于角频率的变化关系。考虑(1.1)的正弦行波解 代入(1-1)得 即 (1-8)上式便是一维标量波动方程的色散关系。由上式得相速度 （1-9）可见，相速与频率无关，称为非色散。非色散意味着对于具有任意调制的包络或脉冲形状的波传播任意距离后波形保持不变。进一步由(1-8)可以得到群速关系 (1-10)这种情况下，群速也是与频率无
6、关。kxtjetxu,222jkcjck ckvpcdkdvg1.2 1.2 数值色散关系数值色散关系(2)(2)上述过程也可用于一维标量波动方程差分近似的数值色散分析。设在离散空间点 ,离散行波解为 ，式中，为存在于有限差分网格中的数值正弦波的波数。一般情况下，不同于连续物理波的波数。正是这种不同导致了数值相速和群速偏离了精确解。进而导致了数值色散误差。将上式代入差分方程(1-6),得 (1-11)重新组合并应用 Euler恒等式，最后得到数值色散关系为 (1-12)nitx,xiktnjninietxuu,ktjxk jxk jtjeeextce22211coscos2xkxtct1.3
7、1.3 数值相速（数值相速（1 1）类似于(1-9)，定义数值相速为 由（1-12）可得（1-13）可见数值相速与频率有关。因此，由FDTD得到的数值波是色散的。取 则数值相速为 。相对误差为-1.27%。如果物理波传播了 距离(100空间格)时，数值模拟波只传播了98.73空间格，相位误差为45.720。取 则 。这时数值相速的相对误差为0.31，减少了4倍。同样，当物理波传播了同样的 时(200空间格)，数值模拟传播了199.378格，相位误差为11.1960，也减少了4倍。误差减少了4倍反映了差分算法是二阶精度的。kvp1cos1cos121ttcxxk1020 xxtc，cvp9873
8、.00102020 xxtc，cvp9969.00101.3 1.3 数值相速（数值相速（2 2）情况情况1 1：非常细网格 根据 ，数值色散关系(1-12)变为即，最后得 ，于是有 。所以，在非常细的网格条件下，差分解逼近精确解。情况情况2 2：魔时间步 (1-12)变为 ,即 。所以，。可见，魔时间步下差分解与精确解相同。0,0 xt0 ,21cos2xxx当112121222xkxtct222kckck ppvv xktcoscosxtckckxkt，ppvv 1.4 1.4 数值群速数值群速 定义数值相速为 （1-14）情况情况1 非常细网格 利用正弦函数的一阶Taylor展开，可得
9、（1-15）所以，群速与相速一样，在细网格条件下趋近精确解。这证明了 当空间步长和时间步长趋于零时，数值解变得精确。情况情况2 魔时间步 将魔时间步条件和波数代入(1-14)，得 （1-16）再次验证了魔时间步下数值解等于精确解。txkxtckddvgsinsin2kctxkxtcvtxg2200ggvcttcctctcvsinsin21.5 1.5 数值稳定性（数值稳定性（1 1）FDTD计算中每一步都是有误差的，随着时间步进，误差会不断积累。如果误差的积累不会造成总误差的增加，就成FDTD法是稳定的，否则成为不稳定的。数值不稳定性会造成计算结果随时间步进无限增加。FDTD法是有条件稳定的，
10、即：时间步必须必须小于一定值以避免数值不稳定性。本节的数值稳定性分析方法是建立在Courant等人几十年前提出的经典方法基础上。这种方法首先把有限差分算法分解为相互分离的时间和空间本征值问题。1.5 1.5 数值稳定性（数值稳定性（1 1）时间本征值问题时间本征值问题 (1-17)差分近似，得 (1-18)定义不变增长因子(1-19)nininumericaluut22ninininiutuuu211211ninininiiuuuuq1.5 1.5 数值稳定性（数值稳定性（2 2）将(1-19)代入(1-18)，有 ，于是 算法稳定性要求 。如果 ，则总有 ，于是 ，满足稳定性要求。这样可得（
11、1-20）这就是稳定的数值差分解所要求的时间本征值谱。1iq012 22iiqtq124222222aattqi012a21ajaqi1122aaqi0422tc1.5 1.5 数值稳定性（数值稳定性（3 3）空间本征值问题空间本征值问题 (1-21)代入中心差分公式，得 (1-22)令 ，Eular公式可得因为 ，所以 (1-23)上式给出了差分网格中任意空间Fourier模的本征值谱。nininumericaluuxc222ninininiuxuuuc21122xik jeuu01cos222xkxc1cosxk0422xc1.5 1.5 数值稳定性（数值稳定性（4 4）稳定性稳定性 为了
12、保证任何空间模式的数值稳定性，(1-23)给出的空间模式的本征值范围必须完全落在(1-20)所给出的时间本征值的稳定范围内，于是即 (1-24)可见，时间步长 必须是有界的。上式称为Courant稳定性条件。有趣的是其上界恰好是魔时间步。22244txccxtt1.6 1.6 激励源的设置激励源的设置 在FDTD模拟电磁波传播时需要设置初始条件和激励源。最简单的源设置方法是“硬源”,即在激励源的位置令 u满足ui=f(n),常用的有 正弦函数 ui=sin(nt+)高斯函数 ui=exp-(n-n0)2/T2 阶跃函数 ui=0 nn1 =(n-n1)/(n2-n1)n1nn2 “硬源”设置简
13、单，但当反射波回到“硬源”位置时，会引起寄生反射，所以，要在这之前“关”掉源。以后会有有关源设置的更详细讨论。1.7 1.7 吸收边界条件吸收边界条件 由于计算机容量所限，计算域必须是有限的。对于理想电壁或磁壁的边界条件的设置是直接的。但如果模拟的是“开”问题，就要设置截断边界。在截断边界上要设置吸收边界条件，使得电磁波可以被完全吸收，模拟波无反射的通过吸收边界。对于一维问题，采用单向波方程 于是利用单向差分近似得到吸收边界条件，详细讨论见后面章节。xuctu结结 论论 1 1 本讲介绍了一维标量波动方程的FDTD求解过程：利用Taylor级数展开方法获取空间/时间导数的二阶中心差分近似，从而
14、得到具有二阶精度的方程数值解的时间步进迭代公式。一般情况下，数值解引入了寄生的数值色散。当空间步长和时间步长非常小时，数值解逼近精确解。当时间步长满足魔时间步条件时，数值解等于精确解。空间步长和时间步长必须满足Courant稳定性条件才能保证数值解的稳定性。习习 题题 1 11.1 利用Taylor级数展开方法分别推导一阶导数 的二阶和四阶精度中心差分近似。1.2 利用数值相速和群速公式分别画出数值相速和群速在 ，和 条件下关于网格空间分辨率 的曲线，并进行相应的讨论。1.3 编写本讲介绍的一维标量波动方程FDTD求解程序。设在网格左边界波源为下列脉冲函数：a.高斯函数 b.矩形函数 在网格右边界 。时间总步数为5000,空间总步为200。在 ，和 条件下每隔1000个时间步画出每个传播脉冲关于位置的分布曲线。xu xtc99.0 xtc10.0 xtc50.0 xtc90.0 xR020)40(expnunnnun40 0400 1000nuxtc90.0 xtcxtc01.1
