有限元-计算结构力学-大作业(共27页).docx

精选优质文档-倾情为你奉上SHANGHAI JIAO TONG UNIVERSITY平面应力问题解的Matlab实现姓 名: heiya168 学 号: 帆哥 班 级： 指导老师： 目录专心-专注-专业1 绪论有限元方法(finite element method)，是求取复杂微分方程近似解的一种非常有效的工具，是现代数字化科技的一种重要基础性原理。将它用于在科学研究中，可成为探究物质客观规律的先进手段。将它应用于工程技术中，可成为工程设计和分析的可靠工具。弹性体在载荷作用下，其基本方程可写成以下的三类方程和两种边界条件。平衡方程应力与外载荷的关系；几何方程应变位移关系；物理方程应力应变关系；力的边界条件；几何边界条件。应用最小位能原理，并利用上述关系，最终建立由刚度方程，节点位移和等效节点载荷所构成的求解方程。带入边界条件求解方程，就可以得出弹性力学问题的一般性解答。本次大作业基于有限元方法的基本原理，使用Matlab这一平台，针对平面应力问题，采用四节点四边形单元编写了求解单元节点位移的程序。主要内容包括：1）介绍有限元的基本原理；2）编程基本思路及流程介绍；3）程序原理及说明； 4）具体算例 这四个部分。2 平面问题的四节点四边形单元2.1 单元的构造（1） 单元的几何和节点描述平面4节点矩形单元如图4-6所示，单元的节点位移共有8个自由度。节点的编号为1、2、3、4，各自的位置坐标为(xi, yi), i=1,2,3,4，各个节点的位移(分别沿x方向和y方向)为(ui,vi), i=1,2,3,4。图2.1 平面4节点矩形单元若采用无量纲坐标 （2.1）则单元4个节点的几何位置为 （2.2） 将所有节点上的位移组成一个列阵，记作；同样，将所有节点上的各个力也组成一个列阵，记作，那么 （2.3） 若该单元承受分布外载，可以将其等效到节点上，也可以表示为如式(2.3)所示的节点力。利用函数插值、几何方程、物理方程以及势能计算公式，可以将单元的所有力学参量用节点位移列阵及相关的插值函数来表示；下面进行具体的推导。（2） 单元位移场的表达从图2-1可以看出，节点条件共有8个，即x方向4个，y方向4个，因此，x和y方向的位移场可以各有4个待定系数，即取以下多项式作为单元的位移场模式 （2.4）它们是具有完全一次项的非完全二次项，以上两式中右端的第四项是考虑到x方向和y方向的对称性而取的，除此外xy项还有个重要特点，就是“双线性”，当x或y不变时，沿y或x方向位移函数呈线性变化，这与前面的线性项最为相容，而或项是二次曲线变化的。因此，未选或项。 由节点条件，在x=xi，y=yi处，有 （2.5）将式(2.5)代入式(2.4)中，可以求解出待定系数a0，a3和b0，b3，然后代回式(4-52)中，经整理后有 （2.6）其中 （2.7）如以无量纲坐标系(2.1)来表达，式(2.7)可以写成 （2.8）将式(2.6)写成矩阵形式，有（2.9）其中，为该单元的形状函数矩阵。（3） 单元应变场的表达由弹性力学平面问题的几何方程(矩阵形式)，有单元应变的表达 （2.10）其中几何矩阵B(x,y)为= （2.11）式(4-59)中的子矩阵为 （2.12）（4） 单元应力场的表达由弹性力学中平面问题的物理方程，可得到单元的应力表达式 （2.13）（5） 单元应力场的表达以上已将单元的三大基本变量用基于节点位移列阵来进行表达，见式(2.9)、式(2.10)及式(2.13)；将其代入单元的势能表达式中，有 （2.14）其中是4节点矩形单元的刚度矩阵。将单元的势能对节点位移取一阶极值，可得到单元的刚度方程 （2.15）2.2 等参变换由前面的单元构造过程可以看出，一个单元的关键就是计算它的刚度矩阵，而由刚度矩阵的构成可知要实现两个坐标系中单元刚度矩阵的变换，必须计算两个坐标系之间的三种映射关系：坐标映射 （2.16）偏导数映射 （2.17）面积映射 （2.18）图2.2矩形单元映射为任意四边形单元（1） 两个坐标系之间的函数映射设如图4-17所示的两个坐标系的坐标映射关系为 （2.19）x和y方向上可以分别写出各包含有4个待定系数的多项式，即 （2.20）其中待定系数a0，a3和b0，b3可由节点映射条件(2.19)来唯一确定。 对照前面4节点矩形单元的单元位移函数式(2.5)，映射函数式(2.20)具有完全相同的形式，同样，将求出的待定系数再代回式(2.20)中，重写该式为 （2.21）其中 （2.22）（2.23）这就可以实现两个坐标系间的映射。 （2）两个坐标系之间的偏导数映射对物理坐标系(x,y)中的任意一个函数(x,y)，求它的偏导数，有 （2.24）则偏导数的变换关系为 （2.25）写成矩阵形式，有 （2.26）其中 （2.27）两个坐标系的偏导数映射关系 （2.28）（3）两个坐标系之间的面积元映射如图2-2所示，在物理坐标系(x, y)中，由d和d所围成的微小平行四边形，其面积为 （2.29）由于d和d在物理坐标系(x, y)中的分量为 （2.30）其中i和j分别为物理坐标系(x, y)中的x方向和y方向的单位向量。由(2.29)式，则有面积微元的变换计算 （2.31）2.3 边界条件的处理置“1”法设边界条件为第r个自由度的位移为零，即；可置整体刚度矩阵中所对应对角元素位置的，而该行和该列的其它元素为零，即，同时也置对应的载荷元素 则，则（2.32）进行以上设置后，这时方程(2.32)应等价于原方程加上了边界条件，下面考察这种等价性，就(2.32)中的第r行，有 （2.33）由于置，则有即为所要求的位移边界条件。而式(5-53)中除第r行外，其它各行在对应于r列的位置上都置了“0”，这相当于考虑了 的影响，除此之外其余各项的影响不变；这恰好就是原方程加上了边界条件的影响。3 有限元分析流程3.1 程序原理和流程该程序的特点如下：问题类型：可用于弹性力学平面应力问题和平面应变问题的分析单元类型：四节点四边形单元 材料性质：单一的均匀的弹性材料节点信息：节点信息文件node.txt单元信息：单元信息文件element.txt载荷信息：等效节点力文件force.txt约束信息：节点约束信息文件constrain.txt结果文件：输出节点位移文件node_displace.txtelement.txtforce.txt程序流程图如下：输入参数constrain.txtnode.txt计算总刚规模function K=Stiffnessfunction z =Assembly形成总刚度阵形成总体载荷向量处理边界约束求解Ka=Pnode_displace.txt输出节点位移图3.1程序流程图3.2 使用的函数Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp)：计算单元的刚度矩阵(具体说明见附录)Assembly(KK,k,i,j,m)：进行单元刚度矩阵的组装（具体说明见附录）3.3 文件管理主程序源文件： FEM.m 输入参数文件： node.txt (节点信息文件) element.txt(单元信息文件) constrain.txt(位移边界条件文件) force.txt(载荷状况文件)输出位移文件：node_displace.txt（节点位移文件）3.4 数据文件格式4 算例开方孔的矩形板拉伸分析4.1 问题的具体参数与载荷 (a) 对边拉伸的带孔矩形薄板 (b) 1/4模型的单元划分 图4-1 受均匀荷载的作用的带孔矩形薄板薄板及有限元分析模型如图4-1(a)所示的矩形带孔薄板对边受均匀荷载的作用，该结构在边界上受正向分布拉力0.3MPa。若取板厚0.02m，弹性模量2.1E11，泊松比=0.3，采用FEM.m程序，和ANSYS分别分析该问题，最后给出各节点位移值。4.2 Matlab程序计算（1）结构的离散化与编号该薄板的荷载和几何形状关于x轴和轴对称，故可只取结构的1/4作为计算模型。将此模型化分为四个全等的直角三角形单元，单元编号和节点编号如图4-1(b)所示。（2）输入节点信息node.txt（3）输入单元信息element.txt（4）输入载荷信息force.txt（5）输入载荷信息constrain.txt（6）运行主函数FEM.m，输入E=2.1E11，NU-0.3，h=0.02。（7）输出各节点位移4.3 ANSYS建模计算选用plane42板单元，建立如下图所示模型，在2、3节点加载水平方向的载荷10000N.图4-2 受均匀荷载的作用的带孔矩形薄板薄板有限元分析模型分析计算所得的最后结果为：4.4 误差分析节点位移2X2Y3X3Y4X4Y2SUM3SUM4SUMMatlab7.65E-062.53E-077.33E-06-1.43E-061.87E-062.025E-077.66E-067.469E-061.88E-06ANSYS7.54E-06-8.00E-077.81E-06-2.58E-061.68E-062.45E-077.58E-068.22E-061.70E-06误差1.02%9.18%9.94%表4-1 Matlab程序与Ansys结果对比表误差范围在10%以内，还是比较大的，导致误差的原因可能有以下几个方面：1）刚度矩阵计算时，使用的积分方法不同，本程序中使用的是两点高斯积分。2）刚度矩阵求解是，使用的算法不同，本程序使用的是高斯消去法。3）约束处理时使用的方法不同，本程序使用的是置一法。5 总结用了两周的时间，完成了这次大作业。在做这次大作业之前，我对于有限元算法的了解还很肤浅，只是模模糊糊的能够记得刚度阵建立的过程，矩阵表示的公式，对于其中具体的物理意义，各个量之间相互关系的表示不太清楚，在做完这次大作业后，对于四边形单元刚度阵的建立过程有了清晰深入的认识和掌握。关于等参变化的认识，是这次大作业中的难点之一。对于和位移差值函数相同的坐标差值函数，如何理解其中“坐标差值”物理意义，在差值完成，刚度矩阵建立后，刚度矩阵的位移模式是局部坐标系中的，还是整体坐标系中的位移，这个问题思考了很久，也和很多同学进行了争论和探讨。最终得到了一致的结论，收获很多。在完成大作业之前，我对于Matlab编程还是一窍不通，在编程的过程中，很多好的想法无法用计算机语言来实现是一件很痛苦的事情，但是在同学的帮助下和参考资料的借鉴下，完成了这次编程。其中，刚度矩阵和组装矩阵是自己独立完成的，主程序借鉴了清华大学曾攀教授编写的有限元分析基础教程中平面三节点单元的主程序。写这个总结的时候，也就是大作业完成的时候，其中还有很多不足的地方有待改进，比如两点高斯积分精确度较差的问题，比如无法直接输入面载荷和体载荷的问题。参考文献1 王勖成.有限单元法M.清华大学出版社.2 曾攀.有限元分析基础教程M.清华大学出版社.3 陈杰.MATLAB 宝典M.电子工业出版社附录function K=Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp)%该函数计算单元的刚度矩阵%输入弹性模量E，泊松比NU，厚度h%输入4个节点i、j、m、p的坐标xi,yi,xj,yj,xm,ym,xp,yp%输出单元刚度矩阵k(8X8)%-syms s t;N1=1/4*(1-s)*(1-t);N2=1/4*(1+s)*(1-t);N3=1/4*(1+s)*(1+t);N4=1/4*(1-s)*(1+t);%坐标变换函数，位移插值函数X=N1*xi+N2*xj+N3*xm+N4*xp; Y=N1*yi+N2*yj+N3*ym+N4*yp;%坐标变换J11=diff(X,s);J12=diff(Y,s);J21=diff(X,t);J22=diff(Y,t);J=J11 J12;J21 J22;%J矩阵的建立H=1/det(J)*J22 -J12 0 0;0 0 -J21 J11;-J21 J11 J22 -J21;%B=HQ,H矩阵建立q11=diff(N1,s);q21=diff(N1,t);q32=diff(N1,s);q42=diff(N1,t);Q1=q11 0;q21 0;0 q32;0 q42;q13=diff(N2,s);q23=diff(N2,t);q34=diff(N2,s);q44=diff(N2,t);Q2=q13 0;q23 0;0 q34;0 q44;q15=diff(N3,s);q25=diff(N3,t);q36=diff(N3,s);q46=diff(N3,t);Q3=q15 0;q25 0;0 q36;0 q46;q17=diff(N4,s);q27=diff(N4,t);q38=diff(N4,s);q48=diff(N4,t);Q4=q17 0;q27 0;0 q38;0 q48;Q=Q1 Q2 Q3 Q4;%Q矩阵的建立B=H*Q;D=(E/(1-NU2)*1 NU 0;NU 1 0;0 0 (1-NU)/2;A=transpose(B)*D*B;k=A*det(J)*h;s=1/3*30.5;t=1/3*30.5;k1=eval(k);s=-1/3*30.5;t=1/3*30.5;k2=eval(k);s=1/3*30.5;t=-1/3*30.5;k3=eval(k);s=-1/3*30.5;t=-1/3*30.5;k4=eval(k);K=k1+k2+k3+k4;function z =Assembly(KK,k,i,j,m,p)%该函数进行单元刚度矩阵的组装%输入单元刚度矩阵k，单元的节点编号i、j、m、p%输出整体刚度矩阵KK%-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;DOF(7)=2*p-1;DOF(8)=2*p;for n1=1:8 for n2=1:8 KK(DOF(n1),DOF(n2)= KK(DOF(n1),DOF(n2)+k(n1,n2); endendz=KK;% FEM.m% main program begin %该程序为计算平面应力问题的主程序%输入节点，单元，载荷，约束信息，调用Stiffness和Assembly%求得总刚，并解得节点位移%clear; %读入节点，单元，载荷，约束信息load ('F:FEMnode.txt')load ('F:FEMelement.txt')load ('F:FEMconstrain.txt') load ('F:FEMforce.txt') %确定节点、单元个数 nnode,ntmp=size(node); nelem,etmp=size(element); nforce,ftmp=size(force); nconstrain,ctmp=size(constrain);%预先设定总体刚度矩阵、节点力向量、节点位移向量 KKG=zeros(2*nnode); FFG=zeros(2*nnode,1); UUG=zeros(2*nnode,1); k=zeros(8,8); % 单元刚度矩阵 %给出相应材料及计算参数 E=input('弹性模量E='); NU=input('泊松比NU='); t=input('板厚='); %单元循环形成总体刚度矩阵 for i=1:nelem K=Stiffness( E,NU,t,node(element(i,1),2),node(element(i,1),3) ,node(element(i,2),2),node(element(i,2),3) ,node(element(i,3),2),node(element(i,3),3),node(element(i,4),2),node(element(i,4),3); KKG =Assembly(KKG,K,element(i,1),element(i,2),element(i,3),element(i,4); end %载荷的处理 for i=1:nforce m=force(i,1); n=force(i,2); FFG(2*(m-1)+n)= force(i,3); end %采用置”1”法处理边界条件 for i=1:nconstrain m=constrain(i,1); n=constrain(i,2); UUG(2*(m-1)+n)=constrain(i,3); KKG(2*(m-1)+n,:)=0; KKG(:,2*(m-1)+n)=0; KKG(2*(m-1)+n,2*(m-1)+n)=1; FFG(2*(m-1)+n)=0; end % 求解节点位移 UUG=KKGFFG; % 输出节点位移值fid=fopen('F:FEMnode_displace.txt','w');fprintf(fid,'n%sn','- NODE DISPLACEMENT -');fprintf(fid,'n%sn','NODE X-disp Y-disp ');for i=1:nnodefprintf(fid,'%d %18.10f %18.10fn',node(i,1),UUG(2*(i-1)+1),UUG(2*(i-1)+2);end