有限元分析课程作业

1、有限元分析课程作业任课教师：徐亚兰学生姓名：陈新杰学 号：班级：1304012时间：2016-01-05、问题描述及分析问题：如图1所示，有一矩形平板，在右侧受到P=10KN/m的分布力，材料常数为：弹性模量E 1 107Pa；泊松比 1/3 ； 板的厚度为t=;试按平面应力问题利用三角形与矩形单元分 别计算各个节点位移及支座反力。P=10KN/m图1平面矩形结构的有限元分析分析：使用两种方案：一、基于 3节点三角形单元的有限元建模，将矩形划分为两个 3节点三角形单元；二、基于4节点矩形单元的有限元建模，使用一个4节点矩形单元。利用MATLABM牛计算由各要求量，再将两种方案的计算结果进行比较

2、、分析、得由结论。二、有限元建模及分析1、基于3节点三角形单元的有限元建模及分析(1)结构的离散化与编号如图2所示，将平面矩形结构分为两个 3节点三角形单元。单元三个节点的编号为1, 2, 4,单元三个节点的编号为3, 4, 2,各个节点的位置坐标为xi,yi ,i 1,2,3,4 ,各个节点的位移（分别沿X方向和y方向）为Ui，Vi ,i 1,2,3,4。图2方案一：使用两个 3节点三角形单元（2）各单元的刚度矩阵及刚度方程 a.单元的几何和节点描述单元有6个节点位移自由度（DOF。将所有节点上的位移组成一个列阵，记作 q；同样，将所有节点上的各个力也组成一个列阵，记作 F（1）,则有q（u

3、1，V1,U2,V2,U4,V4）F （Fx1，Fy1，Fx2，Fy2，Fx4，Fy4）同理，对于单元，有q（U3,V3, U4,V4,U2,V2）F (Fx3,Fy3,Fx4,Fy4,Fx2,Fy2)b.单元的位移场描述对于单元，设位移函数u(x,y) % ax a2y(1-1)v(x,y) b0 bix b2y由节点条件，在x x,y y处，有u(xi,yi) uii 1,2,4v(xi, yi) vi(1-2)将式（1-1）代入节点条件式（1-2）中，可求生式（1-1 )中待定系数，即12A12A12AU1x1U2x2U4x4U2U4x1y1y2y4y1y2y4U1又2 U2x4 u4L

4、gv1 a2v22A1 ,、(a1U1 a2U2 a3U4) 2A1 ,.、(b1U1 b2U2 b3U4) 2A1(G3 C2U2 C3U4)2Aa3v4)(1-6)(1-3)(1-4)(1-5)1-(b!v1 b2v2 b3v4)(1-7)2A1 ,(qv1 C2v2 C3v4)(1-8)2 A在式（1-3） 式（1-8）中1.一 一 一、2(a a? a3)(1-9)1Xy11x2、21x4y4X2y2X4 V4n ； ；2y2 V、(1-10)1 x2qx2 x42 x4 4(1,2,3)上式中的符号(1,2,3)表示下标轮换，如1 2,2 3,3 1同时更换。将单元各节点的位置坐标x

