有限元分析课程作业

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

2、结论。二、有限元建模及分析1、基于3节点三角形单元的有限元建模及分析（1）结构的离散化与编号如图2所示，将平面矩形结构分为两个3节点三角形单元。单元三个节点的编号为1，2，4，单元三个节点的编号为3，4，2，各个节点的位置坐标为，各个节点的位移（分别沿方向和方向）为。Xy1432图2 方案一：使用两个3节点三角形单元（2）各单元的刚度矩阵及刚度方程 a.单元的几何和节点描述单元有6个节点位移自由度（DOF）。将所有节点上的位移组成一个列阵，记作；同样，将所有节点上的各个力也组成一个列阵，记作，则有 同理，对于单元，有b.单元的位移场描述对于单元，设位移函数 (1-1)由节点条件，在处，有 (1

3、-2) 将式(1-1)代入节点条件式(1-2)中，可求出式（1-1）中待定系数，即 (1-3) (1-4) (1-5)(1-6)(1-7)(1-8)在式(1-3)式(1-8)中 (1-9) (1-10)上式中的符号(1,2,3)表示下标轮换，如同时更换。将单元各节点的位置坐标代入得 将系数式(1-3)式(1-8)式代入(1-1)中，重写位移函数，并以节点位移的形式进行表示，有 (1-11)令，则有形状函数矩阵 (1-12)位移函数式(1-11)写成矩阵形式，有 (1-13)对于单元，过程同上，有形状函数矩阵 (1-14)位移函数 (1-15)c.单元的应变场描述对于单元，应变函数 (1-16)

4、其中几何矩阵 (1-17)对于单元，应变方程 (1-18)其中几何矩阵 (1-19)d.单元的应力场描述 (1-20)其中，弹性系数矩阵D(1)为(1-21) (1-22)其中应力函数矩阵S(1)为 (1-23)应力方程为 (1-24)对于单元，过程同上。弹性系数矩阵D(2)为 (1-25)应力函数矩阵S(2)为 (1-26)应力方程为 (1-27)e.单元的势能表达是单元刚度矩阵，即 (1-28)其中薄板厚度。将式(1-17)、式(1-21)代入式(1-29)，得到单元的刚度阵计算得 同理，得到单元的刚度阵为将两个单元按节点位移所对应的位置进行组装，得到总体刚度矩阵为节点力 系统的势能(计算

5、结果在下面呈现)(4)边界条件的处理及方程求解边界条件为。因此，将针对节点2和节点3的位移求解，节点2和节点3对应总体刚度阵KK中的第3行到第6行、第3列到第6列，则需从KK中提出，置给k，然后生成对应的载荷列阵p，再采用高斯消去法进行求解。>> k=KK(3:6,3:6);>> p=500;0;500;0;>> u=kpu =0.0009 0.0001 0.0010 -0.0002 将列排成了行再计算支反力。在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力；先将上面得到的位移结果与边界条件的节点位移进行组合，得到整体的位移列阵U，再代回

