电磁场数值分析.docx
《电磁场数值分析.docx》由会员分享，可在线阅读，更多相关《电磁场数值分析.docx（16页珍藏版）》请在冰点文库上搜索。
电磁场数值分析
《电磁场数值分析》
（作业）
---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;
fori=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;
fori=2:
hy-1;
forj=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吸收边界，观察电磁波的传播过程。
波源为正弦函数：
Ø代码：
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;
fornow=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）;
form=1:
N
forn=1:
N
ifm==n
fi1=fi;
fi2=fi_1;
fi3=fi_2;
z（m,n）=（（k^2*dl^2-2）*fi1+fi2+fi3）;
elseifabs（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','解析值'）
gridon
xlabel（'电长度'）
ylabel（'归一化电流'）
%%方向图
theta=0:
0.01:
2*pi;
abs_f=zeros（1,length（theta））;
forn=1:
1:
N
abs_f=abs_f+I（n）*exp（k*（n*dl-L/2）*cos（theta）*1j）;
end
abs_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））-cos（k*L/2））./sin（theta_2））;
figure
（2）
polar（theta,Far_patten2,'-b'）
holdon
polar（theta_2,Far_theory,'or'）
holdoff
legend（'pocklinton','解析值'）
title（'半波阵子天线E面方向图'）
figure（3）
polar（theta,ones（1,length（theta））,'-b'）
title（'半波阵子天线H面方向图'）
%%半波阵子增益
I_in=I（（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（'%.4fdB',Max_gain_dB）]）
Ø结果：
作业4
基于电场积分方程用MoM分析对称振子天线：
计算振子总长度分别为0.25λ，0.5λ，λ，1.5λ时，振子的输入阻抗和E面方向图。
Ø代码：
lamda=1;
a=0.0001;
me=8.85e-12;
mu=4*pi*（1e-7）;
arg=2*pi*（3e-8）/lamda;
L=0.2*lamda;
k=2*pi/lamda;
N=21;
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）;
form=1:
N
forn=1:
N
ifn==m
Fmnmm=（1/（2*pi*dL））*log（dL/a）-1j*k/（4*pi）;
Fmnee=（1/（2*pi*dL））*log（dL/a）-1j*k/（4*pi）;
Fmnss=（1/（2*pi*dL））*log（dL/a）-1j*k/（4*pi）;
Fmnse=exp（-1j*k*dL）/（4*pi*dL）;
Fmnes=exp（-1j*k*dL）/（4*pi*dL）;
elseifabs（n-m）==1
Fmnmm=exp（-1j*k*dL）/（4*pi*dL）;
Fmnee=exp（-1j*k*dL）/（4*pi*dL）;
Fmnss=exp（-1j*k*dL）/（4*pi*dL）;
ifn>m
Fmnse=exp（-1j*k*2*dL）/（4*pi*2*dL）;
Fmnes=exp（1/（2*pi*dL））*log（dL/a）-1j*k/（4*pi）;
else
Fmnes=exp（-1j*k*2*dL）/（4*pi*2*dL）;
Fmnse=exp（1/（2*pi*dL））*log（dL/a）-1j*k/（4*pi）;
end
else
num=abs（n-m）;
Fmnmm=exp（-1j*k*num*dL）/（4*pi*num*dL）;
Fmnee=exp（-1j*k*num*dL）/（4*pi*num*dL）;
Fmnss=exp（-1j*k*num*dL）/（4*pi*num*dL）;
ifn>m
Fmnse=exp（-1j*k*（num+1）*dL）/（4*pi*（num+1）*dL）;
Fmnes=exp（-1j*k*（num-1）*dL）/（4*pi*（num-1）*dL）;
else
Fmnes=exp（-1j*k*（num+1）*dL）/（4*pi*（num+1）*dL）;
Fmnse=exp（-1j*k*（num-1）*dL）/（4*pi*（num-1）*dL）;
end
end
z（m,n）=1j*arg*mu*dL*dL*Fmnmm+（1/（1j*arg*me））*（Fmnee-Fmnes-Fmnse-Fmnss）;
end
end
V=zeros（N,1）;
fedp=（N+1）/2;
V（fedp）=1;
I=linsolve（z,V）;
Z=V（fedp）/I（fedp）;
theta=0:
pi/100:
2*pi;
ftheta=0;
form=1:
length（theta）
F（m）=0;
forn=1:
N
F（m）=F（m）+I（n）*exp（1j*k*（n-fedp）*dL*cos（theta（m）））;
end
end
F=abs（F）;
F1=F/max（F）;
polar（theta,F1,'b.-'）
Ø结果：
作业5
对课程的建议、自己的收获等。
通过学习这门与数学和计算机编程联系紧密的课程，我发现自己数学理论知识还不是很扎实，而且自己对MATLAB的使用也不是很熟练，希望自己在今后能够加强对这方面知识的学习。
在学习的过程中也锻炼了自己对所学知识的应用，增加了自己的学习兴趣。
升级会员 