有限元编程综述

有限元编程综述——平面四节点等参单元分析姓名：学号：指导老师：日期：有限元编程综述1摘要：有限元法，主要用来解决复杂结构中力与位移的关系。其基本思想是将具有无限个自由度的连续的求解区域离散为具有有限个自由度、且按一定方式(节点)相互连接在一起的离散体(单元)，即将连续体假想划分为数目有限的离散单元，而单元之间只在数目有限的指定点处相互联结，用离散单元的集合体代替原来的连续体。本文主要介绍有限元分析的基本思想及有限元编程的主要流程，并通过一个平面四节点等参单元来具体论述相关内容。关键词：有限元基本原理平面四节点等参单元分析流程1基本原理及概述所谓有限元法（FEA），其基本原理是把连续的几何机构离散成有限个单元，并在每一个单元中设定有限个节点，从而将连续体看作仅在节点处相连接的一组单元的集合体，同时选定场函数的节点值作为基本未知量并在每一单元中假设一个近似插值函数以表示单元中场函数的分布规律，再建立用于求解节点未知量的有限元方程组，从而将一个连续域中的无限自由度问题转化为离散域中的有限自由度问题。求解得到节点值后就可以通过设定的插值函数确定单元上以至个集合体上的场函数。对每个单元，选取适当的插值函数，使得该函数在子域内部、在子域分界面上以及子域与外界面上都满足一定的条件。单元组合体在已知外载荷作用下处于平衡状态时，列出一系列以节点、位移为未知量的线性方程组，利用计算机解出节点位移后，再用弹性力学的有关公式，计算出各单元的应力、应变，当各单元小到一定程度，那么它就代表连续体各处的真实情况。通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、ABAQUS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU提出了有限元编程综述2更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。因此，必须寻找新的编程语言。随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD的实体造型，优化设计等接口。现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。C语言最初是为操作系统、编译器以及文字处理等编程而发明的。随着不断完善，它已应用到其它领域，包括工程应用软件的编程。近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。用C编写的程序，可移植性极好，几乎不用作多少修改，就可在任何一台装有ANSI、C编译器的计算机上运行。C既是高级语言，也是低级语言，也就是说，可用它作数值计算，也可用它对计算机存储进行操作。2编程思想本程序采用C语言编程，编制平面四边形四节点等参元程序，用以求解平面结构问题。程序采用二维等带宽存储整体刚度矩阵，乘大数法引入约束，等带宽高斯消去法求解位移，然后求中间高斯点的应力，最后用绕节点平均法讲单元应力等效到节点上，再将结果写到tecplot文件中。在有限元程序中，变量数据需赋值的可分为节点信息，单元信息，载荷信息等。对于一个节点来说，需以下信息：节点编号（整型），节点坐标（实型），节点已知位移（实型），节点载荷（实型），边界条件（实型）等。同样，对于一个单元来说，需以下信息：单元的节点联接信息（整型），材料信息（弹性模量，泊松比等）（实型）等。在FORTRAN程序中，以上这些变量混合在一起，很难辨认，使程序的可读性不好，如需要进行单元网络的自适应划分，节点及单元的修改将非常困难。在进行C语言编译过程中，采用结构struct使每个节点信息存储在一个结构体有限元编程综述3数组中，提高程序的可读性，使数据结构更趋于合理。3平面四节点等参单元介绍四节点等参单元实际单元与基本单元的映射关系如图3-1所示坐标的映射关系为：其位移模式和坐标的映射有相同的插值函数，形函数为：单元应变矩阵为：xyxyuxuyuvyx上式一般简写为：B其中B的子块矩阵为图3-1有限元编程综述4iiiiiNxNByNNyx由于iN是、的函数，在iB中的x、y要按照复合函数来求导，即iiiiiiNNNxyxxJNNNxyyy从而有：1iiiiNNxJNNy因此，单元应力矩阵为：DB单元刚度矩阵为：TeAKBDBhdxdy其中积分采用三点高斯积分，3311,11111(,)()(,)nipijijiiijifddfWf其中，2nipn（高斯积分点的总数），i和j或iW是加权系数，i和j是单元内的坐标.。对于三点高斯积分，高斯积分点的位置：110.6,5.09.0，220.0,8.09.0，330.6,5.09.0。结构刚度矩阵为：eeKK结构结点荷载列阵为：eePP有限元编程综述5注意，对于上两式中e的理解不是简单的叠加而是按照对应的自由度集成。总刚平衡方程：KP从式上式求出：1KP4有限元分析流程一个典型的有限元分析过程主要包括以下几个步骤：1)读输入数据，定义节点及单元数组。2)由边界条件计算方程个数，赋值荷载列阵。3)读入在带状存储的总刚度矩阵中单元和载荷信息。4)定义总刚度阵数组。5)组装总刚度阵。6)解方程组。其流程图可见下图：输入边界条件(力、位移)形成各荷载工况的节点荷载阵总刚分解回代求出位移及输出计算应变、应力形成单元刚阵单刚向总刚投放坐标变换输入原始参数计算总刚规模形成总刚方程向总节点荷载阵投放形成单元荷载阵调整几何、弹性矩阵调整单元位移列阵有限元编程综述65计算流程图开始输出节点位移计算局部坐标系下的单元刚度矩阵坐标变换矩阵输入控制信息n=nf·npndf=nd·nf计算单元的IS(ndf)数组取出单元在总体坐标系下的节点位移计算局部坐标系下单元节点位移输入原始数据形成LD（n）数组计算局部坐标下的单元刚度矩阵坐标变换矩阵总体坐标下的单元刚度矩阵形成单元的IS(ndf)数组单元刚度矩阵向结构刚度矩阵叠加进行约束处理求解线性方程组求得结构节点位移计算局部坐标系下单元节点力应力输出节点位移和应力结束有限元编程综述76程序变量及函数说明6.1控制信息np：结构节点总数ne：结构离散单元总数nr1，nr2：总的约束的节点数，nr1，x方向；nr2，y方向nd：每个单元的节点数nf：每个节点的自由度数ld：集中力载荷的个数nm：材料的种类nu1，nu2：非零位移边界条件的节点数，nu1，x方向；nu2，y方向u1，u2：非零位移的大小，u1，x方向；u2，y方向n=nf*np：结构的节点位移总数ndf=nd*nf：每个单元的节点自由度数6.2输入的原始数据x(np):节点的x方向坐标y(np):节点的y方向坐标me（nd，ne）：单元节点的总体编号nrr(nr1+nr2)：约束为零的位移所对应的总体位移编号p（n）：载荷向量nu（nu1+nu2）：位移载荷mat（6，nm）：材料参数6.3程序中的其他标识符LD（n）：存放结构刚度阵所以主对角线元素在A（nn）中的序号IS（ndf）：单元节点位移和节点力在总体位移阵列和载荷阵列中对应的序号EK（ndf，ndf）：总体坐标系下的单元刚度矩阵有限元编程综述8A（nn）：架构刚度阵下三角变带宽一维压缩存储的数组nn：数组A的元素个数RSTG(3):高斯积分点的值H(3):高斯积分点的加权系数S（6，ne）：各单元的应力分量XJAC（2,2）：雅阁比矩阵RJAC（2,2）：雅阁比矩阵的逆PN（2,4）：4个节点处形函数对局部坐标的导数DNX（2,4）：4个节点处形函数对整体坐标的导数FUN（4）：形函数的值7计算结果与Abaqus分析结果的比较7.1、中间带圆孔平面应力板的分析宽40m，长50m，圆孔位于板中心，半径为4m，承受水平方向位移载荷，E=200Pa，3.0，取1tm。用abaqus建模离散，并计算。再讲abaqus离散的节点和单元信息拷贝到本程序的输入文件中，用该程序计算，结果输出成tecplot文件，用tecplot可以查看结果。与abaqus的计算结果进行比较，位移和应力云图如下（左边是程序计算结果，右边是abaqus结果，下同）：有限元编程综述97.2、纤维增强复合材料平面应力板的分析宽40m，长50m，增强纤维位于板中心，纤维半径为8m，承受水平方向位移载荷，基体材料E=20000Pa，25.0，纤维材料E=8000Pa，3.0，取1tm。有限元编程综述107.3、半无限大含裂纹板的应力分析宽40m，长60m，裂纹位于板左侧中间位置，裂纹长10m，承受竖直方向位移载荷，E=200Pa，3.0，取1tm。有限元编程综述117.4结果分析比较：通过以上三个算例发现该程序可以用于计算单材料、双材料、带孔、含裂纹等各种平面问题，加载条件可以是加集中力和加位移，因此，该程序的适用范围还是比较广的。以上三个算例的自编程序所得结果与abaqus分析结果进行比较发现，两者的计算结果很接近，而且自编程序对于孔边应力集中和裂尖应力集中都能很好的表达，说明该程序有很好的精度。第三个算例在裂尖处数值上有些区别，但总的分布规律还是很吻合的，主要是因为本程序是用四节点等参单元，对于应力的奇异性表达效果还不是很好。参考文献[1]冷纪桐有限元技术基础.化学工业出版社，2007[2]张洪武刘辉吴敬凯郑永刚平面4节点广义等参单元计算力学学报，2010.6[2]王勖成有限元法清华大学出版社，2003