有限元程序设计报告报告

苏世宏有限元课程设计有限元程序设计报告课程名称:有限元程序设计指导教师：张 亮学 校：重庆大学专 业：工程力学01班姓 名：苏世宏学 号：2015年7月8日有限元程序设计报告一、前言有限元方法（the Finite Element Method）是起源与上个世纪 50、60 年代，基于弹性力 学变分原理的一种近似计算方法，也是当今工程分析中获得最广泛应用的数值计算方法。由 于它的通用性和有效性，受到工程技术界的高度重视。伴随着计算机科学和技术的快速发展， 现已成为计算机辅助设计（CAD）和计算机辅助制造（CAM）的重要组成部分。有限元程 序系统通常包括前处理、有限元程序本体和后处理三部分。前处理包括几何实体模型的建立、 材料参数的赋值、位移边界条件的定义、载荷的定义、分析问题类型的定义、单元类型的选 择和网格的划分等。（分析问题类型如静力分析、动力特性分析、动力响应、温度场分析、 电磁场分析、流体动力学分析等）有限元程序本体是有限元程序系统的核心部分，其功能是 实现各种问题的计算。后处理则是将计算结果用图形、曲线和表格的形式表达。（通常包括 结构的变形图、应力、应变分布云图等）本课程设计则是针对有限元程序本体，参照教学程序（FEATP），编写简单的有限元程 序以计算简单的平面应力、平面应变和轴对称问题，并将其结果与有限元商用软件（ANSYS） 的计算结果，以及问题的理论值进行比较，从而验证程序以及问题模型建立的正确性。1 设计目的1） 通过编写简单的有限元程序熟悉用有限元方法解决实际问题的基本步骤和过程，体会这 种方法的处理手段。2） 在 Visual Fortran 中编写程序，熟悉并巩固 Fortran 语言的语法、算法，学习程序的调试 方法，并体会其在执行某个具体算例时，文件的输入、输出以及程序的执行过程。2 设计内容1）以教学程序（FEATP）为参照，编写程序，计算简单的平面应力（Plane Stress），平面应 变(Plane Strain)问题，验证程序的正确性。2）在具体的算例中，对同一问题，在程序和 ANSYS 中采用不同的单元和网格划分方式， 将其结果与理论值进行对比，体会不同的单元和网格划分对问题解的影响，从而判断模 型的正确性和合理性。3）总结在编写程序和算例中遇到的问题和解决方法，写出自己的心得体会。二、弹性力学平面问题有限元方法的基本公式平面问题1三角形单元第 44 页(1) 位移 ui vi uNiu = = 0N j0Nm0u jv 0Ni0N j0Nm v j ai um vm i= ININ jINm a j ij= NNam mN a e = Na e 其中，形函数 ai ,bi , ci是取决于节点坐标的常数应变矩阵 B 的分块子矩阵是 3 节点单元的应变矩阵是1(2) 形成单元的刚度矩阵和等效节点载荷列阵K e =BT DBdV = BT DBtAVePe = Pe + Pe + Pe + Pefss 0e0 BT(3) 集成结构的刚度矩阵和等效节点载荷列阵K = K e= VDBdVeP = Pfe0+ Ps Psee+ P0+ PFfS s 0e 0F=(Pe + Pe Pee+ Pe ) + P(4) 引入强制边界条件（消除 K 的奇异性）(5) 求解有限元解方程，得到节点位移 aKa = P(6) 计算单元应变和应力 e x yme = ege = Lu = LNa = L NiN jN aemi= Bxy BjB ae = Ba es x s =s yt = De = DBae = Sae xy ijmijm应力矩阵 S = DB = DBBB = SSS S 的分块子矩阵为bin 0 ciSi = DBi =E02n 0bic（i,j,m）i2(1 -n 0 ) A 1 -n 01 -n 02ci2bi 对于平面应力问题E0 = E,n 0 =n对于平面应变问题3 四边形单元 用相同的推导方法可以得到四边形单元位移、应变、应力的有限元表达格式，它具有和三角形单元相似的表达形式，这里就不一一列举了。三、有限元程序设计1程序功能本程序是在教学程序（FEATP）的基础上修改，删减而成的，主要应用于解决各项同性 的弹性力学二维问题，平面应力问题，平面应变问题和轴对称问题。对于原程序中的动态响 应问题（DYNAM）和动力特性分析问题（EIGENVALUE PROBLEM）在本程序中将不涉及。（1）问题类型平面应力问题（MPROB=1）平面应变问题（MPROB=2）轴对称问题（MPROB=3）（2）单元类型36 节点三角形单元（NODE=3 或 6）48 节点四边形单元（NODE=4 或 8）9 节点四边形单元（NODE=9）（3） 求解类型 静力平衡分析：等带宽三角分解法（MOSLV=1）2程序框图（1）程序总体框图输入离散模型数据计算单元刚度阵组集结构刚度矩阵单元循环形成 K计算单元等效结点载荷组集结构结点载荷列阵形成 P引入位移边界条件消除 K 的奇异求解线性方程组求解 Ka=P，得结点位移 a其它辅助计算计算应力、应变等输出结果结束（2）程序调用框图3.输入文件变量名说明输入（1） MND：计算模型的各类单元中最多的节点数； NUMEL：计算模型的单元总数； NUMPT：计算模型的节点总数； MBAND：半带宽（包括主对角元素）。 输入（2）NFIX：有位移约束的节点数； NPC：等效载荷作用的节点数； MPROB：问题类型； MSOLV：分析类型。输入（3）NMATI：材料类型数； GRAV：重力加速度值；若不考虑重力，则输入 0.0； MTYPE：输入控制参数；MTYPE0 输出全部计算结果；MTYPE1 输出除积分点应力以外的全部计算结果；MTYPE2 输出除总体质量矩阵、总体刚度矩阵和总体载荷向量以外的 全部计算结果；MTYPE3 输出除积分点应力以及总体质量矩阵、总体刚度矩阵和总体 载荷向量以外的全部计算结果。输入（4）（4）是节点坐标信息的输入，即（II,(VCOOD(I,J),J=1,2,I=1,NUMPT)。其中II：模型中的节点号，从 1 至 NUMPT 依次按行输入；VCOOD(I,1)：II 节点处的 x 向坐标。 VCOOD(I,2)：II 节点处的 y 向坐标。 输入（5）（5）是单元信息的输入，即（II,(IELEM(I,J),J=1,4+MND),I=1,NUMEL）。其中II：模型中的单元号，从 1 至 NUMEL 依次按行输入。IELEM(I,1)：II 单元的节点数。IELEM(I,2)：II 单元的材料类型号。IELEM(I,3)：II 单元沿 x 方向的高斯积分点数，对于三角形单元，则是 Hammer 积分点数 IELEM(I,4)：II 单元沿 y 方向的高斯积分点数。对于三角形单元，填 1。 IELEM(I,5)：IELEM(I,MND)：依次是 II 单元的局部编号所对应的总体编号。对于 IELEM(I,1) 小于 IELEM（I,MND）的情况，在相应的位置上填 0。输入（6）（6）是位移约束信息的输入，即（II,(IFIXD(I,J),J=1,NF+1),(VFIXD(I,J),J=1,NF),I=1,NFIX）。其中 II：约束信息号，从 1 至 NFIX 依次按行输入。 输入（7）7）是等效节点载荷信息的输入，即（II,(ILOAD(I,J),J=1,NF+1),(VLOAD(I,J),J=1,NF,I=1,NPC)。其中，II：等效节点载荷号。ILOAD(I,1)：第 II 个等效节点载荷作用的节点号。ILOAD(I,2)ILOAD(I,NF+1)：第 II 个等效节点载荷所用节点的自由度开关，1 表示有 载荷作用，0 表示没有载荷作用。VLOAD(I,1)VLOAD(I,NF):第 II 个等效节点载荷作用于节点的自由度方向的载荷值大 小。对于平面问题和轴对称问题分别代表 X,Y 方向的载荷值；对于 Mindlin 板分别代表q x ,q y和 W 方向的载荷值。 输入（8）（8）是材料类型和几何信息的输入，即（II，（VMATI(I,J),J=1,4）,I=1,NMATI） II：材料类型号，从 1 至 NMATI 依次按行输入；VMATI(I,1)：II 号材料的弹性模量（E）。 VMATI(I,2)：II 号材料的泊松比（v）。 VMATI(I,3)：II 号材料的质量密度（dens）。 VMATI(I,4)：II 号材料处板的厚度（th）。四、数值算例1 算例一问题阐述：请采用4节点四边形等参单元对图1所示的无量纲L型框架结构进行有限元分析。材料杨氏模量和泊松比分别为、。(1) 绘制出A点水平位移随均布剪力P取值变化（05000）的关系曲线及结构在3个典型载荷下的变形图；（建议长边采用40个单元，短边采用4个单元对结构进行离散）(2) 结合弹性力学小变形、线弹性假设，谈谈你对有限元分析结果的认识。解：（1）程序编辑和程序调用以及演算过程演示：（i）整体程序（ii）子程序（stineffness，stress，input，output.）（iii）网格的划分：使用ABAQUS对结构进行离散，即用ABAQUS对结构进行网格划分，获得网格节点和单元的信息，然后形成“.txt”文件，保存到相应的文件目录里。网格划分步骤如下：(a)安装ABAQUS并且启动ABAQUS （b）按照ABAQUS建模步骤，一步一步进行。分别为：创建部件 创建材料和截面属性 定义装配件 定义边界条件和载荷 划分网格 提交分析作业 后处理 退出ABAQUS/CAE. （c）利用Ultraedit软件，打开job里面的“.inp”文件，获取节点信息，转换成“.txt”问价输出。（d）利用所得的txt文件，回到MATLAB进行数值计算。（iiii）程序运行结果变形结果数值计算结果（本题只以以此计算的结果为例，其余结果改变外载荷后重新计算即可）本题P在不同取值下，A点位移U的变化如下表所示，其中放大倍数都是1.A点水平位移U(m）施加均布剪力P(KN)6.48E-015001.30E+0010001.94E+0015002.59E+0020003.24E+0025003.89E+0030004.53E+0035005.18E+0040005.83E+0045006.48E+005000绘制出A点水平位移随均布剪力P取值变化（05000）的关系曲线及结构在3个典型载荷下的变形图如下图所示：（2） 结果分析：从程序（FEATP）的运算结果中不难看出，各个节点的位移值（NODEL DISPLCEMENT） 和实际情况符合得很好，得到就是问题的精确解。变形与外在施加呈线性关系。数据分析：自由端的水平位移（节点3,4,47,48,49三点水平位移）相差不大，只有微小的变动，几乎可以忽略。说明节点上的变形几乎是一致的。固定端的水平、竖直位移为零，说明该结构在固定端的数值计算结果与input输入文件里规定的一样，结果可靠。对有限元分析结果的认识：就精度而言，从程序算例结果可以看出，在同样的模型下，位移的值总是与理论值 相等， 而应力的值与理论值则有一定误差。这说明位移的精度要高于应力计算的精 度。而这一点 与有限元法理论是吻合的。（因为他们之间存在间接的一阶倒数关系，经过求导后得到的应变 的精度较位移较低）从应力输出文件中，还不难看出，单元高斯积分点的应力值与理论值相等，而节点处 的应力值 却较理论值有一定误差。这也就验证了有限元理论中所说的：高斯积分点处 的精度最高，而 节点处的精度不理想。数值计算结果与理论值之间误差的减小可以通过细化网格或者提高差值函数阶数完成，一般提高差值函数阶数是最直接的方法，不用重新划网格，节省计算量，提高计算效率。2 算例二问题阐述： 图2所示为一平面悬臂梁结构。结构几何参数为，；材料杨氏模量 和泊松比分别为，；均布载荷。请分别采用 4节点四边形等参单元分析结构在(a)、(b)两种载荷作用下的力学响应。(1) 画出A点的竖向位移随结构总自由度数目变化的曲线，并将有限元分析结果与问题的解析解1进行对比分析。（如果位移误差大于5%，则需通过细化网格来提高有限元解的精度。建议网格划分从疏到密：、和）；(2) 位移解收敛后，在梁中性层的上、下侧任取两个对称位置的高斯积分点及节点，将高斯点应力值和节点应力值分别与解析解进行比较，结合程序分析节点应力精度比高斯积分点应力精度低的原因。图2 悬臂梁结构附杆解析解：悬臂梁解析解：轴线挠度 水平应力分量 （其中，）解：（一）问题一的具体解答（1）理论解：杆解析解：悬臂梁解析解：轴线挠度 水平应力分量 （其中，）(2程序调试和数据设置:(i)主程序不变，将input子程序里面的参数设置和节点、网格信息进行修改，修改以后要与本问题相符。（ii)将外载荷子程序EffecLoad进行外载荷的修改。（2） 各网格情况下的数据和变形结果（i）网格为时变形：在外在a情况下的变形在外在b情况下的变形当上述变形在载荷大小不变，方向相反的变形如下图，可以看到结构的变形范围数据结果：外在在a情况下的数据计算结果载荷在a情况下的水平位移（u）外在在b情况下的数据计算结果载荷在b情况下的扰度（v）结果分析：a情况下的误差（网格划分满足) b情况下的误差分析(网格划分不满足，细化网格重新计算）接下来细化网格重新计算（ii）网格为时由于在a情况下的位移，在10*1网格时已经满足，所以接下来对b情况在不同网阁下的精度进行比较。变形：在外在b情况下的变形a情况下的位移（u)和扰度（v）b情况下的位移（u）和扰度（v)结果分析：a情况下的误差（网格划分满足)b情况下的误差分析(不满足，细化网格重新计算）(iii)网格为时在外在b情况下的变形a情况下的位移（u)和扰度（v）b情况下的位移（u）和扰度（v)结果分析：a情况下的误差（网格划分满足)（网格划分满足)b情况下的误差分析(不满足，细化网格重新计算）(iiii)网格为时在外在b情况下的变形a情况下的位移（u）和扰度（v)b情况下的位移（u）和扰度（v)结果分析：b情况下的误差分析(误差小于允许误差，网格划分满足精度需求，计算完毕）（3） A点的竖向位移随结构总自由度数目变化的曲线（4） 结果分析：从上面每划分一次网格就进行一次运算，将所得数值结果与理论解进行误差分析来看，每进行一次网格细化，计算精度提高一次。a情况在10自由度下误差已在允许范围之内；b情况的误差随自由度逐次增不断精度增提高，最终b情况的误差为4.5%，允许误差为5%，在允许范围之内。由上述结果可以看到，当在a情况下除了水平位移外，竖向位移数值虽然很小，但是数值计算结果显示其存在，然而真实情况确是“无竖直位移（扰度）”，这就是数值计算的误差。在工程实际中，其竖直变形很小，几乎可以忽略不计。同理，b情况下的水平位移也非常小，几乎可以忽略不计。由上面v扰度随自由度变化曲线可以看出，其变化是线型的。（2） 问题2的解答（1） 应力理论解取B,C点作为节点。两点坐标B(-L/2，H/2), C(-L/2，-H/2) B点应力理论解理论应力 C点应力理论解（2） 所有节点应力数值计算结果：（3） 所选取的对称点在网格上位91和99节点，其数值应力为 91节点应力如下 99节点应力如下 对上述数值结果分析可得：对称位置上的数值应力是对称的，与实际情况符合的。误差分析： B点误差分析 C点误差分析 （4） 节点应力精度为何比高斯积分点应力精度低。 在有限元计算中，数值积分点就是高斯积分点，高斯积分点是使单元相对精确积分的，然而节点应力却不是，它只是单元每节点上的应力。结合上述程序来看，应力计算过程中主要有两次近似。第一次是高斯积分的近似（高斯积分是在高斯积分点的积分），虽然精度高，但然而还是一种近似。 第二次是节点应力的近似，是把高斯积分的应力再一次近似为节点应力，通过两次近似后得到节点应力，最后误差累计，造成节点应力精度比高斯积分点应力精度低 。 （a） 高斯积分的近似 for ie=1:Nelem x1=Nod(Ele(ie,1),1);y1=Nod(Ele(ie,1),2); x2=Nod(Ele(ie,2),1);y2=Nod(Ele(ie,2),2); x3=Nod(Ele(ie,3),1);y3=Nod(Ele(ie,3),2); x4=Nod(Ele(ie,4),1);y4=Nod(Ele(ie,4),2); for i=1:Nodes Ue(2*i-1:2*i)=U(2*Ele(ie,i)-1:2*Ele(ie,i); end for IP=1:4 if(IP=1) ksi=ksi2;eta=eta2; else if(IP=2) ksi=ksi1;eta=eta2; else if(IP=3) ksi=ksi1;eta=eta1; else ksi=ksi2;eta=eta1; end end end % J11=(1/4)*(-(1-eta)*x1+(1-eta)*x2+(1+eta)*x3-(1+eta)*x4); J12=(1/4)*(-(1-eta)*y1+(1-eta)*y2+(1+eta)*y3-(1+eta)*y4); J21=(1/4)*(-(1-ksi)*x1-(1+ksi)*x2+(1+ksi)*x3+(1-ksi)*x4); J22=(1/4)*(-(1-ksi)*y1-(1+ksi)*y2+(1+ksi)*y3+(1-ksi)*y4); J=J11,J12;J21,J22; A=(1/det(J)*J22,-J12,0,0;0,0,-J21,J11;-J21,J11,J22,-J12; G1=(1/4)*-(1-eta),0,(1-eta),0,(1+eta),0,-(1+eta),0; G2=(1/4)*-(1-ksi),0,-(1+ksi),0,(1+ksi),0,(1-ksi),0; G3=(1/4)*0,-(1-eta),0,(1-eta),0,(1+eta),0,-(1+eta); G4=(1/4)*0,-(1-ksi),0,-(1+ksi),0,(1+ksi),0,(1-ksi); G=G1;G2;G3;G4; % B=A*G; StrainIP(:,IP,ie)=B*Ue; StressIP(:,IP,ie)=D*StrainIP(:,IP,ie); end %& % Sigma_x(:,1,ie)=T*StressIP(1,:,ie); Sigma_y(:,1,ie)=T*StressIP(2,:,ie); Sigma_xy(:,1,ie)=T*StressIP(3,:,ie); for i=1:4 Sigmax_Nod(1,Ele(ie,i)=Sigmax_Nod(1,Ele(ie,i)+Sigma_x(i,1,ie); Sigmay_Nod(1,Ele(ie,i)=Sigmay_Nod(1,Ele(ie,i)+Sigma_y(i,1,ie); Sigmaxy_Nod(1,Ele(ie,i)=Sigmaxy_Nod(1,Ele(ie,i)+Sigma_xy(i,1,ie); endend%（b）节点应力的近似for iii=1:Nnode n(iii)=size(find(Ele=iii),1); Sigmax_Nod(1,iii)=Sigmax_Nod(1,iii)/n(iii); Sigmay_Nod(1,iii)=Sigmay_Nod(1,iii)/n(iii); Sigmaxy_Nod(1,iii)=Sigmaxy_Nod(1,iii)/n(iii);end Stress_Node=Sigmax_Nod;Sigmay_Nod;Sigmaxy_Nod;end3 算例三问题阐述： 请分别采用3节点三角形单元和4节点四边形等参单元对图3所示的均匀拉伸中心开孔平板进行有限元分析，并验证当圆孔直径远小于平板宽度时，孔边水平应力集中因子为。(1) 给出力学模型图（需反应结构、载荷、边界等信息）；(2) 给出有限元模型（网格剖分图，边界条件；可借助ABAQUS软件生成input文件）；(3) 通过程序计算并提取圆孔附近节点的水平应力分量，得到。图3 均匀拉伸的中心开孔平板解：（1）力学模型：由于结构和载荷都具有对称性，所以只取模型的1/4进行计算分析。（2） 有限元模型：按照ABAQUS建模步骤，一步一步进行。分别为：创建部件、创建材料和截面属性 、定义装配件、定义边界条件和载荷 、划分网格、提交分析作业、后处理、退出ABAQUS/CAE.（网格划分后如下图所示）。 无边界的模型 四节点单元的网格划分 三节点单元的网格划分 获得节点信息（.txt）（3) 四节点单元的应力集中因子K （i）用MATLAB计算得到变形结果如下：(ii)MATLAB计算结果为：与3.0相差还很大，接下来细化网格，再一次计算。然后与ABAQUS的解进行比较。（iii）用细化5倍的网格进行重新划分，网格细化后ABAQUS算法结果显示为：304，此结果与MATLAB计算相符合。，与前面结果相比，细化网格以后的确可以提高精度。(4) 三节点单元的应力集中因子K（i）用ABAQUS计算得到变形结果如下：通过查看可得圆孔水平应力为应力，.与四节点单元相比应力集中因子K小，主要是由于网格细化不够。四节点采用插入点的方式，每条边有40节点；而三节点采用的是长边30节点，短边25节点。所以数值低。加入对三节点细化网格后，结果会不断逼近3.0.（ii）用MATLAB计算结果如下：第5节 点是所需应力点，结果如下：所得结果与ABAQUS结果相差23MPa，.（iii）综上所述：第1、 无论用ABAQUS还是MATLAB对问题进行计算，结果都不同，但结果都十分接近，原因在于上述两种方式都是近似计算，两种方法之间难免都存在误差。第2、 无论用四节点单元还是三节点单元对问题进行离散，只要网格细化足够所得结果都是一样的。这体现了即使离散的方式不同，只要网格细化足够精度还是可以保证的。5、 总结1、 工作仿照教学程序编写了简单的有限元程序（FEATP），然后进行了调试，纠错，运行。进而选择了三个简单的问题进行计算，并将其结果与问题的理论值和 ABAQUS的计算结果进行 了比较，从而验证程序以及模型建立，单元选择的正确性。 进行第一个拉伸问题算例时，我采用 4 节点四 边形单元，得到的结果都还比较理想。这与拉伸问题的常应力分布所得的结果相符合，所以采用四节点单元进行计算。 在进行第二个简支梁问题的算例时，遇到了一些麻烦。首先还是采用四节点单元，程序运行 的结果与理论解有一定的距离，于是猜测有可能是程序本身的问题，也有可能是单元选 择和模型建立的问题；所以采用 ABAQUS进行相同的模拟，得到的结果与程序的运行结 果比较接近，进而判断程序本身并没有错。问题出在模型的建立上。于是鉴于ABAQUS的软件操作的便捷，继