时域有限差分法

时域有限差分法第1页，课件共21页，创作于2023年2月引言（1）1966年,K.S.Yee（美籍香港人）首先提出了Finite-DifferenceTime-DomainMethod,并用于柱形金属柱电磁散射分析。由于当时计算机技术还比较落后，这一方法并未引起重视。1972年,A.Taflovey应用FDTD研究了UHF和微波对人类眼睛的穿透，以了解“微波白内障”的成因。Taflove成功地应用和发展了Yee的FDTD算法。80年代后期，随着高速大容量计算机的普及，FDTD法得到了迅速发展。如今已应用于涉及波动现象的任何领域。至今，FDTD法的研究与应用仍方兴未艾。第2页，课件共21页，创作于2023年2月引言（2）

本课程采用研讨班形式。教师讲授FDTD的基本知识，学生针对某一方向进行较深入的研究。本讲我们考虑描述波动现象的最基本偏微分方程：一维标量波动方程的数值FDTD解，为以后二维、三维Maxwell方程的FDTD分析奠定基础课程内容取自下列的参考书和近年来相关的一些文献

[1]A.Taflove,ComputationalElectrodynamicsTheFinite-DifferenceTime-DomainMethod,ArtechHourse，1995.[2]高本庆，时域有限差分法，国防工业出版社，1995.[3]葛德彪，闫玉波，电磁场时域有限差分法，西电出版社，2002第3页，课件共21页，创作于2023年2月1.1差分近似（1）一维标量波动方程（1-1）上式的解为（1-2）

