作业1
一个二维正方形（边长a=10mm)的静电场区域，电位边界条件如图所示（单位：V)，求区域内的电位分布。要求用超松弛迭代法求解差分方程组进行计算。
clc
clear
close all
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)
surf(v2)
colormap winter
axis([0,11,0,11,0,100])
subplot(1,2,2)
contour(v2,20)
作业2
如图微带线，𝑤∕ℎ=1，ε𝑟 = 9.5。试用有限差分法求其有效介电常数εe及特阻抗Z0。（作为对照：采用保角变换法 解析结果为εe = 6.5，Z0 = 49.5）
clc
clear
close all
hx = 11; % 设置x方向网格节点数 间隔 .5 w
hy = 11; % 设置y方向网格节点数
v1 = zeros(hy, hx); % 设置二维数组，赋初值
v1(hy,:) = zeros(1, hx); % y=a 边界条件
v1(1, :) = zeros(1, hx); % y=0边界条件
v1(:, hx) = zeros(hy, 1); % x=hx边界条件
p = 2; q = 3; e = 1;%9.5; % 相对介电常数
for j = 1:p
v1(q, j) = 1; % 微带线边界条件
end
maxt = 1; t = 0; % 设置误差和最大误差参量
v2 = v1; n = 0;
while maxt>2e-6
n = n+1; % 迭代次数
maxt = 0;
for i = 2:q-1 % 对称轴
v2(i, 1) = (2*v1(i, 2)+v1(i+1, 1)+v2(i-1, 1))/4;
end
for i = q+1:hy-1
v2(i, 1) = (2*v1(i, 2)+v1(i+1, 1)+v2(i-1, 1))/4;
end
for i = 2:q-1 % 均匀区域
for j = 2:hx-1
v2(i, j) = (v1(i, j+1)+v1(i+1, j)+v2(i, j-1)+v2(i-1, j))/4;
end
end
for j = q+1:hy-1
v2(q, j) = (v1(q, j+1)+2/(1+e)*v1(q+1, j)+v2(q, j-1)...
+2/(1+e)*e*v2(q-1, j))/4;
end
for i = q+1:hy-1 % 均匀区域
for j = 2:hx-1
v2(i, j) = (v1(i, j+1)+v1(i+1, j)+v2(i, j-1)+v2(i-1, j))/4;
end
end
for i = 1:hy
for j = 1:hx
t = abs(v2(i, j)-v1(i, j)); % 收敛精度判断
if t>maxt
maxt = t;
end
end
end
v1 = v2;
end
surfl(v1)
colormap jet
作业3
模拟真空中二维TM电磁波的传播，边界设置为一阶Mur吸收边界，观察电磁波的传播过程。波源为正弦函数：
clc
clear
close all
xmesh=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=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);
ez1=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,:)=ez1(2,:)+coef3*(ez(2,:)-ez1(1,:));
ez(xmesh+1,:)=ez1(xmesh,:)+coef3*(ez(xmesh,:)-ez1(xmesh+1,:));
ez(:,1)=ez1(:,2)+coef3*(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.01)
ez1=ez;
end
作业4
基于Pocklington方程用MoM分析半波对称振子天线：
观察天线线径和分段数目分别取不同值对天线阻抗和辐射特性的影响 （半径分别取 0.001λ， 0.0001λ， 0.00001λ，分段数取11,21,31）
clc
clear
close all
% 初始化参数
c=3e8;
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, ':g' ,h,Ithe, 'r' , 'linewidth' ,2)
legend( 'pocklinton' , ' 解析值 ' )
grid on
xlabel( ' 电长度 ' )
ylabel( ' 归一化电流 ' )
% 方向图
theta=0:0.01:2*pi;
abs_f=zeros(1,length(theta));
for n=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)
polarplot(theta,Far_patten2, '-b' )
hold on
polarplot(theta_2,Far_theory, 'or' )
hold off
legend( 'pocklinton' , ' 解析值 ' )
title( ' 半波阵子天线 E 面方向图 ' )
figure(3)
polarplot(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)])
作业5
基于电场积分方程用MoM分析对称振子天线：
计算振子总长度分别为0.25λ ，0.5λ，λ，1.5λ时，振子的输入阻抗和E面方向图。
clc
clear
close all
lamda=1;
c = 3e8;
f = c/lamda; % 频率
w = 2*pi*f; % 角频率
e0 = 8.85e-12; % 介电常数
u0 = 4*pi*1e-7; % 磁导率
a = .001*lamda; % 半径
nn = 1.5;
L = nn*lamda; % 振子长度
k = 2*pi/lamda;
N = 11; % 分段数
dl = L/(N+1);
l = L/2-dl/2;
lz = -1: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;
fi_2 = exp(-k*2*dl*1i)/(2*dl);
for m=1:N
for n=1:N
if n==m
fi11 = fi;
fi12 = fi;
fi13 = fi;
fi2 = fi_1;
fi3 = fi_1;
elseif abs(n-m)==1
fi11 = fi_1;
fi12 = fi_1;
fi13 = fi_1;
if n>m
fi2 = fi_2;
fi3 = fi;
else
fi3 = fi_2;
fi2 = fi;
end
else
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*abs(m-n)*dl*1i)/(abs(m-n)*dl);
if n>m
fi2 = 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(n-m-1)*dl*1i)/(abs(n-m-1)*dl);
fi3 = exp(-k*abs(n-m+1)*dl*1i)/(abs(n-m+1)*dl);
end
end
z(m,n)=1j*w*u0/(4*pi)*dl*dl*fi11+(1/(1i*4*pi*w*e0))*(fi12-fi3-fi2+fi13);
end
end
% 电压矩阵求解
v = zeros(N, 1);
v((N+1)/2) = 1;
% 电流系数矩阵
I = z\v;
% 输入阻抗
zin = 1/I((N+1)/2);
disp(['输入阻抗 = ', num2str(zin)]);
theta = 0:.001:2*pi;
abs_f = zeros(1, length(theta));
for n = 1:1:N
abs_f = abs_f+I(n)*exp(k*(n*dl-L/2)*cos(theta)*1i);
end
abs_f = abs(sin(theta)*dl.*abs_f);
max_f = abs(sum(I)*dl);
far_patten2 = abs_f/max_f(1);
polarplot(theta,far_patten2, 'R--' )
5274
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