6、原整体刚度方程，计算出所有的节点力，按照位置关系找出对应的支反力。>> U=0;0;u;0;0U = 0 0 0.0009 0.0001 0.0010 -0.0002 0 0将列排成了行>> P=KK*UP =-500 -176.4706 500 0 500 0 -500 176.4706将列排成了行所以，节点1的支反力为，节点2的支反力为。根据已求得的位移和支反力计算系统的势能。>> A=0.5*U'*KK*U-P'*UA = -0.4706(5）结果分析上述支反力计算结果满足静力平衡，验证了以上求解过程及MATLAB算法的正确性。2、基于

7、四节点四边形单元的有限元建模及分析（1）结构的离散化与编号如图3所示一个4节点矩形单元，单元的节点位移共有8个自由度（DOF）。节点编号为1,2,3,4，各自的位置坐标为，各个节点的位移（分别沿方向和方向）为。Xy4321图3 方案二：使用一个4节点矩形单元（2）局部坐标系下单元的描述a.单元的几何和节点描述采用无量纲坐标 其中。则单元四个节点的几何位置为 将所有节点上的位移组成一个列阵，记作q；同样，将所有节点上的各个力也组成一个列阵，记作F，则有 b.单元的位移场描述设位移函数为 由节点条件，在处，有 将位移试函数代入节点条件中，求出待定系数和，再代入位移函数中，整理后得 其中 如以无量纲

8、坐标来表达，可写成将带入上式得到形状函数矩阵写成矩阵形式，有 c.单元的应变场描述单元应变为 其中几何矩阵为 d.单元的应力场描述应力表达式为其中，应力函数矩阵。e.单元的势能表达以上已将单元的三大基本变量用基于节点的位移列阵来表示；将其代入单元势能表达式中，有，其中为4节点矩形单元的刚度矩阵，即其中，t为薄板的厚度，上式的各个字块矩阵为 f.单元刚度阵及刚度方程单元刚度阵在上面已经列出。将单元的势能对节点位移取一阶极值，可得到单元的刚度方程(4)边界条件的处理及方程求解 处理方法与3节点三角形单元一致，利用上述求解程序具有的可移植性，简化了求解过程。>> k=K(3:6,3:6)

9、;>> p=500;0;500;0;>> u=kpu =1.0e-03 *0.4815 0.1111 0.4815 -0.1111 将列排成了行再计算支反力。同样注意按照位置关系找出对应的支反力。>> U=0;0;u;0;0U =1.0e-03 * 0 0 0.4815 0.1111 0.4815 -0.1111 0 0将列排成了行>> P=K*UP =-500 -111.1111 500 0 500 0 -500 111.1111将列排成了行所以，节点1的支反力为，节点2的支反力为。根据已求得的位移和支反力计算系统的势能。>> A=

10、0.5*U'*K*U-P'*UA = -0.2407(5）结果分析基于4节点矩形单元计算出的势能小于基于3节点三角形单元计算出的结果，若将该系统分为更多的单元，计算精度也将提高。3.两种方案的比较与分析从以上计算可以看出，用三角形单元计算时，由于形函数是完全一次式，因而其应变场在单元内均为常数；而对于四边形单元，其形函数带有二次式，计算得到的应变场和应力场都是坐标的一次函数，但不是完全的一次函数，对提高精度有一定作用；根据最小势能原理，势能越小，则整体计算精度越高，比较两种单元计算得到的系统势能，可看出，在相同的节点自由度情况下，矩形单元的计算精度要比三角形单元高。三、基于MA

11、TLAB的编程实现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)该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i，j，m，输出整体刚度矩阵KK。Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm

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

13、平面应力，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; 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

14、=2 D=(E/(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;end%function z=Triangle2D3Node_Assemble(KK,k,i,j,m)%该函数进行单元刚度矩阵的组装%输入单元刚度矩阵k%输入单元的节点编号i，j，m%输入整体刚度矩阵KKrf%-DOF(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:6 KK(DOF(n1),DOF(n2)=KK

15、(DOF(n1),DOF(n2)+k(n1,n2); endendz=KK;%function stress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)%该函数计算单元的应力%输入弹性模量E、泊松比NU和厚度t%输入平面问题性质指示参数ID（1位平面应力，2为平面应变），单元的位移列阵u（6*1）%输出单元的应力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=

16、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-NU)/2;elseif ID=2 D=(E/(1+NU)/(1-2*NU)*1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2;endstress=D*B*u;%>> E=1E7;>> NU=1/3;&g

17、t;> t=0.1;>> 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,0,1,1)>> KK=zeros(8,8);>> KK=Triangle2D3No

18、de_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=0.5*U'*KK*U-P'*U（3)计算结果应变CC = 1.0e+05 * 9.3750 5.6250 -7.5000 -1.8750 0 0 -1.8750 -3.7500 5.6250 9.3750 -3.7500 -1.

19、8750 0 0 -1.8750 -7.5000 -7.5000 -3.7500 9.3750 0 -1.8750 -1.8750 0 5.6250 -1.8750 -1.8750 0 9.3750 -3.7500 -7.5000 5.6250 0 0 0 -1.8750 -3.7500 9.3750 5.6250 -7.5000 -1.8750 0 0 -1.8750 -7.5000 5.6250 9.3750 -3.7500 -1.8750 -1.8750 -1.8750 0 5.6250 -7.5000 -3.7500 9.3750 0 -3.7500 -7.5000 5.6250 0

20、-1.8750 -1.8750 0 9.3750位移U = 0 0 0.0009 0.0001 0.0010 -0.0002 0 0节点力P = -500.0000 -176.4706 500.0000 0 500.0000 0 -500.0000 176.4706其中，节点1的支反力为，节点2的支反力为。势能A = -0.4706单元刚度阵KK = 1.0e+05 * 7.5000 3.7500 -5.6250 -1.8750 0 0 -1.8750 -1.8750 3.7500 7.5000 -1.8750 -1.8750 0 0 -1.8750 -5.6250 -5.6250 -1.87

21、50 7.5000 0 -1.8750 -1.8750 0 3.7500 -1.8750 -1.8750 0 7.5000 -1.8750 -5.6250 3.7500 0 0 0 -1.8750 -1.8750 7.5000 3.7500 -5.6250 -1.8750 0 0 -1.8750 -5.6250 3.7500 7.5000 -1.8750 -1.8750 -1.8750 -1.8750 0 3.7500 -5.6250 -1.8750 7.5000 0 -1.8750 -5.6250 3.7500 0 -1.8750 -1.8750 0 7.50002. 基于四节点四边形单元的

22、有限元建模及分析（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)该函数计算单元的应力，输入弹性模量E，泊松比NU，厚度h，4个节点i，j，m，p，平面问题性质指示参数ID（1为平面应力，2为平面应变），单元的位移列阵u(1)，输出单元的应力，由于

23、它为常应力单元，则单元的应力分量为Sx，Sy，Sxy。（2）程序清单%Quad2D4Node%begin%function k=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%输入平面问题性质指示参数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)

24、+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)/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

25、*(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;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

26、;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;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%function stress=Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID)%该函数计算单元的应力%输入弹性模量E、泊松比NU和厚度h%输入平面问题

27、性质指示参数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*(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*

28、(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*(-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

29、 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 0 (1-NU)/2;elseif ID=2 D=(E/(1+NU)/(1-2*NU)*1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2;endstr1=D*B*u;str2=subs(str1,s,t,0,

30、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=0.5*U'*K*U-P'*U>>stress=Quad2D4Node_Stress(E

31、,NU,0,0,1,0,1,1,0,1,10(-3)*0;0;0.4815;0.1111;0.4815;-0.1111;0;0,1)（3)计算结果应变C = 1.0e+06 * 1.2500 0.5625 -0.8750 0.1875 -0.6250 -0.5625 0.2500 -0.1875 0.5625 1.2500 -0.1875 0.2500 -0.5625 -0.6250 0.1875 -0.8750 -0.8750 -0.1875 1.2500 -0.5625 0.2500 0.1875 -0.6250 0.5625 0.1875 0.2500 -0.5625 1.2500 -0.1875 -0.8750 0.5625 -0.6250 -0.6250 -0.5625 0.2500 -0.1875 1.2500 0.5625 -0.8750 0.1875 -0.5625 -0.6250 0.1875 -0.87