《弹力有限元作业.doc》由会员分享，可在线阅读，更多相关《弹力有限元作业.doc（16页珍藏版）》请在三一文库上搜索。
1、.6-4对下图所示的离散结构，试求结点1，2的位移及铰支座3，4，5的反力（按平面应力问题计算）采用MATLAB进行运算程序：一、计算弹性模量E、泊松比NU、厚度t、节点坐标为（xi,yi）、(xj,yj)、(xm,ym)的单元刚度矩阵。p=1表明函数用于平面应力情况。p=2表明函数用于平面应变情况。function y=LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p)A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj)/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai
2、=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=betai 0 betaj 0 betam 0; 0 gammai 0 gammaj 0 gammam; gammai betai gammaj betaj gammam betam/(2*A);if p=1 D=(E/(1-NU*NU)*1 NU 0;NU 1 0; 0 0 (1-NU)/2;精品.elseif p=2 D=(E/(1+NU)/(1-2*NU)*1-NU NU 0; NU 1-NU 0;0 0 (1-2*NU)/2;endy=t*A*B*D*B;二、计算整体刚度矩阵K。function y=LinearTr
3、iangleAssemble(K,k,i,j,m)K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(1,1);K(2*i-1,2*i)=K(2*i-1,2*i)+k(1,2);K(2*i-1,2*j-1)=K(2*i-1,2*j-1)+k(1,3);K(2*i-1,2*j)=K(2*i-1,2*j)+k(1,4);K(2*i-1,2*m-1)=K(2*i-1,2*m-1)+k(1,5);K(2*i-1,2*m)=K(2*i-1,2*m)+k(1,6);K(2*i,2*i-1)=K(2*i,2*i-1)+k(2,1);K(2*i,2*i)=K(2*i,2*i)+k(2,2);K(2
4、*i,2*j-1)=K(2*i,2*j-1)+k(2,3);K(2*i,2*j)=K(2*i,2*j)+k(2,4);K(2*i,2*m-1)=K(2*i,2*m-1)+k(2,5);K(2*i,2*m)=K(2*i,2*m)+k(2,6);K(2*j-1,2*i-1)=K(2*j-1,2*i-1)+k(3,1);K(2*j-1,2*i)=K(2*j-1,2*i)+k(3,2);K(2*j-1,2*j-1)=K(2*j-1,2*j-1)+k(3,3);K(2*j-1,2*j)=K(2*j-1,2*j)+k(3,4);K(2*j-1,2*m-1)=K(2*j-1,2*m-1)+k(3,5);K(
5、2*j-1,2*m)=K(2*j-1,2*m)+k(3,6);K(2*j,2*i-1)=K(2*j,2*i-1)+k(4,1);K(2*j,2*i)=K(2*j,2*i)+k(4,2);K(2*j,2*j-1)=K(2*j,2*j-1)+k(4,3);K(2*j,2*j)=K(2*j,2*j)+k(4,4);K(2*j,2*m-1)=K(2*j,2*m-1)+k(4,5);K(2*j,2*m)=K(2*j,2*m)+k(4,6);K(2*m-1,2*i-1)=K(2*m-1,2*i-1)+k(5,1);K(2*m-1,2*i)=K(2*m-1,2*i)+k(5,2);K(2*m-1,2*j-1
6、)=K(2*m-1,2*j-1)+k(5,3);K(2*m-1,2*j)=K(2*m-1,2*j)+k(5,4);K(2*m-1,2*m-1)=K(2*m-1,2*m-1)+k(5,5);K(2*m-1,2*m)=K(2*m-1,2*m)+k(5,6);K(2*m,2*i-1)=K(2*m,2*i-1)+k(6,1);K(2*m,2*i)=K(2*m,2*i)+k(6,2);精品.K(2*m,2*j-1)=K(2*m,2*j-1)+k(6,3);K(2*m,2*j)=K(2*m,2*j)+k(6,4);K(2*m,2*m-1)=K(2*m,2*m-1)+k(6,5);K(2*m,2*m)=K(
7、2*m,2*m)+k(6,6);y=K;三、计算单元位移矢量为u时的单元应力。p=1表明函数用于平面应力情况。p=2表明函数用于平面应变情况。该函数返回单元应力矢量。function y=LinearTriangleElementStresses(E,NU,xi,yi,yj,xm,ym,p,u)A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj)/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=betai 0 betaj 0 betam 0; 0 gammai 0 g
8、ammaj 0 gammam; gammai betai gammaj betaj gammam betam/(2*A);if p=1 D=(E/(1-NU*NU)*1 NU 0;NU 1 0; 0 0 (1-NU)/2;elseif p=2 D=(E/(1+NU)/(1-2*NU)*1-NU NU 0; NU 1-NU 0;0 0 (1-2*NU)/2;endy=D*B*u假设E=1，F=1，厚度t=1，杆件长度为1，最后算出的位移乘以F/Et，力乘以F，以课本5结点为坐标原点(一) 离散化单元编号结点i结点j结点m123415142324(二) 写出单元刚度矩阵通过调用MATLAB的Lin
9、earTriangleElementStiffness函数，得到两个单元刚度矩阵k1,k2和k3，每个矩阵都是66的。程序调用输出结果为： E=1E =精品. 1 NU=1/6NU = 0.1667 t=1t = 1 k1=LinearTriangleElementStiffness(E,NU,t,0,1,1,2,0,2,1)k1 = 0.2143 0 0 -0.2143 -0.2143 0.2143 0 0.5143 -0.0857 0 0.0857 -0.5143 0 -0.0857 0.5143 0 -0.5143 0.0857 -0.2143 0 0 0.2143 0.2143 -0.
10、2143 -0.2143 0.0857 -0.5143 0.2143 0.7286 -0.3000 0.2143 -0.5143 0.0857 -0.2143 -0.3000 0.7286 k2=LinearTriangleElementStiffness(E,NU,t,1,2,0,1,1,1,1)k2 = 0.2143 0 0 -0.2143 -0.2143 0.2143 0 0.5143 -0.0857 0 0.0857 -0.5143 0 -0.0857 0.5143 0 -0.5143 0.0857 -0.2143 0 0 0.2143 0.2143 -0.2143 -0.2143 0
11、.0857 -0.5143 0.2143 0.7286 -0.3000 0.2143 -0.5143 0.0857 -0.2143 -0.3000 0.7286精品. k3=LinearTriangleElementStiffness(E,NU,t,0,0,1,1,0,1,1)k3 = 0.2143 0 0 -0.2143 -0.2143 0.2143 0 0.5143 -0.0857 0 0.0857 -0.5143 0 -0.0857 0.5143 0 -0.5143 0.0857 -0.2143 0 0 0.2143 0.2143 -0.2143 -0.2143 0.0857 -0.51
12、43 0.2143 0.7286 -0.3000 0.2143 -0.5143 0.0857 -0.2143 -0.3000 0.7286(三) 集成整体刚度矩阵由于结构有5个结点，所以整体刚度矩阵式1010的。程序调用输出结果为：精品. K=zeros(10,10)K = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
13、0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 K=LinearTriangleAssemble(K,k1,4,1,3)K = 0.5143 0 0 0 -0.5143 0.0857 0 -0.0857 0 0 0 0.2143 0 0 0.2143 -0.2143 -0.2143 0 0 0 0 0 0 0 0 0 0 0 0 0精品. 0 0 0 0 0 0 0 0 0 0 -0.5143 0.2143 0 0 0.7286 -0.3000 -0.2143 0.0857 0 0 0.0857 -0.2143 0 0 -0.3000 0.7286 0.2143 -
14、0.5143 0 0 0 -0.2143 0 0 -0.2143 0.2143 0.2143 0 0 0 -0.0857 0 0 0 0.0857 -0.5143 0 0.5143 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 K=LinearTriangleAssemble(K,k2,1,4,2)K = 0.7286 0 -0.2143 0.2143 -0.5143 0.0857 0 -0.3000 0 0 0 0.7286 0.0857 -0.5143 0.2143 -0.2143 -0.3000 0 0 0精品. -0.2143 0.0857 0
15、.7286 -0.3000 0 0 -0.5143 0.2143 0 0 0.2143 -0.5143 -0.3000 0.7286 0 0 0.0857 -0.2143 0 0 -0.5143 0.2143 0 0 0.7286 -0.3000 -0.2143 0.0857 0 0 0.0857 -0.2143 0 0 -0.3000 0.7286 0.2143 -0.5143 0 0 0 -0.3000 -0.5143 0.0857 -0.2143 0.2143 0.7286 0 0 0 -0.3000 0 0.2143 -0.2143 0.0857 -0.5143 0 0.7286 0
16、0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 K=LinearTriangleAssemble(K,k3,5,2,4)K = 0.7286 0 -0.2143 0.2143 -0.5143 0.0857 0 -0.3000 0 0 0 0.7286 0.0857 -0.5143 0.2143 -0.2143 -0.3000 0 0 0精品. -0.2143 0.0857 1.2429 -0.3000 0 0 -1.0286 0.3000 0 -0.0857 0.2143 -0.5143 -0.3000 0.9429 0 0 0.3000 -0.4286 -
17、0.2143 0 -0.5143 0.2143 0 0 0.7286 -0.3000 -0.2143 0.0857 0 0 0.0857 -0.2143 0 0 -0.3000 0.7286 0.2143 -0.5143 0 0 0 -0.3000 -1.0286 0.3000 -0.2143 0.2143 1.4571 -0.3000 -0.2143 0.0857 -0.3000 0 0.3000 -0.4286 0.0857 -0.5143 -0.3000 1.4571 0.2143 -0.5143 0 0 0 -0.2143 0 0 -0.2143 0.2143 0.2143 0 0 0
18、 -0.0857 0 0 0 0.0857 -0.5143 0 0.5143精品.(四) 引入边界条件本题的边界条件如下：U3x=U3y=U4x=U4y=U5x=U5y=U1x=U2x=0F1y=-0.5 F1x=0 (五) 求解 k=K(2:2:4,2:2:4)k = 0.7286 -0.5143 -0.5143 0.9429 f=-0.5;0f = -0.5000 0 u=kfu = -1.1159 -0.6087所以：v1=-1.1159F/Et, v2=-0.6087F/Et(六) 求支座反力 U=0;-1.1159;0;-0.6087;0;0;0;0;0;0U = 0 -1.1159
19、 0 -0.6087 0 0精品. 0 0 0 0 F=K*UF = -0.1304 -0.5000 0.0870 -0.0000 -0.2391 0.2391 0.1522 0.2609 0.1304 0所以，F3x=-0.2391F, F3y=0.2391F F4x=0.1522F, F4y=0.2609F F5x=0.1304F, F5y=03-11挡水墙的密度为，厚度为，如下图1示，水的密度为，试求应力分量。已知：。精品.图1采用ansys软件进行解答解答步骤：1，采用国际单位制2，定义一个四节点平面，并定义其属性，混凝土坝的属性为密度2.8e3，弹性模量3.2e10，泊松比为0.12
20、5，厚度为1.3，网格划分建立四个关键点并对该单元平面作进一步划分：5*25个小单元格（40cmX40cm），并对各小单元进行节点编号、单元编号。4，在底边施加约束，选择线上dof。5，施加荷载。如下图所示：精品.对本题共有两荷载作用于挡水墙，一是挡水墙自重（输入重力加速度9.8），一是水体对其压力。本题将自重和水体压力组合为荷载组共同作用于挡水墙上，其中挡水墙以的厚度计算。模型求解运行ansys求解。求得应力应变图结果如下：、的应力图示如下图示：精品.观察应力图可大致得到如下规律：图：从上往下颜色越来越深，越来越大；图：从中间往两边颜色越来越深，越来越大在两边缘处取得各自的取大值；图：关于轴
21、对称，由上往下越来越大。所计算的结果输出另用txt结果输出。 理论值检验计算建立如图1所示的直角坐标系，经理论推导得到其应力分量理论解如下：；采用matlab进行检验计算结果：精品.程序为：%syms x y sx sy sxy;x=input(输入x坐标)y=input(输入y坐标)%!P2=1000;P1=2500;G=9.8;B=2.5;sx=2*1000*9.8*x3*y/8+3*1000*9.8*x*y/10-4*1000*9.8*x*y3/8-2800*10*x;sy=1000*9.8*x*(2*y3/8-3*y/4-1/2);sxy=-1000*9.8*x2*(3*y2/8-3/8)-1000*9.8*y*(-y3/8+3*y/20-2/80/y);sxsysxy计算结果与理论计算结果相差偏大，可能由于单元格划分过大引起。如有侵权请联系告知删除，感谢你们的配合！精品