有限元分析课程作业.pdf

《有限元分析课程作业.pdf》由会员分享，可在线阅读，更多相关《有限元分析课程作业.pdf（32页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

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

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。1m1my图 2 方案一：使用两个 3 节点三角形单元（2）各单元的刚度矩阵及刚度方程43 a.单元的几何和节点描述单元有 6 个节点位移自由度（DOF） 。将所有节点上的位移组成一个列阵，记作q(1)；同样，将所有节点

3、上的各个力也组成一个列阵，记作F(1)，则有q(1) （u1,v1,u2,v2,u4,v4)F(1) (Fx1,Fy1,Fx2,Fy2,Fx4,Fy4)同理，对于单元，有q(2) （u3,v3,u4,v4,u2,v2)F(2) (Fx3,Fy3,Fx4,Fy4,Fx2,Fy2)b.单元的位移场描述对于单元，设位移函数u(x, y) a0a1xa2yv(x, y) b (1-1)0b1xb2y由节点条件，在x xi, y yi处，有u(xi, yi) uiv(x yi 1,2,4 (1-2)i,i) vi将式(1-1)代入节点条件式(1-2)中，可求出式（中待定系数，即1y1a1u1x02Au2

4、x12y2u2A(a1u1a2u2a3u4) (1-3)4x4y41u1y1ay1112A1u221uy2A(b1u1b2u2b3u4) (1-4)441-1）1x11a21x22A1x4u1u2u41(c1u1c2u2c3u4) (1-5)2A1(a1v1a2v2a3v4)(1-6)2A1b1(b1v1b2v2b3v4)(1-7)2A1b2(c1v1c2v2c3v4)(1-8)2Ab0在式(1-3)式(1-8)中A 11x221x41x1y1y2y41(a1a2a3) (1-9)2y2 x2y4 x4y2y41y2b1 y2 y4 (1-10)1y41x2c1 x2 x41x4(1,2,3)

5、x2a1x4上式中的符号(1,2,3)表示下标轮换， 如12,2 3,31同时更换。将单元各节点的位置坐标x10,y10,x21,y20,x40,y41代入得a11,b1 1,c1 1a2 0,b21,c2 0a3 0,b3 0,c31A12将系数式(1-3)式(1-8)式代入(1-1)中，重写位移函数，并以节点位移的形式进行表示，有u(x, y) N1(x, y)u1 N2(x, y)u2 N4(x, y)u4 (1-11)v(x, y) N1(x, y)v1 N2(x, y)v2 N4(x, y)v4令Ni(26)1(aibixciy),i 1,2,3，则有形状函数矩阵N(x, y)2A(

6、1)NN(x, y) 100N1N200N2N400 1 x y0 x0yN41 x y0 x000y (1-12)位移函数式(1-11)写成矩阵形式，有u1v1(x, y) u(x, y)u02(2u(1)0 x06)v(x, y)1 x yy01 x y0 x0yv2 (1-13)u4v4对于单元，过程同上，有形状函数矩阵(x, y) (2N(2)1 x y01 x01 y0 6)01 x y01 x01 y (1-14)位移函数u01 x01 y(2u(2)(x, y) 6)(x, y)v(x, y)1+x y01 x y01 x0(1-15)c.单元的应变场描述对于单元，应变函数u3v

