梁边界元程序资料

1、6.12.1 梁静态拉压边界元程序BEM1 % * 数据准备 * % = 梁的几何参数 = L=1;D = 0.05;AREA = pi*D*D/4;% = 材料特性参数 = E=200E+9;% = 梁的载荷条件 = NForce=1;Force(1:NForce)=100000;%力的大小，参数：力的序号F_Position(1:NForce)=0.5;%力的位置。参数：力的序号% = 梁的边界条件 = IB(1:4)=0;U(1:4)=0;% 边界编号% _% | POSITION | 受力（N） | 变形(m) |% -% | X=0(首端） | 1 | 3 |% | X=L(末端)

2、| 2 | 4 |% -U(2)=0;U(3)=0;IB(2)=1;IB(3)=1;% = 欲求解的内点坐标 = N=1000; %将梁等分为N等份I=1:N+1;X(I)=L/N*(I-1);% 内点中包含两端边界点% * 数据准备结束 *%*生成力系统矩阵 G 阵*G(1:2,1:2) = 0;G=0 -L/(2*AREA*E);-L/(2*AREA*E),0;%* 生成位移系统矩阵 H 阵 *H(1:2,1:2) = 0;H=-0.5 0.5; 0.5 -0.5;%* 生成力矩阵 P 阵 *P(1:2,1)=0;for I = 1:NForceP(1)=P(1)+Force(I)*F_P

3、osition(I)/(2*AREA*E);P(2)=P(2)+Force(I)*(L-F_Position(I)/(2*AREA*E);end% *转换成FK*X+P的形式*H = -H;K = inv(G)*H;P = inv(G)*P;% * 将已知与末知边界所对应的G，H阵分解成两部分，其中AB阵为未知边界所对应的系数阵;*A=eye(2) -K;AB=;for j=1:4 if IB(j)=0 AB=AB A(:,j); else for i=1:2 P(i)=P(i)-A(i,j)*U(j); end endendResult=inv(AB)*P;ID=0;for j=1:4 if

4、 IB(j)=1 ID=ID+1; U(j)=Result(ID); endend% * 计算内点值 *UI=; %保存内点值for I=1:N+1% X(I)=L/N*I; AA=-X(I)/(2*AREA*E) (X(I)-L)/(2*AREA*E) 0.5 0.5;-0.5/AREA 0.5/AREA 0 0; % PI(1)=0; PI(2)=0; for J=1:NForce PI(1) = PI(1) + (Force(J)*(abs(X(I)-F_Position(J)/(-2*AREA*E); PI(2) = PI(2) + (Force(J)*FUHAO(X(I),F_Pos

5、ition(J)/2); end UI = UI;X(I) (AA*U'+PI')'end% * 结果输出 *message='边界值： 力变形'reshape(U,2,2)'message='内点值：位置坐标 变形应力'UIsubplot(2,1,1),plot(UI(:,1),UI(:,2),'.-')title('图1 挠度') xlabel('长度L / m') ylabel('变形U/ m') grid onsubplot(2,1,2),plot(UI(:,

6、1),UI(:,3),'*-'); title('图2 挠角') xlabel('长度L / m') ylabel('挠角 / rad') grid onreturn;6.12.2 梁纵向振动特性边界元程序BEM2% 梁的纵向振动特性边界元法计算% 输入梁的材料参数E=input('输入梁结构材料弹性模量(GPa)：E=');p=input('输入梁结构材料密度(g/cm3)：p=');E=1.E9*E;p=p*1000.; % 输入梁的几何结构参数type=input(' 输入梁的截面形状

7、(圆形=1，矩形=2),type=');if type=1 D=input('输入梁截面直径(mm)：D='); L=input('输入梁长度(mm)：L='); A=0.25*pi*D*D*1.E-6; L=L*1.E-3;elseif type=2 width=input('输入梁截面的宽度(mm)：width='); high=input('输入梁截面的高度（mm）:high='); L=input('输入梁长度（mm）:L='); A=width*high*1.E-6; L=L*1.E-3;end%

8、 边界条件输入%MM=input('输入约束的数目(< 2)：MM=');for k=1:MM M(k)=input('输入约束的编号(1,2)：=');end% 输入频率扫描计算参数%N=input('输入频率扫描计算点数目，N=');W0=input('输入起始频率(Hz)，W0=');WW=input('输入频率间隔（Hz），WW=');%for k=1:N freq(k)=W0+WW*(k-1); Lmd=(2*pi*freq(k)*(p/E)0.5; SL=sin(Lmd*L); CL=cos(Lm

9、d*L); AES=Lmd*A*E*SL; r(k,1,1)=-CL/AES; r(k,1,2)=-1/AES; r(k,2,1)=-1/AES; r(k,2,2)=-CL/AES; % % % % % % % % % % if MM > 0 if M(1)=1 RR(k)=-r(k,2,1)*r(k,1,2)/r(k,1,1)+r(k,2,2); elseif M(1)=2 RR(k)=r(k,1,1)+r(k,1,2)*(-r(k,2,1)/r(k,2,2); end endend%for kk=1:N if MM=1 xx=abs(RR(kk); elseif MM=0 xx=ab

10、s(r(kk,1,1); end aa(kk)=log(xx);end figure;plot(freq,aa);xlabel('频率（Hz）');ylabel('柔度（dB）');title('梁纵向振动频率特性曲线');for kk=1:N-2 if MM=0 if (r(kk,1,1)<r(kk+1,1,1)&(r(kk+1,1,1)>r(kk+2,1,1) str=num2str(freq(kk+1),' Hz' text(freq(kk+1),aa(kk+1),str); end elseif MM=

11、1 if (RR(kk)<RR(kk+1)&(RR(kk+1)>RR(kk+2) str=num2str(freq(kk+1),' Hz' text(freq(kk+1),aa(kk+1),str); end endend6.12.3 梁的扭转边界元程序BEM3% * 数据准备 * % = 梁的几何参数 = L=1;D = 0.05;%直径IP=pi*D4/32.;AREA = pi*D*D/4;% = 材料特性参数 = GGP=80e+9;% = 梁的载荷条件 = NForce=1;Force(1:NForce)=100000;%力的大小，参数：力的序号F

12、_Position(1:NForce)=0.5;%力的位置。参数：力的序号% = 梁的边界条件 = IB(1:4)=0;U(1:4)=0;% 边界编号% _% | POSITION | 受力（N） | 变形(m) |% -% | X=0(首端） | 1 | 3 |% | X=L(末端) | 2 | 4 |% -U(2)=0;U(3)=0;IB(2)=1;IB(3)=1;% = 欲求解的内点坐标 = N=1000; %将梁等分为N等份I=1:N+1;X(I)=L/N*(I-1);% 内点中包含两端边界点% * 数据准备结束 *%*生成力系统矩阵 G 阵*G(1:2,1:2) = 0;G=0 -L

13、/(2*GGP*IP);-L/(2*GGP*IP),0;%*生成位移系统矩阵 H 阵 *H(1:2,1:2) = 0;H=-0.5 0.5; 0.5 -0.5;%* 生成力矩阵 P 阵 *P(1:2,1)=0;for I = 1:NForceP(1)=P(1)+Force(I)*F_Position(I)/(2*GGP*IP);P(2)=P(2)+Force(I)*(L-F_Position(I)/(2*GGP*IP);end% *转换成FK*X+P的形式*H = -H;K = inv(G)*H;P = inv(G)*P;% *将已知与末知边界所对应的G，H阵分解成两部分，其中AB阵为未知边界

14、所对应的系数阵;A=eye(2) -K;AB=;for j=1:4 if IB(j)=0 AB=AB A(:,j); else for i=1:2 P(i)=P(i)-A(i,j)*U(j); end endendResult=inv(AB)*P;ID=0;for j=1:4 if IB(j)=1 ID=ID+1; U(j)=Result(ID); endend% * 计算内点值 *UI=; %保存内点值for I=1:N+1 AA=-X(I)/(2*GGP*IP) (X(I)-L)/(2*GGP*IP) 0.5 0.5;-0.5 0.5 0 0; % PI(1)=0; PI(2)=0; fo

15、r J=1:NForce PI(1) = PI(1) + (Force(J)*(abs(X(I)-F_Position(J)/(-2*GGP*IP); PI(2) = PI(2) + (Force(J)*FUHAO(X(I),F_Position(J)/2); end UI = UI;X(I) (AA*U'+PI')'end% * 结果输出 *message='边界值： 力变形'reshape(U,2,2)'message='内点值：位置坐标 变形内力'UIsubplot(2,1,1),plot(UI(:,1),UI(:,2),&

16、#39;.-')title('图1 挠度') xlabel('长度L / m') ylabel('变形U/ m') grid onsubplot(2,1,2),plot(UI(:,1),UI(:,3),'*-'); title('图2 挠角') xlabel('长度L / m') ylabel('挠角 / rad') grid onreturn;6.12.4 梁扭转振动边界元程序BEM4% 梁的扭转振动特性边界元法计算% 输入梁的材料参数%G=input('输入梁结构

17、材料剪切模量(GPa)：G=');p=input('输入梁结构材料密度(g/cm3)：p=');G=1.E9*G;p=p*1000.; % 输入梁的几何结构参数%type=input(' 输入梁的截面形状(圆形=1，其它形状=2),type=');if type=1 D=input('输入梁截面直径(mm)：D='); L=input('输入梁长度(mm)：L='); Ip=pi*D4/32*1.E-12; L=L*1.E-3;elseif type=2 Ip=input('输入梁截面的极惯性矩(mm4)：Ip=&

18、#39;); L=input('输入梁长度（mm）:L='); L=L*1.E-3;end% 边界条件输入%MM=input('输入约束的数目(< 2)：MM=');for k=1:MM M(k)=input('输入约束的编号(1,2)：=');end% 输入频率扫描计算参数N=input('输入频率扫描计算点数目，N=');W0=input('输入起始频率(Hz)，W0=');WW=input('输入频率间隔（Hz），WW=');%for k=1:N freq(k)=W0+WW*(k-1);

19、 Lmd=(2*pi*freq(k)*(p/G)0.5; SL=sin(Lmd*L); CL=cos(Lmd*L); AES=Lmd*Ip*G*SL; r(k,1,1)=-CL/AES; r(k,1,2)=-1/AES; r(k,2,1)=-1/AES; r(k,2,2)=-CL/AES; % % % % % % % % % % if MM > 0 if M(1)=1 RR(k)=-r(k,2,1)*r(k,1,2)/r(k,1,1)+r(k,2,2); elseif M(1)=2 RR(k)=r(k,1,1)+r(k,1,2)*(-r(k,2,1)/r(k,2,2); end ende

20、nd%for kk=1:N if MM=1 xx=abs(RR(kk); elseif MM=0 xx=abs(r(kk,1,1); end aa(kk)=log(xx);end figure;plot(freq,aa);xlabel('频率（Hz）');ylabel('柔度（dB）');title('梁扭转振动频率特性曲线');for kk=1:N-2 if MM=0 if (r(kk,1,1)<r(kk+1,1,1)&(r(kk+1,1,1)>r(kk+2,1,1) str=num2str(freq(kk+1),'

21、 Hz' text(freq(kk+1),aa(kk+1),str); end elseif MM=1 if (RR(kk)<RR(kk+1)&(RR(kk+1)>RR(kk+2) str=num2str(freq(kk+1),' Hz' text(freq(kk+1),aa(kk+1),str); end endend6.12.5 梁的弯曲边界元程序BEM5% * 数据准备 * % = 梁的几何参数 = L = 0.3;I=pi*0.054/64.% = 材料特性参数 = E=200E+9;EI=E*I;% = 梁的载荷条件 = NForce=1;

22、Force(1:NForce)=1000;%力的大小，参数：力的序号F_Position(1:NForce)=0.2;%力的位置。参数：力的序号% = 梁的边界条件 = IB(1:8)=0;U(1:8)=0;% 边界编号% _% | POSITION | 力（N）|弯矩(Nm)| 挠度(m)|挠角(rad)|% -% | X=0(首端） | 1 | 2 | 5 | 6 |% | X=L(末端) | 3 | 4 | 7 | 8 |% -U(7)=0;U(4)=0;U(5)=0;U(6)=0;IB(5)=1;IB(6)=1;IB(7)=1;IB(4)=1;% = 欲求解的内点坐标 = N=10;

23、%将梁等分为N等份I=1:N+1;X(I)=L/N*(I-1);% 内点中包含两端边界点% * 数据准备结束 * % *生成力系统矩阵 G 阵*G(1:4,1:4) = 0;G=-L3/(6*EI), 0, 0, -L2/(4*EI); 0, L2/(4*EI), -L3/(6*EI), 0; 0, L/(2*EI), -L2/(4*EI), 0; L2/(4*EI), 0, 0, L/(2*EI);%* 生成位移系统矩阵 H 阵 *H(1:4,1:4) = 0;H=0.5, -L/2 -0.5, 0; -0.5, 0, 0.5 L/2; 0 ,-0.5, 0 ,0.5;0, 0.5, 0,

24、-0.5;%* 生成力矩阵 P 阵 *P(1:4,1)=0;for I = 1:NForceP(1)=P(1)+Force(I)*(2*L3+F_Position(I)3-3*L*F_Position(I)2)/(12*EI);P(2)=P(2)+Force(I)*(2*L3+(L-F_Position(I)3 -3*L*(L-F_Position(I)2) /(12*EI);P(3)=P(3)+Force(I)*( F_Position(I)*(-F_Position(I)+2*L)/(4*EI);P(4)=P(4)+Force(I)*(L-F_Position(I)*(-F_Positio

25、n(I)-L)/(4*EI);end% *转换成FK*X+P的形式*H = -H;K = inv(G)*H;P = inv(G)*P;% *将已知与末知边界所对应的G，H阵分解成两部分，其中AB阵为未知边界所对应的系数阵A=eye(4) -K;AB=;for j=1:8 if IB(j)=0 AB=AB A(:,j); else for i=1:4 P(i)=P(i)-A(i,j)*U(j); end endendResult=inv(AB)*P;ID=0;for j=1:8 if IB(j)=1 ID=ID+1; U(j)=Result(ID); endend% * 计算内点值 *UI=;f

26、or I=1:N+1AA=-(2*L3+X(I)3-3*L*X(I)2)/(12*EI),-(X(I)*(X(I)-2*L)/(4*EI),-(2*L3+(L-X(I)3 -3*L*(L-X(I)2)/(12*EI), (L-X(I)*(L-X(I)-2*L)/(4*EI), -0.5, (X(I)-L)/2., -0.5, X(I)/2;(X(I)*(X(I)-2*L)/(4*EI), (X(I)-L)/(2*EI), -(L-X(I)*(L-X(I)-2*L) /(4*EI), -X(I)/(2*EI), 0, -0.5, 0, -0.5 ; PI(1)=0; PI(2)=0; for J

27、=1:NForce PI(1)=PI(1)+(Force(J)*(2*L3+(abs(F_Position(J)-X(I)3-3*L*(F_Position (J)-X(I)2)/(12*EI);PI(2)=PI(2) + (-Force(J)*(abs(F_Position(J)-X(I) *(abs(F_Position(J)-X(I)-2*L)/(-4*EI) *FUHAO(F_Position(J),X(I); end UI = UI;X(I) (PI'-AA*U')'end% * 结果输出 *message='边界值： 力变形'%Ureshap

28、e(U,4,2)'message='内点值：位置坐标 变形应力'UIsubplot(2,1,1),plot(UI(:,1),UI(:,2),'.-')title('图1 挠度') xlabel('长度L / m') ylabel('变形U/ m') grid onsubplot(2,1,2),plot(UI(:,1),UI(:,3),'*-'); title('图2 挠角') xlabel('长度L / m') ylabel('挠角 / rad'

29、;) grid onreturn;6.12.6 梁弯曲振动边界元程序BEM6% 梁的弯曲振动特性边界元法计算clear all% % 输入计算初始条件% 输入梁的材料参数%E=input('输入梁结构材料弹性模量(GPa)：E=');p=input('输入梁结构材料密度(g/cm3)：p=');E=1.E9*E;p=p*1000.; % 输入梁的几何结构参数%type=input(' 输入梁的截面形状(圆形=1，矩形=2),type=');if type=1 D=input('输入梁截面直径(mm)：D='); L=input(&

30、#39;输入梁长度(mm)：L='); A=0.25*pi*D*D*1.E-6; II=pi*D4/64.*1.E-12; L=L*1.E-3;elseif type=2 width=input('输入梁截面的宽度(mm)：width='); high=input('输入梁截面的高度（mm）:high='); L=input('输入梁长度（mm）:L='); A=width*high*1.E-6; II=width*high3/12*1.E-12; L=L*1.E-3;end% 边界条件输入%MM=input('输入约束的数目(&

31、lt; 4)：MM=');for k=1:MM M(k)=input('输入约束的编号(1,2,3,4)：=');endkkk=0; % 自由度已知元素的数目：kkk if MM=1 for k1=1:4 if M(1)=k1 kkk=kkk+1; MK(kkk)=k1; end endelseif MM=2 for k1=1:4 if M(1)=k1 & M(2)=k1 kkk=kkk+1; MK(kkk)=k1; end endelseif MM=3 for k1=1:4 if M(1)=k1 & M(2)=k1 & M(3)=k1 kkk=

32、kkk+1; MK(kkk)=k1; end endend % 输入频率扫描计算参数N=input('输入频率扫描计算点数目，N=');W0=input('输入起始频率(Hz)，W0=');WW=input('输入频率间隔（Hz），WW=');%for k=1:N freq(k)=W0+WW*(k-1); Lmd=(2*pi*freq(k)0.5*(A*p/E/II)0.25; SL=sin(Lmd*L); CHL=cosh(Lmd*L); CL=cos(Lmd*L); SHL=sinh(Lmd*L); COCH=E*II*(CL*CHL-1);

33、 r(1,1)=(SL*CHL-CL*SHL)/(Lmd3*COCH); r(1,2)=SL*SHL/(Lmd2*COCH); r(1,3)=(SL-SHL)/(Lmd3*COCH); r(1,4)=(CHL-CL)/(Lmd2*COCH); r(2,1)=r(1,2); r(2,2)=(SL*CHL+CL*SHL)/(Lmd*COCH); r(2,3)=(CL-CHL)/(Lmd2*COCH); r(2,4)=(SL+SHL)/(Lmd*COCH); r(3,1)=r(1,3); r(3,2)=r(2,3); r(3,3)=(SL*CHL-CL*SHL)/(Lmd3*COCH); r(3,4

34、)=-SL*SHL/(Lmd2*COCH); r(4,1)=r(1,4); r(4,2)=r(2,4); r(4,3)=r(3,4); r(4,4)=(SL*CHL+CL*SHL)/(Lmd*COCH); % % % % % % % % % % for k1=1:MM for k2=1:MM if MM>1 QR(k1,k2)=r(M(k1),M(k2); else QR(1)=r(M(k1),M(k2); end end % for k2=1:kkk QP(k1,k2)=r(M(k1),MK(k2); end end if MM>0 QRQP=inv(QR)*QP; end %

35、% % % % % % % % % % % for k1=1:kkk for k2=1:kkk if kkk>1 R(k1,k2)=r(MK(k1),MK(k2); else R(1)=r(MK(k1),MK(k2); end end % for k2=1:MM RRP(k1,k2)=r(MK(k1),M(k2); end end % % % % % % % % % % if kkk>1 RRR(k,:,:)=R-RRP*QRQP; elseif kkk=1 RRR(k)=R-RRP*QRQP; else RRR(k,:,:)=r; endend%for kk=1:N if kkk

36、>1 xx=abs(RRR(kk,1,1); elseif kkk=1 xx=abs(RRR(kk); else xx=abs(RRR(kk,1,1); end aa(kk)=log(xx);end figure;plot(freq,aa);xlabel('频率（Hz）');ylabel('柔度（dB）');title(' 梁的弯曲振动频率特性曲线');for kk=1:N-2 if kkk>1 | kkk=0 if (RRR(kk,1,1)<RRR(kk+1,1,1)&(RRR(kk+1,1,1)>RRR(kk+

37、2,1,1) str=num2str(freq(kk+1),' Hz' text(freq(kk+1),aa(kk+1),str); end elseif kkk=1 if (RRR(kk)<RRR(kk+1)&(RRR(kk+1)>RRR(kk+2) str=num2str(freq(kk+1),' Hz' text(freq(kk+1),aa(kk+1),str); end endend%NN=input('输入梁结构内部等份数,NN=');dlt=L/NN;kjj=input('输入要绘制的振型曲线的数目: kj

38、j=');for kj=1:kjj fq=input('输入固有频率（Hz），fq='); Lmd4=(2*pi*fq)0.5*(A*p/E/II)0.25; SL4=sin(Lmd4*L); CHL4=cosh(Lmd4*L); CL4=cos(Lmd4*L); SHL4=sinh(Lmd4*L); COCH4=E*II*(CL4*CHL4-1); % r(1,1)=(SL4*CHL4-CL4*SHL4)/(Lmd43*COCH4); r(1,2)=SL4*SHL4/(Lmd42*COCH4); r(1,3)=(SL4-SHL4)/(Lmd43*COCH4); r(1,4)=(CHL4-CL4)/(Lmd42*COCH4); r(2,1)=r(1,2); r(2,2)=(SL4*CHL4+CL4*SHL4)/(Lmd4*COCH4); r(2,3)=(CL4-CHL4)/(Lmd42*COCH4); r(2,4)=(SL4+SHL4)/(Lmd4*COCH4); r(3,1)=r(1,3); r(3,2)=r(2,3); r(3,3)=(SL4*CHL4-CL4*SHL4)/(Lmd43*COCH4); r(3,4)=-SL4*SHL