电磁仿真实践报告一.doc

电子科技大学物理电子学院

标 准 实 验 报 告

（实验）课程名称 电磁仿真综合实践

电子科技大学教务处制表

电 子 科 技 大 学
实 验 报 告
学生姓名： xxx 学 号： xxxxxxx
指导教师： xxx,xx
实验地点：xxxxx 实验时间：2012年5月—6月
一、实验室名称：xxxxxxxx
二、实验项目名称：电磁仿真综合实践I级
三、实验学时：36学时
四、实验原理：电磁仿真的时域有限差分法
五、实验目的：加强学生的计算机综合应用能力、尤其是运用计算机分析和解决专业问题的能力培养，使学生对独立进行科学研究有初步实践；初步掌握一种纯数值电磁仿真方法——时域有限差分法；初步学会综合应用一种程序设计语言进行科学与工程计算；增强科技报告写作能力，学会相关软件使用。
六、实验内容：
1. 均匀平板传输线传输特性仿真
2. 带挡板的平板传输线传输特性仿真
七、实验器材（设备、元器件）： 电子计算机

八、实验步骤：
1．电磁仿真的时域有限差分法。
数值差分原理:
时空离散及连续取函数样
(1)
在一定体积内和一段时间上对连续电池场的数据取样压缩
差分原理
一阶差分
向前差分 误差O(h) (2)
向后差分 误差O(h) (3)
中心差分 误差O(h2) (4)
用中心差分代替对空间，时间的微分
(5)
(6)

场分量取样节点在空间和时间上采取交替排布，利用电生磁，磁生电的原理

如图1-1所示，Yee单元有以下特点：
（1）与分量在空间交叉放置，相互垂直；每一坐标平面上的分量四周由分量环绕，分量的四周由分量环绕；场分量均与坐标轴方向一致。
（2）每一个Yee元胞有8个节点，12条棱边，6个面。棱边上电场分量近似相等，用棱边的中心节点表示，平面上的磁场分量近似相等，用面的中心节点表示。
（3）每一场分量自身相距一个空间步长，和相距半个空间步长
（4）每一场分量自身相距一个时间步长，和相距半个时间步长
（5）3个空间方向上的时间 步长相等，以保证均匀介质中场量的空间变量与时间变量完全对称。

Maxwell方程FDTD的差分格式:

图8-1 Yee模型
麦克斯韦第一、二方程 （7）
式中，是电流密度，反映电损耗，是磁流密度，单位，反映磁损耗。主要与上式对应。各向同性介质中的本构关系：
（8）
是磁阻率，计算磁损耗的。
以为变量，在直角坐标中，展开麦克斯韦第一、二方程，分别为
（9）
（10）
令代表在直角坐标中的任何一个分量，离散符号取为
（11）
关于时间和空间的一阶偏导数取中心差分近似为
（12）
麦克斯韦方程组可以表示为

蛙跳格式

H0 H1/2 H3/2

E0 E1 E2

解的稳定性:
, （13）
数值色散:
在FDTD网格中,数值波模的传播速度将随频率而改变,即有色散。这种色散由数值网格引起,而非物理上客观存在。
为减小色散,使用中通常取空间步长Δ满足

2．均匀平行板传输线传输特性（一维FDTD）的仿真。
d=0.18m,L=6m,高斯脉冲T=0.5ns,t0=3T,fmax=1Ghz, λmin=0.3m
Δz=λmin/20,Nz=L/Δz
时域波形为高斯波形
数组实现Ex(Nz+1),Hy(Nz)
Yee元胞如图所示。

L
图一维Yee元胞

均匀平面波（TEM波）是一维问题，电磁波沿z轴方向传播，则，场量和介质参数均与x，y无关，即，麦克斯韦方程为
和 （14）
差分格式为
（15）
（16）
截断边界条件：
终端短路 Ex(z=0,t)=F(t)
终端匹配 向z=L传播的单向波方程
（17）
(1) 终端短路；
Ex(Nz)=0
(2) 终端匹配
Ex(Nz+1)=PEx(Nz)+(c*dt-dz)/(c*dt+dz)*(Ex(Nz)-Ex(Nz+1))

