时域有限差分方法林志立

J E
（电极化）
（磁化）
（欧姆定律）
James Clerk Maxwell
(1831–1879)
一切电磁场和电磁波问题均可由以上方程，以及各类 具体的边界条件所决定！
Beihang University
麦克斯韦方程组中的运算符
散度（Divergence) 旋度（Curl)
连续函数的偏微分运算
Beihang University
r
(t)

2
(t)

(s )02 02 02
e0t
sin(
02 02 t)U (t)
频域：D() 0 r () E()
更新方程：
r
()



( s 02
)02 j20 2
b0 Dn b1Dn1 b2 Dn2
% Define the Gaussian source waveform
5.定义场源
time = dt*[0:number_of_time_steps-1].';
Jz_waveform = exp(-((time-2e-10)/5e-11).^2);
source_position_index = round(nx*source_position/domain_size)+1;
(1)电场在时间上取整数倍的 Δt;
t=n *Δt;
(2)磁场在时间上取(整数＋ 1/2)倍的Δt;
t=(n +1/2)*Δt;
Beihang University
麦克斯韦方程的离散化近似
以Hz为例：
Hz 1 ( Ex E y )
t z y x
H
n1/ z
2
(i,
j)

H
n1 z
FDTD的基本思想
－时域和空间域的离散化 －连续偏微分的有限阶近似
设有一连续函数 f ( x) , 现欲求 f '(x) 。 二阶中心差分近似表达式：
f '(x) f (x x) f (x x) 2x
当 x越小时，上式的近似程度 越高。
f (x)
实际上：
f '( x) f ( x x) f ( x x) x2 f '''( x) ...
/
2
(
i,
j)
t

1
(
E
n x
(
i,
j
1)

E
n x
(i,
j)

E
n y
(i
1,
j)

E
n y
(
i,
j) )
z
y
x
H
n1/ z
2
(i,
j)

H
n1 z
/
2
(
i,
j)

t
(
E
n x
(
i,
j
1)

E
n x
(
i,
j)

E
n y
(i
1,
j)

E
n y
(
i,
j) )
z
y
x
上式即为Hz的更新方程，由前一时刻的磁场和前半时刻 的临近空间格点的电场即可求出最新时刻的磁场。
2x
6
Beihang University
FDTD空间域的离散化
(1)空间域的分割离散化
节点
Yee元胞（Δx, Δy, Δz)
Ex分量的空间离散分布图
Hx分量的空间离散分布图
Beihang University
FDTD空间域的离散化
YEE 元胞
例如： Hz 1 ( Ex E y )
矩量法 (Method of Moments) 适合于细线、平面形状结构的电磁场问题
Beihang University
电磁学基本方程
麦克斯韦方程组( Maxwell Equations)
(安培环路定律)
(法拉第感应定律)
(高斯定律－电场) (高斯定律－磁场)
物质本构关系 (Constitutive Relations)
0 (a0 E n a1E n1 a2 E n )
Z-变换
近似， 求系数
r ()

D E

a0 b0
a1Z 1 a2Z 2 b1Z 1 b2Z 2

a0 b0
a1 e j a2 e2 j b1 e j b2 e2 j
,
Beihang University
FDTD的离散参数的稳定性条件
★ 时间步长：Courant 稳定性条件
Ｚ域数值色散方程：（von Neumann method ）
Z1/2 Z 1/2
( c0t
)2r (Z)

x, y,z
4sin2 (k 2
/
2)

0,
为了保持稳定性，该方程的所有解的模必须小于1。
t z y x
Hz为Ex和Ey所环绕。
各电磁场分量在元胞中的位置
K.S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,”IEEE Trans. Antennas Propagat., vol. 14, 1966, pp. 302-307.
z (i, j) y x
球体，物质III 不同的元胞的电磁参量应设
空气
长方体，物质II 长方体，物质I
置为所在空间所代表的介质 的介电常数和磁导率。
场量与介质参数要对应
Beihang University
色散介质的FDTD模拟
以Lorentz 介质为例：
时域：D(t) 0 r (t) * E(t)
磁场磁流部分
Beihang University
编程举例1：一维FDTD问题（续）
% Calculate FDTD updating coefficients Ceze = (2 * eps_r_z * eps_0 - dt * sigma_e_z) ...
./(2 * eps_r_z * eps_0 + dt * sigma_e_z);
% 1D problem space length in meters
dx = 1e-3;
% cell size in meters
dt = 3e-12;
% duration of time step in seconds
number_of_time_steps = 2000; % number of iterations
色散介质的FDTD模拟
模拟Lorentz色散介质的不同方法： {a0, a1, a2,b0,b1,b2}
Error r () r ()
MSE approach :
E n [(b0 Dn b1Dn1 b0 Dn2 ) / 0
(a1E n1 a2 E n2 )] / a0
% Initialize field and material arrays Ceze = zeros(nx+1,1);
3.初始化场量和介质参量阵列
Cezhy = zeros(nx+1,1);
Cezj = zeros(nx+1,1);
Ez = zeros(nx+1,1);
电场电流部分
Jz = zeros(nx+1,1);
1,
j)

H
n1/ x
2
(
i,
j)

H
n1/ x
2
(i,
j
1) )
z
x
y
Beihang University
FDTD的离散参数的选择
★元胞尺寸：边长小于最短波长的1/10, 以减小数值色散。
数值色散方程
c02
sin2 ( k ) 2
( t )2
Beihang University
FDTD空间偏微分的近似
以Hz为例：
Hz 1 ( Ex E y )
t z y x
Ex y

Ex (i,
j
1) Ex (i, 2( y)
j)
2
E y x

Ey(i
1, j) Ey(i, 2( x )
j)
2
类似地，可实现各电磁场分量的 空间偏微分计算。
4.计算更新方程系数
Cezhy = (2 * dt / dx) ... ./(2 * eps_r_z * eps_0 + dt * sigma_e_z);
电场部分
Cezj = (-2 * dt) ... ./(2 * eps_r_z * eps_0 + dt * sigma_e_z);
Chyh = (2 * mu_r_y * mu_0 - dt * sigma_m_y) ... ./(2 * mu_r_y * mu_0 + dt * sigma_m_y);
对于非色散介质，时间步长不能大于以下表达式：
t
1
c
1 x2

1 y2

1 z2
Beihang University
介质电磁参量的设定
Ez (i, j) 1 ( H y H x )
t
z (i, j) x y