matlab计算应力位移,2012年-有限元作业-matlab编程实现有限元求解简单结构位移及应力.doc...

《2012年-有限元作业-matlab编程实现有限元求解简单结构位移及应力.doc》由会员分享，可在线阅读，更多相关《2012年-有限元作业-matlab编程实现有限元求解简单结构位移及应力.doc(8页珍藏版)》请在金锄头文库上搜索。

1、2012年有限元作业对于三角形单元有：1、形函数其中：相乘后得：也可写成：简写为：其中2、应变有弹性力学知道：，可得3、应力4、刚度矩阵题目：如图所示是一平面梁。载荷沿梁的上边均匀分布，其单位长度上的分布载荷P=100N/cm。假定=0，墙梁的厚度t=0.1cm。在不计自重情况下，试求其位移和应力。MATLAB程序：%Finite Element AnalysisE=1;u=0;%对单元1进行分析xy1=0,0;6,6;0,6;A1=0.5*det(1,xy1(1,1),xy1(1,2);1,xy1(2,1),xy1(2,2);1,xy1(3,1),xy1(3,2);B1=zeros(3,6)。

2、;B1(1,1)=xy1(2,2)-xy1(3,2);B1(1,3)=xy1(3,2)-xy1(1,2);B1(1,5)=xy1(1,2)-xy1(2,2);B1(2,2)=xy1(3,1)-xy1(2,1);B1(2,4)=xy1(1,1)-xy1(3,1);B1(2,6)=xy1(2,1)-xy1(1,1);B1(3,1)=xy1(2,2)-xy1(3,2);B1(3,2)=xy1(3,1)-xy1(2,1);B1(3,3)=xy1(3,2)-xy1(1,2);B1(3,4)=xy1(1,1)-xy1(3,1);B1(3,5)=xy1(1,2)-xy1(2,2);B1(3,6)=xy1(2。

3、,1)-xy1(1,1);%对单元2进行分析xy2=0,0;6,0;6,6;A2=0.5*det(1,xy2(1,1),xy2(1,2);1,xy2(2,1),xy2(2,2);1,xy2(3,1),xy2(3,2);B2=zeros(3,6);B2(1,1)=xy2(2,2)-xy2(3,2);B2(1,3)=xy2(3,2)-xy2(1,2);B2(1,5)=xy2(1,2)-xy2(2,2);B2(2,2)=xy2(3,1)-xy2(2,1);B2(2,4)=xy2(1,1)-xy2(3,1);B2(2,6)=xy2(2,1)-xy2(1,1);B2(3,1)=xy2(2,2)-xy2(。

4、3,2);B2(3,2)=xy2(3,1)-xy2(2,1);B2(3,3)=xy2(3,2)-xy2(1,2);B2(3,4)=xy2(1,1)-xy2(3,1);B2(3,5)=xy2(1,2)-xy2(2,2);B2(3,6)=xy2(2,1)-xy2(1,1);%单元1的刚度矩阵K1=zeros(6,6);for i=1:3for j=1:3K1(2*i-1,2*j-1)=(0.1*E)/(4*(1-u2)*A1)*(B1(1,2*i-1)*B1(1,2*j-1)+(1-u)/2*B1(2,2*i)*B1(2,2*j);K1(2*i-1,2*j)=(0.1*E)/(4*(1-u2)*A。

5、1)*(u*B1(1,2*i-1)*B1(1,2*j-1)+(1-u)/2*B1(2,2*i)*B1(1,2*j-1);K1(2*i,2*j-1)=(0.1*E)/(4*(1-u2)*A1)*(u*B1(2,2*i)*B1(1,2*j-1)+(1-u)/2*B1(1,2*i-1)*B1(2,2*j);K1(2*i,2*j)=(0.1*E)/(4*(1-u2)*A1)*(B1(2,2*i)*B1(2,2*j)+(1-u)/2*B1(1,2*i-1)*B1(1,2*j-1);endendK1Ka=zeros(8,8);Ka(1:6,1:6)=K1;temp=Ka(1:1,:);Ka(1:1,:)=。

6、Ka(5:5,:);Ka(5:5,:)=temp;temp=Ka(2:2,:);Ka(2:2,:)=Ka(6:6,:);Ka(6:6,:)=temp;temp=Ka(3:3,:);Ka(3:3,:)=Ka(7:7,:);Ka(7:7,:)=temp;temp=Ka(4:4,:);Ka(4:4,:)=Ka(8:8,:);Ka(8:8,:)=temp;temp=Ka(3:3,:);Ka(3:3,:)=Ka(5:5,:);Ka(5:5,:)=temp;temp=Ka(4:4,:);Ka(4:4,:)=Ka(6:6,:);Ka(6:6,:)=temp;temp=Ka(:,1:1);Ka(:,1:1)=。

7、Ka(:,5:5);Ka(:,5:5)=temp;temp=Ka(:,2:2);Ka(:,2:2)=Ka(:,6:6);Ka(:,6:6)=temp;temp=Ka(:,3:3);Ka(:,3:3)=Ka(:,7:7);Ka(:,7:7)=temp;temp=Ka(:,4:4);Ka(:,4:4)=Ka(:,8:8);Ka(:,8:8)=temp;temp=Ka(:,3:3);Ka(:,3:3)=Ka(:,5:5);Ka(:,5:5)=temp;temp=Ka(:,4:4);Ka(:,4:4)=Ka(:,6:6);Ka(:,6:6)=temp;%Ka%单元2的刚度矩阵K2=zeros(6,6)。

8、;for i=1:3for j=1:3K2(2*i-1,2*j-1)=(0.1*E)/(4*(1-u2)*A2)*(B2(1,2*i-1)*B2(1,2*j-1)+(1-u)/2*B2(2,2*i)*B2(2,2*j);K2(2*i-1,2*j)=(0.1*E)/(4*(1-u2)*A2)*(u*B2(1,2*i-1)*B2(1,2*j-1)+(1-u)/2*B2(2,2*i)*B2(1,2*j-1);K2(2*i,2*j-1)=(0.1*E)/(4*(1-u2)*A2)*(u*B2(2,2*i)*B2(1,2*j-1)+(1-u)/2*B2(1,2*i-1)*B2(2,2*j);K2(2*i。

9、,2*j)=(0.1*E)/(4*(1-u2)*A2)*(B2(2,2*i)*B2(2,2*j)+(1-u)/2*B2(1,2*i-1)*B2(1,2*j-1);endendK2Kb=zeros(8,8);Kb(3:8,3:8)=K2;%Kb%总刚度矩阵K=zeros(8,8);K=Ka+Kb%列方程并求解syms f21 f22 f31 f32 q11 q12 q32 q42yy=K* q11; q12;0; 0; 0 ;q32; 0 ;q42 -0 ;-300; f21; f22 ;f31; f32 ;0; -300 ;f21 f22 f31 f32 q11 q12 q32 q42=solve(eval(yy);q= q11; q12;0; 0; 0 ;q32; 0 ;q42 f=0 ;-300; f21; f22 ;f31; f32 ;0; -300 D=(E/(1-u2)*1 u 0;u 1 0;0 0 (1-u)/2)R1=D*B1;R2=D*B2;stress1=R1*q(3);q(4);q(7);q(8);q(1);q(2)stress2=R2*q(3);q(4);q(5);q(6);q(7);q(8)MATLAB程序运行结果：。