d s=d/3

L
3. 带挡板的平行板传输线传输特性（二维FDTD）的仿真.

d=0.18m,L=6m,高斯脉冲T=0.5ns，t0=3T，fmax=1Ghz， λmin=0.3m
Δz=λmin/20，Nz=L/Δz，Δt=0.01m，Nx=d/Δt
时域波形为高斯波形
数组实现Ex(Nx,Nz+1),Ez(Nx+1,Nz),Hy(Nx,Nz)
由于不连续性的出现，沿X方向场的均匀性招到破坏，沿y方向场的均匀性仍然保持，此时平行板传输线中的场分量有Ex，Ez，Hy。

吸收边界条件为
for m=1:Nx
Ex(m,Nz+1)=PEx(m,Nz)+(c*dt-dz)/(c*dt+dz)*(Ex(m,Nz)-Ex(m,Nz+1));
End
选取观测面分别为L/4,3L/4记录总电压，，
分离出入射电压,反射电压 ,透射电压

在DC-1GHz内计算Nf个频率处的值，频率取样步长Δf=1GHz/Nf，第k个取样点的值为，(k=0,1,2,….,Nf)
二端网络S参数

九、实验数据及结果分析：
1. 均匀平行板传输线（一维FDTD）的仿真（终端短路和匹配）
(1) 流程图

开始

输入各值

否
判断n是否大于Nt

是

切换左边边界条件

蛙跳

结束

(2) 仿真结果及分析
终端匹配时3L/4处电压前600dt时为零，因为波尚未传到3L/4处，又由于终端匹配，波形传到L处被吸收，所以无反射。得到的电压波形如图

终端匹配3L/4处电压时域波形
终端短路时3L/4处电压前600dt时为零，因为波尚未传到3L/4处，又由于终端短路，波形传到L处被反射，反射波形Ex反向因而电压为负。得到的电压波形如图

终端短路3L/4处电压时域波形

3. 带挡板的平行板传输线（二维FDTD）的仿真（终端匹配）。
(1) 流程图
开始

输入各值

判断n是否大于Nt1

否
判断n是否大于Nt1

否

是
是

切换左边边界条件

切换左边边界条件

蛙跳

蛙跳

结束

(2) 仿真结果及分析
由于边界匹配，所以当波传到边界时被吸收，没有反射。而当波经过中间的金属时，大部分波发生透射，但是透过挡板的波最终传至边界时被吸收不会反射回来。刚开始时由于波尚未传至挡板因而反射系数为零，透射系数为一。一段时间后随着频率增大波逐渐被反射，因而透射系数由一逐渐减小，而反射系数由零逐渐增大，如下图所示。

L/4处电压时域波形

3L/4处电压时域波形

S参数图形（S11是’b’，S21是’r’，验证‘g’）
3. 测试题
仿真结果及分析

由于边界短路，且中间所加挡板为金属，所以当波传到挡板及边界时波形发生反射，然而当波传到中间挡板时，由于挡板而发生反射的波很小，而当波传到边界时波全部反射，再传到挡板时又有小部分波发生反射，因而反射系数S11如下图所示刚刚开始十分接近于一，越到后面抖动幅度越大。

L/4处电压时域波形

3L/4处电压时域波形

