平板传输线电磁仿真

1、电磁仿真综合实践Numerical Simulation of Electromagnetic Field作者 张泳 Author Yong Zhang指导老师 王秉中Instructor Bing-Zhong Wang单位 电子科大物理电子学院99777Class 99777 CPE of UESTC设计时间 2002年5月8日6月13日【摘要】仿真是通过建立系统模型，并且对模型进行实际研究，以代替对于实际系统的研究。Matlab 5.3是一种功能强大工程数学软件。本文阐明了怎样利用Matlab 5.3实现一维以及二维电磁仿真实例的建模分析、编程调试、结果分析。【Abstract】Simul

2、ation is an available means in studying real system by Model-building.This paper building a numerical simulation of electromagnetic field model based on the Matlab 5.3,the design is discussed in detail.【关键词】电磁仿真，计算电磁学，电磁场【Key Words】 Electromagnetic Simulation,Computational Electromagnetics,Electroma

3、gnetic Field【设计目的】1 课程性质：属于非计算机专业计算机基础教育第三层次计算机综合应用层次的课程2 任务： 21 加强学生的计算机综合应用能力、尤其是运用计算机分析和解决专业问题的能力培养，使学生对独立进行科学研究有初步实践。 22 初步掌握一种纯数值电磁仿真方法时域有限差分法。 23 初步学会综合应用一种程序设计语言进行科学与工程计算。 24 增强科技报告的写作能力，学会相关软件使用。正文第一部分：电磁仿真的意义现代科学研究的基本模式：“科学实验，理论分析，数值计算”三位一体。计算电磁学：以电磁场理论为基础，以计算机及计算机技术为工具和手段，运用计算科学所提供的各种方法，解决

4、复杂电磁问题的一门应用科学。电磁场工程电磁场理论 电磁规律 实验验证 数学方程计算电磁学 计算数据 建模与仿真 运算手段和方法 优化与设计【电场问题求解方法】解析法电磁场问题的传统解法。强调电磁分析和数学分析通常给出紧凑高效的计算程序。程序的最终拥护只具有很少的弹性，只能改变很少的参数，结构的主要特性已被编进程序，适用于专用程序开发。纯数值法直接以数值的、程序的形式描述电磁场问题。普适性强，用户弹性大。特定问题的边界条件、电气结构、激励等特性可以不编入基本程序，而由用户输入，更好的情况是通过图形界面输入。用户不必具备高深的电磁理论、数学及数值技术方面的知识。适用于通用商业程序开发。受到硬件限制

5、大。 FDTD TLM FEM 惠更斯原理 变分原理麦克斯韦方程数学描述 及 离散化处理结果计算计算程序数值模型 边界条件材料特性 后处理 用户界面 第二部分 时域有限差分法基本原理【差分原理】前项差分 误差 O(h)后项差分 误差 O(h)中心差分 误差 O(h*h) 用中心差分代替对空间、时间坐标的微分 【Yee的差分网格】每个磁场分量被四个电场分量环绕，每个电场分量被四个磁场分量环绕。E和H的空间取样位置相差半个空间步长，E和H时间取样位置相差半个时间步长。Yee的差分算法考虑空间一个无缘区域,其媒质的参数不随时间变化且各项同性,则麦克斯韦旋度方程在直角坐标系中,写成分量式:= (1)

6、(2) (3) (4) (5) (6)【蛙跳格式】H0 H1/2 H3/2E0 E1 E2【解的稳定条件】：； 【数值色散】 在时域有限差分网格中，数值波模的传播速度将随频率改变，即有色散。这种色散由数值网格引起，而非物理上客观存在。为了减小色散，实用中通常取空间步长满足 。第三部分 均匀平行板传输线传输特性仿真 FDTD基本过程【平行板传输线中的主模TEM模】根据平行板传输线中TEM波的特性对MAXWELL旋度方程简化得： 其中 若取时域波形为高斯波形， ，如下图所示。 【一维差分格式】E【网格划分】Ex Hy Ex Hy Ex Hy Ex Hy Ex k=0 k=1/2 k=1 k=3/2

