《有限元作业_三角形单元求解》由会员分享，可在线阅读，更多相关《有限元作业_三角形单元求解（22页珍藏版）》请在金锄头文库上搜索。
1、. . . . 有限元作业年 级2015级学 院机电工程学院专业名称班级学号学生姓名2016年05月 .下载可编辑.如下图所示为一受集中力P作用的结构，弹性模量E为常量，泊松比V=1/6,厚度为I=1。按平面应力问题计算，运用有限元方法，分别采用三角形及四边形单元求解，求节点位移及单元应力（要求三角形单元数量不少于4个，四边形单元不少于2个）图(一)图（二）三角形单元求解图（三）四边形单元求解（1） 如图划分三角形单元，工分成四个分别为（2） 如图分别进行编号1、2、3、4、5、6，并建立坐标系（3） 编程进行求解，得出结果，其中假设力P=2000N调用Triangle2D3Node_Stif
2、fness函数，求出单元刚度矩阵k1 = 1.0e+06 * 7.2857 -3.0000 -2.1429 0.8571 -5.1429 2.1429 -3.0000 7.2857 2.1429 -5.1429 0.8571 -2.1429 -2.1429 2.1429 2.1429 0 0 -2.1429 0.8571 -5.1429 0 5.1429 -0.8571 0 -5.1429 0.8571 0 -0.8571 5.1429 0 2.1429 -2.1429 -2.1429 0 0 2.1429k2 = 1.0e+06 * 5.1429 0 -5.1429 0.8571 0 -0.
3、8571 0 2.1429 2.1429 -2.1429 -2.1429 0 -5.1429 2.1429 7.2857 -3.0000 -2.1429 0.8571 0.8571 -2.1429 -3.0000 7.2857 2.1429 -5.1429 0 -2.1429 -2.1429 2.1429 2.1429 0 -0.8571 0 0.8571 -5.1429 0 5.1429k3 = 1.0e+06 * 2.1429 0 -2.1429 -2.1429 0 2.1429 0 5.1429 -0.8571 -5.1429 0.8571 0 -2.1429 -0.8571 7.285
4、7 3.0000 -5.1429 -2.1429 -2.1429 -5.1429 3.0000 7.2857 -0.8571 -2.1429 0 0.8571 -5.1429 -0.8571 5.1429 0 2.1429 0 -2.1429 -2.1429 0 2.1429k4 = 1.0e+06 * 2.1429 0 -2.1429 -2.1429 0 2.1429 0 5.1429 -0.8571 -5.1429 0.8571 0 -2.1429 -0.8571 7.2857 3.0000 -5.1429 -2.1429 -2.1429 -5.1429 3.0000 7.2857 -0.
5、8571 -2.1429 0 0.8571 -5.1429 -0.8571 5.1429 0 2.1429 0 -2.1429 -2.1429 0 2.1429调用Triangle2D3Node_Assembly函数，求出总体刚度矩阵求出的节点位移U = 0 0 0 0 -0.0004 0.0008 0.0005 0.0010 0.0007 0.0023 -0.0007 0.0026调用Triangle2D3Node_Stress函数，求出应力,S1、S2、S3、中求出的分别为Sx,Sy,SxyS1 = 1.0e+03 * -4.4086 -0.7348 3.5914S2 = 1.0e+03
6、* 4.4086 -0.6405 0.4086S3 = 1.0e+03 * 1.8907 -1.0601 2.1093S4 = 1.0e+03 * -1.8907 2.10931.8907二、（1）如图划分四边形单元，工分成四个分别为（2）如图分别进行编号1、2、3、4、5、6，并建立坐标系（3）编程进行求解，得出结果，其中假设力P=2000N调用 Quad2D4Node_Stiffness函数，求出单元刚度矩阵 调用Quad2D4Node_Assembly函数，求出求出总体刚度矩阵求出节点位移U = 0 0 0 0 0.0012 0.0017 -0.0012 0.0017 0.0016 0.
7、0049 -0.0017 0.0052调用Quad2D4Node_Stress函数，求出单元应力中的的S1、S2、S3分别为Sx,Sy,Sxy应力分量S1 = 1.0e+03 * 0.0000 -0.2478 2.0000S2 = 1.0e+07 * 0.6856 4.1135 -1.7137程序附录一、1、三角形单元总程序：E=1e7;NU=1/6;t=1;ID=1;%调用Triangle2D3Node_Stiffness函数，求出单元刚度矩阵k1=Triangle2D3Node_Stiffness(E,NU,t,0,1,0,0,1,1,ID)k2=Triangle2D3Node_Stiff
8、ness(E,NU,t,0,0,1,0,1,1,ID)k3=Triangle2D3Node_Stiffness(E,NU,t,1,1,1,0,2,0,ID)k4=Triangle2D3Node_Stiffness(E,NU,t,2,0,2,1,1,1,ID)%调用Triangle2D3Node_Assembly函数，求出总体刚度矩阵KK = zeros(12,12);KK=Triangle2D3Node_Assembly(KK,k1,1,2,3);KK=Triangle2D3Node_Assembly(KK,k2,2,4,3);KK=Triangle2D3Node_Assembly(KK,k3
9、,3,4,5);KK=Triangle2D3Node_Assembly(KK,k4,5,6,3)% 边界条件的处理及刚度方程求解k=KK(5:12,5:12) p=0;0;0;0;0;0;0;2000u=kp%支反力的计算U=0;0;0;0;u %为节点位移P=KK*U%调用Triangle2D3Node_Strain函数，求出应变SN1、SN2、SN3中求出的分别为SNx,SNy,SNxyu1=U(1);U(2);U(3);U(4);U(5);U(6);u2=U(3);U(4);U(7);U(8);U(5);U(6);u3=U(5);U(6);U(7);U(8);U(9);U(10);u4=
10、U(9);U(10);U(11);U(12);U(5);U(6); SN1=Triangle2D3Node_Strain(0,1,0,0,1,1,u1) SN2=Triangle2D3Node_Strain(0,0,1,0,1,1,u2) SN3=Triangle2D3Node_Strain(1,1,1,0,2,0,u3) SN4=Triangle2D3Node_Strain(2,0,2,1,1,1,u4) %调用Triangle2D3Node_Stress函数，求出应力,S1、S2、S3、中求出的分别为Sx,Sy,Sxy u1=U(1);U(2);U(3);U(4);U(5);U(6);u2
11、=U(3);U(4);U(7);U(8);U(5);U(6);u3=U(5);U(6);U(7);U(8);U(9);U(10);u4=U(9);U(10);U(11);U(12);U(5);U(6); S1=Triangle2D3Node_Stress(E,NU,0,1,0,0,1,1,u1,ID) S2=Triangle2D3Node_Stress(E,NU,0,0,1,0,1,1,u2,ID) S3=Triangle2D3Node_Stress(E,NU,1,1,1,0,2,0,u3,ID) S4=Triangle2D3Node_Stress(E,NU,2,0,2,1,1,1,u4,ID
12、)2、求刚度矩阵程序function k=Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)%该函数计算单元的刚度矩阵%输入弹性模量E，泊松比NU，厚度t%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)%输出单元刚度矩阵k(6X6)%-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 =
