第1页 / 共8页
第2页 / 共8页
第3页 / 共8页
第4页 / 共8页
第5页 / 共8页
第6页 / 共8页
第7页 / 共8页
第8页 / 共8页
FDTD法模拟TE波的传播.doc
FDTD 法模拟 TE 波的传播
电子院 12 级博士生 雷国伟 学号：2012010199
一、问题提出
模拟二维空间中 TM 波的传播过程。网格空间为 200×200 二维空间，在中
心点（101，101）设置一个线源，令其激发一窄脉冲波。当在该网格空间中执行
TE 波的时域有限差分格式时，就能在网格空间中模拟该脉冲波的传播过程。
二、问题分析
图 2-1 Yee 氏网格及其电磁场分量分布
二维 Yee 网格如图 2-1 所示。运行时时间步长为⊿t=2.5×10-17s，而空间步
⊿s=2v⊿t，满足稳定条件。V 为自由空间中的光速。激励源的脉冲分别为高斯脉
冲和正弦脉冲。
考虑 TE 波的情况。在笛卡儿坐标系中 TE 波的电磁场只有 3 个分量 Ex, Ey,
Hz，媒质的电导率和磁导率分别为σ和σ *，麦克斯韦方程组可以写成（二维）
E
y
x
采用 PML 吸收边界条件，对麦克斯韦方程进行分裂式处理，得 TE 波的时
E
x
y
H
y
H
x
H
t
E
t
E
t
E
E
*
H
x
y
z
x
0
0
z
y
z
0
z
域有限差分方程如下：
E
1
n
xy
(
i
,21
),
kj
CA
)(
(
iEj
n
xy
y
,21
),
kj
CB
)(
j
y
21
H
n
Z
(
i
,21
j
Hk
),21
y
21
n
Z
(
i
,21
j
),21
k
(2-1a)
1
n
xz
(
,21
i
E
),
kj
H
,(
jiEkCA
,21
k
21
)(
,(
ji
n
yz
z
n
x
H
)(
kCB
,(
ji
z
21
n
x
),21
k
)21
z
,21
k
)21
(2-1b)
1
n
yz
,(
ji
E
),21
k
H
,(
jiEkCA
,21
k
21
)(
,(
ji
n
xy
z
n
x
H
)(
kCB
,(
ji
z
21
n
Z
),21
k
)21
z
,21
k
)21
(2-1c)
1
n
yx
,(
ji
E
),21
k
H
1
n
zx
,(
,
E
kji
)21
H
n
z
x
n
yx
)(
(
i
,(
),21
jiEiCA
k
),21
,21
j
Hk
21
x
)(
iCB
x
(
i
21
n
Z
x
n
y
n
zx
)(
(
i
,(
,21
,
)(
)21
kjiEiCA
iCB
x
,
)21
(
kj
i
H
21
21
x
n
y
,21
j
),21
k
(2-1d)
,21
,
kj
)21
(2-1e)
1
n
zy
,(
i
E
,
kj
)21
H
CA
21
n
x
y
(
,(
i
n
zy
)
,(
iEj
,21
j
,
kj
k
CB
x
H
21
n
x
)21
)21
y
(
j
,(
i
)
j
,21
k
)21
(2-1f)
21
n
xy
,(
ji
H
,21
k
DB
y
)21
(
j
DA
y
)21
(
)21
j
H
,(
,1
jiE
n
z
21
n
xy
k
k
,(
,(
,21
)21
ji
)21
,
kjiE
y
n
z
)21
(2-1g)
21
n
xz
,(
,21
ji
H
)21
k
(
kDB
z
21
n
yz
(
,21
,
)21
H
i
kj
(
kDB
z
)21
(
)21
,(
kDA
H
ji
21
n
xz
,(
,21
)1
jiE
k
n
y
z
z
,21
k
,(
jiE
n
y
)21
),21
k
(2-1h)
(
kDA
(
iE
)21
H
,21
n
x
)21
z
(
21
n
yz
,
kj
,
kj
(
iE
n
x
,21
i
)1
z
)21
,21
),
kj
(2-1i)
21
n
yx
(
,21
,
)21
i
H
kj
(
iDB
)21
(
iDA
(
iE
)21
,1
x
n
z
x
H
,
kj
21
n
yx
,
kj
,(
iE
n
z
(
,21
i
)21
x
)21
,
kj
)21
(2-1j)
H
21
n
zx
(
i
,21
j
),21
k
(
iDB
x
)21
(
iDA
n
y
x
(
iE
)21
H
,1
j
j
),21
k
,(
iE
j
),21
k
(2-1k)
,21
n
y
(
i
21
n
zx
),21
k
x
21
n
zy
(
,21
i
H
j
DB
),21
k
j
(
y
DA
y
)21
(
j
(
iE
n
x
)21
H
,21
21
n
zy
j
j
(
iE
n
x
(
,21
i
),1
k
y
),21
k
,21
),
kj
(2-1l)
式中
CA
q
)(
2
)(
t
2
)(
tt
q
q
，
CB
q
)(
2
t
)(
2
q
(2-2)
tt
DA
q
(
)21
2
2
(
(
m
m
)21
)21
t
t
，
DB
q
(
)21
2
t
2
(
m
(2-3)
)21
t
,
i
,
kj
q
,
,
zyx
；
其中，
使用 PML 吸收边界条件时，首先要选定三个参数：PML 层数 N ，电导率
)(R 。图 2-2 为二维空间时 PML 的空间布
分布阶数 m ，PML 表面反射系数
局。
compute zone
perfect metal
PML media
图 2-2 二维情况 PML 的空间布局
三、程序流程图及说明
1. 程序流程图
建立的二维网格空间为 200×200 二维空间，在中心点（101，101）设置一个
线源，令其激发一窄脉冲波。
(1)脉冲波为高斯脉冲
H z
H
0
exp
(
t
2
)
t
0
2
T
H z
H
0
sin(
(
t
t
))
0
(2)脉冲为正弦脉冲
程序的流程图如图 3-1 所示
2. 程序说明
图 3-1 程序流程图
首先需要定义频率、波长、时间步数、真空区域网格数、总网格数、导电率
等参数。再对 TE 波的分量初始化。然后对电磁场迭代计算。最后对结果进行绘
图。
四、结果图
图 4-1 高斯脉冲场源辐射二维效果图
图 4-2 正弦时谐场源辐射二维效果图
五、分析与结论
通过对时域差分法的学习而进行的关于二维 FDTD 的 TE 波仿真，巩固了迭
代式的形式与相应的循环编程技巧。但也存在一定问题，由于在 PML 吸收上的
编程思路不是很规范，所以在吸收效果上与理论上的完全匹配层吸收效果还有
差距，但是已经基本达到了仿真的预期目的。
六、附录
% 二维 FDTD TE 波仿真
clear all;
pi=3.1415;
c=3.0e10;
f=1.0e15;
lambda=c/f;
nmax=400;
del_s=lambda/20;
del_t=0.5*del_s/c;
n=182;
np=9;
N1=n+2*np;
N=N1+1;
M=4;
sigma_max=(M+1)/1.50/pi/del_s; %最大导电率
% TE 波的分量初始化
tic;
figure(1);
axis([0 N 0 N -0.5 0.5]);
%高斯制下光速
%频率
%波长
%时间步数
%每最小波长 20 个采样点
%迭代时间步长
%真空区域网格数
%pml 层数
%总网格数
%采样点数
%导电率渐变指数
Ex=zeros(N1,N);
Ey=zeros(N,N1);
Bz=zeros(N1,N1); %矩阵行为纵向网格数，矩阵列为横向网格数,循环中用 j 表示行数，i
表示列数
%x 方向为横向，采样点为网格的横向边，故行数+1
%y 方向为纵向，采样点为网格的纵向边，故列数+1
Bzx=zeros(N1,N1);
Bzy=zeros(N1,N1);
Bzxx=zeros(nmax,2);
%进入电磁场迭代计算
for tt=1:nmax
for i=1:N1
if i>=np+1&&i<=N1-np
di=0;
elseif i<=np
di=np-i+0.5;
elseif i>=N1-np+1
di=np+i-N1-0.5;
end
sigma_mx=sigma_max*(di/np)^M;
%di 是采样点横向距 PML 内边界的距离
for j=1:N1
if j>=np+1&&j<=N1-np
dj=0;
elseif j<=np
dj=np-j+0.5;
elseif j>=N1-np+1
dj=np+j-N1-0.5;
end
%dj 是采样点纵向距 PML 内边界的距离
sigma_my=sigma_max*(dj/np)^M;
if sigma_mx+sigma_my==0
%真空区
if j==100&&i==100
%可选择高斯脉冲
t=30;
term=(tt-t);
pulse=exp(-4*pi*term^2/20^2);
pulse=sin(2*pi*tt/40); %可选正弦时谐源
c_miu=c*del_t/del_s;
Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j));
Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j));
Bz(i,j)=Bz(i,j)+Eterm1-Eterm2+pulse;%加入脉冲源
%
else
c_miu=c*del_t/del_s;
Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j));
Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j));
Bz(i,j)=Bz(i,j)+Eterm1-Eterm2;
end
else
%PML 区
cpm=(1-2*c*sigma_mx*del_t)/(1+2*c*sigma_mx*del_t);
cqm=c*del_t/(1+2*c*sigma_mx*del_t)/del_s;
Bzx(i,j)=cpm*Bzx(i,j)-cqm*(Ey(i+1,j)-Ey(i,j));
cpm=(1-2*c*sigma_my*del_t)/(1+2*c*sigma_my*del_t);
cqm=c*del_t/(1+2*c*sigma_my*del_t)/del_s;
Bzy(i,j)=cpm*Bzy(i,j)+cqm*(Ex(i,j+1)-Ex(i,j));
Bz(i,j)=Bzx(i,j)+Bzy(i,j);
end
end
end
for i=2:N1
if i>=np+1&&i<=N1-np
di=0;
elseif i<=np
di=np-i+1;
elseif i>=N1-np+1
di=np+i-N1-1;
end
sigma_ex=sigma_max*(di/np)^M;
for j=1:N1
%di 是采样点横向距 PML 内边界的距离
cam=(1-2*c*sigma_ex*del_t)/(1+2*c*sigma_ex*del_t);
cbm=c*del_t/(1+2*c*sigma_ex*del_t)/del_s;
Ey(i,j)=cam*Ey(i,j)-cbm*(Bz(i,j)-Bz(i-1,j));
end
end
for i=1:N1
for j=2:N1
if j>=np+1&&j<=N1-np
dj=0;
elseif j<=np
dj=np-j+1;
elseif j>=N1-np+1
dj=np+j-N1-1;
%dj 是采样点纵向距 PML 内边界的距离
end
sigma_ey=sigma_max*(dj/np)^M;
cam=(1-2*c*sigma_ey*del_t)/(1+2*c*sigma_ey*del_t);
cbm=c*del_t/(1+2*c*sigma_ey*del_t)/del_s;
Ex(i,j)=cam*Ex(i,j)+cbm*(Bz(i,j)-Bz(i,j-1));
%对靠近边界的中央磁场点采样
%可视化处理
%磁场的幅值
end
end
Bzxx(tt,1)=Bz(90,50);
Bzxx(tt,2)=Bz(50,90);
figure(1);
clf;
mesh(Bz);
axis([0 N 0 N -0.5 0.5]);
xlabel('i')
ylabel('j')
drawnow;
end
figure(2);
title('磁场幅值分布图');
surface(Bz);
shading interp;
axis square;
toc
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
