有限元程序课程设计

任务下达日期 2013 年 指导教师
6
月 6 日 （签名）
完成日期 2013 年 07 月 03 日 学 名） 生 （签
1
王顺金：有限元程序课程设计
一、前言
有限单元法是在当今工程分析中获得最广泛应用的数值计算方法，由于其 通用性和有效性，受到工程技术界的高度重视。伴随着计算机科学和技术的快 速发展，现已成为计算机辅助设计和计算机辅助制造的重要组成部分。 由于有限元法是通过计算机实现的，因此有限元程序的编制及相关软件的 研发就变得尤为重要。从二十世纪五十年代以来，有限元软件的发展按目的和 用途可分为专用软件和大型通用商业软件，而且软件往往集成了网格自动划 分，结果分析和显示等前后处理功能，而且随着时间的发展，大型通用商业软 件的功能由线性扩展到非线性，由结构扩展到非结构等等，这一系列强大功能 的实现与运用都要求我们对有限元法的基础理论知识有较为清楚的认识以及对 程序编写的基本能力有较好掌握。 1.元分析程序的理论基础 作为弹性力学微分方程的等效积分形式，虚位移原理和虚应 原理分别是平衡方程与力的边界条件和几何方程与位移边界条件的等效积 分形式。将物理方程引入虚位移原理和虚应力原理可以分别导出最小位能原理 和最小余能原理，它们本质上和等效积分的伽辽金“弱”形式相一致，这是建 立弹性力学有限原方程一般表达格式的理论基础。 对于通过弹性力学变分原理建立的弹性力学问题有限元方法，其未知场变 量是位移，以节点位移为基本未知量，并以最小位能原理为基础建立的有限单 元为位移元。弹性力学平面问题有限元分析表达格式的建立步骤一般为： a．的模型进行单元离散，一般采用三角形单元或四边形单元，不过由于 四边形单元具有更高精度，应用更为普遍，一般而言一次或二次单元已足以满 足精度要求。离散后对单元进行编号，并且给定单元各节点的整体编码以及局 部坐标系下按逆时针局部编码。 b．坐标系下对各单元构造形函得到由单元各节点位移表示的单元位移形 式，进而得到单元刚度矩阵，利用等参元性质和雅阁比矩阵进行组集，建立整 体刚度矩阵。 c．等效结点载荷列阵并组集成结构节点等效载荷列阵 得到单元各节点位移与结构节点等效载荷列阵的线代方程组，求解得到位 移，进而可得应力应变。 根据以上一般步骤，编制相应的计算机程序并采用数值积分方法处理有关 数学计算就可以得到完整的弹性力学问题有限元分析程序。

Bj
Bm

N i 0 x 0 N i 0 （i 1， 2， 3， 4） Ni y N N i i x y
形状函数对局部坐标偏导：
N i 1 N i 1 （ （ ， i 1 i） i 1 i） 4 4
E v v 2 ； 0 1- v 1- v
式中 E 和 v 分别为材料 Young 氏模量和 Poisson 比。 （2）应变矩阵
B Bi
0 x Ni Bi LNi 0 y 0 y x
Pbe P1 P1x

P2 P1y
P3 P2x
P4 P2y
T
P3x
P3y
P4x
P4y

T
体积力产生的wenku.baidu.com效结点载荷 Gauss 积分：
Pbe N T btdxdy
ve

1 1
-1 -1

N T t | J | dd
w p w q N( p , q )T bt | J ( p , q ) |
2
王顺金：有限元程序课程设计
2.问题有限元分析教学程序 本程序可对二维弹性力学问题行有限元分析计算。计算采用的单元形式为 四边形四节点单元或者四边形八节点单元，对于对称和非对称矩阵采用变带宽 存储方法，最终输出结果为单元各节点的位移。
二、平面 4、8 节点有限元公式及计算原理
2.1 基本公式 （1）有限元平衡方程
王顺金：有限元程序课程设计
重庆大学本科学生课程设计任务书
课程设计题目 学院 有限元程序设计 专业 工程力学 年级 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 周星期五完成课程设计报告。
p 1 q 1
np
nq
体积力产生的等效结点载荷子向量：
Pi N i btdxdy
T ve

1 1
-1 -1

