有限单元法编著：丁科
2006年有限单元法
(Finite Element
Method)
有限单元法第１章 绪论有限单元法
1.1 概述问题的求解方法：解析法,数值法,（ 半解析法 ）
解析法,问题的求解可以得到具体的数学表达式,可以计算任意位置的问题的精确解答 。 如悬臂梁求挠度等 。
数值法,得不到具体表达式,只能得到某些离散点处的近似值 。 该方法包括：有限差分,有限元法,边界元法
,变分法等有限单元法有限元法,将结构物看成由有限个划分的单元组成的整体,以单元结点上的值作为整个单元的平均值 。 它是一种化整为零,集零为整,化未知为已知的方法 。 不同的学科,所求解的参数不同 。 在结构力学中,主要有以下三种：
● 位移型：以结点位移为未知量 。
● 力型：以结点力为未知量 。
● 混合型：某些地方以结点位移为未知量,另外一些以结点力为未知量 。
我们主要就,位移型,有限元进行讲解 。
有限单元法有限元方法的实质,将复杂的连续体划分为有限多个简单的单元体,化无限自由度问题为有限自由度问题,将连续场函数的 （ 偏 ） 微分方程的求解问题转化成有限个参数的代数方程组的求解问题 。
有限元方法的基本思想,先化整为零、再积零为整，
也就是把一个连续体人为的分割成有限个单元，即把一个结构看成由若干通过结点项链的单元组成的整体，先进行单元分析，然后再把这些单元组合起来代表原来的结构进行整体分析。
有限单元法
1.2 有限元法的发展
1) 早在公元 3世纪的时候,我国数学家刘徽提出的用割元法求圆周长的方法即是有限元基本思想的体现 。
2) 20世纪 40年代,麦克亨利 （ McHenry）,雷尼柯夫 （
Hrenikoff）,纽马克 （ Newmark） 等首次提出用框架方法求解力学问题,用简单弹性杆排列代替连续体的各个小部分,能够得到连续问题的相当好的解答 。
3) 1943年,柯兰特 （ Courant） 第一次假设饶曲函数在一个划分的三角形单元集合体的每个单元上为简单线性函数 。
4) 1955年，德国斯图加特大学的 J.H,Argyris教授发表了一组能量原理与矩阵分析的论文，奠定了有限元方法的理论基础。
有限单元法
5) 第一个尝试,1956年,特纳 （ Turner）,克拉夫 （ Clough）
等将刚架分析中的位移法扩展到弹性力学平面问题,并用于飞机的结构分析和设计,系统研究了离散杆,梁,三角形的单元刚度表达式,并求得了平面应力问题的正确解答 。
6) 第一次提出,1960年,克拉夫 （ Clough）,建立在虚位移原理或最小势能原理的基础上 。 在处理剖面弹性问题时,
第一次提出并使用,有限元方法,的名称
7) 1972年,Oden出版了第一本处理非线性连续体的专著 。
8) 我国科技工作者：如胡昌海提出了广义变分原理,钱伟长最先研究了拉格朗日乘子法与广义变分原理之间的关系,
冯康研究了有限元方法的精度和收敛性问题,钱令希研究了余能原理等 。
9) 基于变分原理的有限元法得到发展。
有限单元法
1.3 有限元法的分析过程
1,结构物的离散将待分析的结构物从几何上用线或面划分为有限个单元,
其中单元的大小和数目根据计算精度的要求和计算机容量来决定 。 其步骤：
● 建立单元
● 对单元和结点编号
● 准备必需的数据信息
● 建立坐标系有限单元法
x
y
x
y
1 2
3 4
5 6
① ②
③
④ ⑤
⑥
图 1.1
如图 1.1所示,可以将杆系结构分成 6
个单元,这样划分以后,共有 6个结点 。
如图 1.2所示,纵向均匀受拉的带圆孔的薄板,
根据对称性，取其中一部分分析，将其划分为三角形单元。
图 1.1
图 1.2
有限单元法
2,确定单元的位移模式将单元中任意一点的位移近似地表示成单元结点位移的函数,即位移模式或位移函数,用 或 表示,写成：
这里,— 单元中任意一点的位移矩阵,
— 形函数矩阵
— 单元结点位移矩阵 。
位移函数的假设合理与否，直接影响到分析的计算精度、
效率、可靠性。
d d?
e?dN δ
d
N
eδ
有限单元法
3,单元特性分析
（ 1） 几何方程：应变与位移之间的关系这里,— 单元中任意一点的应变矩阵,
— 变形矩阵或应变矩阵,
（ 2） 物理方程：应力与应变之间的关系 （ Hooke定律 ）
这里,— 单元中任意一点的应力矩阵,
— 弹性矩阵,由单元材料的弹性常数确定,(弹性模量 )
— 应力矩阵。
e?ε B δ
ε
B
eeζ D B δ S δ
ζ
D
S
有限单元法
（ 3） 利用虚位移原理或最小势能原理建立单元刚度方程这里,— 单元结点力矩阵
— 单元等效荷载矩阵
— 单元刚度矩阵，
e e e eEk δ F F
eF
eEF
ek
有限单元法
4,建立表示整个结构结点平衡的方程组这里,－整体刚度矩阵
－整体结点位移矩阵
－直接结点荷载
－等效结点荷载
－整体综合结点荷载矩阵
deK Δ P P P
K
Δ
dP
eP
P
有限单元法
1.4 常用软件
（ 1） 软件有限元软件可以分为通用软件和专用软件两类。通用软件适应性广，规格规范，输入方法简单，有比较成熟齐全的单元库，大多提供二次开发的接口。但通用程序功能再强，
对于一些比较专业的问题，尤其是处于研究阶段的内容，往往显得无能为力，因此，针对某些特定领域、特定问题开发的专用软件，在解决专有问题时显得更为有效。目前常用的有限元软件有,ANSYS,MARC,ABQUS,NASTRAN,ADINA、
ALGOR,SAP,STRAND,FEPG等。
有限单元法
（ 2） 程序设计低级语言：机器语言,汇编语言等高级语言,BASIC,C,C++,FORTRAN,JAVA,LISP，
Matlab,Mathmatic等目前所有有限元软件都是用 FORTRAN语言编写的 。 程序设计步骤：
（ 1） 提出问题,拟定解决方案；
（ 2） 构造数学模型；
（ 3） 画出程序流程图；
（ 4） 用选定的算法语言编写程序；
（ 5） 编译调试程序；
（ 6） 试验验证程序；
（ 7） 打包发行；
（ 8）软件维护。
有限单元法第 2章杆系结构的有限元分析有限单元法
2.1 概述所谓 杆件 是指从构造上来说其长度远大于其截面尺寸的一维构件 。 在结构力学上我们通常将承受轴力或扭矩的杆件称为杆,而将承受横向力和弯矩的杆件称为梁 。 在有限单元法中这两种情况的单元分别称为杆单元和梁单元 。 为方便起见,本书都称之为杆单元 。
杆系结构是最简单的一类结构，也是我们在工程上最常见的一类结构。本章以此类结构为基础介绍有限单元法的分析过程。
有限单元法第一步,对结构物进行离散化,划分为有限个单元 。
第二步,对各结点和单元进行编码 。
第三步,建立整体坐标系和各单元的局部坐标系 。
第四步,对已知参数进行准备和整理 。
第五步,对结点位移进行编码,注意前处理法与后处理法的区别 。
第六步,进行单元分析,形成单元刚度矩阵 。
第七步,进行整体分析,形成整体刚度矩阵 。
第八步,引入边界条件 。 边界条件的引入可以使问题具有解的唯一性 。
第九步,求解方程组,计算结构的整体结点位移 。
第十步，求单元内力，对计算成果进行整理、分析，用表格、图示出所需的位移及应力。
有限单元法
1
2
3 4
5
6 1 2
3 4 5
3
1
2
3
( 0 0 0 )
( 0 0 1 )
( 2 3 4 )
4 ( 5 6 7 )
6 ( 1 1 1 2 1 3 )
5 ( 8 9 1 0 )
6
54
21
( 1 2 3 )
2 ( 4 5 6 )
( 7 8 9 )
4 ( 1 0 1 1 1 2 )
6 ( 1 6 1 7 1 8 )5 ( 1 3 1 4 1 5 )
3
1
3
6
54
21
图 2.1 弯曲杆件系统 图 2.2 截面连续变化杆件系统图 2.3 单元位移编码前处理法 后处理法有限单元法
2.2 局部坐标系中的杆单元分析
2.2.1 拉压杆单元
q ( x )?F
i?F j?u i?u j
y
x
i j
图 2.4 拉压杆单元示意图设杆单元长度为,横截面面积为,单元材料的弹性模量为,在局部坐标系中杆端荷载分别为 和,杆端位移分别为 和,单元上的轴向分布荷载为 。
l A
E iF jF
iu ju ()qx
有限单元法
① 单元位移模式 。 用结点位移表示单元上任意截面的位移 。 对拉压杆单元,可以取其位移为一次多项式,即：
由位移的边界条件：
可得系数,为：
这样,任意截面的位移 为：
用矩阵表示为：
其中
u
()u x a b x
(0) iuu? () ju l u?
a b
iau? jiuub
l
( ) (1 ) ijxxu x u ull
i
i i j j i j
j
uu N u N u N N
u
N δ?
1i xN l j xN l?
（ 2-1）
（ 2-2）
有限单元法
② 进行应力,应变分析 。 根据材料力学中应变的定义,
有：
这里为应变矩阵 。 由虎克定律,其应力为：
11
ij
d u d BB
d x d x l l?
N δ δ δ B δ
11
ll
B
EE B δ? （ 2-4）
（ 2-3）
有限单元法
③ 求单元刚度矩阵 。 这里考虑利用虚位移原理求单元刚度矩阵,设杆端 i,j分别产生虚位移,,则由此引起的杆轴任意截面的虚位移为：
对应的虚应变为：
根据虚位移原理虚功方程,有：
将上式整理得：
iu? ju?
Tiiu u uNN δ?
B δ?
0
0
0
()
l
T
d
l
l
TT
W q x d x W
A d x
E A d x
F δ N δ
δ BB δ
外 变
00() TllT T Td q x d x E A d xFN δ δ BB δ
（ 2-5）
（ 2-6）
有限单元法式中,为局部坐标系下单元结点荷载矩阵 。 设：
则可以得到拉压杆单元的单元刚度方程为：
这里 为局部坐标系下的单元刚度矩阵,为局部坐标系下等效结点荷载矩阵,但值得指出的是：分布荷载 中可以包含集中荷载 。 根据定义,可以进一步求得单元刚度矩阵为：
T
d i jFFF?
0 ()
l T
E q x d xFN?
0
l T E A d xk B B?
dEF F k δ
k? EF?
()qx
11
11
EA
l
k? （ 2-10）
（ 2-7）
（ 2-8）
（ 2-9）
有限单元法
2.2.2 扭转杆单元
m ( x )?M
i
M j
y
x
i j
j? i
图 2.4 扭转杆单元示意图设扭转杆单元的长度为,截面惯性矩为,剪切模量为,杆端扭矩分别为,,杆端扭转角分别为,,
单元上的分布荷载集度为,则任意截面的扭转角为：
l
I
G
iM jM i? j?
()mx
( 1 ) ijxxll N δ?
有限单元法式中,为局部坐标系下扭转杆单元的结点位移矩阵 。
由材料力学可知,截面扭矩为：
式中：
我们利用极小势能原理来进行单元分析,杆单元的势能用泛函表示为：
Tijδ?
dM GI GIdx B δ?
11d
d x l l
NB
00
00
1 ( ) ( )
2
1 ( ( ) )
2
llTT
pd
llT T T
d
dM dx m x dx
dx
GI dx m x dx
F δ
δ BB δ NF δ
有限单元法这里 为局部坐标系下扭转杆单元的结点荷载矩阵 。 由极小势能原理,取上述泛函的变分,可得：
或者写为：
设：
可得扭转杆单元的单元刚度方程为：
可以看到,其形式与拉压杆单元的单元刚度方程完全一致 。 同样,由上式可以进一步求得其局部坐标系下得单元刚度矩阵为
：
T
d i jMMF?
0p
00 ()
llT T T
dG I d x m x d xδ B B N F
00( ) ( )
llTT
dG I d x m x d xBB δ NF
0
l T G I d xk B B? 0 ()
l T
E m x d xFN?
dEF F k δ
11
11
GI
l
k?
有限单元法
2.2.3 只计弯曲的杆单元
q ( x )
M
i
M
j
F
y i
F
y j?v i
v j
y
x
m ( x )
i j
j? i
设杆单元的长度为,截面惯性矩为,弹性模量为,杆端剪力为,,杆端弯矩分别为,,杆端横向位移为,,杆端扭转角分别为,,在单元上分布有荷载集度为 的竖向分布荷载和集度为 的分布力偶 。 则结点位移矩阵和结点荷载矩阵分别为：
l I E
yiF yjF iM jM
iv jv i? j?
()qx ()mx
有限单元法取挠曲线方程为 的三次多项式,即单元上任意一点的挠度为：
根据单元的位移边界条件：
时：,
时：,
Ti i j jvvδ?
Td y i i y j jF M F MF?
x
23v a b x c x d x
0x? ivv? idvdx
xl? jvv?
j
dv
dx
2 2 2
3 2 3 2
1 2 3 1
2 1 2 1
i
j
i i j j
i i j j
av
bv
c v v
l l l l
d v v
l l l l
可以得到式中的待定系数：
有限单元法将系数 a,b,c,d代入式,并将挠曲线方程用矩阵形式表示为：
式中 为形函数矩阵,其中：
为平面弯曲单元的形函数。
23
22
3 2 3 2
1 0 0 0
0 1 0 0
3 2 3 1
1
2 1 2 1
i
i
j
j
v
v x x x
vl l l l
l l l l
N δ
1 2 3 4N N N N?N
23
1 23
2
2 2
23
3 23
23
4 2
32
1
2
(1 )
32
xx
N
ll
xx
Nx
ll
xx
N
ll
xx
N
ll
有限单元法根据式 （ 2-19） 确定的单元位移场,可得单元上某一点得曲率为：
截面的弯矩为：
这里：
为平面弯曲杆单元的应变矩阵 。
根据虚位移原理 。 有：
22
22
d v d
d x d x
N δ B δ
TTM E I E I E IB δ δ B
2
2 2 3 2 2 3 2
6 1 2 4 6 6 1 2 2 6d x x x x
d x l l l l l l l l
NB
00
0
( ) ( )
ll T
d
lTT
dW q x d x m x d x W
dx
E I d x
NNF δ
δ BB δ
外 变
有限单元法则平面弯曲杆单元的单元刚度方程为：
其中的单元刚度矩阵可由式 （ 2-23） 求得为：
00( ) ( )
ll
E
dq x d x m x d xdx NFN?
0
l T E I d xk B B?
dEF F k δ
22
3
22
1 2 6 1 2 6
6 4 6 2
1 2 6 1 2 6
6 2 6 4
ll
l l l lEI
lll
l l l l
k?
记：
（ 2-23）
有限单元法
2.2.4 平面一般杆单元
M
i
M
j
F
x i
F
y i
F
x j
F
y j
u i
v i
u j
v j
y
x
j
i
杆单元的长度为,截面面积为,截面惯性矩为,
弹性模量为,单元的,端各有三个力为,,
和,,,其对应的位移为,,和,,,
建立如图 2.7所示的局部坐标系，各物理量的正向如图中所标。
则结点位移矩阵和结点荷载矩阵分别为：
l A I
E i j
xiF yiF iM
xjF yjF jM iu iv i? ju jv j?
有限单元法
Ti i i j j ju v u vδ?
Tx i y i i x j y j jF F M F F MF?
设单元上没有荷载作用,首先考虑轴向力的作用,由于杆端轴力,只引起杆端轴向位移,,根据拉压杆单元的单元刚度方程,有：
xiF xjF iu ju
x i i j
EA EAF u u
ll
x j i j
EA EAF u u
ll
有限单元法其次,杆端弯矩,和杆端剪力,只与杆端的转角位移,和杆端的横向位移,有关系,根据只计弯曲杆单元的单元刚度方程 （ 注意,由于不考虑单元上的荷载作用,
故方程式中的等效结点荷载 等于零 ） 可得：
iM jM yiF yjF
i? j? iv jv
EF?
3 2 3 2
12 6 12 6
y i i i j j
EI EI EI EIF v v
l l l l
3 2 3 2
12 6 12 6
y j i i j j
E I E I E I E IF v v
l l l l
22
6 4 6 2
i i i j j
EI EI EI EIM v v
l l l l
22
6 2 6 4
j i i j j
E I E I E I E IM v v
l l l l
有限单元法这样,上述表达式合并在一起,写成矩阵形式如下：
3 2 3 2
22
3 2 3 2
22
0 0 0 0
12 6 12 6
00
6 4 6 2
00
0 0 0 0
12 6 12 6
00
6 2 6 4
00
ixi
iyi
ii
jxj
jyj
jj
E A E A
ll
E I E I E I E I
uF
l l l l
vF
E I E I E I E I
M l l l l
uF E A E A
ll vF
E I E I E I E I
M
l l l l
E I E I E I E I
l l l l
Fkδ
可以将上式简写为：
有限单元法
2.2.5 空间杆单元
M
x i
M
y j
F
x i
F
y i
F
x j
F
y j
y
x
z
i
M
z i
F
z i
M
y i
M
x j
M
z j?F
z j
j
设局部坐标系的 轴为单元的形心主轴，横截面的两个主轴分别为 轴和 轴（如图所示）。设杆横截面面积为,杆单元长度为,在 平面内抗弯刚度为,在 平面内的抗弯刚度为,杆件的抗扭刚度为 。
x
y z A
l xz yEI xy
zEI EI
有限单元法空间刚架有 6个位移分量和 6个结点力分量,设局部坐标系下它们分别为：
Ti i i x i y i z i j j j x j y j z ju v w u v wδ?
Tx i y i z i x i y i z i x j y j z j x j y j z jF F F M M M F F F M M M?F?
有限单元法
3 2 3 2
3 2 3 2
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0
12 6 12 6
12 6 12 6
z z z z
xi
yi y y y y
zi
xi
yi
zi
xj
yj
zj
xj
yj
zj
EA EA
ll
EI EI EI EI
F
l l l l
F EI EI EI EI
l l l l
F
G I G I
ll
M
M
M
F
F
F
M
M
M
22
22
3 2 3 2
3 2 3 2
22
0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0
6 4 6 2
6 4 6 2
12 6 12 6
12 6 12 6
6 2 6 4
6
y y y y
z z z z
z z z z
y y y y
y y y y
EI EI EI EI
l l l l
EI EI EI EI
l l l l
EA EA
ll
EI EI EI EI
l l l l
EI EI EI EI
l l l l
G I G I
ll
EI EI EI EI
l l l l
22
0 0 0 0 0 0 0
2 6 4
i
i
i
xi
yi
zi
j
j
j
xj
yj
zj
z z z z
u
v
w
u
v
w
EI EI EI EI
l l l l
有限单元法
2.2.6 单元刚度矩阵的性质
① 单元刚度矩阵 为对称矩阵,其元素
② 单元刚度矩阵 中的每个元素代表单位杆端位移引起的杆端力 。 其中的任意元素 的物理意义是第 个杆端位移分量等于 1（ 其余位移分量等于 0） 时,所引起的第 个杆端力的分量值 。
k? ()i j j ik k i j
k?
ijk j
i
有限单元法
③ 一般单元的单元刚度矩阵 是奇异矩阵，它的元素组成的行列式等于零，即 。根据奇异矩阵的性质,没有逆矩阵。也就是说，如果给定杆端位移,根据（ 2-29）
或（ 2-31）式可以求出杆端力 的惟一解，但反过来，如果已知杆端力,则不能根据 来确定杆端位移的惟一解。因为即使在杆端力已知的情况下，由于单元两端无任何约束，因此除出杆端自身变形外，还可以发生任意的刚体位移。举例来说，如果物体处于静止状态，我们可以说其处于平衡状态，但反过来，如果物体处于平衡状态，则我们不能说其一定处于静止。
k?0?k?
k?
δ?
F?
1() dδ kF δ?F?
有限单元法
④ 单元刚度矩阵 具有分块的性质,即可以用子矩阵表示 。 用虚线把 分为四个子矩阵,把 和 各分为两个子矩阵,因此,又可以写为：
这里：
或或或或用子矩阵形式表示单元刚度矩阵和单元刚度方程，可以使其表达的物理意义更加明显。在单元刚度矩阵 中，其任意子矩阵 表示杆端力 和杆端位移 之间的关系。
k?
k?k? F? δ?
i ii ij i
j ji jj j
F k k δ
F k k δ
T
i x i y i iF F MF?
T
i x i y i z i x i y i z iF F F M M MF?
T
j x j y j jF F MF?
T
j x j y j z j x j y j z jF F F M M M?F?
T
i i i iuvδ?
T
i i i i x i y i z iu v wδ?
T
j j j juvδ?
T
j j j j x j y j z ju v wδ?
k?
rsk? rF? sδ?
有限单元法
2.3 杆系结构的整体分析
2.3.1 平面问题坐标转换矩阵图 2.9 平面问题杆端力转换示意图
i
j
o
y
x
M
i
M
j
F
x i
F
y i
F
x j
F
y j
y
x
F
x i
F
x j
F
y i
F
y j
M
j
M
i
i
j
o
y
x
y
x
有限单元法这里 表示由 轴到 轴的角,角度转动的正负由右手定则确定,本书中以顺时针方向转动为正 。 在两个坐标系中,力偶分量保持不变,即有：
同理,对于 端的杆端力,有：
c o s sin
sin c o s
x i x i y i
y i x i y i
F F F
F F F
x x
iiMM?
j
c os sin
sin c os
x j x j y j
y j x j y j
jj
F F F
F F F
MM
有限单元法将这些式子用矩阵形式可表示为：
c o s sin 0 0 0 0
sin c o s 0 0 0 0
0 0 1 0 0 0
0 0 0 c o s sin 0
0 0 0 sin c o s 0
0 0 0 0 0 1
xixi
yiyi
ii
xjxj
yjyj
jj
FF
FF
MM
FF
FF
MM
F T F c os si n 0 0 0 0
si n c os 0 0 0 0
0 0 1 0 0 0
0 0 0 c os si n 0
0 0 0 si n c os 0
0 0 0 0 0 1
T?
上式可以简写成：
这即为两种坐标系下单元杆端力的坐标变换式。其坐标转换矩阵为：
有限单元法从坐标转换矩阵 的表达式可以看出,为正交矩阵,其逆矩阵等于其转置矩阵,即有：
并且有：
式中 为单位矩阵 。
同样的推导,可以得到两种坐标系下的杆端位移之间的转换关系为：
这里 和 分别为局部坐标系和整体坐标系下的杆端位移矩阵,为前面介绍的转换矩阵。
T?
T?
1 TTT
1T T I
I
δ T δ
δ? δ?
T?
有限单元法因此可得：
上式两边同乘以,可以得到：
设可得：
上式即为整体坐标系下的单元刚度方程。
T F k T δ
1?T?
1 TF T k T δ T k T δ
T?k T k T
Fkδ
有限单元法
2.3.2 空间问题坐标转换矩阵
o
F
x i
F
y i
F
x j
F
y j
y
x
x
y
z
i
j?F
z j
F
z i
z
有限单元法设 轴与 x,y,z轴的方向余弦分别为：
，,
则将杆端力,,向 轴投影,可以求得杆端力,
即：
同理可以求得：
x
c o s (,)xxl x x? c o s (,)xyl x y? c o s (,)xzl x z?
xiF yiF ziF
x
xiF
x i x i x x y i x y z i x zF F l F l F l
y i x i y x y i y y z i y zF F l F l F l
z i x i z x y i z y z i z zF F l F l F l
有限单元法用矩阵形式可以表示为：
上式即为结点 i的杆端力在局部坐标系和整体坐标系下的转换关系,其中的矩阵称为关系矩阵 。
与上面的推导类似，同样可以推出以,,表示,,,以及对于结点 j的相对应的转换关系，其中转换关系矩阵都是 。
x i x x x y x z x i
y i y x y y y z y i
z i z x z y z z z i
F l l l F
F l l l F
F l l l F
x x x y x z
y x y y y z
z x z y z z
l l l
l l l
l l l
λ
xiM yiM ziM
xiM yiM ziM
λ
有限单元法综上所述,整体坐标系下单元杆端力矩阵与局部坐标系下单元杆端力矩阵具有如下的关系表达式：
其中的 为：
称为空间坐标系的单元转换矩阵,它是一个正交矩阵,即：
对于杆端位移,同样可推导出在两种坐标系中的转换关系：
这样,可得空间杆件单元在整体坐标系中单元刚度方程为：
其中表示空间单元在整体坐标系中的单元刚度矩阵 。
F T F
T
λ
λ
T
λ
λ
000
0 0 0
0 0 0
000
1 TTT
δ T δ
Fk δ
T?k T k T
有限单元法
2.3.3 整体分析对杆系结构进行单元分析，仅仅是有限元分析中的第一步。
我们的目的是要对整个结构进行分析，研究结构的整体性能。
因此，在对结构的各单元分析完成后，必须将单元分析的结果进行整合，对结构进行整体分析。整体分析的过程实际上是如何将单元分析的结果进行有效组合，建立整体刚度方程并求解结点位移的过程。根据对结点位移的编码方式，可以采用,先处理法,和,后处理法,来建立整体刚度方程。
有限单元法
2.3.3.1 后处理法所谓后处理法,就是由单元刚度矩阵形成整体刚度矩阵,建立刚度方程后再引入支承条件,进而求解结点位移的方法 。 运用这种方法时,假设所有结点位移均为未知量,按照顺序统一进行编码,如图所示的平面杆件单元 。
1 ( 1 2 3 )
2 ( 4 5 6 ) 3 ( 7 8 9 )
4 ( 1 0 1 1 1 2 )
2
1 3
x
y
O
有限单元法结点位移矩阵为：
结点荷载矩阵为：
求出各单元刚度方程后,根据平衡条件和位移连续条件,可以建立整个结构的位移法方程：
1 2 3 4
1 1 1 2 2 2 3 3 3 4 4 4
T
Tu v u v u v u v
δ δ δ δ δ
1 2 3 4
1 1 1 2 2 2 3 3 3 4 4 4
T
T
x y x y x y x yF F M F F M F F M F F M
F F F F F
11
22
33
44
ii ij
ji jj ii ij
ji jj ii ij
ji jj
kkF δ
k k k kF δ
k k k kF δ
kkF δ
00
0
0
00
①①
① ① ② ②
② ② ③ ③
③③
有限单元法简写成：
这里 为结构的整体刚度矩阵,有：
FKδ
K
1 1 1 2 1 3 1 4
2 1 2 2 2 3 2 4
3 1 3 2 3 3 3 4
4 1 4 2 4 3 4 4
K K K K
K K K K
K K K K
K K K K
K
应该注意到，在建立方程的过程中，我们假设所有结点都有位移。因此整个结构在外力作用下，除了发生弹性变形外，还可能发生刚体平动位移，这样各结点位移不能唯一确定。这说明整体刚度矩阵为一奇异矩阵，不能求逆矩阵，即根据整体刚度方程可得到无穷多个解。
有限单元法实际上,在图所示刚架中,结点 1和结点 4均为固定端,其三个位移分量均为 0,即有：
这样,将上述支承条件引入到方程中,对整体刚度方程进行修改,可得：
1 1 1 4 4 4 0u v u v
1 1 1 2 1 3 1 41
2 1 2 2 2 3 2 42 2
3 1 3 2 3 3 3 43 3
4 1 4 2 4 3 4 44
K K K K
K K K K
K K K K
K K K K
F
F δ
F δ
F
0
0
有限单元法对上述方程进行化简，可以得到两组方程：
和
2 2 2 2 3 2
3 3 2 3 3 3
KK
KK
F δ
F δ
1 2 21
4 3 34
K
K
δF
δF
0
0
这样，利用第 1式可以求得结点位移 和,再根据第 2式可以求得支座反力 和 。
2δ 3δ
1F 4F
有限单元法
2.3.3.2 先（前）处理法所谓先处理法,就是先引入支承条件,根据支承条件仅对未知的自由结点位移分量编号,得到的位移矩阵中不包含已知的约束位移分量,即可以直接得到方程求解自由结点位移分量 。
有限单元法如图所示的平面刚架结构 ABCD,由于在 A处和 D处均为固定端,其位移为 0,故位移编码均为 0,在 C处为铰接,故 BC杆在 C端的角位移与 DC杆在 C端的角位移不相同,因此在 C处编两个结点 3和 4,但结点 3和 4的横向位移和竖向位移相同,故采用相同的编号,各结点位移编码如图所示 。
图 2.12 先处理法位移编码示意图
1 ( 0 0 0 )
2 ( 1 2 3 ) 3 ( 4 5 6 )
4 ( 4 5 7 )
A
B
C
D 5 ( 0 0 0 )
2
1 3
x
y
O
有限单元法由单元刚度矩阵集成整体刚度矩阵通常采用刚度集成法。
其计算过程可以分为两步：首先求出各单元的贡献矩阵，然后将它们叠加起来形成整体刚度矩阵。但这样处理在实际中很少采用，因为在编程过程时需先将各单元的贡献矩阵储存起来，
而各单元贡献矩阵的阶数与整体刚度矩阵的阶数相同，因此占用的空间非常巨大，不利于节约资源，并且在实际中有可能耗尽所有资源。
故在实际中并不是采用贡献矩阵法，而是利用各单元的定位数组，采用,边定位，边累加,的方法。
有限单元法所谓单元的定位数组,就是将单元 的始端及末端的位移编码排成一行 （ 始端在前,末端在后 ）,写成如下的形式：
如图中,各单元的定位数组分别为：
这样处理，得到的整体刚度矩阵与叠加所有贡献矩阵得到的结果完全一致，但可以节约大量的存储空间。
1 2 3() nm m m m m
( 0 0 0 1 2 3 )m?①
(1 2 3 4 5 6 )m?②
( 0 0 0 4 5 7 )m?③
有限单元法如图所示的结构有 3个单元,5个结点,7个独立的位移分量,
这样其整体刚度矩阵应为 阶矩阵,即：77?
11 12 13 14 15 16 17
21 22 23 24 25 26 27
31 32 33 34 35 36 37
41 42 43 44 45 46 47
51 52 53 54 55 56 57
61 62 63 64 65 66 67
71 72 73 74 75 76 77
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
有限单元法对于单元,其单元定位数组为,我们将定位数组标记在单元刚度矩阵的上面和右侧,其与单元刚度矩阵一起写成如下形式：
① ( 0 0 0 1 2 3 )m?①
1 1 1 2 1 3 1 4 1 5 1 6
2 1 2 2 2 3 2 4 2 5 2 6
3 1 3 2 3 3 3 4 3 5 3 6
4 1 4 2 4 3 4 4 4 5 4 6
5 1 5 2 5 3 5 4 5 5 5 6
6 1 6 2 6 3 6 4 6 5 6 6
0 0 0 1 2 3
0
0
0
1
2
3
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
4 4 1 1 4 5 1 2 4 6 1 3
5 4 2 1 5 5 2 2 5 6 2 3
6 4 3 1 6 5 3 2 6 6 3 3
k K k K k K
k K k K k K
k K k K k K
① ① ①
① ① ①
① ① ①
有限单元法单元 ② 的刚度矩阵和它的定位数组为：
中各元素在整体刚度矩阵 中的位置为：
1 1 1 2 1 3 1 4 1 5 1 6
2 1 2 2 2 3 2 4 2 5 2 6
3 1 3 2 3 3 3 4 3 5 3 6
4 1 4 2 4 3 4 4 4 5 4 6
5 1 5 2 5 3 5 4 5 5 5 6
6 1 6 2 6 3 6 4 6 5 6 6
1 2 3 4 5 6
1
2
3
4
5
6
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
( 1,2,,6 1,2,,6 )i j i jk K i j②
有限单元法单元 ③ 的刚度矩阵和它的定位数组为：
中各元素在整体刚度矩阵 中的位置为：k
③ K
1 1 1 2 1 3 1 4 1 5 1 6
2 1 2 2 2 3 2 4 2 5 2 6
3 1 3 2 3 3 3 4 3 5 3 6
4 1 4 2 4 3 4 4 4 5 4 6
5 1 5 2 5 3 5 4 5 5 5 6
6 1 6 2 6 3 6 4 6 5 6 6
0 0 0 4 5 7
0
0
0
4
5
7
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
③
③
4 4 4 4 4 5 4 5 4 6 4 7
5 4 5 4 5 5 5 5 5 6 5 7
6 4 7 4 6 5 7 5 6 6 7 7
k K k K k K
k K k K k K
k K k K k K
③③③
③③③
③③③
有限单元法这样,按照以上所讲的定位方法,将,,中的相关元素累加到整体刚度矩阵对应的元素上,可以得到整体刚度矩阵为：
k① k② k③
4 4 1 1 4 5 1 2 4 6 1 3 1 4 1 5 1 6
5 4 2 1 5 5 2 2 5 6 2 3 2 4 2 5 2 6
6 4 3 1 6 5 3 2 6 6 3 3 3 4 3 5 3 6
4 1 4 2 4 3 4 4 4 4 4 5 4 5 4 6 4 6
5 1 5 2 5 3 5 4 5 4 5 5 5 5 5 6 5 6
6
0
0
0
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
① ② ① ② ① ② ② ② ②
① ② ① ② ① ② ② ② ②
① ② ① ② ① ② ② ② ②
② ② ② ② ③ ② ③ ② ③
② ② ② ② ③ ② ③ ② ③
1 6 2 6 3 6 4 6 5 6 6
6 4 6 5 6 6
0
0 0 0 0
k k k k k
k k k
② ② ② ② ② ②
③ ③ ③
有限单元法在实际编程过程中,集成整体刚度矩阵的过程可以概括为：
（ 1） 根据位移编号,建立整体刚度矩阵 Ｋ 的存储空间,
并置 0,若位移个数为 n,则整体刚度矩阵为 阶矩阵 。
（ 2）对每个单元，根据其定位数组，将其元素累加到整体刚度矩阵中对应的元素上，直到处理完所有的单元。
nn?
有限单元法
2.4 等效结点荷载与边界条件的处理
2.4.1 非结点荷载的处理有限单元法实际上,在前面的分析中,我们已经介绍了求等效结点荷载的方法,如 2-7式,2-15式,2-22式分别可用来求不同情况下的等效结点荷载 。
此外,可以这样来考虑：
第一步,在局部坐标系下求单元的固端力 。 对于某个单元,我们假定单元的两端均固定,然后根据静力平衡求得固定端的反力 。
第二步,根据单元固端力求单元的等效结点荷载 。 根据局部坐标系与整体坐标系单元杆端力的变换式,固端内力在两种坐标系下的变换形式可以写成：
因此,整体坐标系下的等效结点荷载矩阵可以由下式计算：
fF?
Tff?F T F
EF?
TEfF T F
有限单元法
2.4.2 边界条件的处理
2.4.2.1 铰结点在杆系结构中,除了刚性结点外,通常会遇到一些杆件通过铰结点与其它杆件联结,如图 2.13所示杆件系统,有 4根杆件汇交于 D点,其中 BD杆在 D端通过铰支座与其它杆件铰接,其余 3根杆为刚性接触 。 对于这样的铰结点,具有如下的性质：
① 铰结点上各杆具有相同的线位移,但截面的转角位移不相同；
② 结点上具有铰接杆端不承受弯矩作用 。 如图 2.13所示结构中,BD杆在 D端的杆端弯矩为 0,只有 CD,ED,GD杆在结点
D上与外弯矩保持平衡 。
有限单元法对于这样的结点,我们在对其进行单元划分时,通常考虑在 D处设置 2个结点 。 按照先处理法,对图示结构进行位移编码,
如图 2.14b所示 。
A
B
C
D
E H
G
F
1 ( 0 0 1 )
2 ( 2 3 4 )
3 ( 5 6 7 )
4 ( 5 6 8 )
5 ( 0 0 0 )
6 ( 9 1 0 1 1 )
7 ( 0 0 0 )
8 ( 1 2 1 3 1 4 )
9 ( 1 5 1 6 1 7 )
图 2.14 铰结点的处理示意图有限单元法
2.4.2.2 弹性支承点在实际工程中，有时会遇到弹性支承的情况（如图 2.15
），这时一般将弹性支座看作是在结构约束点沿约束方向的一个弹簧，弹簧的刚度系数为,在数值上等于使弹簧支座沿约束方向产生单位位移时所需施加的力。
kk
1 2 i - 1 i i + 1 n
k
有限单元法设结构的第 i个结点位移分量 为弹性支座约束,弹簧的刚度系数为,则结构产生位移时所引起的支座反力为：
这里负号表示支座反力方向与约束点位移方向始终相反 。 作用在受约束的结点上,它是结点外力的一部分,由整体刚度方程可知,第个平衡方程应为：
式中 为原有的第 i个结点上的荷载分量 。 进行整理可得：
这样就引入了弹性支承的约束条件 。 根据以上的分析,引入弹性支承的具体做法可以归结为：先解除弹性支承点约束,在 i处给一个结点号,形成总刚度矩阵,然后在总刚度矩阵中将第 i行的主元素加上弹性支承的刚度系数,此时第行变为：
i?
k
RiFk
1 1 2 2i i i i i i n n i RK K K K F F
iF
1 1 2 2 ()i i i i i i n n iK K K k K F
1 2 3i i i i i i nK K K K k K?
有限单元法
2.5 杆系结构分析算例应用有限元位移法对杆系结构进行分析,如果按照,先处理法,的思想,其步骤可以归结为如下：
（ 1） 对单元进行划分,整理原始数据,对单元和结点进行编号,确定每个单元的局部坐标系整个结构的整体坐标系 。
（ 2） 计算局部坐标系中的单元刚度矩阵 。
（ 3） 确定每个单元的坐标转换矩阵,计算整体坐标系下的单元刚度矩阵 。
有限单元法
（ 4） 根据各单元的位移分量编号,形成单元的定位数组
,按照,对号入座,同号叠加,的方法,形成结构的整体刚度矩阵 。
（ 5） 计算总的结点荷载矩阵 。 在这一步里,首先必须将非结点荷载转换成等效结点荷载,再与对应的结点荷载叠加,
形成总的结点荷载矩阵 。
（ 6） 求解结构的整体刚度方程,计算未知的结点位移矩阵 。
（ 7） 计算各单元的杆端力 。
（ 8）对计算结果进行整理，如果需要，计算内力图。
有限单元法例 1
1 ( 0 0 0 )
2 ( 1 2 3 )
3 ( 0 0 4 )
1 0 k N
8 k N6 k N
1 0 k N? m
1
2
k
N
/
m
A
B
C
2,5 m
2
1
x
y
O
2,5 m
如图所示的刚架结构 ABC，
各杆截面尺寸相同，材料性质一样，求刚架的整体刚度矩阵 。
5A B B C l m
20.5Am?
43 1 0E M P a
41
24Im?
K
杆长：
截面面积：
弹性模量：
截面惯性矩：
有限单元法
（ 1） 单元划分,建立局部坐标系和整体坐标系 （ 如图
2.13所示 ）,并对数据进行整理,对单元和结点编号 。
44
44
23
4
3 0 0 1 0 1 0 0 1 0
6 1 2
3 0 1 0 1 2 1 0
E A E I
ll
E I E I
ll
有限单元法
（ 2） 求局部坐标系中的单元刚度矩阵 。 由于单元 ①,
② 的尺寸完全一样,因此其单元刚度矩阵一样,根据 （ 2-30）
式可以得到：
k?
4
300 0 0 300 0 0
0 12 30 0 12 30
0 30 100 0 30 50
10
300 0 0 300 0 0
0 12 30 0 12 30
0 30 50 0 30 100
kk
①②
有限单元法
（ 3） 求整体坐标系中的单元刚度矩阵 。 对于单元 ①,其局部坐标系与整体坐标系的夹角为,故根据 （ 2-40） 式,
计算其坐标转换矩阵为：
单元的定位数组为：
根据 （ 2-46） 式计算单元 ① 在整体坐标系下的单元刚度矩阵为：
90
0 1 0 0 0 0
1 0 0 0 0 0
0 0 1 0 0 0
0 0 0 0 1 0
0 0 0 1 0 0
0 0 0 0 0 1
T
①
(1 2 3 0 0 0 )m?①
4
1 2 3 0 0 0
1 12 0 30 12 0 30
2 0 30 0 0 0 30 0 0
3 30 0 10 0 30 0 50
10
0 12 0 30 12 0 30
0 0 30 0 0 0 30 0 0
0 30 0 50 30 0 10 0
k
①
有限单元法对于单元 ②,由于其局部坐标系与整体坐标系一致,因此两种坐标系下的单元刚度矩阵相同,即有：
（ 4） 形成整体刚度矩阵 。 根据前面介绍的方法建立整体刚度矩阵如下：
4
1 2 3 0 0 4
1 30 0 0 0 30 0 0 0
2 0 12 30 0 12 30
3 0 30 10 0 0 30 50
10
0 30 0 0 0 30 0 0 0
0 0 12 30 0 12 30
4 0 30 50 0 30 10 0
kk
②②
4
3 1 2 0 3 0 0
0 3 1 2 3 0 3 0
10
3 0 3 0 2 0 0 5 0
0 3 0 5 0 1 0 0
K
有限单元法
（ 5） 求总结点荷载 。
首先,求各单元在局部坐标系中的固端内力 。 对于单元 ①
（ 这里应该注意,其 局部坐标系的 轴方向向左 ）,
，,根据表 2.1可得：
对于单元 ②,,,,根据表 2.1可得：
，，
因此有：
fF?
y 1 2 k N /mq
5mal 0b?
0fx i fx jFF①① 3 0 k Nf y i f y jFF①① 2 5 k N mf i f jMM①①
8 kNF? 2,5 mab 5ml?
0fx i fx jFF②② 4 k Nf y i f y jFF②② 5 k N mf i f jMM②②
0 30 25 0 30 25 TfF ①
0 4 5 0 4 5 TfF ②
有限单元法其次,求整体坐标系下各单元的等效结点荷载 。 对于单元 ①,,单元 ②,,因此可以求得各单元的等效结点荷载为：
90 0
3 0 0 0 1 2
3 0 0 2 5 3 0 0 2 5 TEF ①
3 0 01 2 4
0 4 5 0 4 5 TEF ②
第三,求刚架等效结点荷载矩阵 。 按照,同号叠加,的方法有：
第四,求直接作用在结点上的荷载,
最后,求总的结点荷载矩阵
3 0 4 2 0 5 TEF
1 0 6 1 0 0 TdF
2 0 1 0 1 0 5 TdEF F F
有限单元法
（ 6） 求解结构的整体刚度方程,计算未知的结点位移矩阵 。 这里的整体刚度方程为：δ
1
14
1
3
31 2 0 30 0 20
0 31 2 30 30 10
10
30 30 20 0 50 10
0 30 50 10 0 5
u
v
1
1 6
1
3
6,06 54
3,97 29
10
3,58 65
4,39 86
u
v
δ
解得：
有限单元法
（ 7） 计算各单元的杆端内力 。 首先从求出的结点位移中取出各单元在整体坐标系下的杆端位移,有：
δ
δ?1
1
1 6
2
2
2
6,0 6 5 4
3,9 7 2 9
3,5 8 6 5
10
0
0
0
u
v
u
v
δ ①
1
1
1 6
3
3
3
6,0 6 5 4
3,9 7 2 9
3,5 8 6 5
10
0
0
4,3 9 8 6
u
v
u
v
δ ②
然后计算杆端内力。
有限单元法
4
300 0 0 300 0 0
0 12 30 0 12 30
0 30 100 0 30 50
10
300 0 0 300 0 0
0 12 30 0 12 30
0 30 50 0 30 100
6.0654
3.9729
3.5865
0
0 1 0 0 0 0
1 0 0 0 0 0
0 0 1 0 0 0
0 0 0 0 1 0
0
0
0 0 1 0 0
0 0 0 0
0
01
ff
F k F k Tδ δ F
① ① ① ① ① ① ① ①
6
0 11,919 kN
30 28,196 kN
25 19,594 kN m
10
0 11,919 kN
30 31,804 kN
25 28,613 kN m
对于单元①：
有限单元法对于单元②,
4
6
3 0 0 0 0 3 0 0 0 0
0 1 2 3 0 0 1 2 3 0
0 3 0 1 0 0 0 3 0 5 0
10
3 0 0 0 0 3 0 0 0 0
0 1 2 3 0 0 1 2 3 0
0 3 0 5 0 0 3 0 1 0 0
6,0 6 5 4 0 1 8,1 9 6 k
3,9 7 2 9 4
3,5 8 6 5 5
10
00
04
4,3 9 8 6 5
ff
F k F k Fδ δ
② ② ② ② ② ② ②
N
5,9 1 9 k N
9,5 9 4 k N m
1 8,1 9 6 k N
2,0 8 1 k N
0
有限单元法
（ c）剪力图
（ b）弯矩图
（ d）轴力图
（ a）位移曲线图 2.17
1
2 3
2
1
1
2 3
2
8
.6
1
1
9
.5
9
-1
3
.4
0
- 9,5 9
5,2 0
1
2
3
1
.8
0
-2
8
.2
0
5,9 2
- 2,0 8
1
2 3
2
1
1
2
3
-1
1
.9
2
- 1 8,2 0
2
1
有限单元法第 3章弹性力学平面问题的有限元分析有限单元法
3.1 概述在有限元法中,把单元与单元之间设置的相互连接点,称为结点 。 一般用号码 1,2,… 进行结点编号 。 结点可为铰接,
固接或其它形式的连接 。 结点的设置,性质及数目等均视所研究问题的性质,描绘变形状态的需要和计算精度的要求等而定 。
在有限元法中引进结点概念是至关重要的。有了结点，才可将实际连续体看成是仅在结点处相互连接的单元集合组成的离散型结构，从而可使研究的对象转化成可以使用电子计算机计算的教学模型。由单元、结点、结点连线构成的集合称为有限元模型。它是有限元分析与计算的对象。
有限单元法
3.1.1 单元划分类型单元类型：三角形,四边形单元数目：根据计算精度要求来确定结点设置：使单元的的结点编号尽量靠近有限元模型：由单元,结点,结点连线构成的集合有限单元法
3.1.2 位移函数在选择多项式时,为了使有限单元法的计算精度和收敛性得到保障,还需要满足完备性和连续性的要求 。 为了使位移模式尽可能地反映物体中的真实位移形态,它应满足下列条件：
（ 1） 位移模式必须能反映单元的刚体位移；
（ 2） 位移模式必须能反映单元的常量应变；
（ 3） 位移模式应尽可能地反映位移的连续性 。
弹性力学平面问题一般选择多项式函数作为位移函数。
...26524321 yaxyaxayaxaau
.,,26524321 ybxybxbybxbbv
有限单元法
3.2 平面三角形单元三角形单元是一种简单方便、对边界适应性强的单元，
由于以三角形的三个顶点作为结点，因此又成为三结点三角形单元。这种单元的计算精度较低，使用的时候必须进行精细的网格划分，但他仍然是一种常用的单元有限单元法
3.2.1 位移函数的选取 ( x m,y m )
u
i
v
i
( x
i
,y
i
)
( x
j
,y
j
)
i
j
m
x
y
O
u
j
u
m
v
j
v
m
如图 3.4所示的 3结点平面三角形单元,结点 的坐标分别为,,,结点位移分别为,,,,,。 记单元的结点位移向量 和结点力向量 为：
i j m、,
),( ii yx ),( jj yx ),( mm yx iu iv ju
jv mu mv
δ? F?
T
i i j j m mu v u v u vδ?
T
x i y i x j y j x m y mF F F F F F?F?
有限单元法根据完备性和连续性的要求,选取 3结点三角形单元的位移场函数 如下：
（ 3-4）
将 3个结点 上的坐标和位移分别代入式 （ 3-4） 就可以将六个待定系数用结点坐标和结点位移分量表示出来 。
将水平位移分量和结点坐标分别代入 （ 3-4） 中的第一式,得到：
d?
yaxaa
yaxaau
654
321
v
i j m、,
mmm
jjj
iii
yaxaau
yaxaau
yaxaau
321
321
321
有限单元法写成矩阵形式,有：
（ 3-6）
其中,
为三角形单元的面积 。 为了避免出现 的情况,三个结点按逆时针顺序排列 。 由 （ 3-6） 得到：
1
2
3
1
1
1
i i i
j j j
m m m
u x y a
u x y a
u x y a
1
2 de t 1
1
ii
jj
mm
xy
A x y
xy
A 0A?
i j m、,
1
1
2
3
1
1
1
i i i
j j j
m m m
a x y u
a x y u
a x y u
有限单元法将竖向位移和结点坐标代入 （ 3-4） 中的第二式,可以得到：
（ 3-8）
将（ 3-7）、（ 3-8）代回（ 3-4）整理后可得：
1
4
5
6
1
1
1
i i i
j j j
m m m
a x y v
a x y v
a x y v
（ 3-9）
（ 3-10）
其中系数,
（ 下 标轮换 ）
])()()[(2 1 mmmmjjjjiiii uycxbauycxbauycxbaAu
])()()[(2 1 mmmmjjjjiiii vycxbavycxbavycxbaAv
i j m m ja x y x y
i j mb y y
i j mc x x
i j m、,
有限单元法令 （ 3-11）
可得：
（ 3-12）
单元内的位移场函数可以简写成：
（ 3-13）
把 称为形函数矩阵,称为形函数 。 根据形函数的定义,
具有以下性质：
（ 1） 在单元结点上形函数的值为 1或为 0。
（ 2）在单元中的任意一点上，三个形函数之和等于 1。
)(21 ycxbaAN iiii
0 0 0
0 0 0
i
i
i j m j
i j m j
m
m
u
v
N N N uu
N N N vv
u
v
dNδ
N
iN
iN
有限单元法
3.2.2 单元的应变场根据单元的位移场函数式 （ 3-12）,由几何方程可以得到单元的应变场表达式：
（ 3-14）
记为,（ 3-15）
其中,矩阵 称为几何矩阵 。 矩阵 可以表示为分块矩阵的形式
（ 3-16）
这里,
0 0 0
1
0 0 0
2
i
i
i j m
j
i j m
j
i i j j m m
m
m
u
u
vx
b b b
uv
c c c
vyA
c b c b c b
uuv
yx v
ε
ε Bδ?
B B
[]i j m?B B B B
0
1 0 (,,)
2
r
rr
rr
b
c r i j mA
cb
B
有限单元法
3.2.3 单元的应力场由物理方程及式 （ 3-15）,可以得到单元的应力场表达式：
（ 3-18）
其中 为应力矩阵,称为弹性矩阵,对于平面应力问题,
(3-19)
ζ D ε DB δ S δ
S DB D
2
10
10
( 1 )
1
00
2
E
D
有限单元法将应力矩阵表示为分块矩阵的形式,有：
（ 3-20）
其中：
[]i j m?S S S S
2
(,,)
2 ( 1 )
11
22
rr
r r r r
rr
bc
E
b c r i j m
A
cb
S D B
1
1
(,,)
2 ( 1 2 ) 1 1
1 2 1 2
2 1 2 1
rr
r r i r
rr
bc
E
b c r i j m
A
cb
S D B
对于平面应变问题，只需将 换为,换为,
则（ 3-21）变为：
E
21
E
1
有限单元法
3.2.4 单元刚度矩阵
（ 1） 单元的应变能
（ 3-23）
（ 2） 单元上外力的势能
1
2 TUd ε ζ
T T Tvs
AV d d A cf p f p δ F
（ 3-24）
式中,,分别表示单位体积的体积力、单元上的表面力、单元结点上的结点荷载。
vp sp cF?
有限单元法将式（ 3-13）、（ 3-15）、（ 3-18）代入（ 3-23）、（ 3-24）
后相加，得到单元的总势能为：
1
2 T T T T Tv s cAd p d p d Aδ B D B δ δ N N F
δ?
利用最小势能原理，取结点位移 的变分，得到：
0 δδ
由 的任意性，有：?δ
0
i i j j m mu v u v u v
δ
有限单元法考虑到 的对称性，对式（ 3-25）求偏导得到：
TB DB
0T T Tv s cAd d d AB D B δ N p N p F
记：
T d
B D B k
TTv s c
Ap d p d AN N F F
则式（ 3-28）可写为：
k δ F
有限单元法这就是描述单元结点力和结点位移向量之间关系的平衡方程。
其中 称为单元刚度矩阵。在 3结点等厚三角形单元中 和的分量均为常量，则单元刚度矩阵可以表示为：
k? B D
T tA?k B D B?
其中 t,Ａ 分别为单元的厚度和面积。单元刚度矩阵 可以表示为分块矩阵的形式：
k?
ii ij im
ji jj jm
m i m j m m
k k k
k k k k
k k k
有限单元法其中：
,,
,,
,,,r x s x r x s yTr s r s
r y s x r y s y
kk r s i j m
kk
k B D B
对于平面应力问题，其刚度矩阵的显式：
2
11
22,,,
114 ( 1 )
22
r s r s r s r s
rs
r s r s r s r s
b b c c b c c b
Et
r s i j m
A
c b b c c c b b
k
有限单元法
3.2.5 等效结点荷载在式（ 3-28）中括号中的前两个量分别为体积力、
表面力移置到结点上的等效结点力，依次定义为,,即：
、
F?
sF?
T v d?
F N p
Tss
A dAF N p
那么（ 3-30）中的荷载列阵就等于体积力、表面力和单元结点荷载叠加。因此作用在弹性体上的外力，需要移置到相应的结点上成为结点荷载。荷载移置要满足静力等效原则。
有限单元法
3.2.6 整体刚度矩阵得到了单元刚度矩阵后，需要将一系列的将单元组成一个整体结构，然后根据结点载荷平衡的原则进行分析，
得到整体刚度矩阵。整体分析包括以下 4个步骤：
（ 1）建立整体刚度矩阵，
（ 2）根据支承条件修改整体刚度矩阵，
（ 3）解方程组，求出结点的位移，
（ 4）根据结点位移，求出单元的应变和应力。
由单元刚度矩阵得到整体刚度矩阵的基本方法是刚度集成法，即整体刚度矩阵是单元刚度矩阵的集成。
有限单元法
x
y
1
2
3
4 5 6
i
j
m
i
j
m
i
j
m
i
j
m
单元编号 单元结点局部编号 单元结点整体编号
（ 1）
3
1
2
（ 2）
5
2
4
（ 3）
5
3
2
（ 4）
3
5
6
i
j
m
有限单元法
2j（ ）
(2)m
(2)i
2j（ ）
(2)jjk (2)jmk (2)jik
(2)m
(2)mjk (2)
mmk
(2)mik
(2)i
(2)ijk (2)imk (2)iik
整体编号 1 2 3 4 5 6
局部编号 0 0 0
1 0 0 0 0 0 0
2 0 0 0
3 0 0 0 0 0 0
4 0 0 0
5 0 0 0
6 0 0 0 0 0 0
有限单元法有限单元法由上表可以看出，整体刚度矩阵也具有对称性。同时它是一个稀疏矩阵，即其中有大量的零元素，并且非零元素都集中于主对角线附近呈带状。和单元刚度矩阵一样，由于位移函数中包含刚体位移，所以整体刚度矩阵也是一个奇异矩阵。必须要排除刚体位移后，才能变为正定矩阵。
有限单元法
3.2.7 约束条件的处理
（ 1）将方程分组
（ 2）对角线元素改 1，同行同列其他元素改 0
（ 3）对角线元素乘以一个大数有限单元法
（ 2）对角线元素改 1，同行同列其他元素改 0
若第 r个位移分量 为已知值 则将整体刚度矩阵中主对角线元素 改为 1,第 r行、第 r列的其他元素改为 0,
荷载矩阵中第 r行的元素改为,其他元素都减去结点位移的已知值和原 Ｋ 中这行相应列元素的乘积。
ru *ru
rrK
*ru
*
1 1 1 2 1 3 1 1 11
*
2 1 2 2 2 3 2 2 22
*
3 1 3 2 3 3 3 3 33
*
*
1 2 3
0
0
0
0 0 0 1 0
0
n rr
n rr
n rr
r r
n n n n n n n n r r
K K K K u F K u
K K K K u F K u
K K K K u F K u
u u
K K K K u u K u
有限单元法
（ 3）对角线元素乘以一个大数若第 r个位移分量 为已知值 则将整体刚度矩阵中主对角线元素 乘以一个大数（如 ）,第 r行、第 r列的其他元素不变,荷载矩阵中第 r行的元素改为,
其他元素不变。
ru *ru
rrK
* 2010rr rKu?
1 1 1 2 1 3 1 1 1 1
2 1 2 2 2 3 2 2 2 2
3 1 3 2 3 3 3 3 3 3
2 0 * 2 0
1 2 3
1 2 3
1 0 1 0
rn
rn
rn
r r r r r r n r r r r
n n n n r n n n n
K K K K K u F
K K K K K u F
K K K K K u F
K K K K K u K u
K K K K K u u
2010
有限单元法
3.3 平面矩形单元 ( x m,y m )
u
i
v
i ( x
i
,y
i
)
( x
j
,y
j
)
i j
m u
m
v
j
v
m
u
j
( x
k
,y
k
)
k
u
k
v
k
a a
b
b
x,?
y,?
矩形单元也是常用的单元之一，由于采用了比常应变三角形单元更高次数的位移模式，故可以更好地反映弹性体的位移状态和应力状态。
有限单元法如图 3.10所示四结点矩形单元，记单元的结点位移向量和结点力向量 为：
δ?
F?
T
i i j j m m k ku v u v u v u vδ?
T
x i y i x j y j x m y m x k y kF F F F F F F F?F?
为了能推导出简洁的结果，在这里引入无量纲坐标：
a
x
b
y
( x
m
,y
m
)
u
i
v
i ( x
i
,y
i
)
( x
j
,y
j
)
i j
m u m
v
j
v
m
u
j
( x
k
,y
k
)
k
u
k
v
k
a a
b
b
x,?
y,?
有限单元法
3.3.1 单元位移场从图 3.10可以看出，结点条件共有 8个，因此,x和 y方向的位移场可以各有 4个待定系数，可以取以下多项式作为单元的位移场模式：
1 2 3 4
5 6 7 8v
u a a x a y a x y
a a x a y a x y
它们是具有完全一次项的非完全二次项，其中以上两式中右端的第四项 是考虑到 x和 y方向的对称性而取的。
48,a xy a xy
(3-43)
有限单元法由结点条件，在 处，有,(,,,)rrx x y y r i j m k
,(,,,)
,
r r r
r r r
u x y u r i j k m
v x y v
将 (3-44)代入 (3-43)中，可以求解出待定系数,然后再代回 (3-43)中经整理后，有
18aa
(3-44)
0 0 0 0
0 0 0 0
i j m k
i j m k
N N N Nu
N N N Nv
δ N δ
(3-45)
其中,N为单元的形函数矩阵，
有限单元法
1
11
4
1
11
4
1
11
4
1
11
4
i
j
m
k
xy
N
ab
xy
N
ab
xy
N
ab
xy
N
ab
如以无量纲坐标系来表达，则 (3-46)式可以写成
1 1 1 (,,,)4r r rN r i j m k
其中：
,(,,,)rrrrxy r i j m kaa
(3-46)
(3-47)
有限单元法
3.3.2 单元应变场根据单元的位移场函数式（ 3-45）,(3-47)，由几何方程可以得到单元的应变场表达式，
1
uu
b
x
vv
a
y a b
uvuv
ab
yx
ε
记为：
ε Bδ?
(3-48)
(3-49)
有限单元法这里,B矩阵称为几何矩阵。 B矩阵可以表示为分块矩阵的形式
[]i j m kB B B B?B
其中
0
0
00
0
( 1 ) 0
11
0 0 ( 1 ) (,,,)
4
( 1 ) ( 1 )
r
r
r
rr
ri
rr
N
b
b
N
a a r i j k m
ab ab
ab
NN
ab
B
(3-51)
(3-50)
有限单元法
3.3.3 单元应力场由物理方程及式（ 3-49），可以得到单元的应力场表达式，
ζ D ε DB δ S δ
其中 为应力矩阵,D称为弹性矩阵，对于平面应力问题，
S DB
2
10
10
( 1 )
1
00
2
E
D
(3-52)
有限单元法将应力矩阵表示为分块矩阵的形式
[]i j m k?S S S S S
其中：
00
002
00
( 1 ) ( 1 )
( 1 ) ( 1 ) (,,,)
4 ( 1 )
11
( 1 ) ( 1 )
22
rr
r r r r
rr
ba
E
b a r i j k m
ab
ab
S D B
对于平面应变问题，只需将 E换为,换为 。
21
E
1
（ 3-54）
（ 3-53）
有限单元法
3.3.4 单元刚度矩阵和三角形单元一样，可以根据最小势能原理导出结点位移向量和结点力向量之间关系，即单元的刚度矩阵,
可以将其写成分块的形式。
K?
i i i j i m i k
j i j j j m j k
mi mj mm mk
k i k j k m k k
k k k k
k k k k
k
k k k k
k k k k
其中
(,,,,)Tr s r s d r s i j m kk B D B （ 3-56）
（ 3-55）
有限单元法对于平面应力问题，如果单元厚度 t为常数，则得到
（ 3-56）的显式形式：
2
1
1
3 1
211
1
23
,,,,
4( 1 ) 1
1
31
2 11
1
23
r s r s
r s r s
r s r s
rs
r s r s
r s r s
r s r s
b
a
a
bEt
r s i j m k
ba
ab
k
（ 3-57）
有限单元法
3.4 平面问题程序设计有限元方法在分析工程问题时的优越性,在于与计算机的结合 。 通过编写的计算程序,确定所研究问题的信息和数据,问题就可以得到解决 。
前面两节讲述了平面问题的有限元计算原理。本节主要介绍采用三角形单元，计算在给定荷载作用下的弹性力学平面静力问题时，有限元法的程序编制与框图设计。通常有限元法的程序设计和编写分为以下几步：
有限单元法
（ 1） 问题分析使用计算机解决具体问题时，首先要根据所要求解的问题，找出已知的数据和条件，建立数学力学模型，选择计算方案，准备所需要的原始数据。
有限单元法
（ 2） 框图设计计算方案选定以后，
需先画出总框图。平面问题的有限元法的总框图如图 3.11所示。
开 始输 入 原 始 数 据组 成 整 体 刚 度 矩 阵组 成 整 体 荷 载 向 量引 入 边 界 条 件解 刚 度 方 程,输 出 位 移求 应 力,输 出 应 力求 应 变,输 出 应 变求 支 反 力,输 出 支 反 力结 束有限单元法
（ 3） 设计程序根据总框图的规划和有限元法的基本原理,用自己所使用的程序语言把这个解决方案严格地描述出来,也就是编写出程序代码 。 它就是用于求解问题的源程序 。
（ 4） 程序的编译和调试将编写好的源程序输入计算机，进行编译，改正其中的语法等错误。并用实际的输入数据对编好的程序进行调试，
分析所得到的运行结果，进行程序的测试和调整，直至获得预期的结果。
有限单元法第 4章空间轴对称问题的有限元分析有限单元法教学提示,空间轴对称结构是工程中一类特殊的结构 。
对于这类结构,由于其几何形状,约束情况,所受荷载都对称于某一根轴 （ 如图 4.1所示 ）,因此,通常是在柱坐标系中进行分析 。 本章以三角形单元为例介绍了这一类结构的有限元分析过程 。
教学要求,本章要求学生了解 柱坐标系中三角形单元的建立方法，熟悉轴对称结构有限元分析的基本方法。
有限单元法
4.1 概述在工程中有许多结构，如活塞、厚壁容器等，他们的几何形状、约束情况及所受的荷载都对称于空间的某一根轴，
因此在物体中通过该轴的任何平面都是对称面，所有应力、
应变和位移也对称于该轴，这类问题称为空间轴对称问题。
研究轴对称问题时通常采用圆柱坐标系,以 轴为对称轴。
(,,)rz? z
有限单元法
x
ya
b
z
p
如图 4.1所示的受均布内压作用的长圆筒，通过 z轴的一个纵截面就是对称面。由于对称性，轴对问题共有 4个应力分量：
Tr z zrζ
有限单元法其中 表示沿半径方向的正应力，称为径向应力； 表示沿方向的正应力，称为环向应力或切向应力； 表示沿 ｚ 方向的正应力，称为轴向应力； 表示在圆柱面上沿ｚ方向作用的剪应力。
r
z?
zr?
同样，轴对称问题共有 4个应变分量：
Tr z zrε
其中 表示沿半径方向的正应变，称为径向正应变； 表示沿 方向的正应变，称为环向正应变或切向正应变； 表示沿ｚ方向的正应变，称为轴向正应变； 表示沿ｒ和ｚ方向的剪应变。
r
z?
zr?
有限单元法在轴对称问题中，弹性体内任意一点上，不存在切向位移，只存在径向位移 u 和轴向位移 w,两个位移分量表示为，
Tuw?d
有限单元法
4.2 三角形单元轴对称问题分析中所使用的三结点单元，在对称面上是三角形，在整个弹性体中是三棱圆环，各单元中圆环形铰相联接。
有限单元法
4.2.1 位移函数的选取
( r
m
,z
m
)
u
i
w
i
i ( r
i
,z
i
)
j ( r
j
,z
j
)
m
r
z
O
u
j
u
m
w
j
w
m
参照平面问题的三角形单元位移函数，轴对称问题的三结点三角形单元位移函数取为，
1 2 3
4 5 6
u a a r a z
w a a r a z
有限单元法按照平面问题三角形单元的分析过程，将结点坐标和结点位移代入（ 4-4）得到，
1
2
3
1
2
i j m i
i j m j
i j m m
a a a a u
a b b b u
A
a c c c u
4
5
6
1
2
i j m i
i j m j
i j m m
a a a a w
a b b b w
A
a c c c w
其中：
mm
jj
ii
zr
zr
zr
A
1
1
1
2
1
有限单元法
i j m m ja r z z r i j mb z z i m jc r r
j m i i ma r z z r j m ib z z j i mc r r
m i j j ia r z z r
m i jb z z m j ic r r
定义形函数为：
1 ( ) (,,)
2s s s sN a b r c z s i j mA
用矩阵表示的单元位移为：
0 0 0
0 0 0
i
i
i j m j
i j m j
m
m
u
w
N N N uu
N N N ww
u
w
dN δ
有限单元法
4.2.2 单元刚度矩阵轴对称问题的几何方程：
r
z
zr
u
r
u
r
w
z
uw
zr
由（ 4-9）式得，
)(2 1 mmjjii ubububAru
)(2 1 mmjjii ufufufAru
有限单元法其中，
(,,)s s sssa c zf b s i j mrr
)(2 1 mmjjii wcwcwcAzw
)(2 1 mmjjii ucucucAzu
)(2 1 mmjjii wbwbwbArw
用几何矩阵表示单元的应变，
ε Bδ?
[]i j m?B B B B
0
01
(,,)
02
s
s
s
s
ss
b
f
s i j m
cA
cb
B
有限单元法由于 在是坐标 r,z 的函数,分量在单元中不为常量，其它三个应变分量在单元中仍为常量。
sf
由轴对称问题的物理方程，得到弹性矩阵 D，
10
11
10
11( 1 )
( 1 ) ( 1 2 )
10
11
12
000
2( 1 )
E
D
有限单元法由弹性矩阵 D 和几何矩阵 B 可以得到应力矩阵 S，并计算出单元内的应力分量，
[]i j mS D B S S S
其中：
1
11
1
22
( 1 )
(,,)
()2( 1 ) ( 1 2 )
s s s
s s s
ss
s s s
ss
b f A c
A b f A cE
s i j m
A b f cA
A c A b
S D B
式中：
1 1A
2
12
2 (1 )A
有限单元法有了单元应力场和应变场，可以利用虚位移原理或最小势能原理建立单元刚度矩阵 k?
单元刚度矩阵的分块矩阵为，
2 T rdrdzK B D B?
2 (,,,)Tsq s q rd rd z s q i j mK B D B
由于几何矩阵中的元素不是常量，单元刚度矩阵需要通过积分得到，为简化计算，可以用三角形单元形心位置的坐标 代替 B 矩阵中的变量 。
cc zr,
B
,rz
)(31 mjic rrrr )(31 mjic zzzz
有限单元法应变矩阵（ 4-13）变成，
* * * *i j mB B B B
其中,0
01
(,,)
2
0
s
s
s s c
cs
s
ss
b
a
b c z
r s i j m
A
c
cb
*
B
有限单元法单元刚度矩阵的近似表达式为：
2 Tcr* * *K B D B?
单元刚度矩阵的分块矩阵近似表达式为，
1 2 1 2
1 2 2
2
( ) ( )
()
T
s q c s q
s q s q s q s q s q s q s q s q
s q s q s q s q s q
r
b b f f A b f f b A c c A b c f c A c b
A c b c f A b c c c A b b
**
* * * * *
*
K B D B
这里：
* (,,)s s s
ss
a c zf b s i j m
rr
有限单元法
4.2.3 等效结点荷载
（ 1） 集中力集中力的处理很简单，一般直接把集中力作用点取为结点，不需要作特殊处理，就可以直接把集中力加入到结点荷载列阵中去。
（ 2）体积力设单元内单位体积上作用的体积力为,则移置到单元各结点的等效结点力为
2 T v rd rd zF N p?
vp
有限单元法
（ 3）表面力设单元某边上作用的表面力为,则移置到单元各结点的等效结点力为
2 Tsi ss rd sF N p?
sp
轴对称问题分析中，如果直接定义结点荷载，则荷载值是实际弹性体上绕对称轴一周的荷载的累计结果。
有限单元法第 6章平板弯曲问题的有限元分析有限单元法教学提示,对于工程中经常遇到的不规则薄板,经典理论通常显得无能为力,此时必须运用数值方法进行分析 。 本章就 平板弯曲问题的有限元分析进行了介绍,分别运用三角形单元,矩形单元,八结点四边形等参单元等单元划分形式对平板弯曲问题进行介绍,包括位移模式,单元分析,整体分析,等效结点荷载计算等方面 。
教学要求,本章要求学生重点掌握运用三角形单元进行薄板的有限元分析,包括位移模式、单元分析、整体分析、
等效结点荷载计算等。同时要熟悉矩形单元的运用，了解八结点四边形等参单元的分析过程。
有限单元法
6.1 薄板受弯分析的基本方程当薄板上受有一般载荷时,总可以把个荷载分解为两个分量：一个是作用在薄板中面内的中面荷载 （ 也称纵向荷载 ） ；
另一个是垂直于中面的法向荷载 （ 也称横向荷载 ） 。 对于中面荷载,可以认为它们沿薄板的厚度均匀分布,因而它们所引起的位移,应变和应力,可以按平面应力问题进行分析 。 横向荷载将使薄板发生弯曲和翘曲,它们所引起的位移,应变和应力,
应按薄板弯曲问题进行计算 。
当薄板弯曲时,中面所弯成的曲面称为弹性曲面,而中面内各点在垂直于中面方向的位移称为挠度 。
本节只介绍薄板弯曲的小挠度理论,即薄板虽然很薄,但仍具有相当的弯曲刚度,因而它的挠度远小于其厚度 。
有限单元法
6.1.1 基本假设分析薄板弯曲的挠度问题时,和材料力学中分析直梁的弯曲问题时相似 （ 薄板的中面相当于直梁的轴线,薄板的弹性曲面相当于直梁的挠曲线 ）,也采用一些由实践经验得到的基本假设,使问题大大简化,但同时又能在一定程度上反映实际情况 。 这些基本假设是：
（ 1） 薄板的法线变形后没有伸缩 。
（ 2） 变形前的中面法线在变形后仍是弹性曲面的法线 。
（ 3） 薄板中面内各点,没有平等于中面的位称 。
（ 4）忽略挤压应力所引起的变形。
有限单元法
6.1.2 几何方程
t / 2
t / 2
O
x
y
z
取薄板的中面为 xy面,z轴垂直于中面，如图 6.1所示。
由假设（ 1）可知，,再由,可得 。
也就是说，中面法线上的所有各点具有相同的位移,即弹性曲面的挠度。
0x 0z wz (,)w w x y?
w
有限单元法根据假设（ 2），可以推知：薄板的法线（ z方向线段）
与 x方向或 y方向的线段保持垂直，即没有剪应变，也就是：
0,0y z z xvuz y x z
上式也可写成：
,v w u wz y z x
对式（ 6-1）进行积分，注意到 只是 x和 y的函数，不随 z而变，因而得：
w
1
2
(,)
(,)
u z f x y
x
v z f x y
y
（ 6-2）
（ 6-1）
有限单元法假设（ 3）可以表示为,,代入到式（ 6-2）得：
0 0zu 0 0zv
1 (,) 0f x y?
2 (,) 0f x y?
于是式（ 6-2）就简化为：
,wwu z v zxy （ 6-3）
有限单元法现用挠度来表示应变，不难得到：
2
2
2
2
2
2
x
y
xy
wu
xx
vw
z
yy
uv w
yx xy
ε
这就是弯曲薄板的应变与挠度之间的几何方程。
（ 6-4）
有限单元法在小变形情况下,和 分别为弹性曲面在 x和
y方向的曲率 和,而 为弹性曲面在 x和 y方向的扭率 。 这三个参数称作弹性曲面的弯扭变形分量,它们完全确定了薄板内各点的应变分量 。 用矩阵可表示为：
（ 6-5）
将上式代入到式 （ 6-4）,得到：
（ 6-6）
从上式中可以看出，薄板内所有各点的应变分量都可由弹性曲面的弯扭变形求出。因此，有时也把式（ 6-5）称作薄板弯曲问题的几何方程。
22/wx 22/wy
x?y? 22/w x y
xy?
2 2 2
22 2
TT
x y x y
w w w
x y x y
χ
z?ε χ
有限单元法
6.1.3 物理方程假设（ 4）说明：可以忽略挤压应力 引起的变形。因此薄板内各点的应变分量可用应力分量来表示，即：
z?
xy
x E
yx
y E
2 (1 ) xy
xy E
这和薄板平面应力问题中的物理方程相同，由式（ 6-7）解出应力，可得：
（ 6-7）
有限单元法
2 ()1x x y
E
2 ()1y y x
E
2 (1 )x y x y
E
将式（ 6-4）代入上式，得到用挠度表示的应力分量，用矩阵来表示可写成：
2
2
2
2
2
2
x
y
xy
w
x
w
zz
y
w
xy
ζ DD χ
（ 6-8）（ 6-8）
有限单元法式中：
2
10
10
1
1
00
2
E
D （ 6-9）
在薄板的弯曲问题中，由于大多数情况下，都很难使得应力分量在薄板的侧面上（板边上）精确地满足应力边界条件，
而只能使这些应力分量所组成的内力整体地满足边界条件，因此，有必要考察薄板横截面上的内力。
有限单元法
x
t
/
2
t
/
2
yz
O
1
1
y x
x y
y
x
t
/
2
t
/
2
x
yz
O
1
1
从薄板内取出一个平行六面体，它在 x和 y方向上具有单位宽度，在 z方向的高度为 t（图 6.3）。在 x为常量的横截面上，
作用有 和 。由于 和 都和 z成正比，所以它们在薄板全厚度上的代数和分别等于零。只可能分别合成弯矩和扭矩。用表示由 所合成的单位宽度上的弯矩，则得：
x? xy?x?xy?
xM x?
有限单元法
/2
/2
t
xxtM z d z
将式（ 6-8）中的 的表达式代入上式，并对 z积分，得到
x?
3 2 2
2 2 2()1 2 ( 1 )x
E t w wM
xy
xy?与此相似，应力分量 将合成扭矩
/2
/2
t
x y x ytM z d z
将式（ 6-8）中的 表达式代入上式，并对 z积分，得到
xy?
32
1 2 (1 )xy
E t wM
xy?
（ 6-11）
（ 6-10）
有限单元法同样，在 y为常量的横截面上，每单位宽度上的 和 也分别合成如下的弯矩和扭矩
y? yx?
3 2 2/2
2 2 2/2
32/2
/2
()
1 2 ( 1 )
1 2 ( 1 )
t
yy t
t
y x y xt
E t w w
M z d z
yx
E t w
M z d z
xy
在这里可以看到，由剪应力 和 的互等关系，得到扭矩 和 的互等关系。
将式（ 6-10）、（ 6-11）和（ 6-12）中的第一式合并起来，用矩阵表示，则有
yx?xy?
xyM yxM
（ 6-12）
有限单元法
22
22
2
/
/
2/
x
y f f
xy
M w x
M w y
M w x y
M D D χ
其中
3
2
10
10
1 2 ( 1 )
1
00
2
f
Et
D
是薄板弯曲问题的弹性矩阵，它等于平面应力问题中的弹性矩阵乘以 。式（ 6-13）表示了薄板的内力与应变两者之间的关系，因而是薄板弯曲问题中的物理方程。
3 /12t
（ 6-13）
（ 6-14）
有限单元法
t
/
2
t
/
2
z
d
z
x
y
z
x
M
x
y
M
x y
M
y x
M
y
O
x y
y x
1
1
x
y
z
M
x
M
x
M
y
M
y
M
x y
M
x y
M
y x
M
y x
O
弯矩和扭矩,,,的方向及其作用面的位置示于图 6.4a中 。 按右手螺旋法则用双箭头矢量来表示力偶,如图 6.4b所示 （ 图中所示各力偶的方向均为正 ） 。
从式（ 6-13）中解出,,,
再代入式（ 6-8），就得到各应力分量与内力之间的关系式。
xM yM xyM yxM
22/wx 22/wy 22/w x y
图 6.4（ a）
（ b）
有限单元法
3
12TT
x y x y x y x y
z M M M
tζ
从式（ 6-13）中解出,,,
再代入式（ 6-8），就得到各应力分量与内力之间的关系式。
22/wx 22/wy 22/w x y
（ 6-15）
有限单元法
6.2 三角形单元
y
z
x
1
2
3
w
O
x 1
y 1
有限单元法
6.2.1 位移模式问题为保证挠度 w 为坐标的全三次多项式，从帕斯卡三角形可知必须要有 10项，但三角形 3个结点只能有 9个自由度，若舍去三次项中的任一项，显然都无法保证对坐标的不变性，为此 Tocher提出了一种解决方案如下：
2 2 3 2 2 31 2 3 4 5 6 7 8 9()w x y x x y y x x y x y y
但是当三角形的二边分别平行坐标 x轴,y轴时，从式（ 6-16）
无法通过结点的位移条件来确定广义坐标参数,为此必须在离散化时设法避免边界（内部）单元不出现上述现象，分析结果表明此位移模式对应的单元能得到很好的结果。
（ 6-16）
有限单元法另一解决方案是 Adini提出的
2 2 3 2 2 31 2 3 4 5 6 7 8 9w x y x y x x y x y y
此方案由于舍去了二次项 xy，故常扭率 无法得到保证，从而位移偏小，精度很差。这也可以从挠曲函数的
Taylor级数展开来说明，由于不包含完全二次项，故只能有一阶精确度。
2w
xy
有限单元法还有一种 Bell提出的解决方案，他对三角形单元除了
3个角点作为结点外，取形心点挠度作为位移参数，从而取位移模式为全三次多项式利用 10个位移参数条件确定广义坐标参数，通过分析建立起单元刚度、等效结点荷载矩阵，在整体分析之前采取如下措施消去内点自由度。
2 2 3 2 2 31 2 3 4 5 6 7 8 9 1 0w x y x x y y x x y x y y
有限单元法单元刚度方程
（ 6-19）
式中从式 （ 6-19） 的第二个方程可解出：
（ a）
将式 （ a） 代回式 （ 6-19） 第一个方程并进行整理,可得：
,111 1 1 2 1
,1 0102 1 2 2
eq
eq
Fδkk S
Fδkk 0
1 1 1 1 2 2 2 3 3 3 Tx y x y x yw w w wδ
110 22,10 21 1eqδ k F k δ
111 1 1 2 2 2 2 1 1 1,1 1 2 2 2,1 0() e q e qk k k k δ S F k k F
有限单元法设单元刚度矩阵：
等效结点荷载,（ 6-20）
则单元刚度方程为 （ 6-21）
通过静力凝聚最后获得了 9自由度的三角形单元。但是
Zienkiowicz曾指出这样所得到的单元不能保证收敛性，因此
，通常采用另外一种位移模式，即下面介绍的 面积坐标 下的位移模式。
11 1 1 2 2 2 2 1k k k k k?
1,1 1 2 2 2,1 0e q e qF F k k F?
11k δ SF
有限单元法
6.2.2 面积坐标下的位移模式
P
A
2
3 ( x
3
,y
3
)
x
y
O
A
1
A
3
L
3
= 0
L
2
= 0
L
1
= 0
1 ( x
1
,y
1
)
2 ( x
2
,y
2
)
设面积坐标变量为 。
其定义如下：如图 6.6所示，三角形单元中任意一点 P的位置可以由如下三个比值来确定：
1 2 3L L L、,
11 AL A?
22 AL A?
33 AL A?
A为三角形单元的面积，
为小三角形面积
1 2 3A A A、,
有限单元法根据三角形面积公式,可以得到面积坐标与直角坐标的关系为：
式中
1 1 1
1
2 2 2
2
3 3 3
3
2
2
2
a b x c y
L
A
a b x c y
L
A
a b x c y
L
A
1 2 3 3 2a x y x y 1 2 3b y y 1 3 2c x x
2 3 1 1 3a x y x y 2 3 1b y y 2 1 3c x x
3 1 2 2 1a x y x y 3 1 2b y y 3 2 1c x x
有限单元法
Zienkiowicz等采用面积坐标解决了直角坐标下所产生的困难，提出了如下的位移模式 22
1 1 2 2 3 3 4 1 2 1 2 3 5 1 2 1 2 3
222
6 2 3 1 2 3 7 3 2 1 2 3 8 3 1 1 2 3
2
9 1 3 1 2 3
11
( ) ( )
22
1 1 1
( ) ( ) ( )
2 2 2
1
()
2
w L L L L L L L L L L L L L
L L L L L L L L L L L L L L L
L L L L L
利用结点的位移参数条件可以确定广义坐标参数 α,从而建立形函数（当然也可象前面介绍的那样由形函数性质来确定）如下
（ a）
有限单元法
2 2 2 2
1 1 1 2 1 3 1 2 1 3
22
1 3 1 2 1 2 3 2 3 1 1 2 3
22
1 3 1 2 1 2 3 2 3 1 1 2 3
11
( ) ( )
22
11
( ) ( )
22
x
y
N L L L L L L L L L
N b L L L L L b L L L L L
N c L L L L L c L L L L L
式中脚标轮换规则为 1→ 2→ 3→ 1,如果注意到：
1 2 3 1 2 1 2 1 3 1 3 2 3 2 3( 1 ) ( 1 ) ( 1 )L L L L L L L L L L L L L L L
则式（ a）可改写成：
22
1 1 2 2 3 3 4 2 3 5 3 1 6 1 2 7 2 3 3 2
2 2 2 2
8 3 1 1 3 9 1 2 2 1
()
( ) ( )
w L L L L L L L L L L L L L
L L L L L L L L
有限单元法从式 （ c） 出发,可按如下两步法确定形函数：
（ 1） 以,,作为结点自由度 （ 位移参数 ）,求对应它们的形函数；
（ 2）利用关系式
w 1 1wL 2 2wL
11 1 2 1 2
2 1 1 1 1
2
y
x
w w
L c b c bx
wc b c bw
yL
将（ 1）中的,变换成,,然后进行合并整理即可得到对于（,,）结点位移参数的形函数。
1? 2? x? y?
w x?
y?
有限单元法下面具体推导如下：
（ 1） 由,,的位移条件可得
、,(e)
（ 2） 由式 （ c） 求导可得
1 11Lww 2 21Lww 3 31L
ww
11w 22w 33w
2
1 1 3 4 2 5 3 1 6 2 7 2 3 2
1
2 2 2
8 3 1 1 3 9 2 2 1
( ) ( 2 )
( 4 ) ( 2 )
w w w L L L L L L L
L
L L L L L L L
22
2 2 3 4 3 2 5 1 6 1 7 2 3 2 3
2
22
8 1 3 1 9 2 1 1
( ) ( 4 )
( 2 ) ( 2 )
w w w L L L L L L L L
L
L L L L L L
(f)
有限单元法
（ 3）建立如下位移条件
11 11 1 3 5 81L ww
21 1 2 1 3 4 6 7 91L ww
31 13 1 3 5 81L ww
12 2 1 2 3 5 6 8 91L ww
22 2 2 2 3 4 71L ww
32 2 3 2 3 4 71L ww
（ g）
有限单元法并由此求得式：
4 23 22
1 ()
2
5 1 3 1 1
1 ()
2
6 12 11 21 22
1 ()
2
7 3 2 22 23
1 ()
2ww
8 1 3 11 13
1 ()
2ww
9 2 1 1 1 1 2 2 1 2 2
1 ()
2ww
（ h）
有限单元法
（ 4）将式（ e）和式（ h）代回式（ b）并整理，可得：
1 1 1 1 1 1 2 1 2 1 2 2 1 2 1 1 2 2 2 2 3 3 1 3 2 2 2 3 2 1w N w N N N w N N N w N N
式中：
…………
2 2 2 21 1 1 2 1 3 1 2 1 3N L L L L L L L L L
2 2 2 21 1 1 2 3 1 1 2 2 1 3 1 1 31 1 1 1( ) ( )2 2 2 2N L L L L L L L L L L L L
222 1 1 2 1 2 2 111 ()22N L L L L L L
（ j）
（ i）
有限单元法
（ 5）将式（ d）代入式（ i）并再整理即可得：
1 2 3wN δ N N N δ
其中：
( 1,2,3 )i i x i y iN N N iN
1 2 1 1 1 2 1xN b N b N
1 2 1 1 1 2 1yN c N c N （第二脚标轮换） 1→2→3→1
不难验证所得形函数与式（ 6-22）完全相同。但这种推求方法显然比普通广义坐标法直接求 而后得形函数要方便得多。
19
（ 6-23）
有限单元法
6.2.3 单元分析无论是以式（ 6-16）还是用式（ 6-22）、（ 6-23）建立单元的挠度场，有了 w就可用常规方法来进行单元列式。下面以面积坐标的形函数和式（ 6-23）来推导并给出三角形单元的显式单元刚度矩阵等等。
有限单元法注意到面积坐标,设取,为独立坐标，则：1 2 3 1L L L 1L 2L
112
12
2
1
2
Lbbx
ccA
y L
由此不难得到：
2
2
2
2
1
22
2 2 2
2
2 2
12
1
4
2
Lx
y A L
xy LL
T
式中 A是三角形单元中面面积。
（ 6-24）
（ 6-25）
有限单元法
22
1 2 1 2
22
1 2 1 2
1 1 2 2 1 2 2 1
2
2
2 2 2 ( )
b b b b
c c c c
b c b c b c b c
T
因此，可得形变矩阵为：
设：
1
2
2
2
2
22
1 2 3 1 2 32 2 2
2
2 2
12
1
4
2
L
x
w
y A L
xy LL
T N N N δ B B B δ B δ
（ 6-27）
（ 6-26）
有限单元法式中：
1
2
2
2
22
2
2
12
1
( 1,2,3 )
4
ii
L
i
AL
LL
B T N
因为 是面积坐标三次式，因此 是面积坐标一次式。
若记
12(,)i LLN iB
1 2 3
00
00
00
L L L
φ
Φ
（ a）
（ 6-28）
有限单元法则可得
2
1
4iiABTΦ G （ 6-29）
式中
1 2 3 TT T Ti i i iG G G G
22
1 3 2 1 3 2 1
22
6( ) 2 ( ) 2 ( )
62
6( ) 2 ( 2 ) 2 ( 2 )
i i i i i i
i i i i i i
i i i i i i
b z c z
b b b c c c
bc
G
2 1 3 2 1 3
2 1 1
11
62
6( ) 2 ( 2 ) 2 ( 2 )
6( ) ( 2 ) ( 2 )
i i i i i
i i i i i i i
i i i i i i
b b b c c c
bc
bc
G （ d）
（ c）
（ b）
有限单元法
11
22
11
3
22
11
2
11
( 2 3 ) ( 2 3 )
22
26
11
( 3 2 ) ( 3 2 )
22
11
( 2 3 ) ( 2 3 )
22
26
11
( 2 3 ) ( 3 2 )
22
11
( 2 3 ) ( 2 3 )
22
62
1
( 3 2 )
2
i i i i i i
i
i i i i i i
i i i i i i
ii
i i i i i i
i i i i i i
i
i i i
bc
bc
bc
bc
bc
b
G
2
1
( 3 2 )
2
i i i
c
其中,,为 i结点的 3个面积坐标，也即 i结点的。i?
i? i?
（ e）
有限单元法若将单元刚度矩阵写成分块子矩阵形式，则对各向刚性薄板若将单元刚度矩阵写成分块子矩阵形式，则对各向刚性薄板
4
4
10
10
1
00
2
10
10
16
1
00
2
16
T
ij i j
A
T T T
ij
A
TT
ij
A
D dA
D
dA
A
D
dA
A
k B B
G Φ TT Φ G
G Φ H Φ G
（ 6-30）
有限单元法式中：
3
21 2 (1 )
EtD
33ijHH
其中：
2221 1 1 1H b c
221 2 2 1 1 2 1 2 1 2 2 1H H b b c c b c b c
2213 31 1 1 1 2 1 22H H b c b b c c
222 3 3 2 2 2 1 2 1 22H H b c b b c c
2222 2 2 2H b c
222 2 2 23 3 1 1 2 2 1 2 1 2 1 2 2 12 2 2H b c b c b b c c b c b c
有限单元法因为
11 12 13
21 22 23
31 32 33
T
H H H
H H H
H H H
φ φ φ
Φ H Φ φ φ φ
φ φ φ
式中
2
1 1 2 1 3
2
2 1 2 2 3
2
3 1 3 2 3
T
L L L L L
L L L L L
L L L L L
φ φ φ
有限单元法利用
1 2 3
! ! ! 2
( 2 ) !A L L L d A A
并记：
2 1 1
1 2 1
1 1 2
P
则单元刚度子矩阵 为ijk
3 3 3
1 1 2 2 3 33
1 1 1192
T T T
i j i r j r i r j r i r j r
r r r
D
k G P H G G P H G G P H G
（ 6-33）
（ 6-31）
（ 6-32）
有限单元法若单元上受垂直中面的分布荷载 作用，则等效结点荷载为
(,)q x y
33
11
(,) (,)TTe q i i i i
ii
q x y d A q L x L y d A
F N N
当 q为常数时，有式
2 3 2 3 3 1 3 1
1 2 1 2
1 1 1 1 1 1
( ) ( ) ( ) ( )
3 2 4 2 4 3 2 4 2 4
1 1 1
( ) ( )
3 2 4 2 4
eq
T
q b b c c b b c c
b b c c
F
（ 6-35）
（ 6-34）
有限单元法经推导，若记（各向同性时）
2 2 2 2
1 1 2 2 1 2 1 23
2 2 2 2
1 1 2 2 1 2 1 2 3 322
1 1 2 2 1 2 2 1
2 ( )
2 ( ) [ ]
4 8 ( 1 ) ( 1 ) ( 1 ) ( 1 ) ( ) ij
b c b c b b c c
Et b c b c b b c c R
b c b c b c b c
R
3
1
1
3
2
1
3
3
1
r i r
r
i r i r
r
r i r
r
R
R
R
φ G
S φ G
φ G
（ 6-37）
（ 6-36）
有限单元法式中,则薄板内力矩阵可写作：
（ 6-38）
由此式可见,的每一列的物理意义为,单位位移引起板中任一点的内力（点位置取决于,和 ）,。
1 2 3L L L?φ
33
11
[] Ti i i i x i y i
rr
w
MS δ S
iS
1L 2L 3L
有限单元法结束
2006年有限单元法
(Finite Element
Method)
有限单元法第１章 绪论有限单元法
1.1 概述问题的求解方法：解析法,数值法,（ 半解析法 ）
解析法,问题的求解可以得到具体的数学表达式,可以计算任意位置的问题的精确解答 。 如悬臂梁求挠度等 。
数值法,得不到具体表达式,只能得到某些离散点处的近似值 。 该方法包括：有限差分,有限元法,边界元法
,变分法等有限单元法有限元法,将结构物看成由有限个划分的单元组成的整体,以单元结点上的值作为整个单元的平均值 。 它是一种化整为零,集零为整,化未知为已知的方法 。 不同的学科,所求解的参数不同 。 在结构力学中,主要有以下三种：
● 位移型：以结点位移为未知量 。
● 力型：以结点力为未知量 。
● 混合型：某些地方以结点位移为未知量,另外一些以结点力为未知量 。
我们主要就,位移型,有限元进行讲解 。
有限单元法有限元方法的实质,将复杂的连续体划分为有限多个简单的单元体,化无限自由度问题为有限自由度问题,将连续场函数的 （ 偏 ） 微分方程的求解问题转化成有限个参数的代数方程组的求解问题 。
有限元方法的基本思想,先化整为零、再积零为整，
也就是把一个连续体人为的分割成有限个单元，即把一个结构看成由若干通过结点项链的单元组成的整体，先进行单元分析，然后再把这些单元组合起来代表原来的结构进行整体分析。
有限单元法
1.2 有限元法的发展
1) 早在公元 3世纪的时候,我国数学家刘徽提出的用割元法求圆周长的方法即是有限元基本思想的体现 。
2) 20世纪 40年代,麦克亨利 （ McHenry）,雷尼柯夫 （
Hrenikoff）,纽马克 （ Newmark） 等首次提出用框架方法求解力学问题,用简单弹性杆排列代替连续体的各个小部分,能够得到连续问题的相当好的解答 。
3) 1943年,柯兰特 （ Courant） 第一次假设饶曲函数在一个划分的三角形单元集合体的每个单元上为简单线性函数 。
4) 1955年，德国斯图加特大学的 J.H,Argyris教授发表了一组能量原理与矩阵分析的论文，奠定了有限元方法的理论基础。
有限单元法
5) 第一个尝试,1956年,特纳 （ Turner）,克拉夫 （ Clough）
等将刚架分析中的位移法扩展到弹性力学平面问题,并用于飞机的结构分析和设计,系统研究了离散杆,梁,三角形的单元刚度表达式,并求得了平面应力问题的正确解答 。
6) 第一次提出,1960年,克拉夫 （ Clough）,建立在虚位移原理或最小势能原理的基础上 。 在处理剖面弹性问题时,
第一次提出并使用,有限元方法,的名称
7) 1972年,Oden出版了第一本处理非线性连续体的专著 。
8) 我国科技工作者：如胡昌海提出了广义变分原理,钱伟长最先研究了拉格朗日乘子法与广义变分原理之间的关系,
冯康研究了有限元方法的精度和收敛性问题,钱令希研究了余能原理等 。
9) 基于变分原理的有限元法得到发展。
有限单元法
1.3 有限元法的分析过程
1,结构物的离散将待分析的结构物从几何上用线或面划分为有限个单元,
其中单元的大小和数目根据计算精度的要求和计算机容量来决定 。 其步骤：
● 建立单元
● 对单元和结点编号
● 准备必需的数据信息
● 建立坐标系有限单元法
x
y
x
y
1 2
3 4
5 6
① ②
③
④ ⑤
⑥
图 1.1
如图 1.1所示,可以将杆系结构分成 6
个单元,这样划分以后,共有 6个结点 。
如图 1.2所示,纵向均匀受拉的带圆孔的薄板,
根据对称性，取其中一部分分析，将其划分为三角形单元。
图 1.1
图 1.2
有限单元法
2,确定单元的位移模式将单元中任意一点的位移近似地表示成单元结点位移的函数,即位移模式或位移函数,用 或 表示,写成：
这里,— 单元中任意一点的位移矩阵,
— 形函数矩阵
— 单元结点位移矩阵 。
位移函数的假设合理与否，直接影响到分析的计算精度、
效率、可靠性。
d d?
e?dN δ
d
N
eδ
有限单元法
3,单元特性分析
（ 1） 几何方程：应变与位移之间的关系这里,— 单元中任意一点的应变矩阵,
— 变形矩阵或应变矩阵,
（ 2） 物理方程：应力与应变之间的关系 （ Hooke定律 ）
这里,— 单元中任意一点的应力矩阵,
— 弹性矩阵,由单元材料的弹性常数确定,(弹性模量 )
— 应力矩阵。
e?ε B δ
ε
B
eeζ D B δ S δ
ζ
D
S
有限单元法
（ 3） 利用虚位移原理或最小势能原理建立单元刚度方程这里,— 单元结点力矩阵
— 单元等效荷载矩阵
— 单元刚度矩阵，
e e e eEk δ F F
eF
eEF
ek
有限单元法
4,建立表示整个结构结点平衡的方程组这里,－整体刚度矩阵
－整体结点位移矩阵
－直接结点荷载
－等效结点荷载
－整体综合结点荷载矩阵
deK Δ P P P
K
Δ
dP
eP
P
有限单元法
1.4 常用软件
（ 1） 软件有限元软件可以分为通用软件和专用软件两类。通用软件适应性广，规格规范，输入方法简单，有比较成熟齐全的单元库，大多提供二次开发的接口。但通用程序功能再强，
对于一些比较专业的问题，尤其是处于研究阶段的内容，往往显得无能为力，因此，针对某些特定领域、特定问题开发的专用软件，在解决专有问题时显得更为有效。目前常用的有限元软件有,ANSYS,MARC,ABQUS,NASTRAN,ADINA、
ALGOR,SAP,STRAND,FEPG等。
有限单元法
（ 2） 程序设计低级语言：机器语言,汇编语言等高级语言,BASIC,C,C++,FORTRAN,JAVA,LISP，
Matlab,Mathmatic等目前所有有限元软件都是用 FORTRAN语言编写的 。 程序设计步骤：
（ 1） 提出问题,拟定解决方案；
（ 2） 构造数学模型；
（ 3） 画出程序流程图；
（ 4） 用选定的算法语言编写程序；
（ 5） 编译调试程序；
（ 6） 试验验证程序；
（ 7） 打包发行；
（ 8）软件维护。
有限单元法第 2章杆系结构的有限元分析有限单元法
2.1 概述所谓 杆件 是指从构造上来说其长度远大于其截面尺寸的一维构件 。 在结构力学上我们通常将承受轴力或扭矩的杆件称为杆,而将承受横向力和弯矩的杆件称为梁 。 在有限单元法中这两种情况的单元分别称为杆单元和梁单元 。 为方便起见,本书都称之为杆单元 。
杆系结构是最简单的一类结构，也是我们在工程上最常见的一类结构。本章以此类结构为基础介绍有限单元法的分析过程。
有限单元法第一步,对结构物进行离散化,划分为有限个单元 。
第二步,对各结点和单元进行编码 。
第三步,建立整体坐标系和各单元的局部坐标系 。
第四步,对已知参数进行准备和整理 。
第五步,对结点位移进行编码,注意前处理法与后处理法的区别 。
第六步,进行单元分析,形成单元刚度矩阵 。
第七步,进行整体分析,形成整体刚度矩阵 。
第八步,引入边界条件 。 边界条件的引入可以使问题具有解的唯一性 。
第九步,求解方程组,计算结构的整体结点位移 。
第十步，求单元内力，对计算成果进行整理、分析，用表格、图示出所需的位移及应力。
有限单元法
1
2
3 4
5
6 1 2
3 4 5
3
1
2
3
( 0 0 0 )
( 0 0 1 )
( 2 3 4 )
4 ( 5 6 7 )
6 ( 1 1 1 2 1 3 )
5 ( 8 9 1 0 )
6
54
21
( 1 2 3 )
2 ( 4 5 6 )
( 7 8 9 )
4 ( 1 0 1 1 1 2 )
6 ( 1 6 1 7 1 8 )5 ( 1 3 1 4 1 5 )
3
1
3
6
54
21
图 2.1 弯曲杆件系统 图 2.2 截面连续变化杆件系统图 2.3 单元位移编码前处理法 后处理法有限单元法
2.2 局部坐标系中的杆单元分析
2.2.1 拉压杆单元
q ( x )?F
i?F j?u i?u j
y
x
i j
图 2.4 拉压杆单元示意图设杆单元长度为,横截面面积为,单元材料的弹性模量为,在局部坐标系中杆端荷载分别为 和,杆端位移分别为 和,单元上的轴向分布荷载为 。
l A
E iF jF
iu ju ()qx
有限单元法
① 单元位移模式 。 用结点位移表示单元上任意截面的位移 。 对拉压杆单元,可以取其位移为一次多项式,即：
由位移的边界条件：
可得系数,为：
这样,任意截面的位移 为：
用矩阵表示为：
其中
u
()u x a b x
(0) iuu? () ju l u?
a b
iau? jiuub
l
( ) (1 ) ijxxu x u ull
i
i i j j i j
j
uu N u N u N N
u
N δ?
1i xN l j xN l?
（ 2-1）
（ 2-2）
有限单元法
② 进行应力,应变分析 。 根据材料力学中应变的定义,
有：
这里为应变矩阵 。 由虎克定律,其应力为：
11
ij
d u d BB
d x d x l l?
N δ δ δ B δ
11
ll
B
EE B δ? （ 2-4）
（ 2-3）
有限单元法
③ 求单元刚度矩阵 。 这里考虑利用虚位移原理求单元刚度矩阵,设杆端 i,j分别产生虚位移,,则由此引起的杆轴任意截面的虚位移为：
对应的虚应变为：
根据虚位移原理虚功方程,有：
将上式整理得：
iu? ju?
Tiiu u uNN δ?
B δ?
0
0
0
()
l
T
d
l
l
TT
W q x d x W
A d x
E A d x
F δ N δ
δ BB δ
外 变
00() TllT T Td q x d x E A d xFN δ δ BB δ
（ 2-5）
（ 2-6）
有限单元法式中,为局部坐标系下单元结点荷载矩阵 。 设：
则可以得到拉压杆单元的单元刚度方程为：
这里 为局部坐标系下的单元刚度矩阵,为局部坐标系下等效结点荷载矩阵,但值得指出的是：分布荷载 中可以包含集中荷载 。 根据定义,可以进一步求得单元刚度矩阵为：
T
d i jFFF?
0 ()
l T
E q x d xFN?
0
l T E A d xk B B?
dEF F k δ
k? EF?
()qx
11
11
EA
l
k? （ 2-10）
（ 2-7）
（ 2-8）
（ 2-9）
有限单元法
2.2.2 扭转杆单元
m ( x )?M
i
M j
y
x
i j
j? i
图 2.4 扭转杆单元示意图设扭转杆单元的长度为,截面惯性矩为,剪切模量为,杆端扭矩分别为,,杆端扭转角分别为,,
单元上的分布荷载集度为,则任意截面的扭转角为：
l
I
G
iM jM i? j?
()mx
( 1 ) ijxxll N δ?
有限单元法式中,为局部坐标系下扭转杆单元的结点位移矩阵 。
由材料力学可知,截面扭矩为：
式中：
我们利用极小势能原理来进行单元分析,杆单元的势能用泛函表示为：
Tijδ?
dM GI GIdx B δ?
11d
d x l l
NB
00
00
1 ( ) ( )
2
1 ( ( ) )
2
llTT
pd
llT T T
d
dM dx m x dx
dx
GI dx m x dx
F δ
δ BB δ NF δ
有限单元法这里 为局部坐标系下扭转杆单元的结点荷载矩阵 。 由极小势能原理,取上述泛函的变分,可得：
或者写为：
设：
可得扭转杆单元的单元刚度方程为：
可以看到,其形式与拉压杆单元的单元刚度方程完全一致 。 同样,由上式可以进一步求得其局部坐标系下得单元刚度矩阵为
：
T
d i jMMF?
0p
00 ()
llT T T
dG I d x m x d xδ B B N F
00( ) ( )
llTT
dG I d x m x d xBB δ NF
0
l T G I d xk B B? 0 ()
l T
E m x d xFN?
dEF F k δ
11
11
GI
l
k?
有限单元法
2.2.3 只计弯曲的杆单元
q ( x )
M
i
M
j
F
y i
F
y j?v i
v j
y
x
m ( x )
i j
j? i
设杆单元的长度为,截面惯性矩为,弹性模量为,杆端剪力为,,杆端弯矩分别为,,杆端横向位移为,,杆端扭转角分别为,,在单元上分布有荷载集度为 的竖向分布荷载和集度为 的分布力偶 。 则结点位移矩阵和结点荷载矩阵分别为：
l I E
yiF yjF iM jM
iv jv i? j?
()qx ()mx
有限单元法取挠曲线方程为 的三次多项式,即单元上任意一点的挠度为：
根据单元的位移边界条件：
时：,
时：,
Ti i j jvvδ?
Td y i i y j jF M F MF?
x
23v a b x c x d x
0x? ivv? idvdx
xl? jvv?
j
dv
dx
2 2 2
3 2 3 2
1 2 3 1
2 1 2 1
i
j
i i j j
i i j j
av
bv
c v v
l l l l
d v v
l l l l
可以得到式中的待定系数：
有限单元法将系数 a,b,c,d代入式,并将挠曲线方程用矩阵形式表示为：
式中 为形函数矩阵,其中：
为平面弯曲单元的形函数。
23
22
3 2 3 2
1 0 0 0
0 1 0 0
3 2 3 1
1
2 1 2 1
i
i
j
j
v
v x x x
vl l l l
l l l l
N δ
1 2 3 4N N N N?N
23
1 23
2
2 2
23
3 23
23
4 2
32
1
2
(1 )
32
xx
N
ll
xx
Nx
ll
xx
N
ll
xx
N
ll
有限单元法根据式 （ 2-19） 确定的单元位移场,可得单元上某一点得曲率为：
截面的弯矩为：
这里：
为平面弯曲杆单元的应变矩阵 。
根据虚位移原理 。 有：
22
22
d v d
d x d x
N δ B δ
TTM E I E I E IB δ δ B
2
2 2 3 2 2 3 2
6 1 2 4 6 6 1 2 2 6d x x x x
d x l l l l l l l l
NB
00
0
( ) ( )
ll T
d
lTT
dW q x d x m x d x W
dx
E I d x
NNF δ
δ BB δ
外 变
有限单元法则平面弯曲杆单元的单元刚度方程为：
其中的单元刚度矩阵可由式 （ 2-23） 求得为：
00( ) ( )
ll
E
dq x d x m x d xdx NFN?
0
l T E I d xk B B?
dEF F k δ
22
3
22
1 2 6 1 2 6
6 4 6 2
1 2 6 1 2 6
6 2 6 4
ll
l l l lEI
lll
l l l l
k?
记：
（ 2-23）
有限单元法
2.2.4 平面一般杆单元
M
i
M
j
F
x i
F
y i
F
x j
F
y j
u i
v i
u j
v j
y
x
j
i
杆单元的长度为,截面面积为,截面惯性矩为,
弹性模量为,单元的,端各有三个力为,,
和,,,其对应的位移为,,和,,,
建立如图 2.7所示的局部坐标系，各物理量的正向如图中所标。
则结点位移矩阵和结点荷载矩阵分别为：
l A I
E i j
xiF yiF iM
xjF yjF jM iu iv i? ju jv j?
有限单元法
Ti i i j j ju v u vδ?
Tx i y i i x j y j jF F M F F MF?
设单元上没有荷载作用,首先考虑轴向力的作用,由于杆端轴力,只引起杆端轴向位移,,根据拉压杆单元的单元刚度方程,有：
xiF xjF iu ju
x i i j
EA EAF u u
ll
x j i j
EA EAF u u
ll
有限单元法其次,杆端弯矩,和杆端剪力,只与杆端的转角位移,和杆端的横向位移,有关系,根据只计弯曲杆单元的单元刚度方程 （ 注意,由于不考虑单元上的荷载作用,
故方程式中的等效结点荷载 等于零 ） 可得：
iM jM yiF yjF
i? j? iv jv
EF?
3 2 3 2
12 6 12 6
y i i i j j
EI EI EI EIF v v
l l l l
3 2 3 2
12 6 12 6
y j i i j j
E I E I E I E IF v v
l l l l
22
6 4 6 2
i i i j j
EI EI EI EIM v v
l l l l
22
6 2 6 4
j i i j j
E I E I E I E IM v v
l l l l
有限单元法这样,上述表达式合并在一起,写成矩阵形式如下：
3 2 3 2
22
3 2 3 2
22
0 0 0 0
12 6 12 6
00
6 4 6 2
00
0 0 0 0
12 6 12 6
00
6 2 6 4
00
ixi
iyi
ii
jxj
jyj
jj
E A E A
ll
E I E I E I E I
uF
l l l l
vF
E I E I E I E I
M l l l l
uF E A E A
ll vF
E I E I E I E I
M
l l l l
E I E I E I E I
l l l l
Fkδ
可以将上式简写为：
有限单元法
2.2.5 空间杆单元
M
x i
M
y j
F
x i
F
y i
F
x j
F
y j
y
x
z
i
M
z i
F
z i
M
y i
M
x j
M
z j?F
z j
j
设局部坐标系的 轴为单元的形心主轴，横截面的两个主轴分别为 轴和 轴（如图所示）。设杆横截面面积为,杆单元长度为,在 平面内抗弯刚度为,在 平面内的抗弯刚度为,杆件的抗扭刚度为 。
x
y z A
l xz yEI xy
zEI EI
有限单元法空间刚架有 6个位移分量和 6个结点力分量,设局部坐标系下它们分别为：
Ti i i x i y i z i j j j x j y j z ju v w u v wδ?
Tx i y i z i x i y i z i x j y j z j x j y j z jF F F M M M F F F M M M?F?
有限单元法
3 2 3 2
3 2 3 2
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0
12 6 12 6
12 6 12 6
z z z z
xi
yi y y y y
zi
xi
yi
zi
xj
yj
zj
xj
yj
zj
EA EA
ll
EI EI EI EI
F
l l l l
F EI EI EI EI
l l l l
F
G I G I
ll
M
M
M
F
F
F
M
M
M
22
22
3 2 3 2
3 2 3 2
22
0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0
6 4 6 2
6 4 6 2
12 6 12 6
12 6 12 6
6 2 6 4
6
y y y y
z z z z
z z z z
y y y y
y y y y
EI EI EI EI
l l l l
EI EI EI EI
l l l l
EA EA
ll
EI EI EI EI
l l l l
EI EI EI EI
l l l l
G I G I
ll
EI EI EI EI
l l l l
22
0 0 0 0 0 0 0
2 6 4
i
i
i
xi
yi
zi
j
j
j
xj
yj
zj
z z z z
u
v
w
u
v
w
EI EI EI EI
l l l l
有限单元法
2.2.6 单元刚度矩阵的性质
① 单元刚度矩阵 为对称矩阵,其元素
② 单元刚度矩阵 中的每个元素代表单位杆端位移引起的杆端力 。 其中的任意元素 的物理意义是第 个杆端位移分量等于 1（ 其余位移分量等于 0） 时,所引起的第 个杆端力的分量值 。
k? ()i j j ik k i j
k?
ijk j
i
有限单元法
③ 一般单元的单元刚度矩阵 是奇异矩阵，它的元素组成的行列式等于零，即 。根据奇异矩阵的性质,没有逆矩阵。也就是说，如果给定杆端位移,根据（ 2-29）
或（ 2-31）式可以求出杆端力 的惟一解，但反过来，如果已知杆端力,则不能根据 来确定杆端位移的惟一解。因为即使在杆端力已知的情况下，由于单元两端无任何约束，因此除出杆端自身变形外，还可以发生任意的刚体位移。举例来说，如果物体处于静止状态，我们可以说其处于平衡状态，但反过来，如果物体处于平衡状态，则我们不能说其一定处于静止。
k?0?k?
k?
δ?
F?
1() dδ kF δ?F?
有限单元法
④ 单元刚度矩阵 具有分块的性质,即可以用子矩阵表示 。 用虚线把 分为四个子矩阵,把 和 各分为两个子矩阵,因此,又可以写为：
这里：
或或或或用子矩阵形式表示单元刚度矩阵和单元刚度方程，可以使其表达的物理意义更加明显。在单元刚度矩阵 中，其任意子矩阵 表示杆端力 和杆端位移 之间的关系。
k?
k?k? F? δ?
i ii ij i
j ji jj j
F k k δ
F k k δ
T
i x i y i iF F MF?
T
i x i y i z i x i y i z iF F F M M MF?
T
j x j y j jF F MF?
T
j x j y j z j x j y j z jF F F M M M?F?
T
i i i iuvδ?
T
i i i i x i y i z iu v wδ?
T
j j j juvδ?
T
j j j j x j y j z ju v wδ?
k?
rsk? rF? sδ?
有限单元法
2.3 杆系结构的整体分析
2.3.1 平面问题坐标转换矩阵图 2.9 平面问题杆端力转换示意图
i
j
o
y
x
M
i
M
j
F
x i
F
y i
F
x j
F
y j
y
x
F
x i
F
x j
F
y i
F
y j
M
j
M
i
i
j
o
y
x
y
x
有限单元法这里 表示由 轴到 轴的角,角度转动的正负由右手定则确定,本书中以顺时针方向转动为正 。 在两个坐标系中,力偶分量保持不变,即有：
同理,对于 端的杆端力,有：
c o s sin
sin c o s
x i x i y i
y i x i y i
F F F
F F F
x x
iiMM?
j
c os sin
sin c os
x j x j y j
y j x j y j
jj
F F F
F F F
MM
有限单元法将这些式子用矩阵形式可表示为：
c o s sin 0 0 0 0
sin c o s 0 0 0 0
0 0 1 0 0 0
0 0 0 c o s sin 0
0 0 0 sin c o s 0
0 0 0 0 0 1
xixi
yiyi
ii
xjxj
yjyj
jj
FF
FF
MM
FF
FF
MM
F T F c os si n 0 0 0 0
si n c os 0 0 0 0
0 0 1 0 0 0
0 0 0 c os si n 0
0 0 0 si n c os 0
0 0 0 0 0 1
T?
上式可以简写成：
这即为两种坐标系下单元杆端力的坐标变换式。其坐标转换矩阵为：
有限单元法从坐标转换矩阵 的表达式可以看出,为正交矩阵,其逆矩阵等于其转置矩阵,即有：
并且有：
式中 为单位矩阵 。
同样的推导,可以得到两种坐标系下的杆端位移之间的转换关系为：
这里 和 分别为局部坐标系和整体坐标系下的杆端位移矩阵,为前面介绍的转换矩阵。
T?
T?
1 TTT
1T T I
I
δ T δ
δ? δ?
T?
有限单元法因此可得：
上式两边同乘以,可以得到：
设可得：
上式即为整体坐标系下的单元刚度方程。
T F k T δ
1?T?
1 TF T k T δ T k T δ
T?k T k T
Fkδ
有限单元法
2.3.2 空间问题坐标转换矩阵
o
F
x i
F
y i
F
x j
F
y j
y
x
x
y
z
i
j?F
z j
F
z i
z
有限单元法设 轴与 x,y,z轴的方向余弦分别为：
，,
则将杆端力,,向 轴投影,可以求得杆端力,
即：
同理可以求得：
x
c o s (,)xxl x x? c o s (,)xyl x y? c o s (,)xzl x z?
xiF yiF ziF
x
xiF
x i x i x x y i x y z i x zF F l F l F l
y i x i y x y i y y z i y zF F l F l F l
z i x i z x y i z y z i z zF F l F l F l
有限单元法用矩阵形式可以表示为：
上式即为结点 i的杆端力在局部坐标系和整体坐标系下的转换关系,其中的矩阵称为关系矩阵 。
与上面的推导类似，同样可以推出以,,表示,,,以及对于结点 j的相对应的转换关系，其中转换关系矩阵都是 。
x i x x x y x z x i
y i y x y y y z y i
z i z x z y z z z i
F l l l F
F l l l F
F l l l F
x x x y x z
y x y y y z
z x z y z z
l l l
l l l
l l l
λ
xiM yiM ziM
xiM yiM ziM
λ
有限单元法综上所述,整体坐标系下单元杆端力矩阵与局部坐标系下单元杆端力矩阵具有如下的关系表达式：
其中的 为：
称为空间坐标系的单元转换矩阵,它是一个正交矩阵,即：
对于杆端位移,同样可推导出在两种坐标系中的转换关系：
这样,可得空间杆件单元在整体坐标系中单元刚度方程为：
其中表示空间单元在整体坐标系中的单元刚度矩阵 。
F T F
T
λ
λ
T
λ
λ
000
0 0 0
0 0 0
000
1 TTT
δ T δ
Fk δ
T?k T k T
有限单元法
2.3.3 整体分析对杆系结构进行单元分析，仅仅是有限元分析中的第一步。
我们的目的是要对整个结构进行分析，研究结构的整体性能。
因此，在对结构的各单元分析完成后，必须将单元分析的结果进行整合，对结构进行整体分析。整体分析的过程实际上是如何将单元分析的结果进行有效组合，建立整体刚度方程并求解结点位移的过程。根据对结点位移的编码方式，可以采用,先处理法,和,后处理法,来建立整体刚度方程。
有限单元法
2.3.3.1 后处理法所谓后处理法,就是由单元刚度矩阵形成整体刚度矩阵,建立刚度方程后再引入支承条件,进而求解结点位移的方法 。 运用这种方法时,假设所有结点位移均为未知量,按照顺序统一进行编码,如图所示的平面杆件单元 。
1 ( 1 2 3 )
2 ( 4 5 6 ) 3 ( 7 8 9 )
4 ( 1 0 1 1 1 2 )
2
1 3
x
y
O
有限单元法结点位移矩阵为：
结点荷载矩阵为：
求出各单元刚度方程后,根据平衡条件和位移连续条件,可以建立整个结构的位移法方程：
1 2 3 4
1 1 1 2 2 2 3 3 3 4 4 4
T
Tu v u v u v u v
δ δ δ δ δ
1 2 3 4
1 1 1 2 2 2 3 3 3 4 4 4
T
T
x y x y x y x yF F M F F M F F M F F M
F F F F F
11
22
33
44
ii ij
ji jj ii ij
ji jj ii ij
ji jj
kkF δ
k k k kF δ
k k k kF δ
kkF δ
00
0
0
00
①①
① ① ② ②
② ② ③ ③
③③
有限单元法简写成：
这里 为结构的整体刚度矩阵,有：
FKδ
K
1 1 1 2 1 3 1 4
2 1 2 2 2 3 2 4
3 1 3 2 3 3 3 4
4 1 4 2 4 3 4 4
K K K K
K K K K
K K K K
K K K K
K
应该注意到，在建立方程的过程中，我们假设所有结点都有位移。因此整个结构在外力作用下，除了发生弹性变形外，还可能发生刚体平动位移，这样各结点位移不能唯一确定。这说明整体刚度矩阵为一奇异矩阵，不能求逆矩阵，即根据整体刚度方程可得到无穷多个解。
有限单元法实际上,在图所示刚架中,结点 1和结点 4均为固定端,其三个位移分量均为 0,即有：
这样,将上述支承条件引入到方程中,对整体刚度方程进行修改,可得：
1 1 1 4 4 4 0u v u v
1 1 1 2 1 3 1 41
2 1 2 2 2 3 2 42 2
3 1 3 2 3 3 3 43 3
4 1 4 2 4 3 4 44
K K K K
K K K K
K K K K
K K K K
F
F δ
F δ
F
0
0
有限单元法对上述方程进行化简，可以得到两组方程：
和
2 2 2 2 3 2
3 3 2 3 3 3
KK
KK
F δ
F δ
1 2 21
4 3 34
K
K
δF
δF
0
0
这样，利用第 1式可以求得结点位移 和,再根据第 2式可以求得支座反力 和 。
2δ 3δ
1F 4F
有限单元法
2.3.3.2 先（前）处理法所谓先处理法,就是先引入支承条件,根据支承条件仅对未知的自由结点位移分量编号,得到的位移矩阵中不包含已知的约束位移分量,即可以直接得到方程求解自由结点位移分量 。
有限单元法如图所示的平面刚架结构 ABCD,由于在 A处和 D处均为固定端,其位移为 0,故位移编码均为 0,在 C处为铰接,故 BC杆在 C端的角位移与 DC杆在 C端的角位移不相同,因此在 C处编两个结点 3和 4,但结点 3和 4的横向位移和竖向位移相同,故采用相同的编号,各结点位移编码如图所示 。
图 2.12 先处理法位移编码示意图
1 ( 0 0 0 )
2 ( 1 2 3 ) 3 ( 4 5 6 )
4 ( 4 5 7 )
A
B
C
D 5 ( 0 0 0 )
2
1 3
x
y
O
有限单元法由单元刚度矩阵集成整体刚度矩阵通常采用刚度集成法。
其计算过程可以分为两步：首先求出各单元的贡献矩阵，然后将它们叠加起来形成整体刚度矩阵。但这样处理在实际中很少采用，因为在编程过程时需先将各单元的贡献矩阵储存起来，
而各单元贡献矩阵的阶数与整体刚度矩阵的阶数相同，因此占用的空间非常巨大，不利于节约资源，并且在实际中有可能耗尽所有资源。
故在实际中并不是采用贡献矩阵法，而是利用各单元的定位数组，采用,边定位，边累加,的方法。
有限单元法所谓单元的定位数组,就是将单元 的始端及末端的位移编码排成一行 （ 始端在前,末端在后 ）,写成如下的形式：
如图中,各单元的定位数组分别为：
这样处理，得到的整体刚度矩阵与叠加所有贡献矩阵得到的结果完全一致，但可以节约大量的存储空间。
1 2 3() nm m m m m
( 0 0 0 1 2 3 )m?①
(1 2 3 4 5 6 )m?②
( 0 0 0 4 5 7 )m?③
有限单元法如图所示的结构有 3个单元,5个结点,7个独立的位移分量,
这样其整体刚度矩阵应为 阶矩阵,即：77?
11 12 13 14 15 16 17
21 22 23 24 25 26 27
31 32 33 34 35 36 37
41 42 43 44 45 46 47
51 52 53 54 55 56 57
61 62 63 64 65 66 67
71 72 73 74 75 76 77
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
有限单元法对于单元,其单元定位数组为,我们将定位数组标记在单元刚度矩阵的上面和右侧,其与单元刚度矩阵一起写成如下形式：
① ( 0 0 0 1 2 3 )m?①
1 1 1 2 1 3 1 4 1 5 1 6
2 1 2 2 2 3 2 4 2 5 2 6
3 1 3 2 3 3 3 4 3 5 3 6
4 1 4 2 4 3 4 4 4 5 4 6
5 1 5 2 5 3 5 4 5 5 5 6
6 1 6 2 6 3 6 4 6 5 6 6
0 0 0 1 2 3
0
0
0
1
2
3
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
4 4 1 1 4 5 1 2 4 6 1 3
5 4 2 1 5 5 2 2 5 6 2 3
6 4 3 1 6 5 3 2 6 6 3 3
k K k K k K
k K k K k K
k K k K k K
① ① ①
① ① ①
① ① ①
有限单元法单元 ② 的刚度矩阵和它的定位数组为：
中各元素在整体刚度矩阵 中的位置为：
1 1 1 2 1 3 1 4 1 5 1 6
2 1 2 2 2 3 2 4 2 5 2 6
3 1 3 2 3 3 3 4 3 5 3 6
4 1 4 2 4 3 4 4 4 5 4 6
5 1 5 2 5 3 5 4 5 5 5 6
6 1 6 2 6 3 6 4 6 5 6 6
1 2 3 4 5 6
1
2
3
4
5
6
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
( 1,2,,6 1,2,,6 )i j i jk K i j②
有限单元法单元 ③ 的刚度矩阵和它的定位数组为：
中各元素在整体刚度矩阵 中的位置为：k
③ K
1 1 1 2 1 3 1 4 1 5 1 6
2 1 2 2 2 3 2 4 2 5 2 6
3 1 3 2 3 3 3 4 3 5 3 6
4 1 4 2 4 3 4 4 4 5 4 6
5 1 5 2 5 3 5 4 5 5 5 6
6 1 6 2 6 3 6 4 6 5 6 6
0 0 0 4 5 7
0
0
0
4
5
7
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
③
③
4 4 4 4 4 5 4 5 4 6 4 7
5 4 5 4 5 5 5 5 5 6 5 7
6 4 7 4 6 5 7 5 6 6 7 7
k K k K k K
k K k K k K
k K k K k K
③③③
③③③
③③③
有限单元法这样,按照以上所讲的定位方法,将,,中的相关元素累加到整体刚度矩阵对应的元素上,可以得到整体刚度矩阵为：
k① k② k③
4 4 1 1 4 5 1 2 4 6 1 3 1 4 1 5 1 6
5 4 2 1 5 5 2 2 5 6 2 3 2 4 2 5 2 6
6 4 3 1 6 5 3 2 6 6 3 3 3 4 3 5 3 6
4 1 4 2 4 3 4 4 4 4 4 5 4 5 4 6 4 6
5 1 5 2 5 3 5 4 5 4 5 5 5 5 5 6 5 6
6
0
0
0
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
① ② ① ② ① ② ② ② ②
① ② ① ② ① ② ② ② ②
① ② ① ② ① ② ② ② ②
② ② ② ② ③ ② ③ ② ③
② ② ② ② ③ ② ③ ② ③
1 6 2 6 3 6 4 6 5 6 6
6 4 6 5 6 6
0
0 0 0 0
k k k k k
k k k
② ② ② ② ② ②
③ ③ ③
有限单元法在实际编程过程中,集成整体刚度矩阵的过程可以概括为：
（ 1） 根据位移编号,建立整体刚度矩阵 Ｋ 的存储空间,
并置 0,若位移个数为 n,则整体刚度矩阵为 阶矩阵 。
（ 2）对每个单元，根据其定位数组，将其元素累加到整体刚度矩阵中对应的元素上，直到处理完所有的单元。
nn?
有限单元法
2.4 等效结点荷载与边界条件的处理
2.4.1 非结点荷载的处理有限单元法实际上,在前面的分析中,我们已经介绍了求等效结点荷载的方法,如 2-7式,2-15式,2-22式分别可用来求不同情况下的等效结点荷载 。
此外,可以这样来考虑：
第一步,在局部坐标系下求单元的固端力 。 对于某个单元,我们假定单元的两端均固定,然后根据静力平衡求得固定端的反力 。
第二步,根据单元固端力求单元的等效结点荷载 。 根据局部坐标系与整体坐标系单元杆端力的变换式,固端内力在两种坐标系下的变换形式可以写成：
因此,整体坐标系下的等效结点荷载矩阵可以由下式计算：
fF?
Tff?F T F
EF?
TEfF T F
有限单元法
2.4.2 边界条件的处理
2.4.2.1 铰结点在杆系结构中,除了刚性结点外,通常会遇到一些杆件通过铰结点与其它杆件联结,如图 2.13所示杆件系统,有 4根杆件汇交于 D点,其中 BD杆在 D端通过铰支座与其它杆件铰接,其余 3根杆为刚性接触 。 对于这样的铰结点,具有如下的性质：
① 铰结点上各杆具有相同的线位移,但截面的转角位移不相同；
② 结点上具有铰接杆端不承受弯矩作用 。 如图 2.13所示结构中,BD杆在 D端的杆端弯矩为 0,只有 CD,ED,GD杆在结点
D上与外弯矩保持平衡 。
有限单元法对于这样的结点,我们在对其进行单元划分时,通常考虑在 D处设置 2个结点 。 按照先处理法,对图示结构进行位移编码,
如图 2.14b所示 。
A
B
C
D
E H
G
F
1 ( 0 0 1 )
2 ( 2 3 4 )
3 ( 5 6 7 )
4 ( 5 6 8 )
5 ( 0 0 0 )
6 ( 9 1 0 1 1 )
7 ( 0 0 0 )
8 ( 1 2 1 3 1 4 )
9 ( 1 5 1 6 1 7 )
图 2.14 铰结点的处理示意图有限单元法
2.4.2.2 弹性支承点在实际工程中，有时会遇到弹性支承的情况（如图 2.15
），这时一般将弹性支座看作是在结构约束点沿约束方向的一个弹簧，弹簧的刚度系数为,在数值上等于使弹簧支座沿约束方向产生单位位移时所需施加的力。
kk
1 2 i - 1 i i + 1 n
k
有限单元法设结构的第 i个结点位移分量 为弹性支座约束,弹簧的刚度系数为,则结构产生位移时所引起的支座反力为：
这里负号表示支座反力方向与约束点位移方向始终相反 。 作用在受约束的结点上,它是结点外力的一部分,由整体刚度方程可知,第个平衡方程应为：
式中 为原有的第 i个结点上的荷载分量 。 进行整理可得：
这样就引入了弹性支承的约束条件 。 根据以上的分析,引入弹性支承的具体做法可以归结为：先解除弹性支承点约束,在 i处给一个结点号,形成总刚度矩阵,然后在总刚度矩阵中将第 i行的主元素加上弹性支承的刚度系数,此时第行变为：
i?
k
RiFk
1 1 2 2i i i i i i n n i RK K K K F F
iF
1 1 2 2 ()i i i i i i n n iK K K k K F
1 2 3i i i i i i nK K K K k K?
有限单元法
2.5 杆系结构分析算例应用有限元位移法对杆系结构进行分析,如果按照,先处理法,的思想,其步骤可以归结为如下：
（ 1） 对单元进行划分,整理原始数据,对单元和结点进行编号,确定每个单元的局部坐标系整个结构的整体坐标系 。
（ 2） 计算局部坐标系中的单元刚度矩阵 。
（ 3） 确定每个单元的坐标转换矩阵,计算整体坐标系下的单元刚度矩阵 。
有限单元法
（ 4） 根据各单元的位移分量编号,形成单元的定位数组
,按照,对号入座,同号叠加,的方法,形成结构的整体刚度矩阵 。
（ 5） 计算总的结点荷载矩阵 。 在这一步里,首先必须将非结点荷载转换成等效结点荷载,再与对应的结点荷载叠加,
形成总的结点荷载矩阵 。
（ 6） 求解结构的整体刚度方程,计算未知的结点位移矩阵 。
（ 7） 计算各单元的杆端力 。
（ 8）对计算结果进行整理，如果需要，计算内力图。
有限单元法例 1
1 ( 0 0 0 )
2 ( 1 2 3 )
3 ( 0 0 4 )
1 0 k N
8 k N6 k N
1 0 k N? m
1
2
k
N
/
m
A
B
C
2,5 m
2
1
x
y
O
2,5 m
如图所示的刚架结构 ABC，
各杆截面尺寸相同，材料性质一样，求刚架的整体刚度矩阵 。
5A B B C l m
20.5Am?
43 1 0E M P a
41
24Im?
K
杆长：
截面面积：
弹性模量：
截面惯性矩：
有限单元法
（ 1） 单元划分,建立局部坐标系和整体坐标系 （ 如图
2.13所示 ）,并对数据进行整理,对单元和结点编号 。
44
44
23
4
3 0 0 1 0 1 0 0 1 0
6 1 2
3 0 1 0 1 2 1 0
E A E I
ll
E I E I
ll
有限单元法
（ 2） 求局部坐标系中的单元刚度矩阵 。 由于单元 ①,
② 的尺寸完全一样,因此其单元刚度矩阵一样,根据 （ 2-30）
式可以得到：
k?
4
300 0 0 300 0 0
0 12 30 0 12 30
0 30 100 0 30 50
10
300 0 0 300 0 0
0 12 30 0 12 30
0 30 50 0 30 100
kk
①②
有限单元法
（ 3） 求整体坐标系中的单元刚度矩阵 。 对于单元 ①,其局部坐标系与整体坐标系的夹角为,故根据 （ 2-40） 式,
计算其坐标转换矩阵为：
单元的定位数组为：
根据 （ 2-46） 式计算单元 ① 在整体坐标系下的单元刚度矩阵为：
90
0 1 0 0 0 0
1 0 0 0 0 0
0 0 1 0 0 0
0 0 0 0 1 0
0 0 0 1 0 0
0 0 0 0 0 1
T
①
(1 2 3 0 0 0 )m?①
4
1 2 3 0 0 0
1 12 0 30 12 0 30
2 0 30 0 0 0 30 0 0
3 30 0 10 0 30 0 50
10
0 12 0 30 12 0 30
0 0 30 0 0 0 30 0 0
0 30 0 50 30 0 10 0
k
①
有限单元法对于单元 ②,由于其局部坐标系与整体坐标系一致,因此两种坐标系下的单元刚度矩阵相同,即有：
（ 4） 形成整体刚度矩阵 。 根据前面介绍的方法建立整体刚度矩阵如下：
4
1 2 3 0 0 4
1 30 0 0 0 30 0 0 0
2 0 12 30 0 12 30
3 0 30 10 0 0 30 50
10
0 30 0 0 0 30 0 0 0
0 0 12 30 0 12 30
4 0 30 50 0 30 10 0
kk
②②
4
3 1 2 0 3 0 0
0 3 1 2 3 0 3 0
10
3 0 3 0 2 0 0 5 0
0 3 0 5 0 1 0 0
K
有限单元法
（ 5） 求总结点荷载 。
首先,求各单元在局部坐标系中的固端内力 。 对于单元 ①
（ 这里应该注意,其 局部坐标系的 轴方向向左 ）,
，,根据表 2.1可得：
对于单元 ②,,,,根据表 2.1可得：
，，
因此有：
fF?
y 1 2 k N /mq
5mal 0b?
0fx i fx jFF①① 3 0 k Nf y i f y jFF①① 2 5 k N mf i f jMM①①
8 kNF? 2,5 mab 5ml?
0fx i fx jFF②② 4 k Nf y i f y jFF②② 5 k N mf i f jMM②②
0 30 25 0 30 25 TfF ①
0 4 5 0 4 5 TfF ②
有限单元法其次,求整体坐标系下各单元的等效结点荷载 。 对于单元 ①,,单元 ②,,因此可以求得各单元的等效结点荷载为：
90 0
3 0 0 0 1 2
3 0 0 2 5 3 0 0 2 5 TEF ①
3 0 01 2 4
0 4 5 0 4 5 TEF ②
第三,求刚架等效结点荷载矩阵 。 按照,同号叠加,的方法有：
第四,求直接作用在结点上的荷载,
最后,求总的结点荷载矩阵
3 0 4 2 0 5 TEF
1 0 6 1 0 0 TdF
2 0 1 0 1 0 5 TdEF F F
有限单元法
（ 6） 求解结构的整体刚度方程,计算未知的结点位移矩阵 。 这里的整体刚度方程为：δ
1
14
1
3
31 2 0 30 0 20
0 31 2 30 30 10
10
30 30 20 0 50 10
0 30 50 10 0 5
u
v
1
1 6
1
3
6,06 54
3,97 29
10
3,58 65
4,39 86
u
v
δ
解得：
有限单元法
（ 7） 计算各单元的杆端内力 。 首先从求出的结点位移中取出各单元在整体坐标系下的杆端位移,有：
δ
δ?1
1
1 6
2
2
2
6,0 6 5 4
3,9 7 2 9
3,5 8 6 5
10
0
0
0
u
v
u
v
δ ①
1
1
1 6
3
3
3
6,0 6 5 4
3,9 7 2 9
3,5 8 6 5
10
0
0
4,3 9 8 6
u
v
u
v
δ ②
然后计算杆端内力。
有限单元法
4
300 0 0 300 0 0
0 12 30 0 12 30
0 30 100 0 30 50
10
300 0 0 300 0 0
0 12 30 0 12 30
0 30 50 0 30 100
6.0654
3.9729
3.5865
0
0 1 0 0 0 0
1 0 0 0 0 0
0 0 1 0 0 0
0 0 0 0 1 0
0
0
0 0 1 0 0
0 0 0 0
0
01
ff
F k F k Tδ δ F
① ① ① ① ① ① ① ①
6
0 11,919 kN
30 28,196 kN
25 19,594 kN m
10
0 11,919 kN
30 31,804 kN
25 28,613 kN m
对于单元①：
有限单元法对于单元②,
4
6
3 0 0 0 0 3 0 0 0 0
0 1 2 3 0 0 1 2 3 0
0 3 0 1 0 0 0 3 0 5 0
10
3 0 0 0 0 3 0 0 0 0
0 1 2 3 0 0 1 2 3 0
0 3 0 5 0 0 3 0 1 0 0
6,0 6 5 4 0 1 8,1 9 6 k
3,9 7 2 9 4
3,5 8 6 5 5
10
00
04
4,3 9 8 6 5
ff
F k F k Fδ δ
② ② ② ② ② ② ②
N
5,9 1 9 k N
9,5 9 4 k N m
1 8,1 9 6 k N
2,0 8 1 k N
0
有限单元法
（ c）剪力图
（ b）弯矩图
（ d）轴力图
（ a）位移曲线图 2.17
1
2 3
2
1
1
2 3
2
8
.6
1
1
9
.5
9
-1
3
.4
0
- 9,5 9
5,2 0
1
2
3
1
.8
0
-2
8
.2
0
5,9 2
- 2,0 8
1
2 3
2
1
1
2
3
-1
1
.9
2
- 1 8,2 0
2
1
有限单元法第 3章弹性力学平面问题的有限元分析有限单元法
3.1 概述在有限元法中,把单元与单元之间设置的相互连接点,称为结点 。 一般用号码 1,2,… 进行结点编号 。 结点可为铰接,
固接或其它形式的连接 。 结点的设置,性质及数目等均视所研究问题的性质,描绘变形状态的需要和计算精度的要求等而定 。
在有限元法中引进结点概念是至关重要的。有了结点，才可将实际连续体看成是仅在结点处相互连接的单元集合组成的离散型结构，从而可使研究的对象转化成可以使用电子计算机计算的教学模型。由单元、结点、结点连线构成的集合称为有限元模型。它是有限元分析与计算的对象。
有限单元法
3.1.1 单元划分类型单元类型：三角形,四边形单元数目：根据计算精度要求来确定结点设置：使单元的的结点编号尽量靠近有限元模型：由单元,结点,结点连线构成的集合有限单元法
3.1.2 位移函数在选择多项式时,为了使有限单元法的计算精度和收敛性得到保障,还需要满足完备性和连续性的要求 。 为了使位移模式尽可能地反映物体中的真实位移形态,它应满足下列条件：
（ 1） 位移模式必须能反映单元的刚体位移；
（ 2） 位移模式必须能反映单元的常量应变；
（ 3） 位移模式应尽可能地反映位移的连续性 。
弹性力学平面问题一般选择多项式函数作为位移函数。
...26524321 yaxyaxayaxaau
.,,26524321 ybxybxbybxbbv
有限单元法
3.2 平面三角形单元三角形单元是一种简单方便、对边界适应性强的单元，
由于以三角形的三个顶点作为结点，因此又成为三结点三角形单元。这种单元的计算精度较低，使用的时候必须进行精细的网格划分，但他仍然是一种常用的单元有限单元法
3.2.1 位移函数的选取 ( x m,y m )
u
i
v
i
( x
i
,y
i
)
( x
j
,y
j
)
i
j
m
x
y
O
u
j
u
m
v
j
v
m
如图 3.4所示的 3结点平面三角形单元,结点 的坐标分别为,,,结点位移分别为,,,,,。 记单元的结点位移向量 和结点力向量 为：
i j m、,
),( ii yx ),( jj yx ),( mm yx iu iv ju
jv mu mv
δ? F?
T
i i j j m mu v u v u vδ?
T
x i y i x j y j x m y mF F F F F F?F?
有限单元法根据完备性和连续性的要求,选取 3结点三角形单元的位移场函数 如下：
（ 3-4）
将 3个结点 上的坐标和位移分别代入式 （ 3-4） 就可以将六个待定系数用结点坐标和结点位移分量表示出来 。
将水平位移分量和结点坐标分别代入 （ 3-4） 中的第一式,得到：
d?
yaxaa
yaxaau
654
321
v
i j m、,
mmm
jjj
iii
yaxaau
yaxaau
yaxaau
321
321
321
有限单元法写成矩阵形式,有：
（ 3-6）
其中,
为三角形单元的面积 。 为了避免出现 的情况,三个结点按逆时针顺序排列 。 由 （ 3-6） 得到：
1
2
3
1
1
1
i i i
j j j
m m m
u x y a
u x y a
u x y a
1
2 de t 1
1
ii
jj
mm
xy
A x y
xy
A 0A?
i j m、,
1
1
2
3
1
1
1
i i i
j j j
m m m
a x y u
a x y u
a x y u
有限单元法将竖向位移和结点坐标代入 （ 3-4） 中的第二式,可以得到：
（ 3-8）
将（ 3-7）、（ 3-8）代回（ 3-4）整理后可得：
1
4
5
6
1
1
1
i i i
j j j
m m m
a x y v
a x y v
a x y v
（ 3-9）
（ 3-10）
其中系数,
（ 下 标轮换 ）
])()()[(2 1 mmmmjjjjiiii uycxbauycxbauycxbaAu
])()()[(2 1 mmmmjjjjiiii vycxbavycxbavycxbaAv
i j m m ja x y x y
i j mb y y
i j mc x x
i j m、,
有限单元法令 （ 3-11）
可得：
（ 3-12）
单元内的位移场函数可以简写成：
（ 3-13）
把 称为形函数矩阵,称为形函数 。 根据形函数的定义,
具有以下性质：
（ 1） 在单元结点上形函数的值为 1或为 0。
（ 2）在单元中的任意一点上，三个形函数之和等于 1。
)(21 ycxbaAN iiii
0 0 0
0 0 0
i
i
i j m j
i j m j
m
m
u
v
N N N uu
N N N vv
u
v
dNδ
N
iN
iN
有限单元法
3.2.2 单元的应变场根据单元的位移场函数式 （ 3-12）,由几何方程可以得到单元的应变场表达式：
（ 3-14）
记为,（ 3-15）
其中,矩阵 称为几何矩阵 。 矩阵 可以表示为分块矩阵的形式
（ 3-16）
这里,
0 0 0
1
0 0 0
2
i
i
i j m
j
i j m
j
i i j j m m
m
m
u
u
vx
b b b
uv
c c c
vyA
c b c b c b
uuv
yx v
ε
ε Bδ?
B B
[]i j m?B B B B
0
1 0 (,,)
2
r
rr
rr
b
c r i j mA
cb
B
有限单元法
3.2.3 单元的应力场由物理方程及式 （ 3-15）,可以得到单元的应力场表达式：
（ 3-18）
其中 为应力矩阵,称为弹性矩阵,对于平面应力问题,
(3-19)
ζ D ε DB δ S δ
S DB D
2
10
10
( 1 )
1
00
2
E
D
有限单元法将应力矩阵表示为分块矩阵的形式,有：
（ 3-20）
其中：
[]i j m?S S S S
2
(,,)
2 ( 1 )
11
22
rr
r r r r
rr
bc
E
b c r i j m
A
cb
S D B
1
1
(,,)
2 ( 1 2 ) 1 1
1 2 1 2
2 1 2 1
rr
r r i r
rr
bc
E
b c r i j m
A
cb
S D B
对于平面应变问题，只需将 换为,换为,
则（ 3-21）变为：
E
21
E
1
有限单元法
3.2.4 单元刚度矩阵
（ 1） 单元的应变能
（ 3-23）
（ 2） 单元上外力的势能
1
2 TUd ε ζ
T T Tvs
AV d d A cf p f p δ F
（ 3-24）
式中,,分别表示单位体积的体积力、单元上的表面力、单元结点上的结点荷载。
vp sp cF?
有限单元法将式（ 3-13）、（ 3-15）、（ 3-18）代入（ 3-23）、（ 3-24）
后相加，得到单元的总势能为：
1
2 T T T T Tv s cAd p d p d Aδ B D B δ δ N N F
δ?
利用最小势能原理，取结点位移 的变分，得到：
0 δδ
由 的任意性，有：?δ
0
i i j j m mu v u v u v
δ
有限单元法考虑到 的对称性，对式（ 3-25）求偏导得到：
TB DB
0T T Tv s cAd d d AB D B δ N p N p F
记：
T d
B D B k
TTv s c
Ap d p d AN N F F
则式（ 3-28）可写为：
k δ F
有限单元法这就是描述单元结点力和结点位移向量之间关系的平衡方程。
其中 称为单元刚度矩阵。在 3结点等厚三角形单元中 和的分量均为常量，则单元刚度矩阵可以表示为：
k? B D
T tA?k B D B?
其中 t,Ａ 分别为单元的厚度和面积。单元刚度矩阵 可以表示为分块矩阵的形式：
k?
ii ij im
ji jj jm
m i m j m m
k k k
k k k k
k k k
有限单元法其中：
,,
,,
,,,r x s x r x s yTr s r s
r y s x r y s y
kk r s i j m
kk
k B D B
对于平面应力问题，其刚度矩阵的显式：
2
11
22,,,
114 ( 1 )
22
r s r s r s r s
rs
r s r s r s r s
b b c c b c c b
Et
r s i j m
A
c b b c c c b b
k
有限单元法
3.2.5 等效结点荷载在式（ 3-28）中括号中的前两个量分别为体积力、
表面力移置到结点上的等效结点力，依次定义为,,即：
、
F?
sF?
T v d?
F N p
Tss
A dAF N p
那么（ 3-30）中的荷载列阵就等于体积力、表面力和单元结点荷载叠加。因此作用在弹性体上的外力，需要移置到相应的结点上成为结点荷载。荷载移置要满足静力等效原则。
有限单元法
3.2.6 整体刚度矩阵得到了单元刚度矩阵后，需要将一系列的将单元组成一个整体结构，然后根据结点载荷平衡的原则进行分析，
得到整体刚度矩阵。整体分析包括以下 4个步骤：
（ 1）建立整体刚度矩阵，
（ 2）根据支承条件修改整体刚度矩阵，
（ 3）解方程组，求出结点的位移，
（ 4）根据结点位移，求出单元的应变和应力。
由单元刚度矩阵得到整体刚度矩阵的基本方法是刚度集成法，即整体刚度矩阵是单元刚度矩阵的集成。
有限单元法
x
y
1
2
3
4 5 6
i
j
m
i
j
m
i
j
m
i
j
m
单元编号 单元结点局部编号 单元结点整体编号
（ 1）
3
1
2
（ 2）
5
2
4
（ 3）
5
3
2
（ 4）
3
5
6
i
j
m
有限单元法
2j（ ）
(2)m
(2)i
2j（ ）
(2)jjk (2)jmk (2)jik
(2)m
(2)mjk (2)
mmk
(2)mik
(2)i
(2)ijk (2)imk (2)iik
整体编号 1 2 3 4 5 6
局部编号 0 0 0
1 0 0 0 0 0 0
2 0 0 0
3 0 0 0 0 0 0
4 0 0 0
5 0 0 0
6 0 0 0 0 0 0
有限单元法有限单元法由上表可以看出，整体刚度矩阵也具有对称性。同时它是一个稀疏矩阵，即其中有大量的零元素，并且非零元素都集中于主对角线附近呈带状。和单元刚度矩阵一样，由于位移函数中包含刚体位移，所以整体刚度矩阵也是一个奇异矩阵。必须要排除刚体位移后，才能变为正定矩阵。
有限单元法
3.2.7 约束条件的处理
（ 1）将方程分组
（ 2）对角线元素改 1，同行同列其他元素改 0
（ 3）对角线元素乘以一个大数有限单元法
（ 2）对角线元素改 1，同行同列其他元素改 0
若第 r个位移分量 为已知值 则将整体刚度矩阵中主对角线元素 改为 1,第 r行、第 r列的其他元素改为 0,
荷载矩阵中第 r行的元素改为,其他元素都减去结点位移的已知值和原 Ｋ 中这行相应列元素的乘积。
ru *ru
rrK
*ru
*
1 1 1 2 1 3 1 1 11
*
2 1 2 2 2 3 2 2 22
*
3 1 3 2 3 3 3 3 33
*
*
1 2 3
0
0
0
0 0 0 1 0
0
n rr
n rr
n rr
r r
n n n n n n n n r r
K K K K u F K u
K K K K u F K u
K K K K u F K u
u u
K K K K u u K u
有限单元法
（ 3）对角线元素乘以一个大数若第 r个位移分量 为已知值 则将整体刚度矩阵中主对角线元素 乘以一个大数（如 ）,第 r行、第 r列的其他元素不变,荷载矩阵中第 r行的元素改为,
其他元素不变。
ru *ru
rrK
* 2010rr rKu?
1 1 1 2 1 3 1 1 1 1
2 1 2 2 2 3 2 2 2 2
3 1 3 2 3 3 3 3 3 3
2 0 * 2 0
1 2 3
1 2 3
1 0 1 0
rn
rn
rn
r r r r r r n r r r r
n n n n r n n n n
K K K K K u F
K K K K K u F
K K K K K u F
K K K K K u K u
K K K K K u u
2010
有限单元法
3.3 平面矩形单元 ( x m,y m )
u
i
v
i ( x
i
,y
i
)
( x
j
,y
j
)
i j
m u
m
v
j
v
m
u
j
( x
k
,y
k
)
k
u
k
v
k
a a
b
b
x,?
y,?
矩形单元也是常用的单元之一，由于采用了比常应变三角形单元更高次数的位移模式，故可以更好地反映弹性体的位移状态和应力状态。
有限单元法如图 3.10所示四结点矩形单元，记单元的结点位移向量和结点力向量 为：
δ?
F?
T
i i j j m m k ku v u v u v u vδ?
T
x i y i x j y j x m y m x k y kF F F F F F F F?F?
为了能推导出简洁的结果，在这里引入无量纲坐标：
a
x
b
y
( x
m
,y
m
)
u
i
v
i ( x
i
,y
i
)
( x
j
,y
j
)
i j
m u m
v
j
v
m
u
j
( x
k
,y
k
)
k
u
k
v
k
a a
b
b
x,?
y,?
有限单元法
3.3.1 单元位移场从图 3.10可以看出，结点条件共有 8个，因此,x和 y方向的位移场可以各有 4个待定系数，可以取以下多项式作为单元的位移场模式：
1 2 3 4
5 6 7 8v
u a a x a y a x y
a a x a y a x y
它们是具有完全一次项的非完全二次项，其中以上两式中右端的第四项 是考虑到 x和 y方向的对称性而取的。
48,a xy a xy
(3-43)
有限单元法由结点条件，在 处，有,(,,,)rrx x y y r i j m k
,(,,,)
,
r r r
r r r
u x y u r i j k m
v x y v
将 (3-44)代入 (3-43)中，可以求解出待定系数,然后再代回 (3-43)中经整理后，有
18aa
(3-44)
0 0 0 0
0 0 0 0
i j m k
i j m k
N N N Nu
N N N Nv
δ N δ
(3-45)
其中,N为单元的形函数矩阵，
有限单元法
1
11
4
1
11
4
1
11
4
1
11
4
i
j
m
k
xy
N
ab
xy
N
ab
xy
N
ab
xy
N
ab
如以无量纲坐标系来表达，则 (3-46)式可以写成
1 1 1 (,,,)4r r rN r i j m k
其中：
,(,,,)rrrrxy r i j m kaa
(3-46)
(3-47)
有限单元法
3.3.2 单元应变场根据单元的位移场函数式（ 3-45）,(3-47)，由几何方程可以得到单元的应变场表达式，
1
uu
b
x
vv
a
y a b
uvuv
ab
yx
ε
记为：
ε Bδ?
(3-48)
(3-49)
有限单元法这里,B矩阵称为几何矩阵。 B矩阵可以表示为分块矩阵的形式
[]i j m kB B B B?B
其中
0
0
00
0
( 1 ) 0
11
0 0 ( 1 ) (,,,)
4
( 1 ) ( 1 )
r
r
r
rr
ri
rr
N
b
b
N
a a r i j k m
ab ab
ab
NN
ab
B
(3-51)
(3-50)
有限单元法
3.3.3 单元应力场由物理方程及式（ 3-49），可以得到单元的应力场表达式，
ζ D ε DB δ S δ
其中 为应力矩阵,D称为弹性矩阵，对于平面应力问题，
S DB
2
10
10
( 1 )
1
00
2
E
D
(3-52)
有限单元法将应力矩阵表示为分块矩阵的形式
[]i j m k?S S S S S
其中：
00
002
00
( 1 ) ( 1 )
( 1 ) ( 1 ) (,,,)
4 ( 1 )
11
( 1 ) ( 1 )
22
rr
r r r r
rr
ba
E
b a r i j k m
ab
ab
S D B
对于平面应变问题，只需将 E换为,换为 。
21
E
1
（ 3-54）
（ 3-53）
有限单元法
3.3.4 单元刚度矩阵和三角形单元一样，可以根据最小势能原理导出结点位移向量和结点力向量之间关系，即单元的刚度矩阵,
可以将其写成分块的形式。
K?
i i i j i m i k
j i j j j m j k
mi mj mm mk
k i k j k m k k
k k k k
k k k k
k
k k k k
k k k k
其中
(,,,,)Tr s r s d r s i j m kk B D B （ 3-56）
（ 3-55）
有限单元法对于平面应力问题，如果单元厚度 t为常数，则得到
（ 3-56）的显式形式：
2
1
1
3 1
211
1
23
,,,,
4( 1 ) 1
1
31
2 11
1
23
r s r s
r s r s
r s r s
rs
r s r s
r s r s
r s r s
b
a
a
bEt
r s i j m k
ba
ab
k
（ 3-57）
有限单元法
3.4 平面问题程序设计有限元方法在分析工程问题时的优越性,在于与计算机的结合 。 通过编写的计算程序,确定所研究问题的信息和数据,问题就可以得到解决 。
前面两节讲述了平面问题的有限元计算原理。本节主要介绍采用三角形单元，计算在给定荷载作用下的弹性力学平面静力问题时，有限元法的程序编制与框图设计。通常有限元法的程序设计和编写分为以下几步：
有限单元法
（ 1） 问题分析使用计算机解决具体问题时，首先要根据所要求解的问题，找出已知的数据和条件，建立数学力学模型，选择计算方案，准备所需要的原始数据。
有限单元法
（ 2） 框图设计计算方案选定以后，
需先画出总框图。平面问题的有限元法的总框图如图 3.11所示。
开 始输 入 原 始 数 据组 成 整 体 刚 度 矩 阵组 成 整 体 荷 载 向 量引 入 边 界 条 件解 刚 度 方 程,输 出 位 移求 应 力,输 出 应 力求 应 变,输 出 应 变求 支 反 力,输 出 支 反 力结 束有限单元法
（ 3） 设计程序根据总框图的规划和有限元法的基本原理,用自己所使用的程序语言把这个解决方案严格地描述出来,也就是编写出程序代码 。 它就是用于求解问题的源程序 。
（ 4） 程序的编译和调试将编写好的源程序输入计算机，进行编译，改正其中的语法等错误。并用实际的输入数据对编好的程序进行调试，
分析所得到的运行结果，进行程序的测试和调整，直至获得预期的结果。
有限单元法第 4章空间轴对称问题的有限元分析有限单元法教学提示,空间轴对称结构是工程中一类特殊的结构 。
对于这类结构,由于其几何形状,约束情况,所受荷载都对称于某一根轴 （ 如图 4.1所示 ）,因此,通常是在柱坐标系中进行分析 。 本章以三角形单元为例介绍了这一类结构的有限元分析过程 。
教学要求,本章要求学生了解 柱坐标系中三角形单元的建立方法，熟悉轴对称结构有限元分析的基本方法。
有限单元法
4.1 概述在工程中有许多结构，如活塞、厚壁容器等，他们的几何形状、约束情况及所受的荷载都对称于空间的某一根轴，
因此在物体中通过该轴的任何平面都是对称面，所有应力、
应变和位移也对称于该轴，这类问题称为空间轴对称问题。
研究轴对称问题时通常采用圆柱坐标系,以 轴为对称轴。
(,,)rz? z
有限单元法
x
ya
b
z
p
如图 4.1所示的受均布内压作用的长圆筒，通过 z轴的一个纵截面就是对称面。由于对称性，轴对问题共有 4个应力分量：
Tr z zrζ
有限单元法其中 表示沿半径方向的正应力，称为径向应力； 表示沿方向的正应力，称为环向应力或切向应力； 表示沿 ｚ 方向的正应力，称为轴向应力； 表示在圆柱面上沿ｚ方向作用的剪应力。
r
z?
zr?
同样，轴对称问题共有 4个应变分量：
Tr z zrε
其中 表示沿半径方向的正应变，称为径向正应变； 表示沿 方向的正应变，称为环向正应变或切向正应变； 表示沿ｚ方向的正应变，称为轴向正应变； 表示沿ｒ和ｚ方向的剪应变。
r
z?
zr?
有限单元法在轴对称问题中，弹性体内任意一点上，不存在切向位移，只存在径向位移 u 和轴向位移 w,两个位移分量表示为，
Tuw?d
有限单元法
4.2 三角形单元轴对称问题分析中所使用的三结点单元，在对称面上是三角形，在整个弹性体中是三棱圆环，各单元中圆环形铰相联接。
有限单元法
4.2.1 位移函数的选取
( r
m
,z
m
)
u
i
w
i
i ( r
i
,z
i
)
j ( r
j
,z
j
)
m
r
z
O
u
j
u
m
w
j
w
m
参照平面问题的三角形单元位移函数，轴对称问题的三结点三角形单元位移函数取为，
1 2 3
4 5 6
u a a r a z
w a a r a z
有限单元法按照平面问题三角形单元的分析过程，将结点坐标和结点位移代入（ 4-4）得到，
1
2
3
1
2
i j m i
i j m j
i j m m
a a a a u
a b b b u
A
a c c c u
4
5
6
1
2
i j m i
i j m j
i j m m
a a a a w
a b b b w
A
a c c c w
其中：
mm
jj
ii
zr
zr
zr
A
1
1
1
2
1
有限单元法
i j m m ja r z z r i j mb z z i m jc r r
j m i i ma r z z r j m ib z z j i mc r r
m i j j ia r z z r
m i jb z z m j ic r r
定义形函数为：
1 ( ) (,,)
2s s s sN a b r c z s i j mA
用矩阵表示的单元位移为：
0 0 0
0 0 0
i
i
i j m j
i j m j
m
m
u
w
N N N uu
N N N ww
u
w
dN δ
有限单元法
4.2.2 单元刚度矩阵轴对称问题的几何方程：
r
z
zr
u
r
u
r
w
z
uw
zr
由（ 4-9）式得，
)(2 1 mmjjii ubububAru
)(2 1 mmjjii ufufufAru
有限单元法其中，
(,,)s s sssa c zf b s i j mrr
)(2 1 mmjjii wcwcwcAzw
)(2 1 mmjjii ucucucAzu
)(2 1 mmjjii wbwbwbArw
用几何矩阵表示单元的应变，
ε Bδ?
[]i j m?B B B B
0
01
(,,)
02
s
s
s
s
ss
b
f
s i j m
cA
cb
B
有限单元法由于 在是坐标 r,z 的函数,分量在单元中不为常量，其它三个应变分量在单元中仍为常量。
sf
由轴对称问题的物理方程，得到弹性矩阵 D，
10
11
10
11( 1 )
( 1 ) ( 1 2 )
10
11
12
000
2( 1 )
E
D
有限单元法由弹性矩阵 D 和几何矩阵 B 可以得到应力矩阵 S，并计算出单元内的应力分量，
[]i j mS D B S S S
其中：
1
11
1
22
( 1 )
(,,)
()2( 1 ) ( 1 2 )
s s s
s s s
ss
s s s
ss
b f A c
A b f A cE
s i j m
A b f cA
A c A b
S D B
式中：
1 1A
2
12
2 (1 )A
有限单元法有了单元应力场和应变场，可以利用虚位移原理或最小势能原理建立单元刚度矩阵 k?
单元刚度矩阵的分块矩阵为，
2 T rdrdzK B D B?
2 (,,,)Tsq s q rd rd z s q i j mK B D B
由于几何矩阵中的元素不是常量，单元刚度矩阵需要通过积分得到，为简化计算，可以用三角形单元形心位置的坐标 代替 B 矩阵中的变量 。
cc zr,
B
,rz
)(31 mjic rrrr )(31 mjic zzzz
有限单元法应变矩阵（ 4-13）变成，
* * * *i j mB B B B
其中,0
01
(,,)
2
0
s
s
s s c
cs
s
ss
b
a
b c z
r s i j m
A
c
cb
*
B
有限单元法单元刚度矩阵的近似表达式为：
2 Tcr* * *K B D B?
单元刚度矩阵的分块矩阵近似表达式为，
1 2 1 2
1 2 2
2
( ) ( )
()
T
s q c s q
s q s q s q s q s q s q s q s q
s q s q s q s q s q
r
b b f f A b f f b A c c A b c f c A c b
A c b c f A b c c c A b b
**
* * * * *
*
K B D B
这里：
* (,,)s s s
ss
a c zf b s i j m
rr
有限单元法
4.2.3 等效结点荷载
（ 1） 集中力集中力的处理很简单，一般直接把集中力作用点取为结点，不需要作特殊处理，就可以直接把集中力加入到结点荷载列阵中去。
（ 2）体积力设单元内单位体积上作用的体积力为,则移置到单元各结点的等效结点力为
2 T v rd rd zF N p?
vp
有限单元法
（ 3）表面力设单元某边上作用的表面力为,则移置到单元各结点的等效结点力为
2 Tsi ss rd sF N p?
sp
轴对称问题分析中，如果直接定义结点荷载，则荷载值是实际弹性体上绕对称轴一周的荷载的累计结果。
有限单元法第 6章平板弯曲问题的有限元分析有限单元法教学提示,对于工程中经常遇到的不规则薄板,经典理论通常显得无能为力,此时必须运用数值方法进行分析 。 本章就 平板弯曲问题的有限元分析进行了介绍,分别运用三角形单元,矩形单元,八结点四边形等参单元等单元划分形式对平板弯曲问题进行介绍,包括位移模式,单元分析,整体分析,等效结点荷载计算等方面 。
教学要求,本章要求学生重点掌握运用三角形单元进行薄板的有限元分析,包括位移模式、单元分析、整体分析、
等效结点荷载计算等。同时要熟悉矩形单元的运用，了解八结点四边形等参单元的分析过程。
有限单元法
6.1 薄板受弯分析的基本方程当薄板上受有一般载荷时,总可以把个荷载分解为两个分量：一个是作用在薄板中面内的中面荷载 （ 也称纵向荷载 ） ；
另一个是垂直于中面的法向荷载 （ 也称横向荷载 ） 。 对于中面荷载,可以认为它们沿薄板的厚度均匀分布,因而它们所引起的位移,应变和应力,可以按平面应力问题进行分析 。 横向荷载将使薄板发生弯曲和翘曲,它们所引起的位移,应变和应力,
应按薄板弯曲问题进行计算 。
当薄板弯曲时,中面所弯成的曲面称为弹性曲面,而中面内各点在垂直于中面方向的位移称为挠度 。
本节只介绍薄板弯曲的小挠度理论,即薄板虽然很薄,但仍具有相当的弯曲刚度,因而它的挠度远小于其厚度 。
有限单元法
6.1.1 基本假设分析薄板弯曲的挠度问题时,和材料力学中分析直梁的弯曲问题时相似 （ 薄板的中面相当于直梁的轴线,薄板的弹性曲面相当于直梁的挠曲线 ）,也采用一些由实践经验得到的基本假设,使问题大大简化,但同时又能在一定程度上反映实际情况 。 这些基本假设是：
（ 1） 薄板的法线变形后没有伸缩 。
（ 2） 变形前的中面法线在变形后仍是弹性曲面的法线 。
（ 3） 薄板中面内各点,没有平等于中面的位称 。
（ 4）忽略挤压应力所引起的变形。
有限单元法
6.1.2 几何方程
t / 2
t / 2
O
x
y
z
取薄板的中面为 xy面,z轴垂直于中面，如图 6.1所示。
由假设（ 1）可知，,再由,可得 。
也就是说，中面法线上的所有各点具有相同的位移,即弹性曲面的挠度。
0x 0z wz (,)w w x y?
w
有限单元法根据假设（ 2），可以推知：薄板的法线（ z方向线段）
与 x方向或 y方向的线段保持垂直，即没有剪应变，也就是：
0,0y z z xvuz y x z
上式也可写成：
,v w u wz y z x
对式（ 6-1）进行积分，注意到 只是 x和 y的函数，不随 z而变，因而得：
w
1
2
(,)
(,)
u z f x y
x
v z f x y
y
（ 6-2）
（ 6-1）
有限单元法假设（ 3）可以表示为,,代入到式（ 6-2）得：
0 0zu 0 0zv
1 (,) 0f x y?
2 (,) 0f x y?
于是式（ 6-2）就简化为：
,wwu z v zxy （ 6-3）
有限单元法现用挠度来表示应变，不难得到：
2
2
2
2
2
2
x
y
xy
wu
xx
vw
z
yy
uv w
yx xy
ε
这就是弯曲薄板的应变与挠度之间的几何方程。
（ 6-4）
有限单元法在小变形情况下,和 分别为弹性曲面在 x和
y方向的曲率 和,而 为弹性曲面在 x和 y方向的扭率 。 这三个参数称作弹性曲面的弯扭变形分量,它们完全确定了薄板内各点的应变分量 。 用矩阵可表示为：
（ 6-5）
将上式代入到式 （ 6-4）,得到：
（ 6-6）
从上式中可以看出，薄板内所有各点的应变分量都可由弹性曲面的弯扭变形求出。因此，有时也把式（ 6-5）称作薄板弯曲问题的几何方程。
22/wx 22/wy
x?y? 22/w x y
xy?
2 2 2
22 2
TT
x y x y
w w w
x y x y
χ
z?ε χ
有限单元法
6.1.3 物理方程假设（ 4）说明：可以忽略挤压应力 引起的变形。因此薄板内各点的应变分量可用应力分量来表示，即：
z?
xy
x E
yx
y E
2 (1 ) xy
xy E
这和薄板平面应力问题中的物理方程相同，由式（ 6-7）解出应力，可得：
（ 6-7）
有限单元法
2 ()1x x y
E
2 ()1y y x
E
2 (1 )x y x y
E
将式（ 6-4）代入上式，得到用挠度表示的应力分量，用矩阵来表示可写成：
2
2
2
2
2
2
x
y
xy
w
x
w
zz
y
w
xy
ζ DD χ
（ 6-8）（ 6-8）
有限单元法式中：
2
10
10
1
1
00
2
E
D （ 6-9）
在薄板的弯曲问题中，由于大多数情况下，都很难使得应力分量在薄板的侧面上（板边上）精确地满足应力边界条件，
而只能使这些应力分量所组成的内力整体地满足边界条件，因此，有必要考察薄板横截面上的内力。
有限单元法
x
t
/
2
t
/
2
yz
O
1
1
y x
x y
y
x
t
/
2
t
/
2
x
yz
O
1
1
从薄板内取出一个平行六面体，它在 x和 y方向上具有单位宽度，在 z方向的高度为 t（图 6.3）。在 x为常量的横截面上，
作用有 和 。由于 和 都和 z成正比，所以它们在薄板全厚度上的代数和分别等于零。只可能分别合成弯矩和扭矩。用表示由 所合成的单位宽度上的弯矩，则得：
x? xy?x?xy?
xM x?
有限单元法
/2
/2
t
xxtM z d z
将式（ 6-8）中的 的表达式代入上式，并对 z积分，得到
x?
3 2 2
2 2 2()1 2 ( 1 )x
E t w wM
xy
xy?与此相似，应力分量 将合成扭矩
/2
/2
t
x y x ytM z d z
将式（ 6-8）中的 表达式代入上式，并对 z积分，得到
xy?
32
1 2 (1 )xy
E t wM
xy?
（ 6-11）
（ 6-10）
有限单元法同样，在 y为常量的横截面上，每单位宽度上的 和 也分别合成如下的弯矩和扭矩
y? yx?
3 2 2/2
2 2 2/2
32/2
/2
()
1 2 ( 1 )
1 2 ( 1 )
t
yy t
t
y x y xt
E t w w
M z d z
yx
E t w
M z d z
xy
在这里可以看到，由剪应力 和 的互等关系，得到扭矩 和 的互等关系。
将式（ 6-10）、（ 6-11）和（ 6-12）中的第一式合并起来，用矩阵表示，则有
yx?xy?
xyM yxM
（ 6-12）
有限单元法
22
22
2
/
/
2/
x
y f f
xy
M w x
M w y
M w x y
M D D χ
其中
3
2
10
10
1 2 ( 1 )
1
00
2
f
Et
D
是薄板弯曲问题的弹性矩阵，它等于平面应力问题中的弹性矩阵乘以 。式（ 6-13）表示了薄板的内力与应变两者之间的关系，因而是薄板弯曲问题中的物理方程。
3 /12t
（ 6-13）
（ 6-14）
有限单元法
t
/
2
t
/
2
z
d
z
x
y
z
x
M
x
y
M
x y
M
y x
M
y
O
x y
y x
1
1
x
y
z
M
x
M
x
M
y
M
y
M
x y
M
x y
M
y x
M
y x
O
弯矩和扭矩,,,的方向及其作用面的位置示于图 6.4a中 。 按右手螺旋法则用双箭头矢量来表示力偶,如图 6.4b所示 （ 图中所示各力偶的方向均为正 ） 。
从式（ 6-13）中解出,,,
再代入式（ 6-8），就得到各应力分量与内力之间的关系式。
xM yM xyM yxM
22/wx 22/wy 22/w x y
图 6.4（ a）
（ b）
有限单元法
3
12TT
x y x y x y x y
z M M M
tζ
从式（ 6-13）中解出,,,
再代入式（ 6-8），就得到各应力分量与内力之间的关系式。
22/wx 22/wy 22/w x y
（ 6-15）
有限单元法
6.2 三角形单元
y
z
x
1
2
3
w
O
x 1
y 1
有限单元法
6.2.1 位移模式问题为保证挠度 w 为坐标的全三次多项式，从帕斯卡三角形可知必须要有 10项，但三角形 3个结点只能有 9个自由度，若舍去三次项中的任一项，显然都无法保证对坐标的不变性，为此 Tocher提出了一种解决方案如下：
2 2 3 2 2 31 2 3 4 5 6 7 8 9()w x y x x y y x x y x y y
但是当三角形的二边分别平行坐标 x轴,y轴时，从式（ 6-16）
无法通过结点的位移条件来确定广义坐标参数,为此必须在离散化时设法避免边界（内部）单元不出现上述现象，分析结果表明此位移模式对应的单元能得到很好的结果。
（ 6-16）
有限单元法另一解决方案是 Adini提出的
2 2 3 2 2 31 2 3 4 5 6 7 8 9w x y x y x x y x y y
此方案由于舍去了二次项 xy，故常扭率 无法得到保证，从而位移偏小，精度很差。这也可以从挠曲函数的
Taylor级数展开来说明，由于不包含完全二次项，故只能有一阶精确度。
2w
xy
有限单元法还有一种 Bell提出的解决方案，他对三角形单元除了
3个角点作为结点外，取形心点挠度作为位移参数，从而取位移模式为全三次多项式利用 10个位移参数条件确定广义坐标参数，通过分析建立起单元刚度、等效结点荷载矩阵，在整体分析之前采取如下措施消去内点自由度。
2 2 3 2 2 31 2 3 4 5 6 7 8 9 1 0w x y x x y y x x y x y y
有限单元法单元刚度方程
（ 6-19）
式中从式 （ 6-19） 的第二个方程可解出：
（ a）
将式 （ a） 代回式 （ 6-19） 第一个方程并进行整理,可得：
,111 1 1 2 1
,1 0102 1 2 2
eq
eq
Fδkk S
Fδkk 0
1 1 1 1 2 2 2 3 3 3 Tx y x y x yw w w wδ
110 22,10 21 1eqδ k F k δ
111 1 1 2 2 2 2 1 1 1,1 1 2 2 2,1 0() e q e qk k k k δ S F k k F
有限单元法设单元刚度矩阵：
等效结点荷载,（ 6-20）
则单元刚度方程为 （ 6-21）
通过静力凝聚最后获得了 9自由度的三角形单元。但是
Zienkiowicz曾指出这样所得到的单元不能保证收敛性，因此
，通常采用另外一种位移模式，即下面介绍的 面积坐标 下的位移模式。
11 1 1 2 2 2 2 1k k k k k?
1,1 1 2 2 2,1 0e q e qF F k k F?
11k δ SF
有限单元法
6.2.2 面积坐标下的位移模式
P
A
2
3 ( x
3
,y
3
)
x
y
O
A
1
A
3
L
3
= 0
L
2
= 0
L
1
= 0
1 ( x
1
,y
1
)
2 ( x
2
,y
2
)
设面积坐标变量为 。
其定义如下：如图 6.6所示，三角形单元中任意一点 P的位置可以由如下三个比值来确定：
1 2 3L L L、,
11 AL A?
22 AL A?
33 AL A?
A为三角形单元的面积，
为小三角形面积
1 2 3A A A、,
有限单元法根据三角形面积公式,可以得到面积坐标与直角坐标的关系为：
式中
1 1 1
1
2 2 2
2
3 3 3
3
2
2
2
a b x c y
L
A
a b x c y
L
A
a b x c y
L
A
1 2 3 3 2a x y x y 1 2 3b y y 1 3 2c x x
2 3 1 1 3a x y x y 2 3 1b y y 2 1 3c x x
3 1 2 2 1a x y x y 3 1 2b y y 3 2 1c x x
有限单元法
Zienkiowicz等采用面积坐标解决了直角坐标下所产生的困难，提出了如下的位移模式 22
1 1 2 2 3 3 4 1 2 1 2 3 5 1 2 1 2 3
222
6 2 3 1 2 3 7 3 2 1 2 3 8 3 1 1 2 3
2
9 1 3 1 2 3
11
( ) ( )
22
1 1 1
( ) ( ) ( )
2 2 2
1
()
2
w L L L L L L L L L L L L L
L L L L L L L L L L L L L L L
L L L L L
利用结点的位移参数条件可以确定广义坐标参数 α,从而建立形函数（当然也可象前面介绍的那样由形函数性质来确定）如下
（ a）
有限单元法
2 2 2 2
1 1 1 2 1 3 1 2 1 3
22
1 3 1 2 1 2 3 2 3 1 1 2 3
22
1 3 1 2 1 2 3 2 3 1 1 2 3
11
( ) ( )
22
11
( ) ( )
22
x
y
N L L L L L L L L L
N b L L L L L b L L L L L
N c L L L L L c L L L L L
式中脚标轮换规则为 1→ 2→ 3→ 1,如果注意到：
1 2 3 1 2 1 2 1 3 1 3 2 3 2 3( 1 ) ( 1 ) ( 1 )L L L L L L L L L L L L L L L
则式（ a）可改写成：
22
1 1 2 2 3 3 4 2 3 5 3 1 6 1 2 7 2 3 3 2
2 2 2 2
8 3 1 1 3 9 1 2 2 1
()
( ) ( )
w L L L L L L L L L L L L L
L L L L L L L L
有限单元法从式 （ c） 出发,可按如下两步法确定形函数：
（ 1） 以,,作为结点自由度 （ 位移参数 ）,求对应它们的形函数；
（ 2）利用关系式
w 1 1wL 2 2wL
11 1 2 1 2
2 1 1 1 1
2
y
x
w w
L c b c bx
wc b c bw
yL
将（ 1）中的,变换成,,然后进行合并整理即可得到对于（,,）结点位移参数的形函数。
1? 2? x? y?
w x?
y?
有限单元法下面具体推导如下：
（ 1） 由,,的位移条件可得
、,(e)
（ 2） 由式 （ c） 求导可得
1 11Lww 2 21Lww 3 31L
ww
11w 22w 33w
2
1 1 3 4 2 5 3 1 6 2 7 2 3 2
1
2 2 2
8 3 1 1 3 9 2 2 1
( ) ( 2 )
( 4 ) ( 2 )
w w w L L L L L L L
L
L L L L L L L
22
2 2 3 4 3 2 5 1 6 1 7 2 3 2 3
2
22
8 1 3 1 9 2 1 1
( ) ( 4 )
( 2 ) ( 2 )
w w w L L L L L L L L
L
L L L L L L
(f)
有限单元法
（ 3）建立如下位移条件
11 11 1 3 5 81L ww
21 1 2 1 3 4 6 7 91L ww
31 13 1 3 5 81L ww
12 2 1 2 3 5 6 8 91L ww
22 2 2 2 3 4 71L ww
32 2 3 2 3 4 71L ww
（ g）
有限单元法并由此求得式：
4 23 22
1 ()
2
5 1 3 1 1
1 ()
2
6 12 11 21 22
1 ()
2
7 3 2 22 23
1 ()
2ww
8 1 3 11 13
1 ()
2ww
9 2 1 1 1 1 2 2 1 2 2
1 ()
2ww
（ h）
有限单元法
（ 4）将式（ e）和式（ h）代回式（ b）并整理，可得：
1 1 1 1 1 1 2 1 2 1 2 2 1 2 1 1 2 2 2 2 3 3 1 3 2 2 2 3 2 1w N w N N N w N N N w N N
式中：
…………
2 2 2 21 1 1 2 1 3 1 2 1 3N L L L L L L L L L
2 2 2 21 1 1 2 3 1 1 2 2 1 3 1 1 31 1 1 1( ) ( )2 2 2 2N L L L L L L L L L L L L
222 1 1 2 1 2 2 111 ()22N L L L L L L
（ j）
（ i）
有限单元法
（ 5）将式（ d）代入式（ i）并再整理即可得：
1 2 3wN δ N N N δ
其中：
( 1,2,3 )i i x i y iN N N iN
1 2 1 1 1 2 1xN b N b N
1 2 1 1 1 2 1yN c N c N （第二脚标轮换） 1→2→3→1
不难验证所得形函数与式（ 6-22）完全相同。但这种推求方法显然比普通广义坐标法直接求 而后得形函数要方便得多。
19
（ 6-23）
有限单元法
6.2.3 单元分析无论是以式（ 6-16）还是用式（ 6-22）、（ 6-23）建立单元的挠度场，有了 w就可用常规方法来进行单元列式。下面以面积坐标的形函数和式（ 6-23）来推导并给出三角形单元的显式单元刚度矩阵等等。
有限单元法注意到面积坐标,设取,为独立坐标，则：1 2 3 1L L L 1L 2L
112
12
2
1
2
Lbbx
ccA
y L
由此不难得到：
2
2
2
2
1
22
2 2 2
2
2 2
12
1
4
2
Lx
y A L
xy LL
T
式中 A是三角形单元中面面积。
（ 6-24）
（ 6-25）
有限单元法
22
1 2 1 2
22
1 2 1 2
1 1 2 2 1 2 2 1
2
2
2 2 2 ( )
b b b b
c c c c
b c b c b c b c
T
因此，可得形变矩阵为：
设：
1
2
2
2
2
22
1 2 3 1 2 32 2 2
2
2 2
12
1
4
2
L
x
w
y A L
xy LL
T N N N δ B B B δ B δ
（ 6-27）
（ 6-26）
有限单元法式中：
1
2
2
2
22
2
2
12
1
( 1,2,3 )
4
ii
L
i
AL
LL
B T N
因为 是面积坐标三次式，因此 是面积坐标一次式。
若记
12(,)i LLN iB
1 2 3
00
00
00
L L L
φ
Φ
（ a）
（ 6-28）
有限单元法则可得
2
1
4iiABTΦ G （ 6-29）
式中
1 2 3 TT T Ti i i iG G G G
22
1 3 2 1 3 2 1
22
6( ) 2 ( ) 2 ( )
62
6( ) 2 ( 2 ) 2 ( 2 )
i i i i i i
i i i i i i
i i i i i i
b z c z
b b b c c c
bc
G
2 1 3 2 1 3
2 1 1
11
62
6( ) 2 ( 2 ) 2 ( 2 )
6( ) ( 2 ) ( 2 )
i i i i i
i i i i i i i
i i i i i i
b b b c c c
bc
bc
G （ d）
（ c）
（ b）
有限单元法
11
22
11
3
22
11
2
11
( 2 3 ) ( 2 3 )
22
26
11
( 3 2 ) ( 3 2 )
22
11
( 2 3 ) ( 2 3 )
22
26
11
( 2 3 ) ( 3 2 )
22
11
( 2 3 ) ( 2 3 )
22
62
1
( 3 2 )
2
i i i i i i
i
i i i i i i
i i i i i i
ii
i i i i i i
i i i i i i
i
i i i
bc
bc
bc
bc
bc
b
G
2
1
( 3 2 )
2
i i i
c
其中,,为 i结点的 3个面积坐标，也即 i结点的。i?
i? i?
（ e）
有限单元法若将单元刚度矩阵写成分块子矩阵形式，则对各向刚性薄板若将单元刚度矩阵写成分块子矩阵形式，则对各向刚性薄板
4
4
10
10
1
00
2
10
10
16
1
00
2
16
T
ij i j
A
T T T
ij
A
TT
ij
A
D dA
D
dA
A
D
dA
A
k B B
G Φ TT Φ G
G Φ H Φ G
（ 6-30）
有限单元法式中：
3
21 2 (1 )
EtD
33ijHH
其中：
2221 1 1 1H b c
221 2 2 1 1 2 1 2 1 2 2 1H H b b c c b c b c
2213 31 1 1 1 2 1 22H H b c b b c c
222 3 3 2 2 2 1 2 1 22H H b c b b c c
2222 2 2 2H b c
222 2 2 23 3 1 1 2 2 1 2 1 2 1 2 2 12 2 2H b c b c b b c c b c b c
有限单元法因为
11 12 13
21 22 23
31 32 33
T
H H H
H H H
H H H
φ φ φ
Φ H Φ φ φ φ
φ φ φ
式中
2
1 1 2 1 3
2
2 1 2 2 3
2
3 1 3 2 3
T
L L L L L
L L L L L
L L L L L
φ φ φ
有限单元法利用
1 2 3
! ! ! 2
( 2 ) !A L L L d A A
并记：
2 1 1
1 2 1
1 1 2
P
则单元刚度子矩阵 为ijk
3 3 3
1 1 2 2 3 33
1 1 1192
T T T
i j i r j r i r j r i r j r
r r r
D
k G P H G G P H G G P H G
（ 6-33）
（ 6-31）
（ 6-32）
有限单元法若单元上受垂直中面的分布荷载 作用，则等效结点荷载为
(,)q x y
33
11
(,) (,)TTe q i i i i
ii
q x y d A q L x L y d A
F N N
当 q为常数时，有式
2 3 2 3 3 1 3 1
1 2 1 2
1 1 1 1 1 1
( ) ( ) ( ) ( )
3 2 4 2 4 3 2 4 2 4
1 1 1
( ) ( )
3 2 4 2 4
eq
T
q b b c c b b c c
b b c c
F
（ 6-35）
（ 6-34）
有限单元法经推导，若记（各向同性时）
2 2 2 2
1 1 2 2 1 2 1 23
2 2 2 2
1 1 2 2 1 2 1 2 3 322
1 1 2 2 1 2 2 1
2 ( )
2 ( ) [ ]
4 8 ( 1 ) ( 1 ) ( 1 ) ( 1 ) ( ) ij
b c b c b b c c
Et b c b c b b c c R
b c b c b c b c
R
3
1
1
3
2
1
3
3
1
r i r
r
i r i r
r
r i r
r
R
R
R
φ G
S φ G
φ G
（ 6-37）
（ 6-36）
有限单元法式中,则薄板内力矩阵可写作：
（ 6-38）
由此式可见,的每一列的物理意义为,单位位移引起板中任一点的内力（点位置取决于,和 ）,。
1 2 3L L L?φ
33
11
[] Ti i i i x i y i
rr
w
MS δ S
iS
1L 2L 3L
有限单元法结束
