电磁场数值分析

《电磁场数值分析》

（作业）

--- 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)

➢结果：

作业2

模拟真空中二维TM 电磁波的传播，边界设置为一阶Mur 吸收边界，观察电磁波的传播过程。波源为正弦函数：

sin()sin(2)25

z t c E t n t ωπ

==∆

➢ 代码：

xmesh=150;

ymesh=150;

mu0=4*pi*(1.0e-7);

eps0=8.85e-12;

c=3.0e-8;

dx=1.0;

dt=0.7*dx/c;

timestep=200;

ez(1:xmesh+1,1:ymesh+1)=0.0;

hx(1:xmesh+1,1:ymesh)=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);

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);