第三章 空间问题的有限元方法 引言
许多工程实际问题，属于空间问题，由于结构形状或受力的复杂性，用经典弹性理论去求解它们的解析解是不可能的。而有限元法处理此类问题，原则上不存在什么困难，本章将介绍一般空间问题的四面体单元。 一般空间问题的有限元列式
3.2.1 单元位移模式及插值函数
空间问题中，每个单元有四个结点，编码为i,j,m,p 。每个结点有3个位移分量。每个结点的位移可用位移矢量i α表示，即
??
??
?
?????=i i i i w v u α ),,,(p m j i
单元结点的位移向量可表示为
[
]
T
p p p m m m j j j i i i
p m j i e
w v u w v u w v u w v u =??????
????????=ααααα
？
e α为单元结点位移列阵。
假设单元内的位移模式选取一次多项式
z y x u 4321ββββ+++=
z y x v 8765ββββ+++= （3.2.1） z y x w 1211109ββββ+++=
由于四个结点也在单元内，满足位移模式，于是得
i i i i z y x u 4321ββββ+++=
j j j j z y x u 4321ββββ+++= （3.2.2）
m m m m z y x u 4321ββββ+++=
…
p p p p z y x u 4321ββββ+++=
上式是关于4321,,,ββββ的线性方程组。4321,,,ββββ是待定常数，也称为广义坐标。它可由（3.2.2）式求出。上式的系数行列式是
V z y x z y x z y x z y x D p
p
p
m m m j j j i i i 21111==
（3.2.3）
上式中当i,j,m,p 的编号顺序满足右手法则，V 值为正，其大小为四面体体积，因此为了方便单元的编号一般满足右手法则。求得4321,,,ββββ后,回代入位移模式得
p p m m j j i i u N u N u N u N u +++= (3.2.4)
式中
)(61
z d y c x b a V
N i i i i i +++= ),,,(p m j i (3.2.5)
p
p
p
m m m
j j j
i z y x z y x z y x a = -
p p m m
j j
i z y z y z y b 111-=
p
p
m m
j
j
i z x z x z x c 111= ),,,(p m j i （3.2.6）
p
p
m m
j
j
i y x y x y x d 111-= 上式下标),,,(p m j i 轮换，可得j j j j d c b a ,,,,m m m m d c b a ,,,及p p p p d c b a ,,,。 同理,也可得到其它两式,于是得
p p m m j j i i u N u N u N u N u +++=
p p m m j j i i v N v N v N v N v +++= （3.2.7）
p p m m j j i i w N w N w N w N w +++=
其中
、
)(61
z d y c x b a V
N i i i i i +++=
),,,(p m j i (3.2.8)
p m j i N N N N ,,,称为单元的插值函数或形函数，这里它是z y x ,,的一次函数，其中
i i i i d c b a ,,,，j j j j d c b a ,,,,m m m m d c b a ,,,及p p p p d c b a ,,,是常数，由表达式可知，
它完全由单元的大小和方位确定，一旦单元确定了，这些常数也完全确定。
（3.2.7）式的矩阵形式是
??????
??????????????????=??????????=p m j i p p p m
m m j
j j i
i
i
N N N N N N N N N N N N w v u u αααα0
00000
00000
00000
0000 [
]
???????
??????
?=p m j
i p m
j i
IN IN IN IN αααα
[]
e e p m j i N N N N N αα==][][][][ （3.2.9） N 称为插值函数矩阵或形函数矩阵。 3.2.2．应变矩阵和应力矩阵 ⑴应变
(
确定了单元位移后，可以很方便地利用几何方程和物理方程求得单元的应变
和应力。在（1.4.21）式的几何方程中，位移用（）式代入，得到单元应变为
????
?????????
??
????????????
???????+????+????+????????=??????????????????=z u x w y w z v x v y u z w y v x u zx yz xy z y x γγγ
εεεε
[
]
e e p m
j i
B B B B B αα=--=][][][][ （3.2.10）
B 称为应变矩阵。
应变矩阵的分块矩阵i B ][是
????
?
???
??????????=
i i i i i i i i i i b d c d b c d c b V B 00000
000061][ ),,,(p m j i （3.2.11）
可以看出，应变矩阵B 中的元素都是常量，从而单元中的应变都是常量，所以三维线性位移模式的四面体单元是常应变单元。 ⑵ 应力
单元应力可以根据物理方程求得， 其应力应变关系如下：
，
)]
([1
)]([1
)]([1
x y x z z x z y y z y x E
E E σσνσεσσνσεσσνσε+-=+-=+-=
??
?
??
===μτγμτγμτγ///xy xy zx zx yz yz
或
???
???
???+=???
??+-+=+=???
??+-+=+=???
??+-+=
xy xy z z zx zx y y yz yz x x E E E E E E γντεθνννσγντεθνννσγντεθννν
σ)1(2,211)1(2,211)1(2,211
于是应力向量可表示为
e e S DB D ααεσ=== （3.2.12）
式中D 为弹性矩阵，而
>
????
???
????????
????????
?
???????
?-------------+-=)1(2210
0)1(221000)1(2210
1
1111
1111)21)(1()1(νννννν
ν
ν
νν
νν
ννν
ν
νν
νννE D
（3.2.13）
从而可以到，三大物理参量，都可以用单元结点位移向量表示：
e N u α=
[]e T zx yz xy z y x B αγγγεεεε==
[]e e T zx yz xy z y x S DB αατττσσσσ===
由于N ，B ，S 都是已知的矩阵，只要求得e α，则单元内的位移、应变和应力就可以就得，问题是：如何求结点位移向量 3． 单元刚度矩阵和结点载荷向量
对于三维单元，单元刚度矩阵也具有上章所讨论的单元刚度矩阵的一般形式，即
DBV
B DBtdxdy B K T V T e e
==? （3.2.15）
}
写成分块矩阵的形式
?
????????????
?--------=e
pp e pm
e pj
e pi e mp
e mm e mj e mi e jp e jm
e jj
e ji
e ip e im e
ij e ii e
K K K K K K K K K K K K K K K K K （3.2.16） 每个子矩阵为
V B D B K s T r e
rs ][][=
等效结点载荷
?=e
V T e f fdV N P
?=e S T e S TdS N P σ
e S e
f e P P P += （3.2.17）
e P 是单元等效结点载荷（体力和面力引起的等效结点力），e F 是其他单元对该单元的作用力，则单元结点力为e P 与e F 和。 体积力的等效结点载荷：
》
dz dxdy f f f N P P P P e V z y x T e mf e jf e if e f ??????
?????=????
??????= dz dxdy f f f N P P P P e V
z y x i f
e iz e iy e ix e i
f ???????????=????
??????= ),,,(p m j i 面积力的等效结点载荷：
?=e S T e S TdS N P σ
???????????=??????????=e S z y x T S
e m e j e i e S dS
T T T N P P P P σ
??????
?????=??????????=σ
S z y x i S
e iz e iy e ix e
iS dS T T T N P P P P 这里给出两种常见的载荷的等效结点力：
ⅰ）均质单元的自重分配到四个结点的等效结点力，其数值都等于4/gV ρ；ⅱ）设单元的某一边界面上，例如ijm ，受有线性分布载荷，它在m j i ,,三个结点处的强度分别为m j i q q q ,,，则分配到结点i 上的等效结点力的数值为
ijm m j i e i A q q q P ??
?
??++=
212161 ),,(m j i ijm A 为受力面三角形面积。方向与原方向平行。
3.2.4．结构刚度矩阵和结构载荷列阵的集成 由单元分析可得有限元列式为
e e e e F P K +=α （3.2.18）
经叠加，组合，得有限元方程
P K =α 其中
∑==e
N e e K K 1
∑==E
N e e P P 1
式中e K 为扩大后的单元刚度矩阵；e P 为扩大后的单元等效结点载荷；e N 为结构系统的单元数。