电磁散射边界元

电磁散射边界元作业1、边界元的研究背景与意义(BEM)是随着有限元计算方法的发展而在七十年代被提出来的求解偏微分方程初边值问题的数值计算方法。由于它一般只需进行边界剖分，有降低一维进行数值计算的特性，而且所占计算机内存小、计算时间省、与有限元法相棍合能很好地解决工程实际问题等优点，目前这一方法已在各个工程领域获得广泛应用。特别是边界元法中的基本解的性质往往能恰当地反应出无穷远处的边界条件，使得用它来研究无限域问题具有其他数值方法所不能具有的一些独特的优点。2、边界元的研究方法边界元法是在经典的边界积分方程法的基础上吸取了有限元离散技术而发展起来的一种偏微分的数值解法。它把区域内微分方程的边值问题归化为边界上的积分方程，然后在边界上离散求解。其基础在于边界归化，如果边界归化的途径不同，可以从同一边值问题得到几种不同形式的边界积分方程，因而导致不同的边界元方法。目前边界归化一般有以下几种方法:采用Green公式归化的直接边界元法;采用位势理论归化的间接边界元法;还有用Green函数代替基本解的自然边界归化，称为自然边界元法。边界元法也有一定的局限性。首先边界元方法必须事先明确偏微分方程的基本解，这使得它往往局限于求解常系数或某些特殊变系数线性偏微分方程。对于一般的偏微分方程，如果能够采用一些特殊的技术把它们的边值问题或初边值问题化为积分方程，这些积分方程往往是奇异积分方程，甚至是通常意义下不可积的奇异性积分方程，以至于只能在广义函数的意义上去理解这些积分。这一方面在技术上带来了数值积分的困难，另一方面也使对它的数学分析复杂化。3、实例计算1．已知正方形柱的Ⅰ，Ⅲ边界的n,ⅡⅣ边界的，求Ⅰ，Ⅲ边界的和ⅡⅣ边界的n。参考文献：《边界元法基础》上海交大出版社王元淳Page20-24参考资料分析了H，K矩阵元素的求法，其中对角元素20.5(1)2iiiiiiLHGLnLiL为边界元素的长度。非对角元素1214ijiijijhLHdr，111()4iijijLGLndr其中ijr为P（i）点到P（j）点的距离，ijh为P（i）点到含P（j）点边界单元的垂直距离。求解出H，K矩阵后利用1122123412343344HKHKKHKH求出未知边界条件1234MATLAB程序：%BEM.m%本程序用边界元方法求解正方形柱体内电位分布%*********************************初始化***********************************clc;clear;clearall;closeall;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(anal