搜索
电磁散射边界元作业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为常数fori=1:TOTAL;ifiN+1%下侧x(i)=(i-1/2)*step;y(i)=0;elseif(iN&i2*N+1)%右侧x(i)=a;y(i)=(i-N-1/2)*step;elseif(i2*N&i3*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确定fors=1:TOTAL%场点循环fort=1:TOTAL%源点循环if(s==t)%奇异点处理h_st(s,t)=C;elsecurrent_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(t0&tN+1)|(t2*N&t3*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确定fors=1:TOTAL%场点循环fort=1:TOTAL%源点循环if(s==t)%奇异点处理k_st(s,t)=-(log(step/2)-1)*step/(2*pi);elsecurrent_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((t0&tN+1)|(t2*N&t3*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确定foru=1:TOTAL;if((u3*N)&(u4*N+1))%上侧g_u(u)=V_L;elseg_u(u)=0;end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%7.求解剩下电荷电压分布并显示charge_voltage=(H_K1)^(-1)*H_K2*g_u.';disp('下侧电位从左到右:');disp(charge_voltage(1:N));disp('上侧电位从左到右:')disp(charge_voltage(3*N:-1:2*N+1));disp('右侧电荷从上到下:');disp(charge_voltage(2*N:-1:N+1));disp('左侧电荷从上到下:');disp(charge_voltage(3*N+1:4*N));voltage=[charge_voltage(1:N);g_u(N+1:2*N)';charge_voltage(2*N+1:3*N);g_u(3*N+1:4*N)'];charge=[g_u(1:N)';charge_voltage(N+1:2*N);g_u(:,[2*N+1:3*N])';charge_voltage(3*N+1:4*N)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%8.方形内部测试点H1矩阵h_st1确定fort=1:TOTAL%源点循环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(t0&tN+1)|(t2*N&t3*N+1)%上下侧quad=abs(test_y-y(t))./...((test_x-current_x).^2+(y(t)-test_y).^2);h_st1(t)=-(1/(2*pi))*trapz(current_x,quad);else%左右侧quad=abs(test_x-x(t))./...((test_x-x(t)).^2+(current_y-test_y).^2);h_st1(t)=-(1/(2*pi))*trapz(current_y,quad);end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%9.方形内部测试点K1矩阵k_st1确定fort=1:TOTAL%源点循环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((t0&tN+1)|(t2*N&t3*N+1))%上下侧quad=log(((test_x-current_x).^2+...(y(t)-test_y).^2).^(1/2));k_st1(t)=-(1/(2*pi))*trapz(current_x,quad);else%左右侧quad=log(((test_x-x(t)).^2+...(current_y-test_y).^2).^(1/2));k_st1(t)=-(1/(2*pi))*trapz(current_y,quad);end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%10.求解内部测试点电位与解析解resolve=k_st1*charge-h_st1*voltage;%代入离散化泊松公式show=[test_x,test_y];disp('在方柱内部电位值');disp('test_x=test_y=');disp(show);disp('BEM方法为:');disp(resolve);analysis=V_L*(a-test_x)/a;disp('解析解为:')disp(analysis)t2=cputime;t=t2-t1下侧电位从左到右:252.2616150.000047.7384上侧电位从左到右:252.2616150.000047.7384右侧电荷从上到下:-52.9545-48.7611-52.9545左侧电荷从上到下:52.954448.761152.9544在方柱内部电位值test_x=test_y=33BEM方法为:149.9997解析解为:1502．一二维介质圆柱相对介电常数为er，入射波为TM波，电场与柱轴平行，求散射场。参考文献：《有限元法与边界元法》西电出版社曾余庚Page130-133参考资料分析了H，K矩阵元素的求法，其中对角元素非对角元素，其中为P（i）点到P（j）点的距离。MATLAB程序：clearall;closeall;%clc;epsr=input('inputtheer(相对介电常数)?');npoint=input('Howmanypointsalongthecircumference?');radius=input('inputtheradius(unit=wavelength)?');wavelength=1;fori=1:npointpc(i)=radius*sin((2*i-2)*360/(2*npoint)*pi/180)+sqrt(-1)*radius*cos((2*i-2)*360/(2*npoint)*pi/180);endp(npoint)=(pc(npoint)+pc(1))/2;fori=1:1:(npoint-1)p(i)=(pc(i)+pc(i+1))/2;%边界单元的中点坐标endfori=1:1:(npoint-1)dirp(i)=(pc(i+1)-pc(i))/abs(pc(i+1)-pc(i));enddirp(npoint)=(pc(1)-pc(npoint))/abs((pc(1)-pc(npoint)));%边界元素的单位方向矢量Li=abs(pc(1)-pc(2));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%求h1(i,j)，k1(i,j)fori=1:npointforj=1:npo