《梁边界元程序资料》由会员分享，可在线阅读，更多相关《梁边界元程序资料（14页珍藏版）》请在金锄头文库上搜索。
1、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 |
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_Positio
3、n(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)
4、=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(
5、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 梁纵
6、向振动特性边界元程序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)：wid
7、th=); 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(
8、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)