7、 k=2 k=Nz-1/2 k=Nz T=0.5ns，则=1GHz，d=0.18m,L=6m, =0.3m, =0.015m, 400, 【时间步进过程】在初始（n=0）时令所有场两为零。将高斯脉冲激励加在输入口 由FDTD方程算出n+1/2时刻所有场点的磁场 由FDTD方程算出n+1时刻所有场点的电场 令金属边界上的切向电场分量为零，并在截断边界上使用截断边界条件 记录并输出两个观测面处的电压值 n n+1重复上一步骤，直到（n=nt）脉冲波全部传出计算区域。 【截断边界条件】 终端短路： 终端匹配：（吸收边界条件） 【一维仿真实例】终端匹配：入射波大约经过1000步后到达边界，由于采用吸收

8、边界条件，高斯波几乎被完全吸收，通过图可以观测到吸收后只有振幅小于的波形存在，这在误差的允许范围之内。第一条、第二条分别为第一、二观测面的电压。时间步长sZ方向网格数（0015m）Z方向网格数（0015m）电压强度V电场强度V/M电场强度V/M第二条第一条 入射波 被吸收后的波形 观测面电压终端短路：入射波到达终端后由于截断边界条件采取的是短路，所以波形反转，向入射口传输，由于Nt取得足够大（Nt=2000），波形会到达入射口，所以应该在波形到达入射口之前，激励脉冲已完全进入入射口后，将入射口边界条件切换为吸收边界条件。第一条第二条分别为第一、二观测面的电压。在横轴下方的是反射电压，上方为入射

9、电压。第二条电场强度V/M时间步长sZ方向网格数（0015m）Z方向网格数（0015m）第一条第一条电压强度V电场强度V/M 被反射的波形 被吸收后的波形 观测面电压第四部分 带挡板的平行板传输线传输特性仿真X L/2 d s=d/3Z L【二维FDTD方程】由于不连续性的出现，沿X方向场的均匀性遭到破坏，沿Y方向场的均匀性仍然保持，此时平行板传输线中的场分量有Ex,Ez,Hy.【参数设置】：数组EX(NX,NZ+1),EZ(NX+1,NZ),HY(NX,NZ),L/2应选足够长，使得在挡板反射波返回入射口之前，激励脉冲已完全进入入射口，并且入射口边界条件已切换为吸收边界条件。【观测面电压波形

10、】【入射、反射、透射电压频谱】【网络的S参数】现在我们再记录下总电压、就可以分离出入射电压、反射、透射电压。=- = 终端短路的反射系数：终端匹配的S参数： 【终端条件】终端匹配终端短路=0【观测面电压波形】【入射、反射、透射电压频谱】【二维仿真实例】 终端匹配：可以通过对Ex波形的观测知道，波形到达挡板时产生尖端效应，高斯波有一部分被“挤”过挡板，有一部分被挡板反射。此时应该将入射口和终端都设为吸收边界条件。 如下图所示：Z方向网格数(0.015m)X方向网格数(0.015m)电场强度Ex V/MZ方向网格数(0.015m)X方向网格数(0.015m)电场强度Ex V/M X方向网格数(0.

11、015m) 入射波 到达挡板 电场强度Ex V/MZ方向网格数(0.015m)Z方向网格数(0.015m)X方向网格数(0.015m)X方向网格数(0.015m)电场强度Ex V/M 通过挡板 两端都吸收电压幅度(V)时间步长（s）第一个观测面反射电压第二个观测面透射电压第一个观测面入射电压 参考面电压 S参数频率f (0.1GHz)终端匹配时S参数终端短路：由于挡板与短路板之间构成一种谐振腔结构，时域脉冲波在其间来回反射，入射口仍然是吸收边界条件，波完全消失的时间会比终端匹配情况要长。Nt足够大（Nt=20000）如下图所示：电场强度ExV/MX方向网格数(0.015m)z方向网格数(0.0

12、15m) 终端短路、入射口吸收，Ex被终端反射的波形电压幅度V 频率f (0.1GHz) 终端短路时第一参考面电压 反射系数频率f (0.1GHz) 终端短路时反射系数 第五部分 课程总结与讨论短短几周的电磁仿真综合实践很快就结束了，这几周紧张的学习给我留下了深刻的印象。我们已经学习了电磁场与电磁波、微波技术基础等课程，对电磁理论有一定的认识和了解，但那只是完全理论性质的。通过这门课程的学习，我对电磁场等相关知识有了进一步的了解，有了一定感性的认识，感觉电磁场也不是那么和天书一样的难以理解。通过用计算机模拟将电磁现象展示在我们面前，使得我不但加深了对电磁理论的理解，同时也锻炼了我们计算机应用方

