电磁场数值分析期末

1、精选优质文档-倾情为你奉上 电磁场数值分析（期末作业）- 2019学年 -学 院： 学 号： 姓 名： 联系方式： 任课教师： 2019年5月作业1模拟真空中二维TM电磁波的传播，边界设置为一阶Mur吸收边界，观察电磁波的传播过程。波源为正弦函数：Ø 代码：clcclearclose allxmesh =150; ymesh =150;mu0=4*pi*1.0E-7; eps0=8.85E-12;C= 3.0E8;dx=1.0; dt=0.7*dx/C; timestep=150;ez( 1:xmesh+1,1:ymesh+1 ) = 0.0;hx( 1:xmesh+1,1:ymes

2、h ) = 0.0;hy( 1:xmesh,1:ymesh+1 ) = 0.0;coef1 = dt/( mu0 * dx ); coef2 =dt/( eps0 * dx );coef3=(C*dt-dx)/(C*dt+dx);ez1=ez;for now = 1 : timestephx = hx - coef1 * ( ez( :, 2 : ymesh+1 ) - ez( :, 1 : ymesh ) );hy = hy + coef1 * ( ez(2 : xmesh+1, : ) - ez(1 : xmesh, : );ez( 2 : xmesh , 2 : ymesh ) = ez

3、( 2 : xmesh , 2 : ymesh ) - .coef2 * ( hx( 2 : xmesh, 2 : ymesh ) - hx( 2 : xmesh , 1 : ymesh - 1) ) + .coef2 * ( hy( 2 : xmesh ,2 : ymesh ) - hy( 1 : xmesh - 1,2 : ymesh) );ez(1,:)=ez1(2,:)+coef3*(ez(2,:)-ez1(1,:);ez(xmesh+1,:)=ez1(xmesh,:)+coef3*(ez(xmesh,:)-ez1(xmesh+1,:);ez(:,1)=ez1(:,2)+coef3*(

4、ez(:,2)-ez1(:,1);ez(:,ymesh+1)=ez1(:,ymesh)+coef3*(ez(:,ymesh)-ez1(:,ymesh+1);ez( xmesh/2+1, ymesh/2+1) = sin( now * dt * 2 * pi * C / 25.0 );mesh(ez); pause(0.05)ez1=ez;end Ø 结果与分析：第10时间步第100时间步第150时间步作业2基于Pocklington方程用MoM分析半波对称振子天线：观察天线线径和分段数目分别取不同值对天线阻抗和辐射特性的影响（半径分别取 0.001， 0.0001， 0.00001，

5、分段数取11,21,31，可列表说明）Ø 代码：clear all; close all; clc;% 初始化参数c=3e8; % 光速r=1 % 波长f=c/r; % 频率w=2*pi*f; % 角频率e0=8.85e-12; % 介电常数u0=4*pi*1e-7; % 磁导率a=0.00001*r; % 半径L=0.5*r; % 振子长度k=2*pi/r; % 波数N=11; % 分段数（奇数段）dl=L/(N+1); % 每段长度（分母中+1 为两头半段之和）l=L/2-dl/2; % 两头空出半段，满足电流为0的边界条件lz=-l:dl:l;lzs=lz(1:N); % 每一

6、小段的起点坐标lzm=lz(1:N)+dl/2; % 每一小段的中点坐标lze=lz(2:N+1); % 每一小段的终点坐标%阻抗矩阵元素求解fi=log(dl/a)/(2*pi*dl)-k/(4*pi)*1i;fi_1=exp(-k*dl*1i)/(4*pi*dl);fi_2=exp(-k*2*dl*1i)/(8*pi*dl);z=ones(N,N);for m=1:Nfor n=1:Nif m=n fi1=fi;fi2=fi_1;fi3=fi_1;z(m,n)=(k2*dl2-2)*fi1+fi2+fi3);elseif abs(m-n)=1fi1=fi_1;fi2=fi;fi3=fi_2

7、;z(m,n)=(k2*dl2-2)*fi1+fi2+fi3);elsefi1=exp(-k*abs(m-n)*dl*1i)/(4*pi*abs(m-n)*dl);fi2=exp(-k*abs(m+1-n)*dl*1i)/(4*pi*abs(m+1-n)*dl);fi3=exp(-k*abs(n+1-m)*dl*1i)/(4*pi*abs(n+1-m)*dl);z(m,n)=(k2*dl2-2)*fi1+fi2+fi3);end endend%电压矩阵求解V=zeros(N,1);V(N+1)/2)=-1*(1i*w*e0);% 计算电流系数矩阵I=zV;% 计算输入阻抗Z_in=1/I(N+

8、1)/2);disp('输入阻抗 = ',num2str(Z_in);% 计算振子上归一化电流分布I_amp=abs(I); Max=max(I_amp);Iunit2=0;I_amp/Max(1);0; % 两端零电流figure(1);h=0:dl/r:L/r;Ithe=sin(pi*h*r/L); % 半波振子电流解析值plot(h,Iunit2,'b',h,Ithe,'r','linewidth',2);legend('pocklinton','解析值');grid on;xlabel(&#

9、39;电长度L/lambda');ylabel('归一化电流');% 画方向图theta=0:0.01:2*pi;abs_f=zeros(1,length(theta);for n=1:1:Nabs_f=abs_f+I(n)*exp(k*(n*dl-L/2)*cos(theta)*1i);endabs_f=abs(sin(theta)*dl.*abs_f);Max_f=abs(sum(I)*dl);Far_patten2=abs_f/Max_f(1);theta_2=0:0.1:2*pi;Far_theory=abs(cos(k*(L/2)*cos(theta_2)-c

10、os(k*L/2)./sin(theta_2);figure(2);polar(theta,Far_patten2,'-b');hold on;polar(theta_2,Far_theory,'or');hold off;legend('pocklinton','解析值');title('半波振子天线E面方向图');figure(3);polar(theta,ones(1,length(theta),'-b');title('半波振子天线H面方向图');% 半波振子增益I_in=I

11、(N+1)/2);A=(w*u0)2/(4*pi*sqrt(u0/e0)*real(Z_in)*(abs(I_in)2);G_theta=A*abs_f.2;Max_gain=max(G_theta)Max_gain_dB=10*log10(Max_gain);disp('半波振子增益 = ',sprintf('%.4fdBi', Max_gain_dB);Ø 结果与分析： 分段数波长 1121310.001输入阻抗 =81.0029+36.1494iMax_gain =1.6446半波振子增益 = 2.1607dBi输入阻抗 =83.0507+39.

12、2545iMax_gain =1.6503半波振子增益 = 2.1757dBi输入阻抗= 83.9979+41.066iMax_gain =1.6520半波振子增益 = 2.1802dBi0.0001输入阻抗= 77.7374+38.0429iMax_gain =1.6420半波振子增益 = 2.1538dBi输入阻抗= 78.8657+39.7265iMax_gain =1.6469半波振子增益 = 2.1666dBi输入阻抗= 79.3079+40.9582iMax_gain =1.6482半波振子增益 = 2.1700dBi0.00001输入阻抗= 76.3172+39.4969iMax

13、_gain =1.6408半波振子增益 = 2.1506dBi输入阻抗= 77.1417+40.0554iMax_gain =1.6454半波振子增益 = 2.1626dBi输入阻抗= 77.44+40.936iMax_gain =1.6465半波振子增益 = 2.1657dBi作业3基于电场积分方程用MoM分析对称振子天线：计算振子总长度分别为0.25 ，0.5，1.5时，振子的输入阻抗和E面方向图。Ø 代码：clear all; close all; clc;% 初始化参数c=3e8; % 光速r=1; % 波长f=c/r; % 频率w=2*pi*f; % 角频率e0=8.85e

14、-12; % 介电常数u0=4*pi*1e-7; % 磁导率a=0.001*r; % 半径L=1.5*r; % 振子长度k=2*pi/r; % 波数N=11; % 分段数（奇数段）dl=L/(N+1); % 每段长度（分母中+1 为两头半段之和）l=L/2-dl/2; % 两头空出半段，满足电流为0的边界条件lz=-l:dl:l;lzs=lz(1:N); % 每一小段的起点坐标lzm=lz(1:N)+dl/2; % 每一小段的中点坐标lze=lz(2:N+1); % 每一小段的终点坐标%阻抗矩阵元素求解fi=2*log(dl/a)/dl-k*1i;fi_1=exp(-k*dl*1i)/dl;f

15、i_2=exp(-k*2*dl*1i)/(2*dl);for m=1:Nfor n=1:Nif m=n fi11=fi;fi12=fi;fi13=fi;fi2=fi_1;fi3=fi_1;elseif abs(m-n)=1 fi11=fi_1;fi12=fi_1;fi13=fi_1; if n>mfi2=fi_2;fi3=fi; else fi3=fi_2; fi2=fi; endelse fi11=exp(-k*abs(m-n)*dl*1i)/(abs(m-n)*dl); fi12=exp(-k*abs(m-n)*dl*1i)/(abs(m-n)*dl); fi13=exp(-k*ab

16、s(m-n)*dl*1i)/(abs(m-n)*dl); if n>mfi2=exp(-k*(abs(n-m+1)*dl*1i)/(abs(n-m+1)*dl);fi3=exp(-k*(abs(n-m-1)*dl*1i)/(abs(n-m-1)*dl); else fi2=exp(-k*(abs(m-n-1)*dl*1i)/(abs(m-n-1)*dl);fi3=exp(-k*(abs(m-n+1)*dl*1i)/(abs(m-n+1)*dl); end endz(m,n)=1i*w*u0/(4*pi)*dl*dl*fi11+(1/(1i*4*pi*w*e0)*(fi12-fi3-fi2

17、+fi13);endend%电压矩阵求解V=zeros(N,1);V(N+1)/2)=1;% 计算电流系数矩阵I=zV;% 计算输入阻抗Z_in=1/I(N+1)/2);disp('输入阻抗 = ',num2str(Z_in);% 画方向图theta=0:0.01:2*pi;abs_f=zeros(1,length(theta);for n=1:1:Nabs_f=abs_f+I(n)*exp(k*(n*dl-L/2)*cos(theta)*1i);endabs_f=abs(sin(theta)*dl.*abs_f);Max_f=abs(sum(I)*dl);Far_patten2=abs_f/Max_f(1);polar(theta,Far_patten2,'-b');title('振子天线E面方向图');Ø 结果与分析：0.25输入阻抗= = 12.76051-440.6278i