有限元教材-第十章有限元程序设计

第十章有限元程序设计有限元方法作为一门系统的技术，仅学会了它的基本理论是远远不够的，只有形成完整的计算程序，问题才最终得到了解决。完成这样的有限元程序设计是一项工作量很大的工程。本章就是要结合简单的有限元教学程序FEMED，简要介绍有限元程序设计技术。FEMED是专为有限元程序设计教学编制的程序，它不包含复杂的前后处理功能，可进行平面问题及平面桁架的线弹性静力分析，在程序结构上与大型程序类似，具有计算单元的任意扩充功能，在方程的组集和求解上也采用了较为流行的变带宽存储方式。有限元程序大致可分为两类，第一类是专用程序，主要用于研究或教学，一般这类程序规模较小，前后处理功能较弱。用于研究的程序能够解一些特殊的问题，满足研究工作的需要。而教学程序则是为了学生了解有限元的主要结构和设计方法设计的，程序比较简单，FEMED就属于这类程序。第二类是大型通用程序，是大型结构分析的得力工具，目前国际上流行的大约有2000多种。常用的有NASTRAN、MARC、ANSYS、ADINA和ABAQUS等。这类程序一般前后处理功能比较强，有友好的界面，能进行大型计算，但往往无法完成具有特殊要求的计算。通过本章的学习，使读者初步掌握有限元编程的基本方法，具有开发特殊功能的专用程序或为通用程序开发具有特殊功能的计算模块的能力。§10.1有限元程序的基本结构有限元程序一般包括三项基本内容：前处理、结构分析和后处理。早期有限元分析软件的研究重点在于推导新的高效率求解方法和高精度的单元，随着数值分析方法的逐步完善，尤其是计算机内存和运算速度的飞速发展，整个计算系统用于求解运算的时间越来越少，加之求解问题的日益大型化和复杂化，使得数据准备和运算结果的表现问题日益突出。因此目前几乎所有的商业化有限元程序系统都有功能很强的前后处理模块，这直接关系到分析软件的可推广性。它是商用有限元软件不可或缺的部分，但它不是有限元的中心部分，在本书中不作详细介绍。10.1.1前处理在结构分析之前，必须完成的工作就是前处理的任务，它主要包括以下内容：结构造型、选取单元类型，、网格剖分，从而形成节点数、节点编码、节点坐标、单元数、单元节点编码等。另外根据所解题目的不同，还需读入不同的材料参数、边界约束条件及载荷工况等。有限元前处理程序的功能就是为用户提供一种工具，使其能尽可能方便地完成上述工作。而这一问题解决的好坏直接关系到程序使用者需付出的工作量，而网格剖分质量还直接影响计算精度。很多程序都建立了对用户非常友好的图形用户界面，使用户能以可视图形方式直观快速地进行结构造型和网格自动划分，生成有限元分析所需数据。在上述工作中工作量最大的就是结构造型和网格剖分部分，目前越来越倾向于由专用软件完成上述工作，所以当今所有的商业化有限元软件开发商都开发了和著名的CAD软件的接口，从而极大地提高了设计水平和效率。今天，工程师可以在集成的CAD和FEA软件环境中快捷地解决一个在以前无法应付的复杂工程分析问题。在FEMED程序中，只包含了很简单的节点、网格及给定位移的自动生成。所有前处理功能是由INPUT及相关子程序完成的。10.1.2结构分析结构分析部分是有限元程序的主体部分，主要需计算插值函数矩阵[N]，插值函数的导数矩阵[B]，进而进行数值积分得到单元刚度矩阵，组集总刚度矩阵及总载荷向量，解线性代数方程组，得到节点位移。对于动力问题则需计算质量矩阵，求解特征值问题。而对于动力响应问题及非线性问题，则需进行多增量步的计算。在得到节点位移的基础上还需进行应力的计算。在本章中重点介绍的就是结构分析中的静力分析部分，在FEMED程序中主要包含的也就是这部分。通过这部分的学习，希望读者能掌握有限元程序的基本编程技巧。10.1.3后处理在大型分析软件中，程序的后处理功能也是非常重要的，在有限元结构分析完成后，通过强有力的后处理图形功能，可以给出各应力分量、等效应力等的分布云图，结构的变形图，振型图和实时动画等，使得使用者对各物理量在整个结构中的变化情况有一个全面的认识，也可检查网格剖分和所加载荷及约束是否有误。否则面对输出成千上万个计算结果数据，往往使人如坠雾中，不得要领，需花费大量的时间进行数据的整理和解释。由于是简单的教学程序，在FEMED程序中不包含后处理功能，仅仅是将计算结果输出。图10.1给出了一般求解静力学问题的程序框图，FEMED程序的结构基本如图所示，但没有包含带宽优化部分。§10.2数组的半动态分配在有限元程序中需开辟许多不同功能的数组，而这些数组的大小则与所解题目的大小及类型有关，如节点坐标数组、节点位移数组等与节点数相关，而单元节点编号数组等则与单元数相关。此外，上述数组的大小还和自由度维数等相关。对于一个有限元程序而言，应具有求解各类问题及各开始前处理带宽优化，计算对角元地址形成单元刚度阵解方程求应力、节点力后处理装配总刚度阵结构分析部分单元循环单元循环图10.1种规模问题的机动性，因此求解不同的问题，应定义不同大小的数组。FORTRAN语言规定数组的大小必须预先给定，在子程序中可以定义大小可变化的可调数组，但其内存分配必须由该子程序的上一级模块确定。最简单的办法就是每个数组都定义得足够大，使得要求解的问题都能被满足，但这显然会带来内存的巨大浪费，因此，一般采用数组半动态分配的方法来解决上述问题，FEMED也采用了这一方法。在主控程序PCONTR中开辟了两个大数组A(*)和M(*)，分别存储实型量和整型量，将各数组中的数顺序存放在大数组中，在读入控制信息NJ（节点数）、NE（单元数）、NUMMAT（材料数）、NDF（节点自由度数）、NDM（空间维数）、NEN（单元节点数）等之后，就可进行存储单元的分配。在图10.2中以大数组A的存储为例，先存放节点坐标，占用NJ×NDM个存贮单元，再存放材料参数，占用NUMMAT×20个存贮单元……图10.2在上述过程中，最重要的是计算出各数组存放的起始位置又称作地址，如NXC为XC数组在A中的地址，ND为D数组在A中的地址……，这项工作由POINT子程序完成。需要分别计算出大数组A和M中所存数组的所有地址，最后计算出A和M数组需要的总长度，与预设的数组长度作比较，若没有超界计算继续进行，若超界则计算停止，显示错误信息。需要指出的是，在上述数组的半动态分配时，将一系列存储容量不定的一维、二维数组均用地址对应关系分配到事先容量固定的大数组A和M中去，对于大数组A和M来说存储量是不变的，但各数组在大数组中的存储分配却是动态的。我们把这种关系称为数组的半动态分配，这一过程是通过哑实结合完成的。这些数组称为可调数组，即其数组长度是可调的。在子程序的数组说明中，一维数组写成下列形式是等价的，对数组长度没有限制：F(*)，F(1)，F(NJ)，对于二维数组，则只可一个方向自由变动。由于FORTRAN规定二维数组先变行后变列，因此行的长度必须是确定的。，在子程序的数组说明中，下列形式是等价的：XC(2,*)，XC(2,1)，XC(2,NJ)，XC(NDM,*)，(NDM=2)注意，只有对上一模块中已定义过的数组才可作如上的数组说明，否则，则必须准确确定数组的长度。§10.3构造单元刚度矩阵及总刚度矩阵的组集构造单元刚度矩阵是由单元模块完成的，不同的单元调用不同的模块，在FEMED中提供了平面桁架单元（TRUSS）和平面单元（PLANE）两个单元计算模块。读者可根据需要随意开发单元类型，通过ELEMLIB子程序接入程序中。针对不同的单元，在不同的单元模块中构造不同的[B]、[D]矩阵，单元刚度矩阵可表示为：VdVeeTKBDB（10.1）对于常应变单元，上述被积函数为常数，直接用被积函数值乘积分域的大小即可。而对于等参元等非常应变单元，则需对上式做GAUSS积分，最终得到一个二维的单元刚度矩阵。NJ×NDMNUMMAT×20NXL=ND+20*NUMMAT………………NXC=1ND=NXC+NJ*NDM将单元刚度矩阵正确组集为总刚度矩阵，这一工作是由ADDSTF子程序完成的。下面将举例说明总刚度矩阵的组集方法。例1：平面桁架结构如图示，共有5个节点，6个单元。每个单元形成一个4×4的单元刚度矩阵，总刚度阵为10×10。图10.3以1号单元为例，其左节点为2，右节点为4，其4×4的单元刚度矩阵的各个参数应迭加在(10.2)式总刚度矩阵中的图示位置上。需要指出的是，总刚度矩阵的一个位置上往往有多个单元刚度矩阵做贡献，如与3号节点对应的(5,5)，(6,5)，(6,6)位置上就有2、3、4和6号单元都作了贡献。全部单元组集完毕后，总刚度矩阵如式（10.2），形成一对称带状阵，即仅对角元附近的元素不为零，而远离对角线的元素皆为零。解题越大带状性质越明显。00000000444342410033323100222111][称对K（10.2）若储存的是总刚度矩阵的上三角元素，则带宽为每行从对角元，到最后一个非零元素的个数；若储存的是总刚度矩阵的下三角元素，则从第一个非零元素到对角元元素的个数称为这行的带宽。每行的带宽各不相同，且与节点编号的方式有关，在大型商业软件中，一般都包含带宽优化模块，就是为了减少总体带宽。§10.4总刚度矩阵的变带宽存储总刚度阵占用的空间是有限元占用空间的绝大部分，一个软件的解题规模往往是受此限制。所以，一直以来人们在这方面想了很多的办法。很显然，将总刚度阵作为一个二维数组，完全储存起来是最直接和最自然的方法，但这也是占用空间最大、最不合理的方法，一般很少采用。考虑到总刚度阵具有带状稀疏对称的特征，使用半带宽存储无疑是一个最自然的改进，它利用一个二维数组，只存储每行从对角元到这一行的最后一个非零元，如图10.4所示2123451234560000000000000000000000000000000000000称对图10.4由10.4的示意图知，原10×10的矩阵现在可用10×4的矩阵存储，存储空间减少了很多，但若其中有个别行的带宽特别大，则这种方法依然有存贮空间的大量浪费。为解决上述问题，一般多采用变带宽存储又称一维压缩存储的方法，这样刚度阵中大量的零既不需要存储也不参与计算，大大减少了存储空间和计算量，是目前使用较为广泛的方法。在有限元的发展历史中波前法也是一种不可不提的计算方法，当计算机的内存非常有限时，即使采用变带宽存储也无法求解较大规模的问题，波前法曾经起到了不可忽视的作用。由于波前法无需组集总刚，所以对内存要求非常低，曾经被广泛地使用，但它也有非常明显的缺点，内外存交换非常频繁，导致计算速度受到很大的限制。随着计算机技术突飞猛进地发展，计算机的内存有了本质的提高，所以波前法已经很少被采用了。本书重点介绍变带宽存储的方法。变带宽存储就是考虑到总刚度矩阵的对称性及带状稀疏的特点，用一维数组将总刚度矩阵下三角（或上三角）中的非零元素储存起来，如图10.5中虚线以上部分，注意一行中夹在中间的零必须储存。图10.5设A(i,j)为具有带状对称性质的二维数组，若将其存储在一维数组B(k)中，则必须开辟一个存放对角元地址的数组JP(i)，指明每一个对角元在一维数组中的位置。如上图中矩阵的JP(i)应为：1，3，6，9，11，15，18，21，24，27。若第i行，第j列的数据存在一维数组的第k个位置上，则k与i、j之间的对应关系为k=JP(i)-i+j(10.3)第i行的第一个非零元素在mi列，则mi=i-JP(i)+JP(i-1)+1(i1)(10.4)0000000000000000000