13、面的能力，对科学研究有了一定的了解。我有如下几点体会：首先，这次实践使得我对计算机的应用有了进一步的提高。以往我对计算机的认识只是处于一个很简单的层次上，对其应用、原理也只是略知皮毛。但通过这段时间的学习，自己用计算机以及相关软件完成了这次实践活动。虽然在开始的时候对软件以及原理理解不够，但是通过老师的多次讲解，同学之间的互相探讨，渐渐地自己掌握了一定的方法。看着自己编出的程序调试出令人满意的结果，还是相当喜悦的。其次，这次的实践活动让我加深了对专业知识的了解，对电磁理论也不象以前那样见到就躲了。当我用自己编的程序捕捉到了它们，对它们的产生、传播、反射、透射等现象有了更为直接的认识和了解。自己

14、本来对电磁理论有一定的兴趣，苦于没有很好的方法帮助自己更好的理解和学习，或许通过这次实践活动，自己以后能够找到适合自己的学习电磁理论的方法。第三，通过对原理、程序的不断学习、操作，使得我认识到做任何事情都要认真，持之以恒。王老师科学严谨的作风给我留下了深刻的印象，不论是做人还是做科研，王老师都是那么的一丝不苟。这一点让我也不敢有所懈怠。同时我也感觉到自己能力的欠缺，分析问题、解决问题的能力还有待提高。相信在以后我会通过自己的不断努力，提高自己的能力，学习和工作上更加细致认真。第四，这次实践活动，我感觉是我上大学以来收获最大的一门课程。王老师讲课不是传统的灌输式的教育，学生是真正意义上学习的主体

15、。课前需要我们对理论方面的知识多加复习和学习，上课的时候需要我们积极的开动脑筋，如何把理论转化为实际可见，课后也需要我们不断的思考，如何才能使得我们做的软件更加趋于完善。在掌握计算电磁学有关理论知识的同时，更加重要的是我们学到了方法，这一点我觉得才是最为重要的。下面说一下我对这门课程的一些建议吧。我觉得要消除一些学生对这门课程的紧张情绪（开始我就有这样的感觉），适当的增加课时数，毕竟我们目前涉及的真的只是皮毛而已。增加课时数，可以让大家对理论的理解加深，同时也会有更多的同学以后投入到电磁场领域方面的学习和研究。还有上机时间最好不要靠得太近了，毕竟学习有一个消化吸收的过程。这次的科研实践活动是令

16、人愉快的。它让我体会到了快乐，当然也有抓耳挠腮的痛苦；它让我学习到了更多新的知识，新的方法；它让我看到了自己的不足之处，在以后的学习过程中，我一定会加以改正的。最后我还要特别感谢王老师，以及研究生师姐在这次实践活动中给予我的帮助。正是他们一丝不苟的工作态度和严谨的作风深深的教育了我，以后我也会尽我最大的努力去做好每一件事情。第六部分 附录一、参考文献1、计算电磁学（节选），王秉中，电子科技大学，2000年版2、电磁仿真综合实践，王秉中，电子科技大学3、MATLAB基础与应用简明教程，张平 ，北京航空航天大学出版社，2001年1月第一版二、相关程序1、一维终端匹配%定义变量t=0.5E-9;c=

