电磁场数值分析

电磁场数值分析

《电磁场数值分析》

（作业）

--- 2016学年 ---

学院：

学号：

姓名：

联系方式：

任课教师：

2016年6月6日

作业1

一个二维正方形（边长a=10mm)的静电场区域，电位边界条件如图所示（单位：V)，求区域内的电位分布。要求用超松弛迭代法求解差分方程组进行计算。

代码：

hx=11;

hy=11;

v1=zeros(hy,hx);

v1(hy,:)=ones(1,hx)*100;

v1(1,:)=ones(1,hx)*50;

for i=1:hy;

v1(i,1)=0;

v1(i,hx)=100;

end

w=2/(1+sin(pi/(hx-1)));

maxt=1;

t=0;

v2=v1;

n=0;

while(maxt>1e-6)

n=n+1;

maxt=0;

for i=2:hy-1;

for j=2:hx-1;

v2(i,j)=(1-w)*v1(i,j)+w*(v1(i+1,j)+v1(i,j+1)+v2(i-1,j )+v2(i,j-1))/4;

t=abs(v2(i,j)-v1(i,j));

if (t>maxt)

maxt=t;

end

end

end

v1=v2;

end

subplot(1,2,1)

mesh(v2)

axis([0,11,0,11,0,100])

subplot(1,2,2)

contour(v2,20)

结果：

coef1=dt/(mu0*dx);

coef2=dt/(eps0*dx);

coef3=(c*dt-dx)/(c*dt+dx);

ezold=ez;

for now=1:timestep;

hx=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(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,:)=ezold(2,:)+coef3*(ez(2,:)-ezold(1,:));

ez(xmesh+1,:)=ezold(xmesh,:)+coef3*(ez(xmesh,:)-ezold (xmesh+1,:));

ez(:,1)=ezold(:,2)+coef3*(ez(:,2)-ezold(:,1));

ez(:,ymesh+1)=ezold(:,ymesh)+coef3*(ez(:,ymesh)-ezold (:,ymesh+1));

ez(xmesh/2+1,ymesh/2+1)=sin(now*dt*2*pi*c/25.0); mesh(ez)

pause(0.01)

ezold=ez;

end

结果：

作业3

基于Pocklington方程用MoM分析半波对称振子天线：

观察天线线径和分段数目分别取不同值对天线阻抗和辐射特性的影响（半径分别取0.001λ，0.0001λ，0.00001λ，分段数取11,21,31）

代码：

%%初始化参数

c=3e-8;

r=1;

f=c/r;

w=2*pi*f;

e0=8.85e-12;

u0=4*pi*1e-7;

a=0.0001*r;

L=0.5*r;

k=2*pi/r;

N=31;

dl=L/(N+1);

l=L/2-dl/2;

lz=-l:dl:1;

lzs=lz(1:N);

lzm=lz(1:N)+dl/2;

lze=lz(2:N+1);

%%阻抗矩阵元素求解

fi=log(dl/a)/(2*pi*dl)-k/(4*pi)*1j;

fi_1=exp(-k*dl*1j)/(4*pi*dl);

fi_2=exp(-k*2*dl*1j)/(8*pi*dl);

z=ones(N,N);

for m=1:N

for n=1:N

if m==n

fi1=fi;

fi2=fi_1;

fi3=fi_2;

z(m,n)=((k^2*dl^2-2)*fi1+fi2+fi3);

elseif abs(m-n)==1

fi1=fi_1;

fi2=fi;

fi3=fi_2;

z(m,n)=((k^2*dl^2-2)*fi+fi2+fi3);

else

fi1=exp(-k*abs(m-n)*dl*1j)/(4*pi*abs(m-n)*dl);

fi2=exp(-k*abs(m+1-n)*dl*1j)/(4*pi*abs(m+1-n)*dl);

fi3=exp(-k*abs(n+1-m)*dl*1j)/(4*pi*abs(n+1-m)*dl); z(m,n)=((k^2*dl^2-2)*fi+fi2+fi3);

end

end

end

%%电压矩阵求解

V=zeros(N,1);

V((N+1)/2)=-1*(1j*w*e0);

I=z\V;

Z_in=1/I((N+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('电长度')

ylabel('归一化电流')