有限元分析课程作业_第1页
有限元分析课程作业_第2页
有限元分析课程作业_第3页
有限元分析课程作业_第4页
有限元分析课程作业_第5页
已阅读5页,还剩34页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

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-

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论