有限元作业三角形单元求解.docx

1、有限元作业三角形单元求解有限元作业年 级 2015级学 院 机电工程学院 专业名称班级学号学生姓名2016年05月如下图所示为一受集中力P作用的结构，弹性模量E为常量，泊松比V=1/6, 厚度为1=1。按平面应力问题计算，运用有限元方法，分别采用三角形及四边形单元求解，求节点位移及单元应力（要求三角形单元数量不少于4个，四边 形单元不少于2个）图（一）图（二）三角形单元求解图(三)四边形单元求解(1) 如图划分三角形单元，工分成四个分别为 (2)如图分别进行编号1、2、3、4、5、6，并建立坐标系(3)编程进行求解，得出结果，其中假设力P=2000N调用Triangle2D3Node_Stif

2、fness函数，求出单元刚度矩阵 k1 =1.0e+067.2857-3.0000-2.14290.8571-5.14292.1429-3.00007.28572.1429-5.14290.8571-2.1429-2.14292.14292.142900-2.14290.8571-5.142905.1429-0.85710-5.14290.85710-0.85715.142902.1429-2.1429-2.1429002.1429k2 =1.0e+065.14290-5.14290.85710-0.857102.14292.1429-2.1429-2.142905.14292.14297.2

3、857-3.0000-2.14290.85710.8571-2.1429-3.00007.28572.1429-5.14290-2.1429-2.14292.14292.142900.857100.8571-5.142905.1429k3 =1.0e+062.14290-2.1429-2.142902.142905.1429-0.8571-5.14290.85710-2.1429-0.85717.28573.0000-5.1429-2.1429-2.1429-5.14293.00007.2857-0.8571-2.142900.8571-5.1429-0.85715.142902.14290-

4、2.1429-2.142902.1429k4 =1.0e+062.14290-2.1429-2.142902.142905.1429-0.8571-5.14290.85710-2.1429-0.85717.28573.0000-5.1429-2.1429-2.1429-5.14293.00007.2857-0.8571-2.142900.8571-5.1429-0.85715.142902.14290-2.1429-2.142902.1429函数，求出总体刚度矩阵调用 Triangle2D3Node_Assembly0. 7286-0,1000-0.2K39. 0857-0. 5143CL 2

