《有限元分析》课程作业

《有限元分析》课程作业《有限元分析》课程作业《有限元分析》课程作业《有限元分析》课程作业编制仅供参考审核批准生效日期地址：电话：传真:邮编:《有限元分析》课程作业任课教师：徐亚兰学生姓名：陈新杰学号：班级：1304012时间：2016-01-05

一、问题描述及分析问题：如图1所示，有一矩形平板，在右侧受到P=10KN/m的分布力，材料常数为：弹性模量；泊松比；板的厚度为t=；试按平面应力问题利用三角形与矩形单元分别计算各个节点位移及支座反力。P=10KN/m1m1P=10KN/m1m1m图1平面矩形结构的有限元分析分析：使用两种方案：一、基于3节点三角形单元的有限元建模，将矩形划分为两个3节点三角形单元；二、基于4节点矩形单元的有限元建模，使用一个4节点矩形单元。利用MATLAB软件计算出各要求量，再将两种方案的计算结果进行比较、分析、得出结论。二、有限元建模及分析1、基于3节点三角形单元的有限元建模及分析（1）结构的离散化与编号如图2所示，将平面矩形结构分为两个3节点三角形单元。单元①三个节点的编号为1，2，4，单元②三个节点的编号为3，4，2，各个节点的位置坐标为，各个节点的位移（分别沿方向和方向）为。X②①y1432X②①y1432图2方案一：使用两个3节点三角形单元（2）各单元的刚度矩阵及刚度方程a.单元的几何和节点描述单元①有6个节点位移自由度（DOF）。将所有节点上的位移组成一个列阵，记作；同样，将所有节点上的各个力也组成一个列阵，记作，则有同理，对于单元②，有b.单元的位移场描述对于单元①，设位移函数(1-1)由节点条件，在处，有(1-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)其中几何矩阵(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)，得到单元①的刚度阵计算得同理，得到单元②的刚度阵为将两个单元按节点位移所对应的位置进行组装，得到总体刚度矩阵为节点力系统的势能(计算结果在下面呈现)(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=k\pu=[将列排成了行]再计算支反力。在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力；先将上面得到的位移结果与边界条件的节点位移进行组合，得到整体的位移列阵U，再代回原整体刚度方程，计算出所有的节点力，按照位置关系找出对应的支反力。>>U=[0;0;u;0;0]U=0000[将列排成了行]>>P=KK*UP=-50050005000-500[将列排成了行]所以，节点1的支反力为，节点2的支反力为。根据已求得的位移和支反力计算系统的势能。>>A=*U'*KK*U-P'*UA=(5）结果分析上述支反力计算结果满足静力平衡，验证了以上求解过程及MATLAB算法的正确性。2、基于四节点四边形单元的有限元建模及分析（1）结构的离散化与编号如图3所示一个4节点矩形单元，单元的节点位移共有8个自由度（DOF）。节点编号为1,2,3,4，各自的位置坐标为，各个节点的位移（分别沿方向和方向）为。Xy①4321Xy①4321图3方案二：使用一个4节点矩形单元（2）局部坐标系下单元的描述a.单元的几何和节点描述采用无量纲坐标其中。则单元四个节点的几何位置为将所有节点上的位移组成一个列阵，记作q；同样，将所有节点上的各个力也组成一个列阵，记作F，则有b.单元的位移场描述设位移函数为由节点条件，在处，有将位移试函数代入节点条件中，求出待定系数和，再代入位移函数中，整理后得其中如以无量纲坐标来表达，可写成将带入上式得到形状函数矩阵写成矩阵形式，有c.单元的应变场描述单元应变为其中几何矩阵为d.单元的应力场描述应力表达式为其中，应力函数矩阵。e.单元的势能表达以上已将单元的三大基本变量用基于节点的位移列阵来表示；将其代入单元势能表达式中，有，其中为4节点矩形单元的刚度矩阵，即其中，t为薄板的厚度，，上式的各个字块矩阵为f.单元刚度阵及刚度方程单元刚度阵在上面已经列出。将单元的势能对节点位移取一阶极值，可得到单元的刚度方程(4)边界条件的处理及方程求解处理方法与3节点三角形单元一致，利用上述求解程序具有的可移植性，简化了求解过程。>>k=K(3:6,3:6);>>p=[500;0;500;0];>>u=k\pu=*[将列排成了行]再计算支反力。同样注意按照位置关系找出对应的支反力。>>U=[0;0;u;0;0]U=*0000[将列排成了行]>>P=K*UP=-50050005000-500[将列排成了行]所以，节点1的支反力为，节点2的支反力为。根据已求得的位移和支反力计算系统的势能。>>A=*U'*K*U-P'*UA=(5）结果分析基于4节点矩形单元计算出的势能小于基于3节点三角形单元计算出的结果，若将该系统分为更多的单元，计算精度也将提高。3.两种方案的比较与分析从以上计算可以看出，用三角形单元计算时，由于形函数是完全一次式，因而其应变场在单元内均为常数；而对于四边形单元，其形函数带有二次式，计算得到的应变场和应力场都是坐标的一次函数，但不是完全的一次函数，对提高精度有一定作用；根据最小势能原理，势能越小，则整体计算精度越高，比较两种单元计算得到的系统势能，可看出，在相同的节点自由度情况下，矩形单元的计算精度要比三角形单元高。三、基于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)该函数进行单元刚度矩阵的组装，输入单元刚度矩阵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,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=[betai0betaj0betam0;0gammai0gammaj0gammam;gammaibetaigammajbetajgammambetam]/(2*A);ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];endk=t*A*B'*D*B;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionz=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;forn1=1:6forn2=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）%输出单元的应力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=[batai0bataj0betam0;0gammai0gammaj0gammam;gammaibetaigammajbetajgammambetam]/(2*A);ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;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,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=k\p>>U=[0;0;u;0;0]>>P=KK*U>>A=*U'*KK*U-P'*U（3)计算结果①应变CC=+05*0000000000000000②位移U=0000③节点力P=00其中，节点1的支反力为，节点2的支反力为。④势能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，平面问题性质指示参数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)，输出单元的应力，由于它为常应力单元，则单元的应力分量为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%输入平面问题性质指示参数ID（1位平面应力，2为平面应变）%输入单元刚度矩阵k（8*8)%-----------------------------------------symsst;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)/40;0c*(s-1)/4-d*(t-1)/4;c*(-1+s)/4-d*(t-1)/4a*(t-1)/4-b*(s-1)/4];B2=[a*(1-t)/4-b*(-1-s)/40;0c*(-1-s)/4-d*(1-t)/4;c*(-1-s)/4-d*(1-t)/4a*(1-t)/4-b*(-1-s)/4];B3=[a*(t+1)/4-b*(s+1)/40;0c*(s+1)/4-d*(t+1)/4;c*(s+1)/4-d*(t+1)/4a*(t+1)/4-b*(s+1)/4];B4=[a*(-t-1)/4-b*(1-s)/40;0c*(1-s)/4-d*(-t-1)/4;c*(1-s)/4-d*(-t-1)/4a*(-t-1)/4-b*(1-s)/4];Bfirst=[B1B2B3B4];Jfirst=[01-tt-ss-1;t-10s+1-s-t;s-t-s-10t+1;1-ss+t-t-10];J=[xixjxmxp]*[yi;yj;ym;yp]/8;B=Bfirst/J;ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;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%%%%%%%%%%%%%%%%%%%%%%%%%%%%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%--------------------------------------------------symst;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)/40;0c*(s-1)/4-d*(t-1)/4;c*(-1+s)/4-d*(t-1)/4a*(t-1)/4-b*(s-1)/4];B2=[a*(1-t)/4-b*(-1-s)/40;0c*(-1-S)/4-d*(1-t)/4;c*(-1-s)/4-d*(1-t)/4a*(1-t)/4-b*(-1-s)/4];B3=[a*(t+1)/4-b*(s+1)/40;0c*(s+1)/4-d*(t+1)/4;c*(s+1)/4-d*(t+1)/4a*(t+1)/4-b*(s+1)/4];B4=[a*(-t-1)/4-b*(1-s)/40;0c*(1-s)/4-d*(-t-1)/4;c*(1-s)/4-d*(-t-1)/4a*(-t-1)/4-b*(1-s)/4];Bfirst=[B1B2B3B4];Jfirst=[01-tt-ss-1;t-10s+1-s-t;s-t-s-10t+1;1-ss+t-t-10];J=[xixjxmxp]*[yi;yj;ym;yp]/8;B=Bfirst/J;ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;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=Quad