5、1 0,y1 0,x2 1,y2 0,x4 0,y、1代入得a1 1,b11,q 1a2 0,b2 1,C2 0a30h 0,C3 1A 13将系数式(1-3)式(1-8)式代入(1-1)中，重写位移函数，并以节点位移的形式进行表示，有(1-11)u(x,y) N(x, y)。 N2(x,y)u2 N、(x, y)u、v(x,y) N(x,y)v1 N2(x,yM N、(x, y)v、N(x, y)令M(ai bix jy),i 1,2,3,则有形状函数矩阵2A(1)N (x, y)(2 6)N10N20N400N10N20N41 x y 0 x 0 y 001 x y 0 x 0 y(1-1

6、2)位移函数式(1-11)写成矩阵形式，有UiVi u(x, y) 1 x y 0 x 0 y 0U2u (x, y)(2U6)v(x, y)01 x y 0 x 0 yV2U4V4对于单元，过程同上，有形状函数矩阵1 x y 01 x 01 y 0N("y) 01 x y 01 x 01 y(1-13)(1-14)(2)/、 u(x, y) 1+x y11(x, y)(2U6)v(x, y)01U3V301 x 01 y 0 u4x y 01 x 01 y v4位移函数U2V2(1-15)c.单元的应变场描述/ 、(x,y)(3 1)xxyy xy对于单元，应变函数u(x,y) v

7、(x, y)0x1 x y 0 x 0 y 0=0y 01 x y 0 x 0 yUi viU2V2U4V4uiVi-10 10 0 0= 0 -1 0 0 0 1 u2(1-16)V2-1-10 110u4V4/ 、B (x,y) (3"6)-10 10000 -1 0 001-1-10 110(1-17)其中几何矩阵(2)(x,y)(31)u3V3u4V4u2V2对于单元，应变方程一 0 x1 x y 01 x 01 y 00y 01 x y 01 x 01 yu3V310-10003=01000-1u4(1-18)V4110-1-10u2V2其中几何矩阵1(2)B(x,y)0(

8、3 6)/1-10000-100-10-10(1-19)d.单元的应力场描述xx(x,y,Z) yy (3 1) xy001""2xxyyxy=R(1)(3(1)1)(1-20)其中，弹性系数矩阵(1)D) =10012D为107-21-3(1)(3 1)=只（B）9(1)(S)9其中应力函数矩阵s为(1)(1)(1)*D)袅=3753106 10-10-10-1-1应力方程为=3.75106(1-21(1-22)01 =3.750-3106 -1-1-1-3-1(1-23)u1-1300-3100-1011v1- 3 -1 3 0 0 11u2- 1-310032Pav2

9、- 1-101102u4v4-3(1) =3.75 106 -1(3 1)-11(1)3 q =3.75 100 (6 1)(1-25)(1-26)(1-24)对于单元，过程同上弹性系数矩阵D为31(2) 6D =3.75 106 1 3(3 3)00应力函数矩阵S为31-300-1(2)S=13-100-3(3 6)110-1-10应力方程(2) 为3(2)6=3.75 106 1(3 1)1-3-1 00000-1 -1-131-300(2)6-3q=3.75 10613-1000(6 1)110-1-1-1-30u3v3u4Pav4u2v2(1-27)e. 单元的势能表达K (1)是单元

10、刚度矩阵，即K(1)(1)B(1)TD(1)B(1)d(1)B(1)TD(1)B(1)dA t(1-28)其中薄板厚度t 0.1m o将式(1-17)、式(1-21)代入式(1-29),u1v1u2v2u4v47.53.755.6251.8751.8751.875u13.757.51.8751.8751.8755.625v1155.625K(1) 105-1.8755.625001.875u2-1.8751.87501.8751.8750v21.8751.87501.8751.8750u41.8755.6251.875005.625v4同理，得到单元的刚度阵为u3v3u4v4u2v27.53.

11、755.6251.8751.8751.875u13.757.51.8751.8751.8755.625v1255.625K(2)105-1.8755.625001.875u2-1.8751.87501.8751.8750v21.8751.87501.8751.8750u41.8755.6251.875005.625v4得到单元的刚度阵K(1) B(1)TD(1)B(1)tA 3.75 106 0.1计算得-10-10-1 -131001001000101010-10300-101-1-1100000010110将两个单元按节点位移所对应的位置进行组装，得到总体刚度矩阵为KK K(1) K(2)

12、节点力F ( FRx1 FRy1 FPx2 FPy2 FPx3 FPx3 FRx4 F Rx4 )13-0.1 10 10 (FRx1 0 10 10 FRx4 0)0.5 103(FRx1 0 10 10 FRx4 0)N系统的势能U W= ；qTKq-FTq (计算结果在下面呈现) (4)边界条件的处理及方程求解边界条件为U1 V1 u4 v4 0。因此，将针对节点 2和节点 3的位移求解，节点2和节点3对应总体刚度阵 KK中的第3 行到第6行、第3列到第6歹U,则需从KK8 8中提由，置给 k,然后生成对应的载荷列阵 p,再采用高斯消去法进行求解。>> k=KK(3:6,3:

13、6);>> p=500;0;500;0;>> u=kpu=将列排成了行再计算支反力。在得到整个结构的节点位移后，由原整 体刚度方程就可以计算由对应的支反力；先将上面得到的位 移结果与边界条件的节点位移进行组合，得到整体的位移列 阵U(8 1),再代回原整体刚度方程，计算由所有的节点力， 按照位置关系我由对应的支反力。>> U=0;0;u;0;0 将列排成了行>> P=KK*UP =-500 500 0 500 0 -500 将列排成了行所以，节点1 的支反力为FRx1500N, FRy1 -176.4706N ，节点2 的支反力为FRx2500N

14、,FRy2 176.4706N。根据已求得的位移和支反力计算系统的势能。>> A=*U'*KK*U-P'*UA =（5 ）结果分析上述支反力计算结果满足静力平衡，验证了以上求解过程及MATLA徵法的正确性。2、基于四节点四边形单元的有限元建模及分析（ 1）结构的离散化与编号如图 3 所示一个4 节点矩形单元，单元的节点位移共有8 个自由度（DOF） 。节点编号为1,2,3,4 ，各自的位置坐标为X，yi ,i 1,2,3,4 ,各个节点的位移（分别沿x方向和y方向）为ui,vi ,i 1,2,3,4 。4e图3方案二：使用一个4节点矩形单元(2)局部坐标系下单元的描

15、述 a.单元的几何和节点描述采用无量纲坐标=x y a, b其中a 0.5,b 0.5。则单元四个节点的几何位置为1, 11,1,1, 41111q(8 1)F1)b.单元的位移场描述将所有节点上的位移组成一个列阵，记作q；同样，将所有节点上的各个力也组成一个列阵，记作F,则有(U1 V1 U2 V2 U3 V3 U4 V4)T (Fx1 Fy1 Fx2 Fy2 Fx3 Fy3 Fx4 Fy4)T设位移函数为u(x,y) a ax v(x, y)凤 b)x由节点条件，在u(xi,y。 Uiiv(xi,yj via2 y a3xyb2 y b3xyx x,y y处，有1,2,3,4a0，a，a2

16、，%将位移试函数代入节点条件中，求生待定系数 和b0,b,b2,b3,再代入位移函数中，整理后得u(x,y) Ni(x, y)ui N2(x,y)U2 N3(x, y)U3 N4(x, y)u4 v(x, y) Ni(x, y)vi N?(x, yM N3(x,y)v3 N4(x, yM其中1Ni(x, y) - 1 2x 1 2y41N2(x,y)1 2x 1 2y41M(x,y) - 1 2x 1 2y41N4(x,y) 2y 0=-01 2x 1 2x 1 2y4如以无量纲坐标来表达，可写成1Ni -11 ,i 1,2,3,441带入上式J等 1 1, 1 1， 21, 2 1，31,

17、31， 4 1, 4得到形状函数矩阵1N11141N21-141N3113 41N4 1 1-14写成矩阵形式，有UiVlH(x，y)u(x, y)v(x, y)U2N10N20N30N40v20N10N20N30N4u3N qlU4V4C.单元的应变场描述单元应变为xx(31)(x，y)- Nq 且q(3 1)(3 2) 1 1 2x 1 2y 8 8 1(3 8) 8 1xy其中几何矩阵B(x,y)为0N40xNi0N20N30N40y 0N10N20N30y x(1 2y)001 2x1 2x(1 2y)(1 2y)00(1 2x)(1 2x)(1 2y)1 2y 00(1 2x)(1

18、2x)1 2yd.单元的应力场描述应力表达式为DD息q(sq(3 1)(3 3)(3 1)(3 3)(3 8) (8 1)(3 8)(8 1)其中，应力函数矩阵S DB。e.单元的势能表达以上已将单元的三大基本变量U,用基于节点的位移列阵q来表示；将其代入单元势能表达式中，有=lqTKq FTq ,其中K为4节点矩形单元的刚度矩阵，即kiiK A BT DBdA tk21k22symk31k32k33k41k42k43k44其中，t为薄板的厚度，t 0.1m,上式的各个字块矩阵为krs。BTDBsr,s 1,2,3,4(2 2)(2 3) (3 3) (3 2)f.单元刚度阵及刚度方程单元刚度

19、阵在上面已经列由。将单元的势能对节点位移q取一阶极值，可得到单元的刚度方程风9）（F）（4）边界条件的处理及方程求解处理方法与3节点三角形单元一致，利用上述求解程序具有的可移植性，简化了求解过程。>> k=K(3:6,3:6);>> p=500;0;500;0;>> u=kp将列排成了行再计算支反力。同样注意按照位置关系我由对应的支反力。>> U=0;0;u;0;0U = *0 00 0 将列排成了行>> P=K*UP =-500 500 0 500 0 -500 将列排成了行所以，节点1的支反力为Frxi500N,FRyi -111

20、.1111N ,节点2 的支反力为FRx2500N,FRy2 111.1111N 。根据已求得的位移和支反力计算系统的势能。>> A=*U'*K*U-P'*UA =(5 )结果分析基于 4 节点矩形单元计算出的势能小于基于3 节点三角形单元计算出的结果，若将该系统分为更多的单元，计算精度也将提高。3. 两种方案的比较与分析从以上计算可以看出，用三角形单元计算时，由于形函数是完全一次式，因而其应变场在单元内均为常数；而对于四边形单元，其形函数带有二次式，计算得到的应变场和应力场都是坐标的一次函数，但不是完全的一次函数，对提高精度有一定作用；根据最小势能原理，势能越小，

21、则整体计算精度越高，比较两种单元计算得到的系统势能，可看出，在相同的节点自由度情况下，矩形单元的计算精度要比三角 形单元高。三、基于MATLAB勺编程实现1. 基于 3 节点三角形单元的有限元编程实现（ 1）程序编写说明Triangle2D3Node_Stiffness（E,NU,t,xi,yi,xj,yj,xm ,ym,ID）该函数计算单元的刚度矩阵，输入模量 E,泊松比NU, 厚度t ,三个节点i , j , mi平面问题性质指示参数（1为平 面应力，2 为平面应变），输出单元刚度矩阵k（ 6 6） 。Triangle2D3Node_Assemble（KK,k,i,j,m）该函数进行单元刚

22、度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i , j , mi输由整体刚度矩阵 K左Triangle2D3Node_Stress（E,NU,xi,yi,xj,yj,xm,ym,u,ID）该函数计算单元的应力，输入弹性模量 E,泊松比NU, 厚度t ,三个节点i , j , mi平面问题性质指示参数（1为平 面应力，2 为平面应变），单元的位移列阵u（6 1），输出单元的应力，由于它为常应力单元，则单元的应力分量为Sx，Sy， Sxy。（ 2）程序清单%Triangle2D3Node%begin% functionk=Triangle2D3Node_Stiffness(E,NU,t,xi,

23、yi,xj,yj,xm,ym,ID)%该函数计算单元的刚度矩阵%俞入弹T生模量E、泊松比NU和厚度t%俞入3个节点i , j , m的坐标xi,yi,xj,yj,xm,ym%输入平面问题性质指示参数ID( 1 位平面应力，2为平面应变)%输入单元刚度矩阵k( 6*6)%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;

24、gammai betai gammaj betaj gammambetam/(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;00 (1-2*NU)/2;endk=t*A*B'*D*B;end%function z=Triangle2D3Node_Assemble(KK,k,i,j,m)%该函数进行单元刚度矩阵的组装%输入单元刚度矩阵k%输入单元的节点编号i ， j ， m%输入整体刚度矩阵KKrf%DOF(1)=2*i-1

25、;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;DOF(5)=2*m-1;DOF(6)=2*m;for n1=1:6for n2=1:6KK(DOF(n1),DOF(n2)=KK(DOF(n1),DOF(n2)+k(n1,n2);endendz=KK;% functionstress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)%该函数计算单元的应力%俞入弹T生模量E、泊松比NU和厚度t%输入平面问题性质指示参数ID( 1 位平面应力，2为平面应变) ，单元的位移列阵u( 6*1 )%输出单元的应力stress

26、( 3*1 ) ， 由于它为常应力单元，则单元的应力分量Sx， Sy， Sxy%A=(xi*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=batai 0 bataj 0 betam 0;0 gammai 0 gammaj 0 gammam;gammai betai gammaj betaj gammambetam/(2*A);if ID=1D=(E/(1-NU*NU)*1 NU 0;NU 1 0;0 0 (1-NU)/2;elseif ID=

27、2D=(E/(1+NU)/(1-2*NU)*1-NU NU 0;NU 1-NU 0;00 (1-2*NU)/2;endstress=D*B*u;%> > E=1E7;> > NU=1/3;> > t=;>> c=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,0,1,2)> > CC=zeros(8,8);>> CC=Triangle2D3Node_Assemble(KK,k1,1,2,4) ；>> CC=Triangle2D3Node_Assemble(KK,k1,3,4,2

28、)>> k1=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,0,1,1)>> KK=zeros(8,8);>> KK=Triangle2D3Node_Assemble(KK,k1,1,2,4)>> KK=Triangle2D3Node_Assemble(KK,k1,3,4,2)>> k=KK(3:6,3:6);>> p=500;0;500;0;>> u=kp>>U=0;0;u;0;0>>P=KK*U>> A=*U'*KK*U-P&#

29、39;*U( 3) 计算结果应变CC =+05 *0000000000000000位移U =00节点力P =00其中，节点1 的支反力为FRx1500N, FRy1 -176.4706N ，节点2 的支反力为FRx2500N,FRy2 176.4706N。势能A =单元刚度阵KK =+05 *00000000000000002. 基于四节点四边形单元的有限元建模及分析 （ 1）程序编写说明Quad2D4Node_Stiffness（E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID）该函数计算单元的刚度矩阵，输入模量E,泊松比NU,厚度h, 4个节点i , j , m p,平面

30、问题性质指示参数ID （ 为平面应力，2 为平面应变），输出单元刚度矩阵k（ 8 8） 。Quad2D4Node_Stress（E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID）该函数计算单元的应力，输入弹性模量E,泊松比NU, 厚度h, 4个节点i , j , mi p,平面问题性质指示参数ID (1 为平面应力，2为平面应变),单元的位移列阵u( 8 1),输 出单元的应力，由于它为常应力单元，则单元的应力分量为Sx， Sy， Sxy。( 2)程序清单%Quad2D4Node%begin%functionk=Quad2D4Node_Stiffness(E,NU,h,xi,

31、yi,xj,yj,xm,ym,xp,yp,ID)%该函数计算单元的刚度矩阵%俞入弹T生模量E、泊松比NU和厚度h%输入4 个节点 i ， j ， m， p 的坐标xi,yi,xj,yj,xm,ym,xp,yp%输入平面问题性质指示参数ID( 1 位平面应力，2为平面应变)%输入单元刚度矩阵k( 8*8)%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)+

32、xj*(-1-s)+xm*(1+s)+xp*(1-s)/4;B1=a*(t-1)/4-b*(s-1)/4 0;0 c*(s-1)/4-d*(t-1)/4;c*(-1+s)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4;B2=a*(1-t)/4-b*(-1-s)/4 0;0 c*(-1-s)/4-d*(1-t)/4;c*(-1-s)/4-d*(1-t)/4 a*(1-t)/4-b*(-1-s)/4;B3=a*(t+1)/4-b*(s+1)/4 0;0 c*(s+1)/4-d*(t+1)/4;c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4;B4=a

33、*(-t-1)/4-b*(1-s)/4 0;0 c*(1-s)/4-d*(-t-1)/4;c*(1-s)/4-d*(-t-1)/4 a*(-t-1)/4-b*(1-s)/4;Bfirst=B1 B2 B3 B4;Jfirst=0 1-t t-s s-1;t-1 0 s+1 -s-t;s-t -s-1 0 t+1;1-s s+t -t-1 0;J=xi xj xm xp*yi;yj;ym;yp/8;B=Bfirst/J;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

34、0;NU 1-NU 0;0 0 (1-2*NU)/2;endBD=J*transpose(B)*D*B;r=int(int(BD,t,-1,1),s,-1,1);z=h*r;k=double(z);end% functionstress=Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID)%该函数计算单元的应力%俞入弹T生模量E、泊松比NU和厚度h%输入平面问题性质指示参数ID( 1 位平面应力，2为平面应变)%输入单元的位移列阵u( 8*1 )%输出单元的应力stress ( 3*1 ) ， 由于它为常应力单元，则单元的应力分量Sx， S

35、y， Sxy% sym s t;a=(yi*(s-1)+yi*(-1-s)+ym*(1+S)+yp*(1-s)/4;b=(yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t)/4c=(xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t)/4d=(xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s)/4B1=a*(t-1)/4-b*(s-1)/4 0;0 c*(s-1)/4-d*(t-1)/4;c*(-1+s)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4;B2=a*(1-t)/4-b*(-1-s)/4 0;0 c*(-1-S)/4-d*(1-t)/4;c*(-1-