S参数图形（S11是’b’）
十、实验结论：
均匀平板传输线，电磁波沿z轴方向传播只存在Y方向的磁场和X方向的电场,当终端匹配时波传至终端被吸收,当终端短路时, 波传至终端反射回来且Ex反向。
带挡板的平行板传输线，由于不连续性的出现，沿X方向场的均匀性招到破坏，沿y方向场的均匀性仍然保持，此时平行板传输线中的场分量有Ex，Ez，Hy。
当终端匹配时波传至终端被吸收,当终端短路时, 波传至终端反射回来且Ex反向。
十一、总结及心得体会：（作文1篇，1000字以上）
通过此次实验,以及对电磁仿真的学习才体会到自己对电磁场与波以及微波技术基础学习上的不足。以前学习电磁场与波以及微波技术基础时只是理论上的理解，最终的目的只是通过考试、学会做题，然而对其实际的应用却并不十分了解，也不懂得如何在实际中去设计微波器件。总之学过之后却无法真正用于实际。而通过此次实验，才真正学会了如何将学到的知识用于实际，解决问题。
通过本次的仿真实验我的收获很大。以前学习过matlab的知识以及对其的应用但是在此次实验中却发现许多的知识都忘了，一些命令也不会运用，而通过此次运用matlab进行电磁仿真，不仅让我学会如何运用matlab进行电磁仿真而且熟悉并加深了matlab知识，提高了自己的计算机综合应用能力。其次，通过此次实验熟悉并加深了对以前学过的专业知识的理解，如均匀平板传输线的传输特性，终端匹配及短路时波形如何传播；带挡板的平板传输线传输特性终端匹配及短路时波形如何传播等等。另外此次运用电脑仿真经过自己的亲自实践不仅掌握了更多的知识，而且也提高对专业知识学习的兴趣。
通过此次实验也提高了自己分析问题和解决问题的能力以及独立工作的能力。在此次实验过程中计算均匀平板传输线的传输特性时，匹配边界要如何添加,短路边界如何添加；带挡板的平行板传输线又要怎样实现边界匹配等等，通过自己不断的思考,不断完善得到最终结果。使得自己解决问题的能力有了很大的提高。另外，在此次实验过程中通过自己的努力，运用自己所学的知识来解决问题，不仅以前的知识得到巩固，并且也发现更容易学到知识，理解问题。
总之觉得学校开这门课程是十分必要的，这种类似于做项目的实践经验在我们平时的学习中是很难得到的。

十二、对本实验过程及方法、手段的改进建议：
报告评分：
指导教师签字：

附件

18

附件1：一维FDTD终端匹配仿真源代码
c=3e8; % ×ÔÓÉ¿Õ¼ä¹âËÙ
mu0=4*pi*1e-7; % ×ÔÓÉ¿Õ¼ä´Åµ¼ÂÊ
eps0=8.85e-12; % ×ÔÓÉ¿Õ¼ä½éµç³£Êý
d=0.18;l=6;T=0.5e-9;
t0=3*T;fmax=1e9;
bc=0.3;dz=bc/20;
Nz=l/dz;dt=dz/(2*c);
Nt=6*T/dt+50;
Nt1=6*T/dt+4*Nz;
Ex=zeros(1,Nz+1);
Hy=zeros(1,Nz);
for n=1:Nt1
if n<Nt
Ex(1)=exp(-(n*dt-t0)^2/T^2); % ¼Ó¼¤Àø
else Ex(1)=PEx(2)+(c*dt-dz)/(c*dt+dz)*(Ex(1)-Ex(2));
end
for k=1:Nz
Hy(k)=Hy(k)+((dt/mu0)*(Ex(k)-Ex(k+1))/dz);
end
PEx=Ex;
for k=2:Nz
Ex(k)=Ex(k)+((dt/eps0)*(Hy(k-1)-Hy(k))/dz);
end
Ex(Nz+1)=PEx(Nz)+(c*dt-dz)/(c*dt+dz)*(Ex(Nz)-Ex(Nz+1));
% PExÊDZÈExÔçÒ»¸öʱ¼ä²½Ê±¿ÌµÄµç³¡Öµ
Ver(n)=d*Ex(Nz*3/4);
figure(1);plot(Ex)
axis([0 400 -1 1])
figure(2);plot(Hy)
pause(0.0001)
axis([0 400 -0.01 0.01])