7、30u4 yv4u2v21u0 xxxx(1)(x, y) vu(x, y1)yyy0(3y)y)xyv(x,uvyxyxu1x0v1=01 x0 x0y0u2yy01 x y0 x0yv2yxu4v4u1-101000v1=0-10001u2 (1-16)-1-10110v2u4v4其中几何矩阵-101000(1)(x, y) 0-10001 (1-17)(3B6)-1-10110对于单元，应变方程x0(2)(x,1 x01 x01 y(3y) 1)0yy01 x y01 x0yxu30v31 yu4v4u2v2u3v310-1000u4 (1-18)=01000-1110-1-10v4u2

8、v2其中几何矩阵10-1000 (2)(x, y) 01000-1 (1-19)B(36)110-1-10d.单元的应力场描述1xxE(1)(x, y,z) yy12(31)xy 0(1)0 xx(1)(1)10yy=D (1-20)(33)(31)10 xy2其中，弹性系数矩阵 D 为1E(1)=D21(33) 010710110=23111-03 20131031060=3.7510130001113 20(1-21)(1)(31)=D(1)(33)(36)Bq(1)(1)(61)(36)Sq (1-22)(1)(61)(1)其中应力函数矩阵 S 为(1)(36)S(1)310 -1010

9、00-3-13001(1)(1)D B=3.751061300-10001=3.75106-1-31003(33)(36)001-1-10110-1-10110 (1-23)应力方程(1)为u1v1-3-13001-3-13001(1)(1)u266=3.7510-1-31003q=3.7510-1-31003Pa(31)-1-10110(61)-1-10110v2u4v4 (1-24)对于单元，过程同上。弹性系数矩阵 D 为310(2)=3.75106130 (1-25)D(33)001(2)应力函数矩阵 S 为31-300-1(2)=13-100-3 (1-26)S(36)110-1-10

10、(2)应力方程(2)为u3v331-300-131-300-1(2)(2)u466=3.751013-100-3q=3.751013-100-3Pa(31)110-1-10(61)110-1-10v4u2v2 (1-27)e.单元的势能表达K(1)是单元刚度矩阵，即K(1)(1)B(1)TD(1)B(1)d A(1)B(1)TD(1)B(1)dAt (1-28)其中薄板厚度t 0.1m。将式(1-17)、式(1-21)代入式(1-29)，得到单元的刚度阵-10 1(1)T(1)(1)6 BD B tA 3.7510 0.1000-1-1-1310 -10100000 1300-10001010

11、01-1-1011001100K(1)计算得u1v1u2v2u4v4u1v1u2v2u4v4K(1)3.755.6251.8751.8751.8757.53.757.51.8751.8751.8755.6255.625-1.8755.625001.875 510 -1.8751.87501.8751.87501.8751.87501.8751.87501.8755.6251.875005.625同理，得到单元的刚度阵为u3v3u4v4u2v2u1v1u2v2u4v4K(2)3.755.6251.8751.8751.8757.53.757.51.8751.8751.8755.6255.625-1

12、.8755.625001.875 510 -1.8751.87501.8751.87501.8751.87501.8751.87501.8755.6251.875005.625将两个单元按节点位移所对应的位置进行组装，得到总体刚度矩阵为KK K(1) K(2)节点力F (FRx1FRy1FPx2FPy2FPx3FPx3FRx4FRx4)T10.110103(FRx10 1 0 1 0 FRx40)2 0.5103(FRx10 1 0 1 0 FRx40)N系统的势能1 U W=qTKq-FTq(计算结果在下面呈现)2(4)边界条件的处理及方程求解边界条件为u1v1u4v40。因此，将针对节点

13、2 和节点3 的位移求解，节点 2 和节点 3 对应总体刚度阵 KK 中的第 3行到第 6 行、第3 列到第 6 列，则需从KK88中提出，置给k， 然后生成对应的载荷列阵 p， 再采用高斯消去法进行求解。 k=KK(3:6,3:6); p=500;0;500;0; u=kpu = 将列排成了行再计算支反力。在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力；先将上面得到的位移结果与边界条件的节点位移进行组合，得到整体的位移列阵 U(81)，再代回原整体刚度方程，计算出所有的节点力，按照位置关系找出对应的支反力。 U=0;0;u;0;0U = 0 0 0 0将列排成了行 P=

14、KK*UP =-500 500 0 500 0 -500将列排成了行所以，节点1 的支反力为FRx1 500N,FRy1 -176.4706N，节点2 的支反力为FRx2 500N,FRy2176.4706N。根据已求得的位移和支反力计算系统的势能。 A=*U*KK*U-P*UA =(5）结果分析上述支反力计算结果满足静力平衡，验证了以上求解过程及 MATLAB 算法的正确性。2、基于四节点四边形单元的有限元建模及分析（1）结构的离散化与编号如图 3 所示一个 4 节点矩形单元，单元的节点位移共有8 个自由度（DOF） 。节点编号为 1,2,3,4，各自的位置坐标为xi, yi,i 1,2,3

15、,4，各个节点的位移（分别沿x方向和y方向）为ui,vi,i 1,2,3,4。y4图 3 方案二：使用一个 4 节点矩形单元3（2）局部坐标系下单元的描述a.单元的几何和节点描述采用无量纲坐标=,xayb其中a 0.5,b 0.5。则单元四个节点的几何位置为11,112 1,213 1,3 141,4 1将所有节点上的位移组成一个列阵，记作 q；同样，将所有节点上的各个力也组成一个列阵，记作 F，则有q(81)(81) (u1v1u2x1v2Fx2u3Fy2v3u4Fx3Fy3v4)TFx4Fy4)TF (FFy1 b.单元的位移场描述设位移函数为u(x, y) a0a1xa2ya3xyv(x

16、, y) b0b1xb2yb3xy由节点条件，在x xi, y yi处，有u(xi, yi) uii 1,2,3,4v(xi, yi) vi将位移试函数代入节点条件中，求出待定系数a0,a1,a2,a3和b0,b1,b2,b3，再代入位移函数中，整理后得u(x, y) N1(x, y)u1 N2(x, y)u2 N3(x, y)u3 N4(x, y)u4v(x, y) N1(x, y)v1 N2(x, y)v2 N3(x, y)v3 N4(x, y)v4其中112x12y41N2(x,y) 12x12y41N3(x,y) 12x12y41N4(x,y) 12x12y4N1(x,y) 如以无量纲

17、坐标来表达，可写成Ni11i1i,i 1,2,3,44将11,11，2 1,21，3 1,3 1，41,4 1带入上式得到形状函数矩阵11141N21-141N311-41N41-1-4N1写成矩阵形式，有u1v1u20 v2qN4u3N2881v3u4v4u(x, y)N1(x, y) 0uv(x, y)(21)0N1N200N2N300N3N40c.单元的应变场描述单元应变为xx(x, y) yyNqBq(31)(32)2881(38)81xy其中几何矩阵B(x, y)为x(x, y) 0NB(38)(32)28y0 N1y0 x0N1N200N2N300N3N400 N40(12y)0(

18、12y)012y012y1=012x012x0(12x)0(12x)212x(12y)(12x)(12y)(12x)12y12x12y d.单元的应力场描述应力表达式为DD BqSq(31)(33)(31)(33)(38)(81)(38)(81)其中，应力函数矩阵S DB。e.单元的势能表达以上已将单元的三大基本变量u,用基于节点的位移列阵q来表示； 将其代入单元势能表达式中， 有=1qTKqFTq，2其中K为 4 节点矩形单元的刚度矩阵，即k11kkTK B DBdAt 2122Ak31k32k41k42sym.k33k43k44其中，t 为薄板的厚度，t 0.1m，上式的各个字块矩阵为Tk

19、rs 0.1BrDBs(22)(23)(33)(32)r,s 1,2,3,4f.单元刚度阵及刚度方程单元刚度阵在上面已经列出。将单元的势能对节点位移q取一阶极值，可得到单元的刚度方程(88)(81)KqF(81)(4)边界条件的处理及方程求解处理方法与 3 节点三角形单元一致，利用上述求解程序具有的可移植性，简化了求解过程。 k=K(3:6,3:6); p=500;0;500;0; u=kpu = * 将列排成了行再计算支反力。同样注意按照位置关系找出对应的支反力。 U=0;0;u;0;0U = *0 0 0 0将列排成了行 P=K*UP =-500 500 0 500 0 -500将列排成了

20、行所以，节点 1 的支反力为FRx1 500N,FRy1 -111.1111N，节点2 的支反力为FRx2 500N,FRy2111.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，m，平面问题性质指示参数（1 为平面应力，2 为平面应变） ，输出单元刚度矩阵 k（66） 。Triangle2D3Node_Assemble(KK,k,i,j,m)

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

23、,xi,yi,xj,yj,xm,ym,ID)%该函数计算单元的刚度矩阵%输入弹性模量 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 gamma

24、j 0 gammam; gammai betai gammaj betaj gammam betam/(2*A);if ID=1 D=(E/(1-NU*NU)*1 NU 0;NU 1 0;0 0 (1-NU)/2;elseif ID=2 D=(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%-DO

25、F(1)=2*i-1;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:6 for 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)%该函数计算单元的应力%输入弹性模量 E、泊松比 NU 和厚度 t%输入平面问题性质指示参数 ID（1 位平面应力，2 为平面应变） ，单元的位移列阵 u（6*1）%

26、输出单元的应力 stress（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 gammam betam/(2*A);if ID=1 D=(E/(1-NU*NU)*1 NU 0;NU 1 0;0 0 (1-N

27、U)/2;elseif ID=2 D=(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)k1=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,

28、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=kpU=0;0;u;0;0P=KK*U A=*U*KK*U-P*U（3)计算结果应变CC = +05 * 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0位移U = 0 0 0 0节点力P = 0 0其中，节点1 的支反力为FRx1 500N,FRy1 -176.4706N，节点2 的支反力为FRx2 500N

29、,FRy2176.4706N。势能A =单元刚度阵KK = +05 * 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 02. 基于四节点四边形单元的有限元建模及分析（1）程序编写说明Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID)该函数计算单元的刚度矩阵，输入模量 E，泊松比 NU，厚度 h，4 个节点 i，j，m，p，平面问题性质指示参数ID（1为平面应力，2 为平面应变） ，输出单元刚度矩阵 k（88） 。Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID)该

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

