梁边界元程序资料

6.12.1 梁静态拉压边界元程序BEM1 % * 数据准备 * % = 梁的几何参数 = L=1;D = 0.05;AREA = pi*D*D/4;% = 材料特性参数 = E=200E+9;% = 梁的载荷条件 = NForce=1;Force(1:NForce)=;%力的大小，参数：力的序号F_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/(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_Position(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 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_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),.-)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.2 梁纵向振动特性边界元程序BEM2% 梁的纵向振动特性边界元法计算% 输入梁的材料参数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(输入梁长度(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% 边界条件输入%MM=input(输入约束的数目( 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=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+2,1,1) str=num2str(freq(kk+1), Hz; text(freq(kk+1),aa(kk+1),str); end elseif MM=1 if (RR(kk)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)=;%力的大小，参数：力的序号F_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/(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阵为未知边界所对应的系数阵;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; for 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),.-)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(输入梁结构材料剪切模量(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=); L=input(输入梁长度（mm）:L=); L=L*1.E-3;end% 边界条件输入%MM=input(输入约束的数目( 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=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+2,1,1) str=num2str(freq(kk+1), Hz; text(freq(kk+1),aa(kk+1),str); end elseif MM=1 if (RR(kk)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;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; %将梁等分为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, -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_Position(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=;for 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=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=边界值： 力变形%Ureshape(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) 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(输入梁长度(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(输入约束的数目(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 MM0 QRQP=inv(QR)*QP; end % % % % % % % % % % % % for k1=1:kkk for k2=1:kkk if kkk1 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 kkk1 RRR(k,:,:)=R-RRP*QRQP; elseif kkk=1 RRR(k)=R-RRP*QRQP; else RRR(k,:,:)=r; endend%for kk=1:N if kkk1 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 kkk1 | kkk=0 if (RRR(kk,1,1)RRR(kk+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+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(输入要绘制的振型曲线的数目: kjj=);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*SHL4/(Lmd42*COCH4); r(4,1)=r(1,4); r(4,2)=r(2,4); r(4,3)=r(3,4); r(4,4)=(SL4*CHL4+CL4*SHL4)/(Lmd4*COCH4); % % % % % % % % % % % for k1=1:MM for k2=1:MM if MM1 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