空间问题的四面体单元
第三章 轴对称、三维和高次单元
§ 3-2空间问题的四面体单元
空间问题的有限单元法，和平面问题及轴对称问题的有限单元法的原理和分析过程完 全相同。由于空间问题应采用三维坐标系，因此单元的自由度、刚度矩阵的元素个数，方 程组内方程个数等要较平面问题和轴对称问题多，所以空间问题的规模一般比轴对称问题 和平面问题大得多。它要求计算机的内存大，且计算时间长，费用高。这些问题都给三维 有限单元法的具体运用带来许多困难。
和平面问题一样，空间有限单元法采用单元 也是多种多样的，其中最简单的是四节点四面体 单元。采用四面体单元和线性位移模式来处理空 间问题，可以看作平面问题中三角形单元的推广。
在采用四面体单元离散化后的空间结构物 中，一系列不相互重叠的四面体之间仅在节点处 以空间铰相互连接。四节点四面体单元仅在四个 顶点处取为节点，其编号为i,j,m,p 。每个单元的 计算简图如图3-7所示。
在位移法中，取节点位移为基本未知量，四 节点四面体单元共有十二个自由度 （位移分量），
其节点位移列阵为
U i V i W i (i,j,m)
相应的节点力列阵为
U i V
i
w i U j V j w j
U m
T
W m U p V p W p
其子矩阵
图3-7空间四面体单元
F i F j F m F p
其子矩阵
F i U i V i w
一、单元法位移函数
结构中各点的位移是坐标 X 、
y 、z 的函数。 当单元足够小时， 单元内各点的位移可用 简单的线性多项式来近似描述， 即
u
1 2 X
3y 4Z
v
5
6 X
7
y 8
Z
(3-49)
w
0 10X
ny
12Z
曰
2，…，
12是
卜二个待定系数，它们可由单元的节点位移和坐标确定。假定节 点 i,j,m,p 的坐标分别为(x i y i Z i )、
、(x
j y j z j )
、(X m
将它们代入 (3-49)式的第一式可得各个节点在
X 方向的位移
U i
1 2X i 3Y i
4
Z
u j
1 2X j
3Y j
4Z j
U m 1 2
X m 3Y m 4 Z m
U p
1
2
X p
3
Y p
4 Z p
解上述线性方程组，可得到
1 ,
2 ,
3 ,
4 , 再代入
U
6V
[(a i bX
cy d i Z)U i (a j
b j x
(a m b m X C m y
d m z)U m
(a p b p X C
(3-50)
y d p Z )U p ] 1 X i Y i Z i 1 X j y j Z j 1 X m y m Z m
1
X P
Y P
Z P
(3-52)
(3-50)式，得
y m Z m
)、(X p y p Z p )，
5y 3)5 (3-51)
式中1 ,
其中V 为四面体ijmp 的体积，a,b i ,…，c p ,d P 为系数。
为了使四面体的体积 v 不致为负值，单元四个节点的标号 在右手坐标系中，要使得右手螺旋在按照
3-1中单元那样。
综合表达式(3-51)、(3-54)及(3-55)，可以将位移分量表示成为
T
e
f u v w
[ N ]
IN i IN j IN m IN p
e
(3-56) 其中1是三阶的单位矩阵， [N]为形函数矩阵，而各个形函数为
N i
(a i b i x c i y d i z)/6V
(i,m) (3-57) N j (a b i x c y d i z)/6V
(j, p)
1 ,
5 ,
6代表刚性移动
U 0 , V 0 , W o ;系
6个系数反映了刚性转动
W x , W y , W z 和常量
剪应变。这就是说，12个系数充分反映了单元的刚体位移和常量应变。同时，可以证明： 由于位移模式是线性的， 两个相邻单元的共同边界在变形过程中
，始终是相互贴合的，使 得离
散的模型变形中保持为连续体。这样，选用的位移函数满足收敛的充分必要条件，保 证了有限单元法解答收敛于精确解。
a i
(i,j,m,p)
(3-53)
C i,j,m,p
i T j T m 的转向转动时, 必须按照一定的顺序：
向 p 的方向前进，象图
用同样方法, 可以得出其余二个位移分量:
1
6V
(a m 1 6V
(a m [(a i b i x cy d i Z)v (a j b j X C j y d j Z)V j
b m X C m Y
d m Z)V m (a p b p X C p y d p Z)V p ]
[(a i b i x qy d j Z )w (a j b j X c I Y d j Z)W j
b m X C m Y
d m Z)W m (a p b p X C p y d p Z)W p ]
和平面问题相似，(3-49)式中的系数
12代表常量的正应变；其余 (3-54)
(3-55)
二、载荷移置
空间问题的单元载荷移置和平面问题一样，也是根据静力等效原则，将不作用在节点上的集中力、体力、面力移置成作用在节点上的等效节点载荷。其通用公式的形式和平面问题也是一样的，只不过多出一维空间分量。
1.集中力
T
设单元上某点(x,y,z)作用有集中力P P x P y P,
则仍然得到等效节点载荷
R [N]T P (3-58)
这里 e
R [X i Y 乙X j Y j Z
j X m Y m Z m
X p Y p Z p]
2.分布体力
单元上作用有分布体力P [X Y Z]T，则
R e[N]T PdV (3-59)
其中dV是单元中的微分体积，对于直角坐标糸上式为
R e[N]T p dxdydz (3-60) e
3.分布面力
单元的某一边界面S上作用有一分布面力P X Y Z T
R e [N]T P dA
其中dA是边界面S上的微分面积。
4.常见载荷的移置
上列公式是空间问题载荷移置的通用公式。 对于四节点四 面体单元，由于其采用线性位移模式， 采用直接计算虚功的方 法求出节点载荷比较简单。下面介绍常见的二种载荷的移置。
⑴重力
四面体单元的自重为 W 作用在质心C 处(如图3-8)。为 求得节点载荷
X,Y i ,Z i ,可分别假想发生 u * 1 , V * 1或 w ； 1的虚位移。
在U i* 1或V ； 1时，整个单元上各点的均没有 z 方向上
的虚位移，重力 W 不做功，所以 X=Y i =O 。
*
1 “
W
w
c
-，
Z
i 4
4
对于其余三个节点可得同样结论，于是有
e R i
0 0
T
W
(i,j,m,
p)
4
即，对于四节点四面体单元承受的重力载荷，只需要把共
(2)界面压力
设四面体的一个边界面
ijm 上受有一线性分布的压力
所得各节点载荷的方向和分布力的方向相同，要求各节点载荷分量还需乘上相应的方向余 弦。
由上述面力移置结果，可求出任意线性分布面的等效节点载荷。如在 ijm 面受有线性
分布面力在各点强度分别为
q i , q j ，q m ,时，在i 节点的等效载荷为
作用于
1/4。于
ijm 面上的d 点， 二是可得
d 点到ij
边和im 边的距离分别为 m 至U ij 及j 到im 边的距离的
T
T
e
P
P P c
1
,1 1
R
i
-0
q ijm 1 -
0 (3-62)
2 4 4
6 j
2 2
q i ,0,0。很容易看出，该力向p 点移置的等效节点力为零。
由水力学知，总压力P 1
q i ijm ,
3
当W i 1
时，jmp 面上各点的虚位移为零，即
*
1
W b 0，又因bc —bi ,所以有
4 (3-61)
1
移置到每个节点上即可。
4
P ,共在三个节点上的强度分别为
图3-8重力移置
111
P (q i q j q m) jm (i，j，m) (3-63)
6 2 2
三、应力应变矩阵
空间问题几何方程为
y
T
z x y u v w
x y z
u v
y x
u z w
y
w u
x z
z
将四面体单兀之位移表达式
(3-52)、(3-54) 和(3-55) 代入几何方程，即得单兀应变。用节点位移可表示为 e
[B]
E 3i B j B m B p e
(3-64)
式中应变矩阵子矩阵为6X 3矩阵:
b i 0 0
0 C i 0
1 0 0 d i
[B i] c (i，j，m,p) (3-65)
6V c i b i 0
0 d i q
d i 0 b i
由上式可以看出，每一个单元的应变矩阵是一个常量矩阵；因此，采用线性位移模式的四面体单元是常应变单元。这与平面问题中的三角形单元是一样的。而与平面问题的不同之处仅在于应变矩阵的阶数不同。
将表达式(3-16)代入空间问题的物理方程，即可得出用单元节点位移表示的单元应力:
e e
[D] D[B] [S] (3-66)
式中弹性矩阵［D］为
应力矩阵[D]
1 1
0 0
1 2
2(1 )
0 0 0
0 0 0
称
1 2
2(1 )
1 2
2(1 )
[S]
A1
S i S j S m S p (3-67
)
A2
b i
Ab
[S i
]
E(1 ) Ab
6(1 )(1 2 )V A2 G
A2d i
A?A?
C i A?
A?d i
(i,j,m,p)
A2b i 0
A2d i
A2G
0 A2b i
［S］是常量矩阵，所以，四面体单元是
显然，式(3-68)中各元素均为常量，应力矩阵常应
力单元。
四、单元刚度矩阵
空间问题的单元刚度由虚功方程导出。假设该单元发生某虚位移，相应节点虚位移为
e
。此时相应的虚应变为
将上式及式(3-66)代入虚功方程，有
* e T e *
e
T
e
( )F ([B]
)[D][B] dxdydz
v
通过与平面问题一样的处理，并注意到矩阵 [B]中的元素为常量，可以得到
F e [B]T [D][B]dxdydz e [B]T [D][B]
e
V [K]e e (3-69)
v
式中，[K]e
为单元刚度矩阵:
[K]e [B]T [D][B]dxdydz [B]T [D][B]V
(3-70)
e
将式(3-64)和(3-68)式代入，可以得出
其中，[K rs ]e
为3X 3阶方阵:
(r,s=i,j,m,p) (3-72)
有了单元节点力和节点位移之间的关系之后，通过分析每个节点的平衡条件可得到
这个矩阵形式的方程实际上代表了关于
r 节点三个坐标轴方向的力平衡方程式。将关
[B]
K H
K ij K im K ip [K]e
心
K jj
K jm ? K
mi K
mj
K
mm K mp
K pi K pj
K
pm
K pp
(3-71)
[K rs ]e
E(1 ) 36(1 )(1 )V
bR s A 2(g d r d s ) AGb s A 2b r C s AdR s
AbA
Ab r C s AC r C s C r C s A 2(bQ s d r d s )
Adi s
Acd
Ab r d s A z d r b s AC r d s A z d r C s d r d s A 2(b r b s C r C s )
[K rs ]e s
e s i,j ,m, p
R r
e
于结构物所有节点的线性方程式集合起来，可以得到
[K]
式中代表整个结构的节点的位移，是所求之基本未知量；R 代表整个结构的节点载荷；
[ K ]为整体刚度矩阵，其是由每个单元刚度矩阵升阶后组集得到，即
NE
e
[K] [K]e
e1
其为3NP阶方阵。显然，对每一个子矩阵，应有
NE
[K rs] [K rs ]
e1
和平面问题一样，[K]是对称、带状、稀疏矩阵，在消除刚体位移之后，它是正定的。
由平衡方程组可以解出节点位移，随后即可求得所需节点和单元应力。
五、形成四面体的对角线划分方法
在实际计算中，用一系列的四面体来组合成一个空间物体，这个形象是很难想象的。
但是如果先用一系列较为直观的六面体（图3-9）来划分弹性体，然后由计算机来将这些六面体及三棱柱划分为若干个四面体，则要方便得多。同时减少许多准备及输入工作，也为将来结果分析带来方便。
现在介绍一种适合计算机进行自动划分四面体的方法——对角线划分法。
1. 将六面体划分为四面体的方法
通过连接六面体上一些四边形的对角线，可以把一个六面体划分为五个或六个四面体。为叙述方便，先将六面体的八个角点进行局部编号，编号原则是先顶面后底面，对于顶面或底面的节点来说，则是先前后后，从左到右（见图3-9）排列。
(1)
将一个六面体划分为五个四面体
这种方法是先过六面体的一些四边形的对角线，从六面体的四个角上切下四个四面体，
最后剩下中心的一个四面体，共得五个六面体单元。选择被切下的角点不同，有二种不同的划分结果，如图3-10(a )和(b )所示。我们分别称之为A5型划分和B5型划分。A5型划
分所得五个四面体为1246, 1347, 1467, 1567, 4678;而B5型划分则得到1235，2348, 2358 , 2568, 3578五个四面体。
以上二种划分方法的共同特点是，六面体二对面四边形的剖分对角线是交叉的。这就
使得如果一个六面体按A5型划分，那么与之相邻的各个六面体必定要按B5型划分。
图3-10 六面体划分为五个四面体
(a) A5 型剖分；(b)B5 型剖分
(2) 将一个六面体划分成六个四面体
六面体和三棱柱
将六面体划分成六个四面体有很多种划分方法。这里介绍两种，如图3-11所示。它们
的共同特点是，六面体上两对面四边形的剖分对角线是“平行”的。所不同的是在A6型剖分中取大对角线36作为划分线，而在B6型中则是取大对角线45作为划分线。
为清楚起见，可将A6型划分理解为先将六面体沿2367分成两个三棱柱，再将每个三
棱柱分成三个四面体，分别得到1235，2356，3567和2346，3467, 4678六个四面体(见图
3-12。当然，A6型划分也可看成先将六面体沿3456面剖分，得到两个不同于前的三棱柱，
但最后得到的六个四面体是相同的(图3-13)。对于B6型划分，六面体先以折面2457为分
界面拆分成两个三棱柱，如图3-14所示。于是可见，每一个“三棱柱”被划分为三个四面
体，它们分别是1235，2345，3456和2456，4567，4678。同时，也不难证明，若以3456 为分界面按B6型划分将六面体拆成的二个“三棱柱”虽与前面的不同，但是划分成的六个
四面体和前面得到的完全相同。
(a) (b)
图3-11 六面体划分为六个四面体
(a) A6 型剖分； (b)B6 型剖分
图3-12 A6型划分，以折面2376为两个“三棱柱”的分界面
我们看到，A6型和B6型划分，由于其相对四边形的对角线“平行” ，而剖分大对角线
35和45并不在六面体表面上，其相邻的六面体可以全部采用
A6型划分或B6型划分，两
种划分也可以交替使用。一个六面体划分为六个四面体，各四面体的体积大小一般较为均 匀；但是在相等的
六面体数目下，
A6型和B6型划分所产生的四面体单元的总数，要比
A5
型和B5型产生的多六分之一。
此外，在A6型与B6型划分中，如果离散体的节点整体编号是按本节开头所述，从上 到下，从左到右连续进行的；同时每个六面体八个节点的整体编号的大小次序与其局部编 号的大小次序相一致(由小到大)的话，那么在划分中，对底面上任一节点，与它构成四面 体的三个节点中的最小号码，不会比其正上方那个节点的号码更小。这是由于在
A6型及
B6型划分中，注意到在连各个四边形对角线时，不使节点编号之差较大的三个节点出现在 同一个三角形中的
结果。
例如图3-11中，对于1357四边形，我们连接了 35对角线而使1、
7两节点分别属于两个三角形中。这样的划分能获得一个带宽较窄的刚度矩阵。
(3) 编号推算
如果将六面体的八个顶点的节点整个编号置于数组 D[ 1:8 ]中，而将前述图中的局
部编号1?8理解为数组
D[ 1:8 ]的下标时，于是上述问题就转化为：要在有八个元素
图3-13 A6 型划分，以折面3456为两个“三棱柱”的分界面
图3-14 B6 型划分拆成两个“三棱柱”
的数组中按一定规律，每次取四个元素构成一个四面体单元的节点编号问题。对于A6 型和A5 型划分所
得四面体顶点编号的规律性进行一些分析之后，可以导出下列公式，分别表示按预定规律划分成的四面体的
各节点编号：
D[m(1 3(l 1)J (1 m)(l J)]
D[m(3l J 1) (1 m)(1 l J(J 1)/2)]
D[m(2(1 l J) J(l J) (1 m)(2 l J(5 J)/2)] (3-73)
D[m(5 l J(3 J) /2 (1 m)(4 l J)]
m 0,1)
(l 1,2; J 0,1,
2;
通过直接代入数字检验，知道m=0对应着A6型划分所得的六个四面体，m=1则对应着
A5型划分所得的五个四面体，(此时l=J=2形成的四面体应舍去)。例如m=0的情况，当1=1 ，J=0，1，2时，得到D[1]D[2]D[3]D[5] ，D[2]D[3]D[5]D[6] ，D[3]D[5] D[6]D[7] 三个四面
体，而l=2，J= 0, 1, 2 时，得到D[2]D[3]D[4]D[6] ，D[3]D[4]D[6]D[7] ，D[4]D[6]D[7]D[8]
三个四面体，与前述结果一致。
对于B6 型及B5 型划分,同样可以导出一个相似的计算公式。但是也可以利用(3-73) 式，只需将原来D〔1：8〕中元素位置作一定更动。更动的办法是，对于B6型划分，将8、6、4、2位置
的元素置于1、2、3、4位置上,将7、5、3、1 位置的元素置于5, 6, 7, 8 位置；对于B5型划分，将5，6，7, 8位置的元素分别和1, 2, 3, 4位置的元素交换即可。
以上讨论的是四面体八个顶点编号自动产生的情形，如果八个顶点的整体编号是外部输入的,则在输入前作位
置的更换,不必增添更换位置的附加程序。
2. 三棱柱划分为四面体的方法
在弹性体的实际分割中, 边角位置常出现六节点的三棱柱。由前述知道, 任一三棱柱, 可以看成是某个
六面体的一半,并能划分成三个四面体。为直接使用公式(3-73) ,只需要
把图3-9 所示之三棱柱视情况添加两个节点, 使之构成图3-11 所示六面体, 利用式(3-73) 时,置数组D[1：8 ]中两个添加点对应元素为零,并规定,所得四面体的四个节点编号中若有一个号码是零,就表示为空单元,
不进行单元编号,从而只留下三个四面体单元。具体地说,如图3-9 所示的三棱柱,可以认为是图3-12 中六面体以折面2376为分界面拆开的两个三棱柱中的一个——三棱柱234678。于是将数组D[1：8]中的第一
个和第五个元
素置零,其他元素则存入相应的节点整体编号。
实际计算中还可能遇到五节点的五面体,同样可以通过增添三个节点,使之构成一个六面体。此时应在
D[1 ：8] 的适当位置补上三个零,使用公式(3-73) 后,便得两个非空的四面体单元节点号数组。
最后,应当指出, 在我
们的划分中,都没有要求六面体某个四边形是平面的,即没有要求四个顶点在同一个
平面内。因为在划分中出现的四边形，都是用连对角线的办法一分为二的，因此，任一四边形都可以是由二
个折面构成的，也就是说，八个角点的位置可以
是空间任意安置的。这对于将弹性体划分为一系列六面体、三棱柱来说是便利的。