有限元程序课程设计.doc

王顺金：有限元程序课程设计 重庆大学本科学生课程设计任务书课程设计题目有限元程序设计学院资源及环境科学学院专业工程力学年级2010级已知参数和设计要求：1.独立完成有限元程序设计。2.独立选择计算算例，并能通过算例判断程序的正确性。3.独立完成程序设计报告，报告内容包括理论公式、程序框图、程序本体、计算算例，算例结果分析、结论等。学生应完成的工作：1. 复习掌握有限单元法的基本原理。2. 掌握弹性力学平面问题3节点三角形单元或4节点等参单元有限元方法的计算流程，以及单元刚度矩阵、等效节点载荷、节点应变、节点应力和高斯积分等的计算公式。3. 用Fortran语言编写弹性力学平面问题3节点三角形单元或4节点等参单元的有限元程序。4. 在Visual Fortran 程序集成开发环境中完成有限元程序的编辑和调试工作。5. 利用编写的有限元程序，计算算例，分析计算结果。6. 撰写课程设计报告。目前资料收集情况（含指定参考资料）：1.王勖成，有限单元法，北京：高等教育出版社，2002。2.O.C. Zienkiewicz, R. L. Taylor, Finite Element Method, 5th Eition, McGraw-Hall Book Company Limited, 2000。3.张汝清，董明，结构计算程序设计，重庆：重庆大学出版社，1988。课程设计的工作计划：1. 第1周星期一上午：教师讲解程序设计方法，程序设计要求和任务安排。2. 第1周星期一至星期二完成程序框图设计。3. 第1周星期三至第2周星期四完成程序设计。4. 第2周星期五完成课程设计报告。任务下达日期 2013 年 6 月 6 日完成日期 2013 年 07 月 03 日指导教师 （签名）学 生 （签名）一、前言有限单元法是在当今工程分析中获得最广泛应用的数值计算方法，由于其通用性和有效性，受到工程技术界的高度重视。伴随着计算机科学和技术的快速发展，现已成为计算机辅助设计和计算机辅助制造的重要组成部分。由于有限元法是通过计算机实现的，因此有限元程序的编制及相关软件的研发就变得尤为重要。从二十世纪五十年代以来，有限元软件的发展按目的和用途可分为专用软件和大型通用商业软件，而且软件往往集成了网格自动划分，结果分析和显示等前后处理功能，而且随着时间的发展，大型通用商业软件的功能由线性扩展到非线性，由结构扩展到非结构等等，这一系列强大功能的实现与运用都要求我们对有限元法的基础理论知识有较为清楚的认识以及对程序编写的基本能力有较好掌握。1.元分析程序的理论基础作为弹性力学微分方程的等效积分形式，虚位移原理和虚应原理分别是平衡方程与力的边界条件和几何方程与位移边界条件的等效积分形式。将物理方程引入虚位移原理和虚应力原理可以分别导出最小位能原理和最小余能原理，它们本质上和等效积分的伽辽金“弱”形式相一致，这是建立弹性力学有限原方程一般表达格式的理论基础。对于通过弹性力学变分原理建立的弹性力学问题有限元方法，其未知场变量是位移，以节点位移为基本未知量，并以最小位能原理为基础建立的有限单元为位移元。弹性力学平面问题有限元分析表达格式的建立步骤一般为：a的模型进行单元离散，一般采用三角形单元或四边形单元，不过由于四边形单元具有更高精度，应用更为普遍，一般而言一次或二次单元已足以满足精度要求。离散后对单元进行编号，并且给定单元各节点的整体编码以及局部坐标系下按逆时针局部编码。b坐标系下对各单元构造形函得到由单元各节点位移表示的单元位移形式，进而得到单元刚度矩阵，利用等参元性质和雅阁比矩阵进行组集，建立整体刚度矩阵。c等效结点载荷列阵并组集成结构节点等效载荷列阵得到单元各节点位移与结构节点等效载荷列阵的线代方程组，求解得到位移，进而可得应力应变。 根据以上一般步骤，编制相应的计算机程序并采用数值积分方法处理有关数学计算就可以得到完整的弹性力学问题有限元分析程序。2. 问题有限元分析教学程序本程序可对二维弹性力学问题行有限元分析计算。计算采用的单元形式为四边形四节点单元或者四边形八节点单元，对于对称和非对称矩阵采用变带宽存储方法，最终输出结果为单元各节点的位移。二、平面4、8节点有限元公式及计算原理2.1基本公式（1） 有限元平衡方程（2） 单元刚度矩阵（3）单元等效结点载荷2.2单元位移插值及坐标变换（1） 通过Serendipity四边形单元格式构造插值函数。对于4节点单元，插值函数为： 对于8节点单元，插值函数为： 其中： （i=1,2,3,4）（2） 位移插值 （3） 坐标变换 2.3单元刚度的计算（1）弹性矩阵平面应力问题：；平面应变问题：；式中E和v分别为材料Young氏模量和Poisson比。（2） 应变矩阵形状函数对局部坐标偏导：，形状函数对整体坐标偏导：Jacobi矩阵：Jacobi矩阵的逆：其中：（3） 单元刚度矩阵计算单元刚度矩阵Gauss积分：单元刚度矩阵子矩阵： （i=1，2）2.4等效结点载荷（1） 集中载荷直接施加在结点上（2） 体积力产生的等效结点载荷向量体积力产生的等效结点载荷Gauss积分：体积力产生的等效结点载荷子向量：其中：若体积力作用方向于整体坐标系的y轴正方向角度，则：2.5单元应变和应力计算（1） 单元应变单元中任意点的应变：单元中任意点的应变：其中可作为Gauss积分点或单元结点的整体坐标，其于局部坐标间的关系为：,(2) 单元应力三、程序结构（1）有限元程序系统的组成及分析过程我们本次设计主要针对二维弹性力学静力学问题。有限元程序系统通常包括处理、有限元程序本体和后处理三部分。有限元软件系统的组成及分析过程如下图所示：CAD模型、材料参数、边界条件、分析类型定义、网格划分等。有限元分析本体程序前 处 理后 处 理核心：各种计算方法的实现。以图形、曲线和表格等方式表达、分析计算结果。（2）有限元计算程序框图 输入离散模型数据型数据计算单元刚度阵组集结构刚度矩阵计算单元等效结点载荷组集结构结点载荷列阵引入位移边界条件求解线性方程组其它辅助计算输出结果结 束单元循环形成 K形成 P消除 K 的奇异求解 Ka=P得结点位移 a计算应力、应变等（3）有限元程序主要编程技巧a.结构模块化主程序：组织调用子程序，过程清晰，不涉及复杂计算。子模块：实现特定功能，可用多个子程序组成。子程序：实现简单功能，便于调用和扩展。b.覆盖技术在主程序中定义一个整数组和一个实型数组，分别存放整型和实型边界数组。设计动态数组表。动态数组覆盖技术。（4）有限元程序主程序打开输入数据，调用PMESH，调用PROFIL计算K的对角元地址和方程号，分配数组大小，调用ASSEM组集整体K，调用PLOAD组集载荷列向量，调用DATRI对刚度矩阵进行三角分解，调用DASOL解方程，调用PRTDIS输出结点位移，调用TRESSS计算GAUSS点和结点应变和应力。（5）变量说明NUMNP：节点总数NUMEL：单元总数NUMMAT：材料种类NLOAD：载荷数NDM：每个节点坐标方向数NDF：节点自由度NEN：每个单元节点数。ID（K，J）：节点J的第K个自由度约束编号X（K，J）：节点J的第K个坐标IX（K，J）：J单元中K节点的全局节点编号F（K，J）：结点J第K方向上的集中载荷F（NDF，NLOAD）：集中载荷X（NDM，NUMNP）：结点坐标IX（NEN，NUMEL）：单元定义F(K,J)：节点J在K方向上集中载荷大小ALFL=TRUE：非对称矩阵集合ALFL=FALSE：对称矩阵集合EE：杨氏弹性模量XNU：泊松比ITYPE：问题类型：单元刚度矩阵：载荷和内力向量AD：对角元素AU：上三角元素AL：下三角元素JP：上（下）三角元素每行（列）最后元素LD：一个单元中各自由度方程数JD：列高ID：边界条件，=0：无约束；=1：固定约束。（6）数据输入格式整数：I5格式小数：F10.4格式（7）程序本体（见附录）4、 程序验证算例算例：平面板受轴向拉伸长3m宽0.4m的平面板左端铰接固定，右端受到水平方向向右、大小为10N的均布载荷，模型如下图所示：F=10NF=10NF=10N1284519其中结构参数：F=10N,L=3m,h=0.4m 材料参数：E=200GPa,=0.31、建立坐标系，将板划分成6个如上图所示的四边形四节点单元，进行单元编号，均布载荷处理为右端三节点每个节点受水平向右大小10N的力，竖直方向不受力。1.1程序input文件中输入：example(NEN=4) 12 6 1 3 2 2 4 1 1 1 0.0000 0.0000 2 0 0 1.0000 0.0000 3 0 0 2.0000 0.0000 4 0 0 3.0000 0.0000 5 1 1 0.0000 0.2000 6 0 0 1.0000 0.2000 7 0 0 2.0000 0.2000 8 0 0 3.0000 0.2000 9 1 1 0.0000 0.4000 10 0 0 1.0000 0.4000 11 0 0 2.0000 0.4000 12 0 0 3.0000 0.4000 1 1 2 6 5 2 2 3 7 6 3 3 4 8 7 4 5 6 10 9 5 6 7 11 10 6 7 8 12 11 200.0E+9 0.3000 1 4 10.0000 0.0000 8 10.0000 0.0000 12 10.0000 0.00001.2调试程序得到结果为：nodal displacements 1 0.0000E+00 0.0000E+00 2 0.3638E-09 0.2813E-10 3 0.7369E-09 0.2030E-10 4 0.1126E-08 0.2437E-10 5 0.0000E+00 0.0000E+00 6 0.3609E-09 -0.2981E-23 7 0.7429E-09 -0.1095E-22 8 0.1104E-08 -0.2132E-22 9 0.0000E+00 0.0000E+00 10 0.3638E-09 -0.2813E-10 11 0.7369E-09 -0.2030E-10 12 0.1126E-08 -0.2437E-102、建立坐标系，将板划分成6个如上图所示的四边形八节点单元，进行单元编号，均布载荷处理为右端三节点每个节点受水平向右大小10N的力，竖直方向不受力。其中结构参数：F=10N,L=3m,h=0.4m 材料参数：E=200GPa,=0.32.1程序input文件中输入：example(NEN=8) 29 6 1 3 2 2 8 1 1 1 0.0000 0.0000 2 0 0 1.0000 0.0000 3 0 0 2.0000 0.0000 4 0 0 3.0000 0.0000 5 1 1 0.0000 0.2000 6 0 0 1.0000 0.2000 7 0 0 2.0000 0.2000 8 0 0 3.0000 0.2000 9 1 1 0.0000 0.4000 10 0 0 1.0000 0.4000 11 0 0 2.0000 0.4000 12 0 0 3.0000 0.4000 13 0 0 0.5000 0.0000 14 0 0 1.5000 0.0000 15 0 0 2.5000 0.0000 16 1 1 0.0000 0.1000 17 0 0 1.0000 0.1000 18 0 0 2.0000 0.1000 19 0 0 3.0000 0.1000 20 0 0 0.5000 0.2000 21 0 0 1.5000 0.2000 22 0 0 2.5000 0.2000 23 1 1 0.0000 0.3000 24 0 0 1.0000 0.3000 25 0 0 2.0000 0.3000 26 0 0 3.0000 0.3000 27 0 0 0.5000 0.4000 28 0 0 1.5000 0.4000 29 0 0 2.5000 0.4000 1 1 2 6 5 13 17 20 16 2 2 3 7 6 14 18 21 17 3 3 4 8 7 15 19 22 18 4 5 6 10 9 20 24 27 23 5 6 7 11 10 21 25 28 24 6 7 8 12 11 22 26 29 25 200.0E+9 0.3000 1 4 10.0000 0.0000 8 10.0000 0.0000 12 10.0000 0.00002.2调试程序得到结果为：nodal displacements 1 0.0000E+00 0.0000E+00 2 0.3734E-09 0.2492E-10 3 0.7549E-09 0.4768E-10 4 0.1185E-08 0.6785E-10 5 0.0000E+00 0.0000E+00 6 0.3770E-09 -0.1562E-22 7 0.7371E-09 -0.3614E-22 8 0.1086E-08 -0.7600E-22 9 0.0000E+00 0.0000E+00 10 0.3734E-09 -0.2492E-10 11 0.7549E-09 -0.4768E-10 12 0.1185E-08 -0.6785E-10 13 0.1849E-09 0.2673E-10 14 0.5567E-09 0.1616E-10 15 0.9277E-09 0.1568E-11 16 0.0000E+00 0.0000E+00 17 0.3737E-09 0.1232E-10 18 0.7517E-09 0.2430E-10 19 0.1114E-08 0.3177E-10 20 0.1810E-09 -0.5215E-23 21 0.5559E-09 -0.2565E-22 22 0.9462E-09 -0.5238E-22 23 0.0000E+00 0.0000E+00 24 0.3737E-09 -0.1232E-10 25 0.7517E-09 -0.2430E-10 26 0.1114E-08 -0.3177E-10 27 0.1849E-09 -0.2673E-10 28 0.5567E-09 -0.1616E-10 29 0.9277E-09 -0.1568E-113.数据分析其理论计算过程如下，对于平面应力问题有： 其中 ： 代入计算得 任一距原点水平距离为x的点其位移为：理论解与数值解对比图如下：可见接近自由端处两曲线很接近，当X=3m时，u=15x10E-9m，对应输出文件中的节点3、7、11的输出值和3、7、11点有限元解较为接近，误差为很小。同时可以发现靠近固定端处理论解与数值解误差变大。分别对x=2处的理论解和有限元解进行比较发现误差较小。考虑到实际情况固定端的铰支约束条件对实际位移的影响以及离散单元数目较少等可能造成的误差，可以认为该程序是正确的。另外数值由可以看出同的对称性，符合基本力学理论结果。程序可用于平面问题的求解，可获得较好的位移分布。程序适用于平面应力、平面应变状态结构4、8节点单元有限元分析，可求出结构的位移、应力以及应变。五、算例1、算例：选择超静定梁中点受集中载荷作用的结构，如下图所示：其中结构参数：F=10kN,L=4m,h=0.4m 材料参数：E=200GPa,=0.32、单元网格划分及规定对上述结构采用 101,102三种单元网格划分形式，分别采用四结点和八结点四边形单元。 规定单元编号先沿 y 正方向，再沿 x 正方向。为了方便比较，所有编号位置固定，例如第12结点，不管是101单元四结点，还是101单元八结点，其位置不变。其不同之处在于八结点增加的结点编号在四结点的基础上插入，而不是替换原有的结点编号。特殊点如下所示：11611217221233116281712322(1).10X1单元四结点input文件中输入：WSJ_10X1(NEN=4) 22 10 1 1 2 2 4 1 1 1 0.0000 0.0000 2 0 0 0.4000 0.0000 3 0 0 0.8000 0.0000 4 0 0 1.2000 0.0000 5 0 0 1.6000 0.0000 6 0 0 2.0000 0.0000 7 0 0 2.4000 0.0000 8 0 0 2.8000 0.0000 9 0 0 3.2000 0.0000 10 0 0 3.6000 0.0000 11 1 1 4.0000 0.0000 12 1 1 4.0000 0.4000 13 0 0 3.6000 0.4000 14 0 0 3.2000 0.4000 15 0 0 2.8000 0.4000 16 0 0 2.4000 0.4000 17 0 0 2.0000 0.4000 18 0 0 1.6000 0.4000 19 0 0 1.2000 0.4000 20 0 0 0.8000 0.4000 21 0 0 0.4000 0.4000 22 1 1 0.0000 0.4000 1 1 2 21 22 2 2 3 20 21 3 3 4 19 20 4 4 5 18 19 5 5 6 17 18 6 6 7 16 17 7 7 8 15 16 8 8 9 14 15 9 9 10 13 14 10 10 11 12 13 200.0E+9 0.3000 1 17 0.0000 -10.0E+310X1单元四结点input文件中输出：nodal displacements 1 0.0000E+00 0.0000E+00 2 -0.2029E-06 -0.2673E-06 3 -0.3048E-06 -0.8380E-06 4 -0.3055E-06 -0.1509E-05 5 -0.2059E-06 -0.2083E-05 6 0.2994E-21 -0.2330E-05 7 0.2059E-06 -0.2083E-05 8 0.3055E-06 -0.1509E-05 9 0.3048E-06 -0.8380E-06 10 0.2029E-06 -0.2673E-06 11 0.0000E+00 0.0000E+00 12 0.0000E+00 0.0000E+00 13 -0.2015E-06 -0.2671E-06 14 -0.3018E-06 -0.8376E-06 15 -0.3012E-06 -0.1510E-05 16 -0.1986E-06 -0.2077E-05 17 -0.1853E-21 -0.2365E-05 18 0.1986E-06 -0.2077E-05 19 0.3012E-06 -0.1510E-05 20 0.3018E-06 -0.8376E-06 21 0.2015E-06 -0.2671E-06 22 0.0000E+00 0.0000E+00(2).10X1单元八结点input文件中输入：WSJ_10X1(NEN=8) 53 10 1 1 2 2 8 1 1 1 0.0000 0.0000 2 0 0 0.4000 0.0000 3 0 0 0.8000 0.0000 4 0 0 1.2000 0.0000 5 0 0 1.6000 0.0000 6 0 0 2.0000 0.0000 7 0 0 2.4000 0.0000 8 0 0 2.8000 0.0000 9 0 0 3.2000 0.0000 10 0 0 3.6000 0.0000 11 1 1 4.0000 0.0000 12 1 1 4.0000 0.4000 13 0 0 3.6000 0.4000 14 0 0 3.2000 0.4000 15 0 0 2.8000 0.4000 16 0 0 2.4000 0.4000 17 0 0 2.0000 0.4000 18 0 0 1.6000 0.4000 19 0 0 1.2000 0.4000 20 0 0 0.8000 0.4000 21 0 0 0.4000 0.4000 22 1 1 0.0000 0.4000 23 0 0 0.2000 0.0000 24 0 0 0.6000 0.0000 25 0 0 1.0000 0.0000 26 0 0 1.4000 0.0000 27 0 0 1.8000 0.0000 28 0 0 2.2000 0.0000 29 0 0 2.6000 0.0000 30 0 0 3.0000 0.0000 31 0 0 3.4000 0.0000 32 0 0 3.8000 0.0000 33 1 1 4.0000 0.2000 34 0 0 3.6000 0.2000 35 0 0 3.2000 0.2000 36 0 0 2.8000 0.2000 37 0 0 2.4000 0.2000 38 0 0 2.0000 0.2000 39 0 0 1.6000 0.2000 40 0 0 1.2000 0.2000 41 0 0 0.8000 0.2000 42 0 0 0.4000 0.2000 43 1 1 0.0000 0.2000 44 0 0 0.2000 0.4000 45 0 0 0.6000 0.4000 46 0 0 1.0000 0.4000 47 0 0 1.4000 0.4000 48 0 0 1.8000 0.4000 49 0 0 2.2000 0.4000 50 0 0 2.6000 0.4000 51 0 0 3.0000 0.4000 52 0 0 3.4000 0.4000 53 0 0 3.8000 0.4000 1 1 2 21 22 23 42 44 43 2 2 3 20 21 24 41 45 42 3 3 4 19 20 25 40 46 41 4 4 5 18 19 26 39 47 40 5 5 6 17 18 27 38 48 39 6 6 7 16 17 28 37 49 38 7 7 8 15 16 29 36 50 37 8 8 9 14 15 30 35 51 36 9 9 10 13 14 31 34 52 35 10 10 11 12 13 32 33 53 34 200.0E+9 0.3000 1 17 0.0000 -10.0E+310X1单元八结点input文件中输出：nodal displacements 1 0.0000E+00 0.0000E+00 2 -0.2896E-06 -0.3882E-06 3 -0.4446E-06 -0.1198E-05 4 -0.4464E-06 -0.2169E-05 5 -0.2954E-06 -0.2984E-05 6 -0.1067E-19 -0.3336E-05 7 0.2954E-06 -0.2984E-05 8 0.4464E-06 -0.2169E-05 9 0.4446E-06 -0.1198E-05 10 0.2896E-06 -0.3882E-06 11 0.0000E+00 0.0000E+00 12 0.0000E+00 0.0000E+00 13 -0.2881E-06 -0.3880E-06 14 -0.4419E-06 -0.1198E-05 15 -0.4415E-06 -0.2170E-05 16 -0.2928E-06 -0.2992E-05 17 0.1160E-19 -0.3406E-05 18 0.2928E-06 -0.2992E-05 19 0.4415E-06 -0.2170E-05 20 0.4419E-06 -0.1198E-05 21 0.2881E-06 -0.3880E-06 22 0.0000E+00 0.0000E+00 23 -0.1582E-06 -0.1219E-06 24 -0.3874E-06 -0.7546E-06 25 -0.4638E-06 -0.1684E-05 26 -0.3906E-06 -0.2617E-05 27 -0.1666E-06 -0.3254E-05 28 0.1666E-06 -0.3254E-05 29 0.3906E-06 -0.2617E-05 30 0.4638E-06 -0.1684E-05 31 0.3874E-06 -0.7546E-06 32 0.1582E-06 -0.1219E-06 33 0.0000E+00 0.0000E+00 34 0.7195E-09 -0.3635E-06 35 0.1556E-08 -0.1195E-05 36 0.2019E-08 -0.2174E-05 37 0.3347E-08 -0.3008E-05 38 0.7962E-21 -0.3388E-05 39 -0.3347E-08 -0.3008E-05 40 -0.2019E-08 -0.2174E-05 41 -0.1556E-08 -0.1195E-05 42 -0.7195E-09 -0.3635E-06 43 0.0000E+00 0.0000E+00 44 0.1574E-06 -0.1217E-06 45 0.3852E-06 -0.7543E-06 46 0.4599E-06 -0.1684E-05 47 0.3863E-06 -0.2615E-05 48 0.1608E-06 -0.3252E-05 49 -0.1608E-06 -0.3252E-05 50 -0.3863E-06 -0.2615E-05 51 -0.4599E-06 -0.1684E-05 52 -0.3852E-06 -0.7543E-06 53 -0.1574E-06 -0.1217E-06(3).10X2单元四结点input文件中输入：WSJ_10X2(NEN=4) 33 20 1 1 2 2 4 1 1 1 0.0000 0.0000 2 0 0 0.4000 0.0000 3 0 0 0.8000 0.0000 4 0 0 1.2000 0.0000 5 0 0 1.6000 0.0000 6 0 0 2.0000 0.0000 7 0 0 2.4000 0.0000 8 0 0 2.8000 0.0000 9 0 0 3.2000 0.0000 10 0 0 3.6000 0.0000 11 1 1 4.0000 0.0000 12 1 1 4.0000 0.4000 13 0 0 3.6000 0.4000 14 0 0 3.2000 0.4000 15 0 0 2.8000 0.4000 16 0 0 2.4000 0.4000 17 0 0 2.0000 0.4000 18 0 0 1.6000 0.4000 19 0 0 1.2000 0.4000 20 0 0 0.8000 0.4000 21 0 0 0.4000 0.4000 22 1 1 0.0000 0.4000 23 1 1 0.0000 0.2000 24 0 0 0.4000 0.2000 25 0 0 0.8000 0.2000 26 0 0 1.2000 0.2000 27 0 0 1.6000 0.2000 28 0 0 2.0000 0.2000 29 0 0 2.4000 0.2000 30 0 0 2.8000 0.2000 31 0 0 3.2000 0.2000 32 0 0 3.6000 0.2000 33 1 1 4.0000 0.2000 1 1 2 24 23 2 2 3 25 24 3 3 4 2