5、143000000-0. 30000. 72即0.2 M30 0S57-0.214.30Qg000-0, 21430-2U30- 72S60Q-0, 3G00-0=3143D- 085?Q000CL 0057-0-5143D-a.3ooo00.2143Th 2 hBa00Q-0/51430adB570-0. 3QOO1. 45710-0.42册00a. 3ooo-0-5113-0.03570, 2143-0. 2143-0. 3WD001,45710-1.02W(L 300(10-0-2143-0,214300-0,5343-a伽0I.4E714-S5I43-0.51430000山 0857-

6、Ou 2143a-|a 02BS0-ft. OSS-0.21430000000CL 3000-0-5143-0.08570.728S0-0-2113-0.2130000仇 30000-0.2H300.7266-0.0S57-0r5H300007 5143-0,3143Q0-Ou 2143-0,065：Q. 12B8乳 300000D0仇 D857-CL 2U300-0. 2143-0.SU30.3M0D.72A6求出的节点位移000-0.00040.00080.00050.00100.00070.0023-0.00070.0026中求出的分别为调用 Triangle2D3Node_Stress

7、 函数，求出应力，S1、S2、S3、Sx,Sy,SxyS1 =1.0e+03 *-4.4086-0.73483.5914S2 =1.0e+034.4086-0.64050.408653=1.0e+03 *1.8907-1.06012.109354=1.0e+03 *-1.89072.10931.8907二 二、(1)如图划分四边形单元，工分成四个分别为(2)如图分别进行编号1、2、3、4、5、6,并建立坐标系(3)编程进行求解,得出结果，其中假设力P=2000N调用Quad2D4Node_Stiffness函数，求出单元刚度矩阵4. 8S71-1.5H00.2857-0.6429-2.4286

8、1.5000-2. 71430.6429-1.50001 85710. 6429-2.T143L5QOO-2.4286-0. 642d山 28670= 28570. M294.8571L50QQ-2.7143-CL 6429-Z 4286-L50W-0.6429-2.71431.50004. 85710.61290.2857-L 5000-2.4286-2.42861. 6000-2.71430.石竝94. 3571-1-5000CL 2857-0.6428L5000-2,42B6-0. C4290,2B57IrSOOO4B57i0, 9429一 Z 7143-2.7143-0. 6t2942

9、86-L 50000. 28B70.64294. 86711.5000Q.6429Du 2857-L5000-N 4 2鯛-0.6i2S-2- 7U3L 5Q004.S57Jk2 二L Oe+06 *4 8S711.5Q0C-2.7143-0. 6429-2.4286-1.50000, 28S7Ou 4291B 5000札 35710.61200. 2B57-LSOOO-2.4236-Du 6429-2/7143-2.71430. 04294.S571-i.noo0.2857-0. 6429-2, 4286LSOOO-0.64290. 2867-L.EOOO4.86710.6429*2.714

10、31.6000-2.4286-2. 4236-L 50000.235：4.85711-5000-2-7143-0.625-1.5000-2. 4286-0.6429-2.71431,50004 S57I0. S4290l 28570.2S5T-0. 642B-2.1286L 500S-2.71430.64294 8571-LSOOOQ.429-2. 7U3hSQOQ-N 42S6-Q+64290-2857-扎 5Q0Q4.85711. Ofl+Ofl *4.8571-1=5000CL 2S57山 439-2. 42BSK5000-2.711130.04290000-1.54004.8571CL

11、 6429-2,7143LGflOO2 42B6-0.61290,26570000D.2S6T0.6i24 8B71. 5000-2.7143-0.6429-2.42M00DQ-0.6129-2.71431. 50004. 8S71(L 64290.2557-1.5000-2.428fi0000-2.283.5000-瓷 71430. 64209-. 71430O571也-2. 7H3-CL 開2B- 42&6-1. 50 DOK500Q-2.4280-CL B42-9CL 285709.7U30-5.428A0. 6429CL 2857-L 5000-2u 42BG-2.71437*429-

12、2.42SC-3, 5000OU 571409.71430監 428CL50G0-2.71435 6429D仙0.2B67-I.S004-2. 4S860-6.42B609.71431.5000-2a 42&67、6429也 2857000-2. 71430.6429-2.428flh &DQO4,9571-1. EOOfl0, 2357-a. 6429000-0U S42B0.2S51.5000-2.-I, 5000L B571OL 6429-2-7143a0fl0-i 42B6-K5G00-2.7U3-D. A420. 2B57CL B42B4. 8571K50D0Q000-LEOOO临

13、428E0. C42B0.39E77 9429临 7143L 50004.8571求出节点位移0000.00120.0017 -0.00120.00170.00160.0049-0.00170.0052S3分别为调用Quad2D4Node_Stress函数，求出单元应力中的的 S1、S2、Sx,Sy,Sxy应力分量S1 =1.0e+03 *0.0000 -0.24782.0000S2 =1.0e+070.6856 4.1135-1.7137程序附录、1、三角形单元总程序：E=1e7;NU=1/6;t=1;ID=1;%调用Triangle2D3Node_Stiffness函数，求出单元刚度矩阵k

14、仁 Tria ngle2D3Node_Stiffness(E,NU,t,0,1,0,0,1,1,ID)k2=Tria ngle2D3Node_Stiffness(E,NU,t,0,0,1,0,1,1,ID)k3=Tria ngle2D3Node_Stiffness(E,NU,t,1,1,1,0,2,0,ID)k4=Tria ngle2D3Node_Stiffness(E,NU,t,2,0,2,1,1,1,ID)%调用Triangle2D3Node_Assembly 函数，求出总体刚度矩阵KK = zeros(12,12);KK=Tria ngle2D3Node_Assembly(KK,k1,1

15、,2,3);KK=Tria ngle2D3Node_Assembly(KK,k2,2,4,3);KK=Tria ngle2D3Node_Assembly(KK,k3,3,4,5);KK=Tria ngle2D3Node_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);

16、U ;U (5) ;U (6);u2=U (3);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);SN1= Tria ngle2D3Node_Strai n(0,1,0,0,1,1,u1)SN2=Tria ngle2D3Node_Strai n(0,0,1,0,1,1,u2)SN3=Tria ngle2D3Node_Strai n(1,1,1,0,2,0,u3)SN4=Tria ngle2D3Node_Strai n(2,0,2,1,1,1,u4

17、)%调用Triangle2D3Node_Stress函数，求出应力，S1、S2、S3、中求出的分别为 Sx,Sy,Sxyu仁U(1);U (2);U (3);U ;U (5) ;U (6);u2=U (3);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= Tria ngle2D3Node_Stress(E,NU,0,1,0,0,1,1,u1,ID)S2=Tria ngle2D3Node_Stress(E,NU,0,0,1,0,1,1,u2

18、,ID)S3=Tria ngle2D3Node_Stress(E,NU,1,1,1,0,2,0,u3,ID)S4=Tria ngle2D3Node_Stress(E,NU,2,0,2,1,1,1,u4,ID)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*(

19、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 ID = 1D = (E/(1-NU*NU)*1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2;elseif ID = 2D = (E

20、/(1+NU)/(1-2*NU)*1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2; endk= t*A*B*D*B;3、求整体刚度矩阵function z = Triangle2D3Node_Assembly(KK,k,i,j,m)%该函数进行单元刚度矩阵的组装%输入单元刚度矩阵k%输入单元的节点编号I、j、m%输出整体刚度矩阵KK% DOF(1)=2*i-1;DOF(2)=2*i;DOF(3)=2*j-1;DOF=2*j;DOF(5)=2*m-1;DOF(6)=2*m;for n仁 1:6for n2=1:6KK(DO F(n 1),DO F(n 2)= KK(

21、DO F(n 1),DO F(n 2)+k (n 1, n2);endendz=KK;4、求应变程序function strain=Triangle2D3Node_Strain(xi,yi,xj,yj,xm,ym,u)%该函数计算单元的应变%输入三个节点 i、j、m 的坐标 xi,yi,xj,yj,xm,ym%输入单元的位移列阵u(6X1)%输出单元的应力strain(3X1)，由于它为常应变单元，则单元的应变分量为SNx,SNy,SNz% A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj)/2;betai = yj-ym;betaj = ym-yi;betam

22、 = 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);strain = B*u;5、求应力程序function stress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)%该函数计算单元的应力%输入弹性模量E,泊松比NU ,厚度t%输入三个节点 i、j、m 的坐标 xi,yi,xj,

23、yj,xm,ym%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)，单元的位移列阵 u(6X1)%输出单元的应力stress(3X1)，由于它为常应力单元，则单元的应力分量为Sx,Sy,Sxy% 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 be

24、tai gammaj betaj gammam betam/(2*A);if ID = 1D = (E/(1-NU*NU)*1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2;elseif ID = 2D = (E/(1+NU)/(1-2*NU)*1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2; end stress = D*B*u;1、四边形单元总程序：E=1e7;NU=1/6;h=1;ID=1;%调用Quad2D4Node_Stiffness函数,求出单元刚度矩阵k1= Quad2D4Node_Stiffness(E,NU,h,0,1,0,0,1,0

25、,1,1,ID)k2= Quad2D4Node_Stiffness(E,NU,h,1,0,2,0,2,1,1,1,ID)%调用Quad2D4Node_Assembly函数,求出求出总体刚度矩阵KK=zeros(12,12);KK= Quad2D4Node_Assembly(KK,k1,1,2,3,4);KK= Quad2D4Node_Assembly(KK,k2,3,5,6,4)%边界条件的处理及刚度方程求解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%调用Quad2D4Node_Stress函数，

26、求出单元应力中的的 S1、S2、S3分别为Sx,Sy,Sxy应力分量u1=U(1);U (2) ;U (3) ;U ;U( 5);U (6) ;U (7);U(8);u2=U (5);U (6) ;U(9);U(10);U(11);U(12);U(7);(8);S1= Quad2D4Node_Stress(E,NU,0,1,0,0,1,0,1,1,u1,ID)S2= Quad2D4Node_Stress(E,NU,1,0,2,0,2,1,1,1,u2,ID)2、求刚度矩阵程序function k= Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,x

27、p,yp,ID)%该函数计算单元的刚度矩阵%输入弹性模量E，泊松比NU，厚度h%输入 4 个节点 i、j、m、p 的坐标 xi,yi,xj,yj,xm,ym,xp,yp%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)%输出单元刚度矩阵k(8X8)% syms s t;a = (yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s)/4;b = (yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t)/4;c = (xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t)/4;d = (xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*