有限元程序设计.doc

前言有限单元法是在当今工程分析中获得最广泛应用的数值计算方法。由于它的通用性和有效性, 受到工程技术界的高度重视。伴随着计算机科学与技术的快速发展, 现已成为计算机辅助设计与计算机辅助制造的重要组成部分。由于有限元法是通过计算机实现的, 因此它的计算机程序的研制和开发是其理论和方法应用于生产和科研实际的前提和基础。同时所研制和开发的计算机程序又是有限元理论的和方法的研究必要平台。本次程序设计是将 Fortran 程序设计与有限元理论结合。根据有限元理论知识,进行设计程序,从而获得简单平面问题的计算方法。十、平面 4、8 节点有限元公式及计算原理(1 )通过 Serendipity 四边形单元格式构造插值函数。对于 4 节点单元,插值函数为: )1 )(1(4 1???? iiiN???( i=1,2,3,4 ) 对于 8 节点单元,插值函数为: 65112 12 1?NNNN??? 65222 12 1?NNNN??? 76332 12 1?NNNN??? 87442 1?NNNN???)1 )(1(2 1 25?????N)1 )(1(2 1 26?????N )1 )(1(2 1 27?????N)1 )(1(2 1 28?????N)1 )(1(4 1????? iiiN???( i=1,2,3,4 ) (2 )通过????????????????????????)1(2 2100 011 01 1)21 )(1( )1(?????????ED ,?????????????????????????????????????????????????????????????x Ny N y N x NN Nxy y x LN B ii i ii iii0 00 00 0?? mjiBBBB?计算出 D和B 矩阵。(3 )运用高斯积分公式求积分公式近似数值积分解: ???? ni ii baFHdF 1)()(???其中权函数??dlH ba nii???)( )1((4 )单元矩阵的变换: ?????????????????????????????????y N x NJN N i ii i??,因此?????????????????????????yx yxJ 所以面积微元可表示成????ddJdddA||???。(5 )单元刚度矩阵与载荷列向量可表示为: ???????? ni Ti TeddJ DBt BH DBtdxdy BK 1|| ,??????????? ni T ni TiS TTeddJ TtNddJ ftNH Ttds N ftdxdy NP 11||||????将其代入整体刚度矩阵与载荷列向量可得整体刚度矩阵 K 和载荷列向量 P。(6 )模型泛函总势能为: pa Ka a TTp??2 1?。根据最小势能原理 0???a p?可得: P Ka ?运用上式即可求得求得各节点位移。十一、程序结构(1) 有限元程序系统的组成及分析过程工程力学系前处理几何模型、材料参数、边界条件、分析类型定义、网格划分等。核心:各种计算方法的实现。有限元分析本体程序(2 )程序框图(3 )变量说明 NUMNP : 节点数; NUMEL : 单元数; NUMMAT : 材料种类; NLOAD : 载荷数; NDM : 每个节点坐标方向数; NDF :节点自由度; NEN :每个单元节点数。 ID(K,J) :节点 J 的第 K 个自由度约束编号; X(K,J) :节点 J 的第 K 个坐标。 IX(K,J):J 单元中 K 节点的全局节点编号。 EE :杨氏弹性模量; XNU :泊松比; ITYPE :问题类型。 F(K,J) :节点 J在K 方向上集中载荷大小。 ALFL=TRUE :非对称矩阵集合; ALFL=FALSE :对称矩阵集合:S:单元刚度矩阵; P:载荷和内力向量; AD :对角元素; AU :上三角元素; AL :下三角元素; JP :上(下)三角元素每行(列)最后元素; LD :一个单元中各自由度方程数。 JD :列高; ID :边界条件。(4 )数据输入格式整数: I5 格式;小数: 格式。输入离散模型数据型数据计算单元刚度阵组集结构刚度矩阵计算单元等效结点载荷组集结构结点载荷列阵引入位移边界条件求解线性方程组其它辅助计算输出结果结束单元循环形成 K形成 P消除 K的奇异求解 Ka =P,得结点位移 a 计算应力、应变等以图形、曲线和表格等方式表达、分析计算结果。后处理(5 )程序本体 C---------------------------------------- C C..... FEA2DP---A FINITE ELEMENT ANALYSIS PROGRAM FOR C 2D ELA