有限元程序课程设计

王顺金：有限元程序课程设计1重庆大学本科学生课程设计任务书课程设计题目有限元程序设计学院资源及环境科学学院专业工程力学年级2010级已知参数和设计要求：1.独立完成有限元程序设计。2.独立选择计算算例，并能通过算例判断程序的正确性。3.独立完成程序设计报告，报告内容包括理论公式、程序框图、程序本体、计算算例，算例结果分析、结论等。学生应完成的工作：1.复习掌握有限单元法的基本原理。2.掌握弹性力学平面问题3节点三角形单元或4节点等参单元有限元方法的计算流程，以及单元刚度矩阵、等效节点载荷、节点应变、节点应力和高斯积分等的计算公式。3.用Fortran语言编写弹性力学平面问题3节点三角形单元或4节点等参单元的有限元程序。4.在VisualFortran程序集成开发环境中完成有限元程序的编辑和调试工作。5.利用编写的有限元程序，计算算例，分析计算结果。6.撰写课程设计报告。目前资料收集情况（含指定参考资料）：1.王勖成，有限单元法，北京：高等教育出版社，2002。2.O.C.Zienkiewicz,R.L.Taylor,FiniteElementMethod,5thEition,McGraw-HallBookCompanyLimited,2000。3.张汝清，董明，结构计算程序设计，重庆：重庆大学出版社，1988。课程设计的工作计划：1.第1周星期一上午：教师讲解程序设计方法，程序设计要求和任务安排。2.第1周星期一至星期二完成程序框图设计。3.第1周星期三至第2周星期四完成程序设计。4.第2周星期五完成课程设计报告。任务下达日期2013年6月6日完成日期2013年07月03日指导教师（签名）学生（签名）王顺金：有限元程序课程设计2一、前言有限单元法是在当今工程分析中获得最广泛应用的数值计算方法，由于其通用性和有效性，受到工程技术界的高度重视。伴随着计算机科学和技术的快速发展，现已成为计算机辅助设计和计算机辅助制造的重要组成部分。由于有限元法是通过计算机实现的，因此有限元程序的编制及相关软件的研发就变得尤为重要。从二十世纪五十年代以来，有限元软件的发展按目的和用途可分为专用软件和大型通用商业软件，而且软件往往集成了网格自动划分，结果分析和显示等前后处理功能，而且随着时间的发展，大型通用商业软件的功能由线性扩展到非线性，由结构扩展到非结构等等，这一系列强大功能的实现与运用都要求我们对有限元法的基础理论知识有较为清楚的认识以及对程序编写的基本能力有较好掌握。1.元分析程序的理论基础作为弹性力学微分方程的等效积分形式，虚位移原理和虚应原理分别是平衡方程与力的边界条件和几何方程与位移边界条件的等效积分形式。将物理方程引入虚位移原理和虚应力原理可以分别导出最小位能原理和最小余能原理，它们本质上和等效积分的伽辽金“弱”形式相一致，这是建立弹性力学有限原方程一般表达格式的理论基础。对于通过弹性力学变分原理建立的弹性力学问题有限元方法，其未知场变量是位移，以节点位移为基本未知量，并以最小位能原理为基础建立的有限单元为位移元。弹性力学平面问题有限元分析表达格式的建立步骤一般为：a．的模型进行单元离散，一般采用三角形单元或四边形单元，不过由于四边形单元具有更高精度，应用更为普遍，一般而言一次或二次单元已足以满足精度要求。离散后对单元进行编号，并且给定单元各节点的整体编码以及局部坐标系下按逆时针局部编码。b．坐标系下对各单元构造形函得到由单元各节点位移表示的单元位移形式，进而得到单元刚度矩阵，利用等参元性质和雅阁比矩阵进行组集，建立整体刚度矩阵。c．等效结点载荷列阵并组集成结构节点等效载荷列阵得到单元各节点位移与结构节点等效载荷列阵的线代方程组，求解得到位移，进而可得应力应变。根据以上一般步骤，编制相应的计算机程序并采用数值积分方法处理有关数学计算就可以得到完整的弹性力学问题有限元分析程序。王顺金：有限元程序课程设计32.问题有限元分析教学程序本程序可对二维弹性力学问题行有限元分析计算。计算采用的单元形式为四边形四节点单元或者四边形八节点单元，对于对称和非对称矩阵采用变带宽存储方法，最终输出结果为单元各节点的位移。二、平面4、8节点有限元公式及计算原理2.1基本公式（1）有限元平衡方程PKa（2）单元刚度矩阵eniTiTeddJDBtBHDBtdxdyBK1||（3）单元等效结点载荷niTniTiSTTeddJTtNddJftNHTtdsNftdxdyNPee11||||2.2单元位移插值及坐标变换（1）通过Serendipity四边形单元格式构造插值函数。对于4节点单元，插值函数为：)1,2,3,4(i)1)(1(41iiiN对于8节点单元，插值函数为：65112121ˆNNNN65222121ˆNNNN76332121ˆNNNN874421ˆNNNN)1)(1(2125N)1)(1(2126N)1)(1(2127N)1)(1(2128N其中：)1)(1(41iiiN（i=1,2,3,4）（2）位移插值i41iiuNui41iivNv（3）坐标变换王顺金：有限元程序课程设计4i41iixNxi41iiyNy2.3单元刚度的计算（1）弹性矩阵210001v0v11000200ED平面应力问题：EE0；vv0平面应变问题：20v-1EE；v-1vv0式中E和v分别为材料Young氏模量和Poisson比。（2）应变矩阵mjiBBBB），，，（4321i000000xNyNyNxNNNxyyxLNBiiiiiiii形状函数对局部坐标偏导：）（iii141N，）（iii141N形状函数对整体坐标偏导：）（）（iiii1-1-141141JJNNyNxNiiiiJacobi矩阵：王顺金：有限元程序课程设计541iiii41iiii41iiii41iiii41iii41iii41iii41iiiy)1(41x)1(41y)1(41x)1(41yNxNyNxNyxyxJJacobi矩阵的逆：41iiii41iiii41iiii41iiii1-x)1(41x)1(41y)1(41y)1(41|J|1x--y|J|1xyJ其中：)y)1(41)(x)1(41(-)y)1(41)(x)1(41(|J|41iiii41iiii41iiii41iiii（3）单元刚度矩阵计算8887868584838281787776757473727168676665646362615857565554535251484746454443424138373635343332312827262524232221181716151413121122211211ekkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkKKKKK单元刚度矩阵Gauss积分：王顺金：有限元程序课程设计6pqn1pn1qqp11-11-v|),(|),(),(ww||dvqpqpTqpTTeJtDBBddJDBtBDBBKe单元刚度矩阵子矩阵：pqn1pn1qjiqp11-11-jivjiij|),(|),(),(ww||dvqpqpTqpTTJtDBBddJtDBBDBBKe（i=1，2）2.4等效结点载荷（1）集中载荷直接施加在结点上（2）体积力产生的等效结点载荷向量T4y4x3y3x2y2x1y1xT4321ebPPPPPPPPPPPPP体积力产生的等效结点载荷Gauss积分：pqn1pn1qqp11-11-TvTeb|),(|b),(Nww||NbtdxdyNPqpTqpJtddJte体积力产生的等效结点载荷子向量：pqn1pn1qiqp11-11-TivTii|),(|b),(Nww||NbtdxdyNPqpTqpJtddJte王顺金：有限元程序课程设计7其中：),(00),(),(NiqpqpqpNN若体积力作用方向于整体坐标系的y轴正方向角度，则：|),(|cossin),(11qpnpnqqpiqpiyixJtNwwPPpq2.5单元应变和应力计算（1）单元应变eBa单元中任意点的应变：43214321xyyx),(),(),(aaaaBBBByxyxyx单元中任意点的应变：43214321xyyx),(),(),(),(),(),(),(aaaaBBBByxyxyxqqqqqqqqqqqqqq其中),(qqyx可作为Gauss积分点或单元结点的整体坐标，其于局部坐标间的关系为：i41iiqx),(Nxqq,i41iiqy),(Nyqq(2)单元应力eDBay)(x,y)(x,y)(x,Dy)(x,),(),(xyyxxyyxyxyx三、程序结构（1）有限元程序系统的组成及分析过程我们本次设计主要针对二维弹性力学静力学问题。有限元程序系统通常包括处理、有限元程序本体和后处理三部分。王顺金：有限元程序课程设计8有限元软件系统的组成及分析过程如下图所示：CAD模型、材料参数、边界条件、分析类型定义、网格划分等。有限元分析本体程序前处理后处理核心：各种计算方法的实现。以图形、曲线和表格等方式表达、分析计算结果。（2）有限元计算程序框图输入离散模型数据型数据计算单元刚度阵组集结构刚度矩阵计算单元等效结点载荷组集结构结点载荷列阵引入位移边界条件求解线性方程组其它辅助计算输出结果结束单元循环形成K形成P消除K的奇异求解Ka=P得结点位移a计算应力、应变等王顺金：有限元程序课程设计9（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个坐标I