朱月恒 电磁场与电磁波实验.doc

电磁场与电磁波实验报告同实验者：冯泽涵 041010113 魏仁奎 041010119实验 二维静电场边值问题的数值仿真实验对于静电场边值问题，某些情况下用解析法求解很困难，甚至不可能。在实际求解过程中，可以将待求解电磁场问题化为离散的数值问题，即求解待求量在一些离散点上的数值。本实验根据二维静电场所满足的微分方程和边界条件，利用有限差分法来求解无源区域内一系列离散点上的电位，从而观测在给定边界条件下无源区域静电场的电位分布。1、 实验目的 1、学习用有限差分法求解静电场的方法。 2、观察边界条件对场的影响。2、 实验原理 1、方程的离散有限差分法对于解决任何偏微分方程都是一种可行的方法，因为所有的电磁场问题都可以表示成标量或矢量偏微分方程，所以可以利用有限差分法解决各种介 质中电磁场的空间分布。有限差分法把要解决的区域划分为有限个离散点，并将偏微分方程用一组差分方程代替。因此，这种方法是近似的。但是，如果我们把离散点取得足够密，就能够把误差减小到可以接受的程度。 2、边界条件的处理边界条件在电磁场问题的求解中起着至关重要的作用，因为不同电磁场问题可以满足相同的偏微分方程，它们之所以有不同的解就是因为满足的边界条件不同。为求得区域内各节点上近似解，还需将原问题中的边界条件化成离散的边界条件，即确定出在边界节点处的值。 3、差分方程的求解3、 实验内容 1、自行给定区域尺寸 a、b 和，取 x=a/ 20 ， y=b /10 ，在 y=b 的边界上令f ( x)=10 ，计算各内部节点处的电位值，并绘制出区域内的电位分布（可用Matlab 中的Mesh或Surf 函数）。 程序结果： 区域内电位分布： 2、其它条件不变，在y=b 的边界上令f(x )=10sin(x)，计算并绘制出区域内的电位分布。观察两种情况下电位分布的差别，理解边界条件对场分布的决定作用。程序结果：区域内电位分布： 3、对上述两种情况，分别以列表的形式给出 x= 10 x 和 y= 5 y 处各节点的电位与分离变量法（无穷级数截断到第100 项）所得结果的比较。分离变量法在y= 5 y的结果： 分离变量法在x= 10 x 的结果：程序代码：（其中a=20,b=10）a=20;B=10;pi=3.14159;dx=a/20;dy=B/10;M=a/dx-1;N=B/dx-1;m=M;x=1:a/20:a;%f=10*sin(pi*x/a);%f=x;A=zeros(N*M,N*M); % 事先把系数矩阵各元素置0，在重新赋值前保持为0b=zeros(N*M,1); % 事先把已知列向量各元素置0，在重新赋值前保持为0for j=1:Nfor i=1:Mk=(j-1)*M+i;if j=1 & j=N&i=1 & i=MA(k,k-1)=1.0/dx2;A(k,k)=-2.0*(1.0/dx2+1.0/dy2);A(k,k+1)=1.0/dx2;A(k,k-m)=1.0/dy2;A(k,k+m)=1.0/dy2;elseif i=1&j=1&j=NA(k,k)=-2.0*(1.0/dx2+1.0/dy2);A(k,k+1)=1.0/dx2;A(k,k-m)=1.0/dy2;A(k,k+m)=1.0/dy2;elseif i=M&j=1&j=NA(k,k-1)=1.0/dx2;A(k,k)=-2.0*(1.0/dx2+1.0/dy2);A(k,k-m)=1.0/dy2;A(k,k+m)=1.0/dy2;elseif j=1&i=1&i=MA(k,k-1)=1.0/dx2;A(k,k)=-2.0*(1.0/dx2+1.0/dy2);A(k,k+1)=1.0/dx2;A(k,k+m)=1.0/dy2;elseif j=N&i=M&i=1A(k,k-1)=1.0/dx2;A(k,k)=-2.0*(1.0/dx2+1.0/dy2);A(k,k+1)=1.0/dx2;A(k,k-m)=1.0/dy2;b(k)=(-1)*10/dy2; % f(i*dx)是上边界给定的已知电位elseif i=1&j=1A(k,k)=-2.0*(1.0/dx2+1.0/dy2);A(k,k+1)=1.0/dx2;A(k,k+m)=1.0/dy2;elseif i=1&j=NA(k,k)=-2.0*(1.0/dx2+1.0/dy2);A(k,k+1)=1.0/dx2;A(k,k-m)=1.0/dy2;b(k)=(-1)*10/dy2;elseif i=M&j=1A(k,k-1)=1.0/dx2;A(k,k)=-2.0*(1.0/dx2+1.0/dy2);A(k,k+m)=1.0/dy2;elseif i=M&j=NA(k,k-1)=1.0/dx2;A(k,k)=-2.0*(1.0/dx2+1.0/dy2);A(k,k-m)=1.0/dy2;b(k)=(-1)*10/dy2;endendendV=Ab;%surf(V);U=zeros(N,M);for i=1:M for j=1:N U(j,i)=V(j-1)*M+i); endend %subplot(1,2,1);mesh(U); %三维图xlabel(X);ylabel(Y);zlabel(U);axis(0,20,0,20,0,10);title(电位U的分布曲线);%subplot(1,2,2),contour(U,18); %等电位线% hold on;U1=0;for m=1:100U1=U1+40/pi*1/(m*(exp(m*pi/20*10)/2-exp(-1)*m*pi/20*10)/2)*sin(m*pi/20*10)*(exp(m*pi/20*5)/2-exp(-1)*m*pi/20*5)/2);end实验 矩形波导中场结构模拟实验一、实验目的要求： 1、通过实验编程及图像动态演示，形象具体的了解电磁波在波导中传播特性。 2、通过编写Matlab 程序，加深矩形波导中电磁波公式推导以及单模电磁波在矩形波导中的传播理解。 二、实验内容： 电磁场本身比较复杂和抽象，是涉及空间和时间的多维矢量场，需要具有较强的空间想象能力来理解它。 1、实验原理： 矩形波导是截面形状为矩形的金属波导管，如图一所示。波导内壁面位置坐标设为：x=0 和x=a ；y=0 和y=b 。波导中填充介电常数为 、磁导率为、电导率为 的媒质，通常波导内填充理想介质（ =0）。由于波导内没有自由电荷和传导电流，所以传播的电磁波是正弦电磁波。理想导电壁矩形波导中不可能传输TEM模，只能传输TE模或TM模。对于矩形波导中TEMN模的电场强度E、磁场强度H 场分量表达式为：其中： 为微波角频率；m和n 值可以取 0 或正整数，代表不同的TE波场结构模式，称为TE模，波导中可有无穷多个TE模式；kc为临界波束，为相位常数，波导中的一个重要参数为截止频率fc，有 当工作频率低于截止频率fc时，电磁场衰减很快，不可能传播很远，所以波导呈现高通滤波器的特性，只有工作频率高于截止频率fc时电磁波才能通过。具有最低截止频率的式，成为最低模式，也称为主模，其他模式都成为高次模式。在矩形波导内传输的所有模型中，TE10 模为主模。 2、实验步骤： 设置矩形波导宽边a=22.86mm，窄边b =10.16mm，波导内媒质为空气，当工作频率f 为9.84GHz 时，波导中只能传输TE10 模。 利用Matlab 显示矩形波导TE10模的电磁场分布的程序设计过程： （1）根据已知参数 m，n ，a ，b 和f 编程计算kc ， 和 角频率等参数。Matlab 中代码实现： a=22.86*1e-3; b=10.16*1e-3; f=9.84*1e9; m=1; n=0; miu=4*pi*1e-7; eps=8.854*1e-12; %E=2.71828; kc=(m*pi/a)2+(n*pi/b)2)0.5; w=2*pi*f; beta=(miu*eps*w2-kc2)0.5; （2）根据式 1-6 定义的各场强变量，以电场强度、磁场强度各分量为因变量，以时间t 为自变量。Matlab 中代码实现： ngrid=20; x=0:a/ngrid:a;y=0:b/2:b; z=0:0.04/ngrid:0.04;%定义 x ，y，z 坐标空间矩阵 % 公式表示 for p=0:ngrid%执行循环p 赋初值0 ，循环步长为 1 ，总步长 ngrid for q=0:2 for r=0:ngrid% 三层循环，赋值 ex 、ey 、ez 、hx 、hy、hz 空间上的数值 ex(p,q,r)=j*(w*miu/kc2)*(n*pi/b)*cos(m*pi/a)*x(p)*sin(n*pi/b)*y(q)*exp(j*(w*t)-(beta*z(r); ey(p,q,r)=j*(w*miu/kc2).*(m*pi/a).*sin(m*pi/a).*x(p).*cos(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); ez(p,q,r)=0; hx(p,q,r)=j*(beta/kc2).*(m*pi/a).*sin(m*pi/a).*x(p).*cos(n*pi/b)*y(q).*exp(j*(w*t)-(beta.*z(r); hy(p,q,r)=j*(beta/kc2).*(m*pi/a).*cos(m*pi/a).*x(p).*sin(n*pi/b)*y(q).*exp(j*(w*t)-(beta.*z(r); hz(p,q,r)=cos(m*pi/a).*x(p).*cos(n*pi/b ).*y(q).*exp(j* (w*t)-(beta.*z(r); end end end（3）通过绘图函数 quiver （）绘制二维场强分布图。 Matlab 中代码实现： 需要将三维的数组转换成二维形式，下面是参考代码，具体取值请自己编写 ex2(:,:)=ex(:,:,1); ey2(:,:)=ey(:,:,1); ez2(:,:)=ez(:,:,1); hx2(:,:)=hx(:,:,1); hy2(:,:)=hy(:,:,1); hz2(:,:)=hz(:,:,1); % 二维实现： X,Y=meshgrid(0:a/ngrid:a,0:b/2:b); % 生成 xoy 坐标面的二维网格个点矩阵 quiver(X,Y,abs(hx2),abs(hy2); %绘制磁场强度分布图，这里 abs 取矢量模值 % 也可以取幅值或者相位，分别为不加函数和angle（）函数实现 hold on; %图像保留，将电场强度绘制到同一张图上 quiver(X,Y,abs(ex),abs(ey); %绘制电场强度分布图 依次可以绘制出xoz、yoz 面的电磁场分布图。（4）通过绘图函数 quiver3（）绘制三维场强分布图。 Matlab中代码实现： % 三维实现 X,Y,Z=meshgrid(0:a/ngrid:a,0 :b/2:b,0:0.04/ngrid:0.04); quiver3(X,Y,Z,abs(hx),abs(hz),abs(hy);% 绘制三维磁场分布图 hold on; quiver3(X,Y,Z,abs(ex),abs(ez),abs(ey) ,r);% r 为标记红颜色 3、扩展练习： 扩展练习1 ，改变初始变量 m、n 、a 、b 、f 的值，编写多种模式在矩形波导中传输的程序，通过程序实现可观测到2 种及2 种以上 TE模式的电磁波在矩形波导中传输，具体编程以上程序给出了扩展。 扩展练习2 ，运用 for 循环对式1-6 中的时间 t 进行循环取值，用三维绘图绘制场分布图，利用 getframe（）函数捕获当前画面，产生一个数据向量，创建一个帧动画矩阵。程序每循环一次将绘制一张该时刻的电磁场分布图，当时间变量 t 大于预设值时跳出循环，再利用播放动画函数 movie（）将各帧图连续播放，形成波导内部三维空间电磁场的动态分布图。 Matlab 中代码提示： for N=1：20 F（N）=getframe； end movie（F）；% 将此段程序中的循环内加入三维绘图程序即可。 四、实验讨论及实验程序及数据：1、简要描述矩形波导中传播电磁波的原理。答：矩形波导中传播电磁波的原理: 电磁波入射到金属导体的表面时，透射系数为零，因此矩形波导将电磁波能量局限在矩形金属内部，电磁波在金属管内壁发生斜入射或反射，而沿管轴向形成导波。2、简单画出程序实验中绘制出的二维电场磁场分布图各一张（yoz 面的电场线，xoz面的磁场线）。 yoz面的电场线: xoz面的磁场线：附程序：a=22.86*1e-3;b=10.16*1e-3;f=9.84*1e9;m=1;n=0;miu=4*pi*1e-7;eps=8.854*1e-12;%E=2.71828;kc=(m*pi/a)2+(n*pi/b)2)0.5;w=2*pi*f;t=2;beta=(miu*eps*w2-kc2)0.5;ngrid=20;x=0:a/ngrid:a;y=0:b/2:b;z=0:0.04/ngrid:0.04;for p=1:21 for q=1:3 for r=1:21 ex(p,q,r)=j*(w*miu/kc2).*(n*pi/b).*cos(m*pi/a).*x(p).*sin(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); ey(p,q,r)=j.*(w*miu/kc2).*(m*pi/a).*sin(m*pi/a).*x(p).*cos(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); ez(p,q,r)=0; hx(p,q,r)=j.*(beta/kc2).*(m*pi/a).*sin(m*pi/a).*x(p).*cos(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); hy(p,q,r)=j.*(beta/kc2).*(n*pi/b).*cos(m*pi/a).*x(p).*sin(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); hz(p,q,r)=cos(m*pi/a).*x(p).*cos(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); end endendfor p=1:21 for q=1:3 k=ex(p,q,1);ex2(p,q)=k; k=ey(p,q,1);ey2(p,q)=k; k=ez(p,q,1);ez2(p,q)=k; k=hx(p,q,1);hx2(p,q)=k,k=hy(p,q,1);hy2(p,q)=k,k=hz(p,q,1);hz2(p,q)=k; endendX,Y=meshgrid(0:a/ngrid:a,0:b/2:b);quiver(X,Y,abs(hx),abs(hy);hold on;quiver(X,Y,abs(ex),abs(ey),r);3、简单画出程序实验中绘制出的三维电磁场分布图（用实线表示电场线，虚线表示磁场线）。 附程序：a=22.86*1e-3;b=10.16*1e-3;f=9.84*1e9;m=1;n=0;miu=4*pi*1e-7;eps=8.854*1e-12;%E=2.71828;kc=(m*pi/a)2+(n*pi/b)2)0.5;w=2*pi*f;t=2;beta=(miu*eps*w2-kc2)0.5;ngrid=20;x=0:a/ngrid:a;y=0:b/2:b;z=0:0.04/ngrid:0.04;for p=1:21 for q=1:3 for r=1:21 ex(p,q,r)=j*(w*miu/kc2).*(n*pi/b).*cos(m*pi/a).*x(p).*sin(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); ey(p,q,r)=j.*(w*miu/kc2).*(m*pi/a).*sin(m*pi/a).*x(p).*cos(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); ez(p,q,r)=0; hx(p,q,r)=j.*(beta/kc2).*(m*pi/a).*sin(m*pi/a).*x(p).*cos(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); hy(p,q,r)=j.*(beta/kc2).*(n*pi/b).*cos(m*pi/a).*x(p).*sin(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); hz(p,q,r)=cos(m*pi/a).*x(p).*cos(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); end endendfor p=1:21 for q=1:3 k=ex(p,q,1);ex2(p,q)=k; k=ey(p,q,1);ey2(p,q)=k; k=ez(p,q,1);ez2(p,q)=k; k=hx(p,q,1);hx2(p,q)=k,k=hy(p,q,1);hy2(p,q)=k,k=hz(p,q,1);hz2(p,q)=k; endendX,Y,Z=meshgrid(0:a/ngrid:a,0:b/2:b,0:0.04/ngrid:0.04);quiver3(X,Y,Z,abs(hz),abs(hx),abs(hy),:r);xlabel(传播方向);ylabel(a);zlabel(b);hold on;quiver3(X,Y,Z,abs(ez),abs(ex),abs(ey),-g);动画：附程序：a=22.86*1e-3;b=10.16*1e-3;f=9.84*1e9;m=1;n=0;miu=4*pi*1e-7;eps=8.854*1e-12;%E=2.71828;kc=(m*pi/a)2+(n*pi/b)2)0.5;w=2*pi*f;beta=(miu*eps*w2-kc2)0.5;ngrid=20;x=0:a/ngrid:a;y=0:b/2:b;z=0:0.04/ngrid:0.04;for t=1:20for p=1:21 for q=1:3 for r=1:21 ex(p,q,r)=j*(w*miu/kc2).*(n*pi/b).*cos(m*pi/a).*x(p).*sin(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); ey(p,q,r)=j.*(w*miu/kc2).*(m*pi/a).*sin(m*pi/a).*x(p).*cos(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); ez(p,q,r)=0; hx(p,q,r)=j.*(beta/kc2).*(m*pi/a).*sin(m*pi/a).*x(p).*cos(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); hy(p,q,r)=j.*(beta/kc2).*(n*pi/b).*cos(m*pi/a).*x(p).*sin(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); hz(p,q,r)=cos(m*pi/a).*x(p).*cos(n*pi/b).*y(q).*exp(j*(w*t)-(beta.*z(r); end endendfor p=1:21 for q=1:3 k=ex(p,q,1);ex2(p,q)=k; k=ey(p,q,1);ey2(p,q)=k; k=ez(p,q,1);ez2(p,q)=k; k=hx(p,q,1);hx2(p,q)=k,k=hy(p,q,1);hy2(p,q)=k,k=hz(p,q,1);hz2(p,q)=k; endendX,Y=meshgrid(0:a/ngrid:a,0:b/2:b);quiver(X,Y,abs(ex),abs(ez);X,Y,Z=meshgrid(0:a/ngrid:a,0:b/2:b,0:0.04/ngrid:0.04);quiver3(X,Y,Z,abs(hz),abs(hx),abs(hy);xlabel(传播方向);ylabel(a);zlabel(b);%hold on;%quiver3(X,Y,Z,abs(ez),abs(ex),abs(ey),r);F(t)=getframe;endmovie(F);4、通过实验，学会简单的 matlab 计算绘图等简单编程，根据电流场公式，试讨论编写电流场分布的二维三维绘图及动态仿真。二维：a=22.86*1e-3; b=10.16*1e-3; f=9.84*1e9; m=2; n=1;j=sqrt(-1);miu=4*pi*1e-7; eps=8.854*1e-12; %E=2.71828;kc=(m*pi/a)2+(n*pi/b)2)0.5;w=2*pi*f;beta=(miu*eps*w2-kc2)0.5;ngrid=20;x=0:a/ngrid:a;y=(0:b/2:b); z=(0:0.04/ngrid:0.04);%定义x，y，z坐标空间矩阵%公式表示ex=zeros(ngrid+1,3,ngrid+1);ey=zeros(ngrid+1,3,ngrid+1);ez=zeros(ngrid+1,3,ngrid+1);hx=zeros(ngrid+1,3,ngrid+1);hy=zeros(ngrid+1,3,ngrid+1);hz=zeros(ngrid+1,3,ngrid+1);%t=1;t=0;for p=1:ngrid+1%执行循环p赋初值0，循环步长为1，总步长ngridfor q=1:3for r=1:ngrid+1%三层循环，赋值ex、ey、ez、hx、hy、hz空间上的数值ex(p,q,r)=j*(w*miu/kc2)*(n*pi/b)*cos(m*pi/a)*x(p)*sin(n*pi/b)*y(q)*exp(j*(w*t)-(beta*z(r); ey(p,q,r)=j*(w*miu/kc2)*(m*pi/a)*sin(m*pi/a)*x(p)*cos(n*pi/b).*y(q)*exp(j*(w*t)-(beta*z(r);ez(p,q,r)=0;hx(p,q,r)=j*(beta/kc2)*(m*pi/a)*sin(m*pi/a)*x(p)*cos(n*pi/b)*y(q)*exp(j*(w*t)-(beta*z(r);hy(p,q,r)=j*(beta/kc2)*(m*pi/a)*cos(m*pi/a)*x(p)*sin(n*pi/b)*y(q)*exp(j*(w*t)-(beta*z(r);hz(p,q,r)=cos(m*pi/a)*x(p)*cos(n*pi/b)*y(q)*exp(j*(w*t)-(beta*z(r);endendendex2=zeros(3,ngrid+1);ey2=zeros(3,ngrid+1);ez2=zeros(3,ngrid+1);for q1=1:3for r1=1:ngrid+1ex2(q1,r1)=j*(w*miu/kc2)*(n*pi/b)*cos(m*pi/a)*x(1)*sin(n*pi/b)*y(q1)*exp(j*(w*t)-(beta*z(r1);ey2(q1,r1)=j*(w*miu/kc2)*(m*pi/a)*sin(m*pi/a)*x(1)*cos(n*pi/b).*y(q1)*exp(j*(w*t)-(beta*z(r1);ez2=0;endendhx2=zeros(ngrid+1,ngrid+1);hy2=zeros(ngrid+1,ngrid+1);hz2=zeros(ngrid+1,ngrid+1);for p2=1:ngrid+1 for r2=1:ngrid+1 hx2(p2,r2)=j*(beta/kc2)*(m*pi/a)*sin(m*pi/a)*x(p)*cos(n*pi/b)*y(1)*exp(j*(w*t)-(beta*z(r); hy2(p2,r2)=j*(beta/kc2)*(m*pi/a)*cos(m*pi/a)*x(p)*sin(n*pi/b)*y(1)*exp(j*(w*t)-(beta*z(r); hz2(p2,r2)=cos(m*pi/a)*x(p)*cos(n*pi/b)*y(1)*exp(j*(w*t)-(beta*z(r); endend %二维实现：u,v=meshgrid(1:b/2:b,1:a/ngrid:a);%生成xoy坐标面的二维网格个点矩阵subplot(3,1,1);quiver(u,v,abs(hx2),abs(hy2),y);%绘制磁场强度分布图，这里abs取矢量模值%也可以取幅值或者相位，分别为不加函数和angle（）函数实现title(磁场强度H在xoy平面分布图);xlabel(X);ylabel(Y);hold on;%图像保留，将电场强度绘制到同一张图上A,B=meshgrid(1:a/ngrid:a,1:b/2:b);subplot(3,1,2);quiver(A,B,abs(ex2),abs(ey2),g);title(电场强度E在yoz平面分布图);xlabel(Y);ylabel(Z);%三维实现ssX,Y,Z=meshgrid(0:a/ngrid:a,0:b/2