31、面问题性质指示参数 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)+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)

32、/4-b*(s-1)/4;B2=a*(1-t)/4-b*(-1-s)/4 0;0c*(-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*(-t-1)/4-b*(1-s)/4 0;0c*(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 B

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

34、%functionstress=Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID)%该函数计算单元的应力%输入弹性模量 E、泊松比 NU 和厚度 h%输入平面问题性质指示参数 ID（1 位平面应力，2 为平面应变）%输入单元的位移列阵 u（8*1）%输出单元的应力 stress（3*1） ，由于它为常应力单元，则单元的应力分量 Sx，Sy，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

35、*(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;0c*(-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

36、)/4; c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4;B4=a*(-t-1)/4-b*(1-s)/4 0;0c*(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=1 D=(E/(1-NU*NU)*1 NU 0;NU 1 0;0

37、0 (1-NU)/2;elseif ID=2 D=(E/(1+NU)/(1-2*NU)*1-NU NU 0;NU 1-NU 0;00 (1-2*NU)/2;endstr1=D*B*u;str2=subs(str1,s,t,0,0);stress=double(str2);end=Quad2D4Node_Stiffness(E,NU,h,0,0,1,0,1,1,0,1，)K=Quad2D4Node_Stiffness(E,NU,h,0,0,1,0,1,1,0,1,1) k=K(3:6,3:6); p=500;0;500;0; u=kp U=0;0;u;0;0 P=K*U A=*U*K*U-P*U

38、stress=Quad2D4Node_Stress(E,NU,0,0,1,0,1,1,0,1,10(-3)*0;0;0;0,1)（3)计算结果应变C = +06 *单元刚度阵K =+06 *位移U = * 0 0 0 0节点力P =其中，节点 1 的支反力为FRx1 500N,FRy1 -111.1111N，节点2 的支反力为FRx2 500N,FRy2111.1111N。势能A =应力stress = +04 * 0四、结束语通过本学期有限元方法的学习，我了解在力学计算过程中，将结构“化整为零” ，即将结构分成有限个单元进行处理的方法，十分巧妙。由于存在几何、物理方程的变换，就不可避免涉及到矩阵变换， 再到具体实现的软件， 如 MATLAB、ANSYS 等，将对不同模型的算法高度集成与软件本身，减少了工程人员的负担。在接下来的学习过程中，应该加深对有限元基本原理的认识，灵活运用，做到熟练掌握不同的有限元分析软件。