end
figure(3); plot((1:Nt1)*dt,Ver)
XLabel('Time (s)','FontSize',15,'FontWeight','b');
YLabel('Ver (V)','FontSize',15,'FontWeight','b');
axis([0 4.5e-8 -0.2 0.2]

附件2：二维FDTD终端匹配仿真源代码

clear
clc
c=3e8; % ×ÔÓÉ¿Õ¼ä¹âËÙ
mu0=4*pi*1e-7; % ×ÔÓÉ¿Õ¼ä´Åµ¼ÂÊ
eps0=8.85e-12; % ×ÔÓÉ¿Õ¼ä½éµç³£Êý
d=0.18;l=6;T=0.5e-9;
t0=3*T;fmax=1e9;
bc=0.3;dz=bc/20;
Nx=d/dz;
Nz=l/dz;dt=dz/(2*c);
Nt1=6*T/dt+100;
Nt=6*T/dt+800+100;
Ex=zeros(Nx,Nz+1);
Ez=zeros(Nx+1,Nz);
Hy=zeros(Nx,Nz);
Ex1=zeros(Nx,Nz+1);
Ez1=zeros(Nx+1,Nz);
Hy1=zeros(Nx,Nz);
Vt1=zeros(1,Nt);
Vt2=zeros(1,Nt);
PEx=zeros(Nx,Nz+1);
PEx1=zeros(Nx,Nz+1);
for n=1:Nt
if n<=Nt1
Ex(:,1)=exp(-(n*dt-t0)^2/T^2); % ¼Ó¼¤Àø
else
for m=1:Nx Ex(m,1)=PEx(m,2)+(c*dt-dz)/(c*dt+dz)*(Ex(m,2)-Ex(m,1));
end
end
for i=1:Nx
for k=1:Nz Hy(i,k)=Hy(i,k)+(dt/mu0)*((Ez(i+1,k)-Ez(i,k))+(Ex(i,k)-Ex(i,k+1)))/dz;
end
end
PEx=Ex;
for i=1:Nx
for k=2:Nz
Ex(i,k)= Ex(i,k)+(dt/eps0)*(Hy(i,k-1)-Hy(i,k))/dz;
end
end
for i=2:Nx
for k=1:Nz Ez(i,k)=Ez(i,k)+((dt/eps0)*(Hy(i,k)-Hy(i-1,k))/dz);
end
end
Ez(1,:)=0;Ez(Nx+1,:)=0;
Ex(1:Nx/3,Nz/2+1)=0;Ex(2*Nx/3:Nx,Nz/2+1)=0;
for m=1:Nx
Ex(m,Nz+1)=PEx(m,Nz)+(c*dt-dz)/(c*dt+dz)*(Ex(m,Nz)-Ex(m,Nz+1));
end
if n<=Nt1
Ex1(:,1)=exp(-(n*dt-t0)^2/T^2); % ¼Ó¼¤Àø
else
for m=1:Nx Ex1(m,1)=PEx1(m,2)+(c*dt-dz)/(c*dt+dz)*(Ex1(m,2)-Ex1(m,1));
end
end
for i=1:Nx
for k=1:Nz Hy1(i,k)=Hy1(i,k)+(dt/mu0)*((Ez1(i+1,k)-Ez1(i,k))+(Ex1(i,k)-Ex1(i,k+1)))/dz;
end
end
PEx1=Ex1;
for i=1:Nx
for k=2:Nz
Ex1(i,k)= Ex1(i,k)+(dt/eps0)*(Hy1(i,k-1)-Hy1(i,k))/dz;
end
end
for i=2:Nx
for k=1:Nz Ez1(i,k)=Ez1(i,k)+((dt/eps0)*(Hy1(i,k)-Hy1(i-1,k))/dz);
end
end
Ez1(1,:)=0;Ez1(Nx+1,:)=0;
for m=1:Nx Ex1(m,Nz+1)=PEx1(m,Nz)+(c*dt-dz)/(c*dt+dz)*(Ex1(m,Nz)-Ex1(m,Nz+1));
end
Vt1(n)=d*sum(Ex(:,Nz/4));
Vt2(n)=d*sum(Ex(:,3*Nz/4));
Vt3(n)=d*sum(Ex1(:,Nz/4));
Vt4(n)=d*sum(Ex1(:,3*Nz/4));
figure(1)
if mod(n,5)==0
mesh(Ex);
pause(0.001)
end
end
Nf=Nt;df=fmax/Nf;
for k=1:Nf V1f(k)=dt*sum(Vt1(1:Nt).*exp(-j*2*pi*(k*fmax/Nf)*(1:Nt)*dt)); Vin(k)=dt*sum(Vt3(1:Nt).*exp(-j*2*pi*(k*fmax/Nf)*(1:Nt)*dt)); V2f(k)=dt*sum(Vt2(1:Nt).*exp(-j*2*pi*(k*fmax/Nf)*(1:Nt)*dt)); V2n(k)=dt*sum(Vt4(1:Nt).*exp(-j*2*pi*(k*fmax/Nf)*(1:Nt)*dt));
end
S11=(V1f-Vin)./Vin;
S21=V2f./Vin;
figure(2)
plot((1:Nt)*dt,Vt1)
XLabel('Time (s)','FontSize',15,'FontWeight','b');
YLabel('Vt1 (V)','FontSize',15,'FontWeight','b');
axis([0 0.28e-7 -3.5 3.5])
figure(3)
plot((1:Nt)*dt,Vt2)
XLabel('Time (s)','FontSize',15,'FontWeight','b');
YLabel('Vt2 (V)','FontSize',15,'FontWeight','b');
axis([0 0.28e-7 -3.5 3.5])
figure(4)
plot((1:Nf)*df,abs(S11))
XLabel('f (Hz)','FontSize',15,'FontWeight','b');
YLabel('S²ÎÊý','FontSize',15,'FontWeight','b')
hold on
plot((1:Nf)*df,abs(S21),'r')
plot((1:Nf)*df,(abs(S11)).^2+(abs(S21)).^2,'g')
hold off

附件3：测试题源代码
clear
clc
c=3e8; % ×ÔÓÉ¿Õ¼ä¹âËÙ
mu0=4*pi*1e-7; % ×ÔÓÉ¿Õ¼ä´Åµ¼ÂÊ
eps0=8.85e-12; % ×ÔÓÉ¿Õ¼ä½éµç³£Êý
d=0.18;l=6;T=0.5e-9;
t0=3*T;fmax=1e9;
bc=0.3;dz=bc/20;
dx=d/18; Nx=18;
Nz=l/dz;dt=dz/(2*c);
Nt1=6*T/dt+100;
Nt=6*T/dt+8*800+100;
Ex=zeros(Nx,Nz+1);
Ez=zeros(Nx+1,Nz);
Hy=zeros(Nx,Nz);
Ex1=zeros(Nx,Nz+1);
Ez1=zeros(Nx+1,Nz);
Hy1=zeros(Nx,Nz);
Vt1=zeros(1,Nt);
Vt2=zeros(1,Nt);
PEx=zeros(Nx,Nz+1);
PEx1=zeros(Nx,Nz+1);
for n=1:Nt
if n<=Nt1
Ex(:,1)=exp(-(n*dt-t0)^2/T^2); % ¼Ó¼¤Àø
else
for m=1:Nx Ex(m,1)=PEx(m,2)+(c*dt-dz)/(c*dt+dz)*(Ex(m,2)-Ex(m,1));
end
end
for i=1:Nx
for k=1:Nz Hy(i,k)=Hy(i,k)+(dt/mu0)*((Ez(i+1,k)-Ez(i,k))+(Ex(i,k)-Ex(i,k+1)))/dz;
end
end
PEx=Ex;
for i=1:Nx
for k=2:Nz
Ex(i,k)= Ex(i,k)+(dt/eps0)*(Hy(i,k-1)-Hy(i,k))/dz;
end
end
for i=2:Nx
for k=1:Nz Ez(i,k)=Ez(i,k)+((dt/eps0)*(Hy(i,k)-Hy(i-1,k))/dx);
end
end
Ez(1,:)=0;Ez(Nx+1,:)=0;
Ex(9:10,Nz/2+1)=0; Ex(7:12,209)=0;
Ex(7:9,204)=0;Ex(10:12,204)=0;
Ez(9,201:204)=0;Ez(10,201:204)=0;
Ez(7,204:209)=0;Ez(12,204:209)=0;
if n<=Nt1
Ex1(:,1)=exp(-(n*dt-t0)^2/T^2); % ¼Ó¼¤Àø
else
for m=1:Nx Ex1(m,1)=PEx1(m,2)+(c*dt-dz)/(c*dt+dz)*(Ex1(m,2)-Ex1(m,1));
end
end
for i=1:Nx
for k=1:Nz Hy1(i,k)=Hy1(i,k)+(dt/mu0)*((Ez1(i+1,k)-Ez1(i,k))+(Ex1(i,k)-Ex1(i,k+1)))/dz;
end
end
PEx1=Ex1;
for i=1:Nx
for k=2:Nz
Ex1(i,k)= Ex1(i,k)+(dt/eps0)*(Hy1(i,k-1)-Hy1(i,k))/dz;
end
end
for i=2:Nx
for k=1:Nz Ez1(i,k)=Ez1(i,k)+((dt/eps0)*(Hy1(i,k)-Hy1(i-1,k))/dz);
end
end
Ez1(1,:)=0;Ez1(Nx+1,:)=0;
for m=1:Nx Ex1(m,Nz+1)=PEx1(m,Nz)+(c*dt-dz)/(c*dt+dz)*(Ex1(m,Nz)-Ex1(m,Nz+1));
end
Vt1(n)=d*sum(Ex(:,Nz/4));
Vt2(n)=d*sum(Ex(:,3*Nz/4));
Vt3(n)=d*sum(Ex1(:,Nz/4));
Vt4(n)=d*sum(Ex1(:,3*Nz/4));
Vtr(n)=Vt1(n)-Vt3(n);
figure(1)
if mod(n,30)==0
mesh(Ex);
pause(0.001)
end
end
Nf=Nt;df=fmax/Nf;
for k=1:Nf V1f(k)=dt*sum(Vtr(1:Nt).*exp(-j*2*pi*(k*fmax/Nf)*(1:Nt)*dt)); Vin(k)=dt*sum(Vt3(1:Nt).*exp(-j*2*pi*(k*fmax/Nf)*(1:Nt)*dt)); V2f(k)=dt*sum(Vt2(1:Nt).*exp(-j*2*pi*(k*fmax/Nf)*(1:Nt)*dt)); V2n(k)=dt*sum(Vt4(1:Nt).*exp(-j*2*pi*(k*fmax/Nf)*(1:Nt)*dt));
end
figure(2)
plot((1:Nt)*dt,Vt1)
XLabel('Time (s)','FontSize',15,'FontWeight','b');
YLabel('Vt1 (V)','FontSize',15,'FontWeight','b');
axis([0 1.7e-7 -3.5 3.5])
figure(3)
plot((1:Nt)*dt,Vt2)
XLabel('Time (s)','FontSize',15,'FontWeight','b');
YLabel('Vt2 (V)','FontSize',15,'FontWeight','b');
axis([0 1.7e-7 -3.5 3.5])
S11=V1f./Vin;
S21=V2f./Vin;
figure(4)
plot((1:Nf)*df,abs(S11))
XLabel('S²ÎÊý','FontSize',15,'FontWeight','b');
YLabel('f (Hz)','FontSize',15,'FontWeight','b');
axis([0 1e9 0.95 1.02])