采用Taylor展开（1-3）第4页，课件共21页，创作于2023年2月1.1差分近似（2）于是，有(1-4）同理，有(1-5）上式称为二阶偏导数的二阶中心差分格式。将它们代入(1-1)，得（1-6）忽略高次项，便可得到求解的差分迭代公式。第5页，课件共21页，创作于2023年2月1.1差分近似（3）NoYesn=0在所有空间点给uin,uin-1(i=1:imax)赋初值n=n+1由(1-6)在所有空间点求uin+1(i=1:imax)结束n>nmax?图1.1一维波动方程FDTD流程图第6页，课件共21页，创作于2023年2月

1.1差分近似（4）

应当注意，在一般情况下(1-6)对时间或空间具有二阶精度。但对于的特殊情况，根据解(1-2)，可以证明

于是

所以，（1-6）中的两个剩余项抵消，得到了精确的数值差分公式

（1-7）正因为有这样的奇妙特性，为“魔时间步”(Magictimestep).第7页，课件共21页，创作于2023年2月1.2数值色散关系(1)色散关系定义为行波的波长随频率的变化关系。为方便起见，色散关系也常表示为行波的波数关于角频率的变化关系。考虑(1.1)的正弦行波解代入(1-1)得即

(1-8)上式便是一维标量波动方程的色散关系。

由上式得相速度（1-9）可见，相速与频率无关，称为非色散。非色散意味着对于具有任意调制的包络或脉冲形状的波传播任意距离后波形保持不变。进一步由(1-8)可以得到群速关系

(1-10)这种情况下，群速也是与频率无关。第8页，课件共21页，创作于2023年2月1.2数值色散关系(2)

上述过程也可用于一维标量波动方程差分近似的数值色散分析。设在离散空间点,离散行波解为，式中，为存在于有限差分网格中的数值正弦波的波数。一般情况下，不同于连续物理波的波数。正是这种不同导致了数值相速和群速偏离了精确解。进而导致了数值色散误差。将上式代入差分方程(1-6),得

(1-11)重新组合并应用Euler恒等式，最后得到数值色散关系为

(1-12)第9页，课件共21页，创作于2023年2月1.3数值相速（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倍反映了差分算法是二阶精度的。第10页，课件共21页，创作于2023年2月1.3数值相速（2）情况1：非常细网格根据，数值色散关系(1-12)变为即，，最后得，于是有。所以，在非常细的网格条件下，差分解逼近精确解。情况2：魔时间步

(1-12)变为,即。所以，。可见，魔时间步下差分解与精确解相同。第11页，课件共21页，创作于2023年2月1.4数值群速定义数值相速为（1-14）情况1

非常细网格利用正弦函数的一阶Taylor展开，可得（1-15）所以，群速与相速一样，在细网格条件下趋近精确解。这证明了当空间步长和时间步长趋于零时，数值解变得精确。

情况2

魔时间步将魔时间步条件和波数代入(1-14)，得（1-16）再次验证了魔时间步下数值解等于精确解。第12页，课件共21页，创作于2023年2月1.5数值稳定性（1）FDTD计算中每一步都是有误差的，随着时间步进，误差会不断积累。如果误差的积累不会造成总误差的增加，就成FDTD法是稳定的，否则成为不稳定的。数值不稳定性会造成计算结果随时间步进无限增加。FDTD法是有条件稳定的，即：时间步必须必须小于一定值以避免数值不稳定性。本节的数值稳定性分析方法是建立在Courant等人几十年前提出的经典方法基础上。这种方法首先把有限差分算法分解为相互分离的时间和空间本征值问题。第13页，课件共21页，创作于2023年2月1.5数值稳定性（1）时间本征值问题

(1-17)差分近似，得(1-18)定义不变增长因子(1-19)第14页，课件共21页，创作于2023年2月1.5数值稳定性（2）将(1-19)代入(1-18)，有，于是

算法稳定性要求。如果，则总有，于是，满足稳定性要求。这样可得（1-20）这就是稳定的数值差分解所要求的时间本征值谱。第15页，课件共21页，创作于2023年2月1.5数值稳定性（3）空间本征值问题

(1-21)代入中心差分公式，得

(1-22)令，Eular公式可得因为，所以(1-23)上式给出了差分网格中任意空间Fourier模的本征值谱。第16页，课件共21页，创作于2023年2月1.5数值稳定性（4）稳定性为了保证任何空间模式的数值稳定性，(1-23)给出的空间模式的本征值范围必须完全落在(1-20)所给出的时间本征值的稳定范围内，于是即(1-24)可见，时间步长必须是有界的。上式称为Courant稳定性条件。有趣的是其上界恰好是魔时间步。第17页，课件共21页，创作于2023年2月1.6激励源的设置在FDTD模拟电磁波传播时需要设置初始条件和激励源。最简单的源设置方法是“硬源”,即在激励源的位置令u满足ui=f(n),常用的有正弦函数ui=sin(nt+)高斯函数ui=exp[-(n-n0)2/T2]阶跃函数ui=0n<n1=(n-n1)/(n2-n1)n1<n<n2=1n>n2“硬源”设置简单，但当反射波回到“硬源”位置时，会引起寄生反射，所以，要在这之前“关”掉源。以后会有有关源设置的更详细讨论。第18页，课件共21页，创作于2023年2月1.7吸收边界条件由于计算机容量所限，计算域必须是有限的。对于理想电壁或磁壁的边界条件的设置是直接的。但如果模拟的是“开”问题，就要设置截断边界。在截断边界上要设置吸收边界条件，使得电磁波可以被完全吸收，模拟波无反射的通过吸收边界。对于一维问题，采用单向波方程于是利用单向差分近似得到吸收边界条件，详细讨论见后面章节。第19页，课件共21页，创作于2023年2月结论1

本讲介绍了一维标量波动方程的FDTD求解过程：利用Taylor级数展开方法获取空间/时间导数的二阶中心差分近似，从而得到具有二阶精度的方程数值解的时间步进迭代公式。一般情况下，数值解引入了寄生的数值色散。当空间步长和时间步长非常小时，数值解逼近精确解。当时间步长满足魔时间步条件时，数值解等于精确解。空间步长和时间步长必须满足Courant稳定性条件才能保证数值解的稳定性。第20页，课件共21页，创作于2023年2月习题11.1利用Taylor级数展开方法分别推导一阶导数的二阶和四阶精度中心差分近似。1.2利用数值相速和群速公式分别画出数值相速和群速在，，和条件下关于网