《有限元程序课程设计》由会员分享，可在线阅读，更多相关《有限元程序课程设计（51页珍藏版）》请在金锄头文库上搜索。
1、王顺金：有限元程序课程设计 重庆大学本科学生课程设计任务书课程设计题目有限元程序设计学院资源及环境科学学院专业工程力学年级2010级已知参数和设计要求：1.独立完成有限元程序设计。2.独立选择计算算例，并能通过算例判断程序的正确性。3.独立完成程序设计报告，报告内容包括理论公式、程序框图、程序本体、计算算例，算例结果分析、结论等。学生应完成的工作：1. 复习掌握有限单元法的基本原理。2. 掌握弹性力学平面问题3节点三角形单元或4节点等参单元有限元方法的计算流程，以及单元刚度矩阵、等效节点载荷、节点应变、节点应力和高斯积分等的计算公式。3. 用Fortran语言编写弹性力学平面问题3节点三角形单
2、元或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. 第
3、1周星期一上午：教师讲解程序设计方法，程序设计要求和任务安排。2. 第1周星期一至星期二完成程序框图设计。3. 第1周星期三至第2周星期四完成程序设计。4. 第2周星期五完成课程设计报告。任务下达日期 2013 年 6 月 6 日完成日期 2013 年 07 月 03 日指导教师 （签名）学 生 （签名）一、前言有限单元法是在当今工程分析中获得最广泛应用的数值计算方法，由于其通用性和有效性，受到工程技术界的高度重视。伴随着计算机科学和技术的快速发展，现已成为计算机辅助设计和计算机辅助制造的重要组成部分。由于有限元法是通过计算机实现的，因此有限元程序的编制及相关软件的研发就变得尤为重要。从二十世
4、纪五十年代以来，有限元软件的发展按目的和用途可分为专用软件和大型通用商业软件，而且软件往往集成了网格自动划分，结果分析和显示等前后处理功能，而且随着时间的发展，大型通用商业软件的功能由线性扩展到非线性，由结构扩展到非结构等等，这一系列强大功能的实现与运用都要求我们对有限元法的基础理论知识有较为清楚的认识以及对程序编写的基本能力有较好掌握。1.元分析程序的理论基础作为弹性力学微分方程的等效积分形式，虚位移原理和虚应原理分别是平衡方程与力的边界条件和几何方程与位移边界条件的等效积分形式。将物理方程引入虚位移原理和虚应力原理可以分别导出最小位能原理和最小余能原理，它们本质上和等效积分的伽辽金“弱”形
5、式相一致，这是建立弹性力学有限原方程一般表达格式的理论基础。对于通过弹性力学变分原理建立的弹性力学问题有限元方法，其未知场变量是位移，以节点位移为基本未知量，并以最小位能原理为基础建立的有限单元为位移元。弹性力学平面问题有限元分析表达格式的建立步骤一般为：a的模型进行单元离散，一般采用三角形单元或四边形单元，不过由于四边形单元具有更高精度，应用更为普遍，一般而言一次或二次单元已足以满足精度要求。离散后对单元进行编号，并且给定单元各节点的整体编码以及局部坐标系下按逆时针局部编码。b坐标系下对各单元构造形函得到由单元各节点位移表示的单元位移形式，进而得到单元刚度矩阵，利用等参元性质和雅阁比矩阵进行
6、组集，建立整体刚度矩阵。c等效结点载荷列阵并组集成结构节点等效载荷列阵得到单元各节点位移与结构节点等效载荷列阵的线代方程组，求解得到位移，进而可得应力应变。 根据以上一般步骤，编制相应的计算机程序并采用数值积分方法处理有关数学计算就可以得到完整的弹性力学问题有限元分析程序。2. 问题有限元分析教学程序本程序可对二维弹性力学问题行有限元分析计算。计算采用的单元形式为四边形四节点单元或者四边形八节点单元，对于对称和非对称矩阵采用变带宽存储方法，最终输出结果为单元各节点的位移。二、平面4、8节点有限元公式及计算原理2.1基本公式（1） 有限元平衡方程（2） 单元刚度矩阵（3）单元等效结点载荷2.2单
7、元位移插值及坐标变换（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） 体积力产生的
8、等效结点载荷向量体积力产生的等效结点载荷Gauss积分：体积力产生的等效结点载荷子向量：其中：若体积力作用方向于整体坐标系的y轴正方向角度，则：2.5单元应变和应力计算（1） 单元应变单元中任意点的应变：单元中任意点的应变：其中可作为Gauss积分点或单元结点的整体坐标，其于局部坐标间的关系为：,(2) 单元应力三、程序结构（1）有限元程序系统的组成及分析过程我们本次设计主要针对二维弹性力学静力学问题。有限元程序系统通常包括处理、有限元程序本体和后处理三部分。有限元软件系统的组成及分析过程如下图所示：CAD模型、材料参数、边界条件、分析类型定义、网格划分等。有限元分析本体程序前 处 理后 处
9、理核心：各种计算方法的实现。以图形、曲线和表格等方式表达、分析计算结果。（2）有限元计算程序框图 输入离散模型数据型数据计算单元刚度阵组集结构刚度矩阵计算单元等效结点载荷组集结构结点载荷列阵引入位移边界条件求解线性方程组其它辅助计算输出结果结 束单元循环形成 K形成 P消除 K 的奇异求解 Ka=P得结点位移 a计算应力、应变等（3）有限元程序主要编程技巧a.结构模块化主程序：组织调用子程序，过程清晰，不涉及复杂计算。子模块：实现特定功能，可用多个子程序组成。子程序：实现简单功能，便于调用和扩展。b.覆盖技术在主程序中定义一个整数组和一个实型数组，分别存放整型和实型边界数组。设计动态数组表。动
10、态数组覆盖技术。（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节点的全局节点
11、编号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）程序本体（见附录）四、 程序验
12、证算例算例：平面板受轴向拉伸长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
