有限单元法的编程实现，内附Fortran源代码.docx

作业2：有限单元法的编程实现 2.1有限元概述 有限元分析的基本概念是将求解域离散为若干子域，并通过它们边界上的节点相互联结成为组合体，对每一单元假定一个合适的近似解，然后推导求解这个域总的满足条件，从而得到问题的解。这个解只是近似解，因为实际问题被较简单的问题所代替。由于大多数实际问题难以得到准确解，而有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的工程分析手段。 对于不同物理性质和数学模型的问题，有限元求解法的基本步骤是相同的，只是具体公式推导和运算求解不同。有限元求解问题的基本步骤通常为： 2.1.1确定求解域及求解域的离散化，即子域的剖分 根据实际问题近似确定求解域的物理性质和几何区域。将求解域近似为具有不同有限大小和形状且彼此相连的有限个单元组成的离散域，习惯上称为有限元网络划分，求解域的离散化是有限元法的核心技术之一。选择单元的形式，确定单元数，节点数及单元及节点的编号。 2.1.2单元分析 一个具体的物理问题通常可以用一组包含问题状态变量边界条件的微分方程式表示，为适合有限元求解，通常将微分方程化为等价的泛函形式。 作为弹性力学微分方程的等效积分的形式，虚位移原理与虚应力原理分别是平衡方程与力的边界条件和几何方程与位移边界条件的等效积分形式。在导出它们的的过程中都未涉及到物理方程所以它们适用于线弹性、非线性线弹性及弹塑性的问题。现有的有限元计算多采用以位移为未知量的形式，可用虚位移原理来描述其平衡方程，其矩阵形式为： (2.1) 我们需要对单元构造一个适合的近似解，即推导有限单元的格式，其中包括选择合理的单元坐标系，建立单元试函数，以某种方法给出单元各状态变量的离散关系，从而形成单元矩阵（结构力学中称刚度阵或柔度阵）。 这里以平面问题为例，单元位移与节点位移的关系表示为： (2.2) 应变与节点位移的关系为： (2.3) 应力与节点位移的关系为： (2.4) 单元刚度矩阵的表达式为： (2.5) 单元等效节点载荷矩阵表示为 (2.6) 2.1.3 整体分析 将单元总装形成离散域的总矩阵方程（联合方程组），反映对近似求解域的离散域的要求，即单元函数的连续性要满足一定的连续条件。总装是在相邻单元结点进行，状态变量及其导数（可能的话）连续性建立在结点处。 结构整体刚度矩阵和结构节点载荷矩阵表示为： (2.7) 其中G的作用就是对单元刚度矩阵的组装。 2.1.4联立方程组求解 有限元求解方程的最终的求解格式为： (2.8) 有限元法最终导致联立方程组，其未知量为节点位移。联立方程组的求解可用直接法、选代法和随机法。求解结果是单元结点处状态变量的近似值。存在的主要问题在于利用系数矩阵的对称、正定、稀疏性减少内存空间，增加效率。 2.1.5后处理 后处理则是采集处理分析结果，使用户能简便提取信息，了解计算结果。 2.2有限元方法的编程实现 为加深自己对有限元方法的理解和应用，本人从教材出发，结合已掌握的编程知识，编写了简单的有限元程序。本FEM程序针对弹性平面问题，能计算结点约束、自重和外部荷载的边界条件下的位移、应力大小，并通过商业软件ANSYS的后处理功能，能将计算结果通过图形和数值的方式显示出来。 程序使用的语言为Fortran 95，程序代码详见附件一。 程序使用说明： ①利用ANSYS对研究的结构进行离散，得到模型的单元、结点、约束、荷载信息，用ANSYS导出的文件分别为“ELIST.lis”，“NLIST.lis”，“DLIST.lis”，“FLIST.lis”。 ②平面问题的材料参数弹性模量E、泊松比V在程序中修改。 ③程序计算结果为两个文件，“Model.mac”文件为模型文件，“Output.mac”为位移、应力结果文件。 2.3 FEM程序的验证及应用 2.3.1中心圆孔位移分析验证 一薄板中心有一圆孔，薄板在右边边界受到均布压力的作用，左边和下边分界受到约束，划分三角形单元后如图2-1所示。 计算参数为小孔直径为2m，矩形薄板边长为10m，弹性模量E=106 Pa，泊松比取V=0，右边的均布荷载大小为1N/m。 用编写的程序和ANSYS分别进行计算，将结果进行对比，两种计算得到的位移用云图如图2-2和图2-3所示。 图2-1 中心圆孔网格图 （a）FEM （b）ANSYS 图2-2 X方向的位移云图(负号表示向左，cm) （a）FEM