N i t | J | dd
T
w p w q N i ( p , q )T bt | J ( p , q ) |
p 1 q 1
np
nq
6
王顺金：有限元程序课程设计
形状函数对整体坐标偏导：
N i N i 1 （ i 1 i） x -1 -1 4 N J N J 1 i i （ 1 ） i i 4 y
Ka P
（2）单元刚度矩阵
K
e

e

B T DBtdxdy
H
i 1
n
n
i
B T DBt | J | d d
（3）单元等效结点载荷
P
e



N
e
T
ftdxdy

S

N T Ttds
e
H iN
i 1
T
ft | J | d d

n
N T Tt | J | d d
4 1 (1 i )i y i 1 4 i 14 1 |J| (1 i )i x i i 1 4
其中：
4 4 4 4 1 1 1 1 | J | ( (1 i ) i x i )( (1 i ) i y i ) - ( (1 i ) i x i )( (1 i ) i y i ) i 1 4 i 1 4 i 1 4 i 1 4
i 1
2.2 单元位移插值及坐标变换 （1）通过 Serendipity 四边形单元格式构造插值函数。对于 4 节点单元，插 值函数为：
Ni
1 (1 i )(1 i ) (i 1,2,3,4 ) 4
对于 8 节点单元，插值函数为：
ˆ 1N 1N N1 N 1 5 6 2 2 ˆ 1N 1N N3 N 3 6 7 2 2 ˆ 1N N N4 N 4 7 8 2 1 N 6 (1 2 )(1 ) 2 1 N 7 (1 2 )(1 ) 2 ˆ 1N 1N N2 N 2 5 6 2 2
（3）单元刚度矩阵计算
K K e 11 K 21
k11 k12 k13 k14 k15 k16 k17 k18 k 21 k 22 k 23 k 24 k 25 k 26 k 27 k 28 k 31 k 32 k 33 k 34 k 35 k 36 k 37 k 38 K12 k 41 k 42 k 43 k 44 k 45 k 46 k 47 k 48 K 22 k 51 k 52 k 53 k 54 k 55 k 56 k 57 k 58 k 61 k 62 k 63 k 64 k 65 k 66 k 67 k 68 k 71 k 72 k 73 k 74 k 75 k 76 k 77 k 78 k 81 k 82 k 83 k 84 k 85 k 86 k 87 k 88
单元刚度矩阵 Gauss 积分：
5
王顺金：有限元程序课程设计
K e B T DBdv
ve

1 1
-1 -1

B T DBt | J | dd
w p w q B ( p , q )T DB( p , q )t | J ( p , q ) |
p 1 q 1
np
nq
其中：
0 N ( p , q ) N i ( p , q ) 0 N ( p , q )
若体积力作用方向于整体坐标系的 y 轴正方向 角度，则：
单元刚度矩阵子矩阵：
K ij Bi DBjdv
T ve

1 1
-1 -1

Bi DBjt | J | dd
T
（i=1，2）
w p w q Bi ( p , q )T DBj ( p , q )t | J ( p , q ) |
p 1 q 1
np
nq
2.4 等效结点载荷 （1）集中载荷 直接施加在结点上 （2）体积力产生的等效结点载荷向量
x Ni x i
i 1
4
y N i yi
i 1
4
2.3 单元刚度的计算 （1）弹性矩阵
1 E0 D v 2 0 1 0 0
平面应力问题： E 0
v0 1 0
0 0 1 0 2
E ； v0
v
平面应变问题： E 0
Jacobi 矩阵：
4
王顺金：有限元程序课程设计
x J x
y y
4 N i x i 1 i 4 N i x i i 1
yi i 1 4 N i y i i 1 4 4 1 1 4 (1 i ) i x i 4 (1 i ) i y i 1 i 1 i 4 4 1 1 (1 ) x ( 1 ) y i i i i i i 4 4 i 1 i 1

4
N i
Jacobi 矩阵的逆：
y 1 J -1 | J | - x y x
4 1 (1 i ) i y i i 1 4 4 1 (1 i ) i x i 4 i 1
N5
1 (1 2 )(1 ) 2
N8
1 (1 2 )(1 ) 2
其中： N i （2）位移插值
1 (1 i )(1 i ) 4
（i=1,2,3,4）
u Ni u i
i 1
4
v N i vi
i 1
4
（3）坐标变换
3
王顺金：有限元程序课程设计