时域有限差分法

~ 2 d c t sin k x ~ v g ~ x sint dk


1.5 数值稳定性（1）
• FDTD计算中每一步都是有误差的，随着时间步进，误 差会不断积累。如果误差的积累不会造成总误差的增 加，就成FDTD法是稳定的，否则成为不稳定的。数值 不稳定性会造成计算结果随时间步进无限增加。 • FDTD法是有条件稳定的，即：时间步必须必须小于一 定值以避免数值不稳定性。 • 本节的数值稳定性分析方法是建立在Courant等人几十 年前提出的经典方法基础上。这种方法首先把有限差 分算法分解为相互分离的时间和空间本征值问题。

x 2
4c 2

t 2
(1-24)
4
即
t
x c
可见，时间步长 t 必须是有界的。上式称为Courant稳 定性条件。有趣的是其上界恰好是魔时间步。
1.6 激励源的设置
在FDTD模拟电磁波传播时需要设置初始条件和激励 源。最简单的源设置方法是“硬源”, 即在激励源的位置 令 u满足ui=f(n), 常用的有 正弦函数 ui=sin(nt+) 高斯函数 ui=exp[-(n-n0)2/T2] 阶跃函数 ui= 0 n<n1 = ( n-n1)/(n2-n1) n1<n<n2 =1 n>n2 “硬源”设置简单，但当反射波回到“硬源”位置时， 会引起寄生反射，所以，要在这之前“关”掉源。 以后会有有关源设置的更详细讨论。
~ cosk x 1
~ 4c 2 0 因为 cos k x 1 ，所以 2 x

(1-23)
上式给出了差分网格中任意空间Fourier模的本征值谱。
1.5 数值稳定性（4）
•稳定性 为了保证任何空间模式的数值稳定性，(1-23)给出的 空间模式的本征值范围必须完全落在(1-20)所给出的时间 本征值的稳定范围内，于是
引言（2）
• 本课程采用研讨班形式。教师讲授 FDTD 的基本知识， 学生针对某一方向进行较深入的研究。 • 本讲我们考虑描述波动现象的最基本偏微分方程：一 维 标 量 波 动 方 程 的 数 值 FDTD 解 ， 为 以 后 二 维 、 三 维 Maxwell方程的FDTD分析奠定基础 • 课程内容取自下列的参考书和近年来相关的一些文献
uin
~ j nt k ix
将上式代入差分方程(1-6),得
2
~ ~ ct jk jt x jk x e 2e 2 e jt e x 重新组合并应用 Euler恒等式，最后得到数值色散关系为



(1-11)
~ ct cost cos k x 1 1 x
色散关系定义为行波的波长随频率的变化关系。为方便起见， 色散关系也常表示为行波的波数关于角频率的变化关系。 考虑(1.1)的正弦行波解 ux, t e j t kx 代入(1-1)得
k 即 c 上式便是一维标量波动方程的色散关系。 由上式得相速度
vp k c
1.5 数值稳定性（2）
将(1-19)代入(1-18)，有 qi2 [2 t 2 ]qi 1 0，于是
qi 2 t
2
2 t
2
2
4
2
a a2 1
算法稳定性要求 qi 1 。如果 a 2 1 0 ，则总有 qi a j 1 a 2 ，于是 qi a 2 1 a 2 1 ，满足稳定性 要求。这样可得 2 （1-20） 这就是稳定的数值差分解所要求的时间本征值谱。
j 2 c 2 jk 2
(1-8)
（1-9） 可见，相速与频率无关，称为非色散。非色散意味着对于具有任意 调制的包络或脉冲形状的波传播任意距离后波形保持不变。进一步 由(1-8)可以得到群速关系 d vg c (1-10) dk 这种情况下，群速也是与频率无关。
1.2 数值色散关系(2)
xi ,t n
（1-3）
1.1 差分近似（2）
于是，有 同理，有
2u x
2 xi ,t n

u in1 2u in u in1 x
2
O x 2