17、3E8;fmax=1/(2*t);lammin=c/fmax;epslon=8.854E-12;mu=4E-7*p;y0=120*pi;dz=lammin/20;dt=dz/(2*c);t0=3*t;l=6;d=0.18;nz=400;nt=1000;k1=101;k2=301;cm=dt/(mu*dz);ce=dt/(epslon*dz);cmur=(dt*c-dz)/(dt*c+dz);ex(1:nz+1)=0;hy(1:nz)=0;vref1(1:nt)=0;vref2(1:nt)=0;%程序for n=1:nt ex(1)=exp(-(n*dt-t0)2/t2); for i=1:nz

18、 hy(i)=hy(i)+cm*(ex(i)-ex(i+1); end exm=ex(nz); for j=2:nz ex(j)=ex(j)+ce*(hy(j-1)-hy(j); end ex(nz+1)=exm+cmur*(ex(nz)-ex(nz+1);%一维匹配终端吸收条件 vref1(n)=ex(k1)*d; vref2(n)=ex(k2)*d; %plot(1:nz+1,ex);pause(0.001);%电场强度监控 x=1:nt; plot(x,vref1(x),'r',x,vref2(x),'g');pause(0.001);%观测面电压图 %g

19、ridend%一维终端匹配2、一维终端短路%定义变量y0=120*pi;dz=lammin/20;dt=dz/(2*c);t0=3*t;l=6;d=0.18;nz=400;nt=2000;k1=101;k2=301;cm=dt/(mu*dz);ce=dt/(epslon*dz);cmur=(dt*c-dz)/(dt*c+dz);ex(1:nz+1)=0;hy(1:nz)=0;vref1(1:nt)=0;vref2(1:nt)=0;for n=1:nt for i=1:nz hy(i)=hy(i)+cm*(ex(i)-ex(i+1); end if n<=800 ex(1)=exp(-(n

20、*dt-t0)2/t2); for j=2:nz ex(j)=ex(j)+ce*(hy(j-1)-hy(j); endelse exm=ex(2); for j=2:nz ex(j)=ex(j)+ce*(hy(j-1)-hy(j); end ex(nz+1)=0; %终端短路边界条件 ex(1)=exm+cmur*(ex(2)-ex(1); %高斯波完全进入后，入口改为吸收边界条件end vref1(n)=ex(k1)*d; vref2(n)=ex(k2)*d; x=1:nt;plot(x,vref1(x),'r',x,vref2(x),'g');%两个观测面的

21、电压波形 pause(0.01);% plot(1:nz+1,ex);pause(0.001);%实时监控Ex波形 %gridend%一维终端短路3、二维终端匹配%定义变量t=0.5E-9;c=3E8;fmax=1/(2*t);lammin=c/fmax;epslon=8.854E-12;mu=4E-7*pi;y0=120*pi;dz=lammin/20;dx=lammin/20;dt=dz/(2*c);t0=3*t;l=6;d=0.18;nz=400;nt=1000;k1=101;k2=301;cm=dt/(mu*dz);ce=dt/(epslon*dz);cf=dt/(epslon*dx)

22、;cmur=(dt*c-dz)/(dt*c+dz);nx=12;ez(1:nx+1,1:nz)=0;ex(1:nx,1:nz+1)=0;hy(1:nx,1:nz)=0;ver1in(1:nt)=0;ver1re(1:nt)=0;ver2tr(1:nt)=0;exb(1:nx,1)=0;exe(1:nx,1)=0;nf=10;df=1e8;jj=sqrt(-1);gver1in(1:nf)=0;gver1re(1:nf)=0;gver2tr(1:nf)=0;s11(1:nf)=0;s21(1:nf)=0;for n=1:nt %时间循环 for i=1:nx %二维FDTD过程 for k=1:

23、nz hy(i,k)=hy(i,k)+cm*(ez(i+1,k)-ez(i,k)*dz/dx+ex(i,k)-ex(i,k+1); end end if n<=300 ex(1:12,1)=exp(-(n*dt-t0)2/t2); for i=1:nx for k=2:nz ex(i,k)=ex(i,k)+ce*(hy(i,k-1)-hy(i,k); end end ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0; %设置挡板 for i=2:nx for k=1:nz ez(i,k)=ez(i,k)+cf*(hy

24、(i,k)-hy(i-1,k); end end ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0; else %高斯波完全进入，入射端改为吸收边界条件 exb(1:nx,1)=ex(1:nx,2);exe(1:nx,1)=ex(1:nx,nz); for i=1:nx for k=2:nz ex(i,k)=ex(i,k)+ce*(hy(i,k-1)-hy(i,k); end end for i=2:nx for k=1:nz ez(i,k)=ez(i,k)+cf*(hy(i,k)-hy(i-1,k); end endex

25、(1:nx,1)=exb(1:nx,1)+cmur*(ex(1:nx,2)-ex(1:nx,1);%入口边界条件ex(1:nx,nz+1)=exe(1:nx,1)+cmur*(ex(1:nx,nz)-ex(1:nx,nz+1);%吸收边界条件 ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0;endfor i=1:nx %分离出入射电压、反射电压、透射电压 ver2tr(n)=ver2tr(n)+ex(i,301)*d/12; end if n<=500 for i=1:nx ver1in(n)=ver1in(n)+

26、ex(i,101)*d/12; end else for i=1:nx ver1re(n)=ver1re(n)+ex(i,101)*d/12; end end %x=1:nt%plot(x,ver1in(x),'r',x,ver1re(x),'b',x,ver2tr(x),'m');%各个电压波形surf(ex); %监控Ex波形pause(0.01)endfor f=1:nf for n=1:nt %傅氏变换 gver1in(f)=gver1in(f)+dt*ver1in(n)*exp(-2*pi*f*df*n*dt*jj); gver1re(

27、f)=gver1re(f)+dt*ver1re(n)*exp(-2*pi*f*df*n*dt*jj); gver2tr(f)=gver2tr(f)+dt*ver2tr(n)*exp(-2*pi*f*df*n*dt*jj); endendfor f=1:nf %S参数 s11(f)=gver1re(f)/gver1in(f); s21(f)=gver2tr(f)/gver1in(f);endfor i=1:nf s1(i)=abs(s11(i); s2(i)=abs(s21(i);end%x=1:nf;%plot(x,s1,'g',x,s2,'r');%S参数%g

28、rid %二维终端匹配4、二维终端短路%定义变量t=0.5E-9;c=3E8;fmax=1/(2*t);lammin=c/fmax;epslon=8.854E-12;mu=4E-7*pi; y0=120*pi;dz=lammin/20;dx=lammin/20;dt=dz/(2*c); t0=3*t;l=6;d=0.18;nz=400; nt=20000;k1=101;k2=301; cm=dt/(mu*dz);ce=dt/(epslon*dz); cf=dt/(epslon*dx);cmur=(dt*c-dz)/(dt*c+dz);nx=12; ez(1:nx+1,1:nz)=0;ex(1:

29、nx,1:nz+1)=0;hy(1:nx,1:nz)=0; vref1(1:nt)=0;vref2(1:nt)=0;ver1in(1:nt)=0; ver1re(1:nt)=0;ver2tr(1:nt)=0;exb(1:nx,1)=0;exe(1:nx,1)=0; nf=10;df=1e8;jj=sqrt(-1); gver1in(1:nf)=0;gver1re(1:nf)=0;gver2tr(1:nf)=0; s11(1:nf)=0;s21(1:nf)=0;for n=1:nt for i=1:nx %二维FDTD过程 for k=1:nz hy(i,k)=hy(i,k)+cm*(ez(i+

30、1,k)-ez(i,k)*dz/dx+ex(i,k)-ex(i,k+1); end end if n<=300 ex(1:12,1)=exp(-(n*dt-t0)2/t2); for i=1:nx for k=2:nz ex(i,k)=ex(i,k)+ce*(hy(i,k-1)-hy(i,k); end end ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0; %设置挡板 for i=2:nx for k=1:nz ez(i,k)=ez(i,k)+cf*(hy(i,k)-hy(i-1,k); end end ez(1,1:nz)=0; ez(nx+1,1:nz)=0; ex(1:4,201)=0; ex(9:12,201)=0; else %高斯波完全进入，入射端改为吸收边界条件 ex