电磁散射第二次作业

电磁散射边界元作业 10级电磁场专业

1．已知正方形柱的Ⅰ，Ⅲ边界的,ⅡⅣ边界的，求Ⅰ，Ⅲ边界的和

ⅡⅣ边界的。

参考文献：《边界元法基础》上海交大出版社王元淳 Page20-24

参考资料分析了H，K矩阵元素的求法，其中对角元素

为边界元素的长度。非对角元素，其中为P（i）点到P（j）点的距离，

为P（i）点到含P（j）点边界单元的垂直距离。求解出H，K矩阵后利用

求出未知边界条件

MATLAB程序：

% BEM.m

% 本程序用边界元方法求解正方形柱体内电位分布

clear;clc;

t1=cputime; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 1.常数定义

a=6; % 正方形长

N=3; % 每边分段数

step=a/N; % 每段长度

TOTAL=N*4; % 共剖分成TOTAL段

C=1/2; % 常数定义

NN=100; % 积分离散精度

V_L=300; % 已知电压矩阵

test_x=a/2; % 方形内部任意一点X坐标

test_y=a/2; % 方形内部任意一点Y坐标%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 2.坐标定位,计算各段中点对应的坐标

% 以方柱左下角为坐标原点建立坐标系

% 方柱左右两边X为常数,方柱上下两边Y为常数

for i=1:TOTAL;

if i

x(i)=(i-1/2)*step;

y(i)=0;

elseif (i>N & i<2*N+1) % 右侧

x(i)=a;

y(i)=(i-N-1/2)*step;

elseif (i>2*N & i<3*N+1) % 上侧

x(i)=a-(i-2*N-1/2)*step;

y(i)=a;

else % 左侧

x(i)=0;

y(i)=a-(i-3*N-1/2)*step;

end;

end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 3.H矩阵h_st确定

for s=1:TOTAL % 场点循环

for t=1:TOTAL % 源点循环

if (s==t) % 奇异点处理

h_st(s,t)=C;

else

current_x=linspace(x(t)-step/2,x(t)+step/2,NN); % X积分变量离散

current_y=linspace(y(t)-step/2,y(t)+step/2,NN); % Y积分变量离散

if (t>0 & t2*N & t<3*N+1) % 上下侧

quad=abs(y(s)-y(t))./...

((x(s)-current_x).^2+(y(t)-y(s)).^2);

h_st(s,t)=-(1/(2*pi))*trapz(current_x,quad);

else % 左右侧

quad=abs(x(s)-x(t))./...

((x(s)-x(t)).^2+(current_y-y(s)).^2);

h_st(s,t)=-(1/(2*pi))*trapz(current_y,quad);

end;

end;

end;

end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 4.K矩阵k_st确定

for s=1:TOTAL % 场点循环

for t=1:TOTAL % 源点循环

if (s==t) % 奇异点处理

k_st(s,t)=-(log(step/2)-1)*step/(2*pi);

else

current_x=linspace(x(t)-step/2,x(t)+step/2,NN); % X积分变量离散

current_y=linspace(y(t)-step/2,y(t)+step/2,NN); % Y积分变量离散

if ((t>0 & t2*N & t<3*N+1)) % 上下侧

quad=log(((x(s)-current_x).^2 + ...

(y(t)-y(s)).^2).^(1/2));

k_st(s,t)=-(1/(2*pi))*trapz(current_x,quad);

else % 左右侧

quad=log( ( (x(s)-x(t)).^2 + ...

(current_y-y(s)).^2 ).^(1/2) );

k_st(s,t)=-(1/(2*pi))*trapz(current_y,quad);

end;

end;

end;

end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 5.矩阵整序

% 上下侧电荷分布已知

% 左右侧电压分布已知

H_K1=[h_st(:,[1:N]),-k_st(:,[N+1:2*N]),h_st(:,[2*N+1:3*N]),-k_st(:,

[3*N+1:4*N])];

H_K2=[k_st(:,[1:N]),-h_st(:,[N+1:2*N]),k_st(:,[2*N+1:3*N]),-h_st(:,

[3*N+1:4*N])]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 6.已知部分电压和电荷矩阵g_u确定

for u=1:TOTAL;

if ( (u>3*N) & (u<4*N+1) ) % 上侧

g_u(u)=V_L;

else

g_u(u)=0;

end;

end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 7.求解剩下电荷电压分布并显示

charge_voltage=(H_K1)^(-1)*H_K2*g_u.';