(1-4） (1-5）
2u t
2 xi ,t n

u in 1 2u in u in 1 t
2
O t 2
时域有限差分法
第1讲 一维标量波动方程
引言（1）
• 1966年,K.S. Yee（美籍香港人）首先提出了FiniteDifference Time-Domain Method,并用于柱形金属柱 电磁散射分析。由于当时计算机技术还比较落后，这 一方法并未引起重视。
• 1972年,A.Taflovey应用FDTD研究了UHF和微波对人类 眼睛的穿透，以了解“微波白内障”的成因。Taflove 成功地应用和发展了Yee的FDTD算法。 • 80年代后期，随着高速大容量计算机的普及，FDTD法 得到了迅速发展。如今已应用于涉及波动现象的任何 领域。至今，FDTD法的研究与应用仍方兴未艾。
1.1 差分近似（1）
一维标量波动方程 上式的解为 采用Taylor 展开
2u t
2
c2
2u x
2
（1-1） （1-2）
u F ( x ct ) G( x ct )
u x 2 2 u u ( xi x, t n ) u xi ,tn x xi ,t n x 2! x 2 x 3 3 u x 4 4 u ,t 3 xi ,t n 3! x 4! x 4 1 n

t
4c
2
0
1.5 数值稳定性（3）
• 空间本征值问题
c2 2 x 2
numerical
u in u in
(源自文库-21)
代入中心差分公式，得
c
2
u in1 2u in u in1
令 u u0 e
~ jk ix
x2
u in
2c 2
(1-22)
，Eular公式可得 x 2
上述过程也可用于一维标量波动方程差分近似的数值色散分析。
， 设在离散空间点 xi , t n ,离散行波解为 uxi , t n e ~ 为存在于有限差分网格中的数值正弦波的波数。一般情况 式中，k 下，不同于连续物理波的波数。正是这种不同导致了数值相速和群 速偏离了精确解。进而导致了数值色散误差。
0 x ~ • 取 ct 则数值相速为v p 0.9873c 。相对误差为 ，x 2 10 -1.27%。如果物理波传播了 100 距离(100空间格)时，数值模拟波只
传播了98.73空间格，相位误差为45.720。
0 x ~ 0.9969c 。这时数值相速的相对 ，x • 取 ct 则v p 2 20 误差为0.31，减少了4倍。同样，当物理波传播了同样的 100 时
应当注意，在一般情况下(1-6)对时间或空间具有二阶精度。但对 于 ct x 1 的特殊情况，根据解(1-2)，可以证明
4u x
于是
2
4 ,t n
c4
4u t 4
xi ,
2 c x c 2 Ox
4u x
4
12
,t n
t 2
2

(1-12)
1.3 数值相速（1）
~ ~ k 类似于(1-9)，定义数值相速为 v p
由（1-12）可得 （1-13）
2 ~ 1 x 1 k cos 1 cost 1 x ct
可见数值相速与频率有关。因此，由FDTD得到的数值波是色散的。
[1]A.Taflove,Computational Electrodynamics The FiniteDifference Time-Domain Method, Artech Hourse，1995. [2]高本庆，时域有限差分法，国防工业出版社，1995. [3]葛德彪，闫玉波，电磁场时域有限差分法，西电出版社，2002
2





忽略高次项，便可得到求解的差分迭代公式。
1.1 差分近似（3）
n=0 在所有空间点给uin, uin-1(i=1:imax)赋初值
n=n+1
由(1-6)在所有空间点求uin+1(i=1:imax)
No
n>nmax? Yes 结束
图1.1 一维波动方程FDTD流程图
1.1 差分近似（4）
于是利用单向差分近似得到吸收边界条件，详细讨论 见后面章节。
结 论 1
本讲介绍了一维标量波动方程的FDTD求解过程： • 利用Taylor级数展开方法获取空间/时间导数的二阶 中心差分近似，从而得到具有二阶精度的方程数值解的 时间步进迭代公式。 • 一般情况下，数值解引入了寄生的数值色散。当空 间步长和时间步长非常小时，数值解逼近精确解。当时 间步长满足魔时间步条件时，数值解等于精确解。 • 空间步长和时间步长必须满足Courant稳定性条件才能 保证数值解的稳定性。
定义数值相速为 （1-14） • 情况1 非常细网格 利用正弦函数的一阶Taylor展开，可得 ~ 2 2~ （1-15） c t k x c k ~ v g x0 x t t 0 所以，群速与相速一样，在细网格条件下趋近精确解。这证明了 当空间步长和时间步长趋于零时，数值解变得精确。 • 情况2 魔时间步 将魔时间步条件和波数代入(1-14)，得 2 （1-16） c t sin c ct ~ vg c v g ct sint 再次验证了魔时间步下数值解等于精确解。
4u
4
12 t
xi ,
O t 2


所以，（1-6）中的两个剩余项抵消，得到了精确的数值差分公式 （1-7） ct x 1为“魔时间步”(Magic time 正因为有这样的奇妙特性， step).
uin1 uin1 uin1 uin1
1.2 数值色散关系(1)
(200空间格)，数值模拟传播了199.378格，相位误差为11.1960，也 减少了4倍。误差减少了4倍反映了差分算法是二阶精度的。
1.3 数值相速（2）
• 情况1：非常细网格 t 0, x 0
2 cos x 1 x 2 , 当 x 0 ，数值色散关系(1-12)变为 根据 ~ 2 2 2 t k x ct 1 1 1 1 2 2 x ~ ~ v 。 2 2 ~2 即， c k ，最后得 k k ，于是有 v p p c

所以，在非常细的网格条件下，差分解逼近精确解。
• 情况2： 魔时间步 ct x
(1-12)变为 cost cos k x ,即 t k x ，k
~ v v 所以， p p

~
~
~

c
k 。
。可见，魔时间步下差分解与精确解相同。
1.4 数值群速

上式称为二阶偏导数的二阶中心差分格式。 将它们代入(1-1)，得
u in 1 ct n n n n n 1 u i 1 2u i u i 1 2u i u i x （1-6） 2 2 2 2 t c O x O t
1.5 数值稳定性（1）
• 时间本征值问题
2 t 2
n n u u numerical i i
(1-17)
差分近似，得
u in 1 2u in u in 1
t 2
u in 1 u in
u in
(1-18)
定义不变增长因子
qi

u in u in 1
(1-19)
1.7 吸收边界条件
由于计算机容量所限，计算域必须是有限的。对于 理想电壁或磁壁的边界条件的设置是直接的。但如果模 拟的是“开”问题，就要设置截断边界。在截断边界上 要设置吸收边界条件，使得电磁波可以被完全吸收，模 拟波无反射的通过吸收边界。 对于一维问题，采用单向波方程
u u c t x