有限元程序课程设计

重庆大学本科学生课程设计任务书课程设计题目有限元程序设计学院资源及环境科学学院专业工程力学年级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, 5thEition, 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）有限元平衡方程 P Ka =（2）单元刚度矩阵ηξ⎰∑Ω===eni T i Ted d J DBt B H DBtdxdy BK 1||（3）单元等效结点载荷∑∑⎰⎰==Ω+=+=ni T ni Ti S TTed d J Tt N d d J ft NH Ttds Nftdxdy NP e e11||||ηξηξσ2.2单元位移插值及坐标变换（1）通过Serendipity 四边形单元格式构造插值函数。

对于4节点单元，插值函数为：)1,2,3,4(i )1)(1(41=++=ηηξξi i i N对于8节点单元，插值函数为：65112121ˆN N N N --= 65222121ˆN N N N --= 76332121ˆN N N N --=874421ˆN N N N --= )1)(1(2125ηξ--=N )1)(1(2126ξη+-=N )1)(1(2127ηξ+-=N )1)(1(2128ξη--=N其中： )1)(1(41ηηξξi i i N ++= （i=1,2,3,4）（2）位移插值i 41i i u N u ∑== i 41i i v N v ∑==（3）坐标变换i 41i i x N x ∑== i 41i i y N y ∑==2.3单元刚度的计算 （1）弹性矩阵⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡--=210001v 0v 11000200ννE D 平面应力问题：E E 0=；v v 0=平面应变问题：20v -1E E =；v-1vv 0=式中E 和v 分别为材料Young 氏模量和Poisson 比。

（2）应变矩阵[]mjiB B B B =），，，（4321i 000000=⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂=⎥⎦⎤⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂==x N y N yN x N N N x y y x LN B ii i i i ii i形状函数对局部坐标偏导：）（ηηξξi i i 141N +=∂∂，）（ξξηηi i i 141N +=∂∂形状函数对整体坐标偏导：⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡++=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∂∂∂∂=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∂∂∂∂）（）（ξξηηηξηξi i i i 1-1-141141J J N N y N x N i i i iJacobi 矩阵：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡++++=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂=∑∑∑∑∑∑∑∑========41i i i i 41i i i i 41i i i i 41i i i i 41i i i 41i ii 41i i i 41i i iy )1(41x )1(41y )1(41x )1(41y N x N y N x N ηξξηξξξηηξηηηηξξηηξξy x y x JJacobi 矩阵的逆：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡++-+-+=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂=∑∑∑∑====41i i i i 41i i i i 41i i i i 41i ii i 1-x )1(41x )1(41y )1(41y )1(41|J |1x --y |J |1ξηηηξξξηηηξξξηξηxy J其中：)y )1(41)(x )1(41(-)y )1(41)(x )1(41(|J |41i i i i 41i i i i 41i i i i 41i i i i ∑∑∑∑====++++=ξηηηξξηξξξηη（3）单元刚度矩阵计算⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎦⎤⎢⎣⎡=8887868584838281787776757473727168676665646362615857565554535251484746454443424138373635343332312827262524232221181716151413121122211211e k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k K K K K K单元刚度矩阵Gauss 积分：∑∑⎰⎰⎰=====p qn 1p n 1q q p 11-11-v |),(|),(),(w w ||dv q p q p T q p T T e J t DB B d d J DBt B DB B K eηξηξηξηξ单元刚度矩阵子矩阵：∑∑⎰⎰⎰=====p qn 1p n 1q j i q p 11-11-j i v j i ij |),(|),(),(w w ||dv q p q p T q p TTJ t DB B d d J t DB B DB B K eηξηξηξηξ （i=1，2）2.4等效结点载荷 （1）集中载荷直接施加在结点上（2）体积力产生的等效结点载荷向量[][]T4y4x 3y 3x 2y2x1y1xT4321e b P P P P P P P P P P P P P ==体积力产生的等效结点载荷Gauss 积分：∑∑⎰⎰⎰=====p qn 1p n 1q q p 11-11-T v T e b |),(|b ),(N w w ||N btdxdyN P q p T q p J t d d J t eηξηξηξ体积力产生的等效结点载荷子向量：∑∑⎰⎰⎰=====p qn 1p n 1q i q p 11-11-Ti v Ti i |),(|b ),(N w w ||N btdxdyN P q p T q p J t d d J t eηξηξηξ其中：⎥⎦⎤⎢⎣⎡=),(00),(),(N i q p q p q p N N ηξηξηξ若体积力作用方向于整体坐标系的y 轴正方向θ角度，则：|),(|cos sin ),(11q p n p n q q p i q p iy ix J t N w w P P p q ηξθθηξ∑∑==⎥⎦⎤⎢⎣⎡-=⎥⎦⎤⎢⎣⎡2.5单元应变和应力计算（1）单元应变e Ba =ε单元中任意点的应变：[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡43214321xy y x ),(),(),(a a a a B B B B y x y x y x γεε 单元中任意点的应变：[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡43214321xy y x ),(),(),(),(),(),(),(a a a a B B B B y x y x y x q q q q q q q q q q q q q q ηξηξηξηξγεε 其中),(q q y x 可作为Gauss 积分点或单元结点的整体坐标，其于局部坐标间的关系为：i 41i i q x ),(N x ∑==q q ηξ,i 41i i q y ),(N y ∑==q q ηξ(2)单元应力e DBa =σ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡y)(x,y)(x,y)(x,D y)(x,),(),(xy y x xy y x y x y x γεετσσ三、程序结构（1）有限元程序系统的组成及分析过程我们本次设计主要针对二维弹性力学静力学问题。