弹力有限元作业

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=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=LinearTriangleAssemble(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*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(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)=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(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 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=D*B*u假设E=1，F=1，厚度t=1，杆件长度为1，最后算出的位移乘以F/Et，力乘以F，以课本5结点为坐标原点(一) 离散化单元编号结点i结点j结点m123415142324(二) 写出单元刚度矩阵通过调用MATLAB的LinearTriangleElementStiffness函数，得到两个单元刚度矩阵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.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.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.5143 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 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 -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.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 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 -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 -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 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.125，厚度为1.3，网格划分建立四个关键点并对该单元平面作进一步