03-02 空间问题的四面体单元解析

第三章 轴对称、三维和高次单元

§3-2 空间问题的四面体单元

空间问题的有限单元法，和平面问题及轴对称问题的有限单元法的原理和分析过程完全相同。由于空间问题应采用三维坐标系，因此单元的自由度、刚度矩阵的元素个数，方程组内方程个数等要较平面问题和轴对称问题多，所以空间问题的规模一般比轴对称问题和平面问题大得多。它要求计算机的内存大，且计算时间长，费用高。这些问题都给三维有限单元法的具体运用带来许多困难。 和平面问题一样，空间有限单元法采用单元也是多种多样的，其中最简单的是四节点四面体单元。采用四面体单元和线性位移模式来处理空

间问题，可以看作平面问题中三角形单元的推广。 在采用四面体单元离散化后的空间结构物中，一系列不相互重叠的四面体之间仅在节点处以空间铰相互连接。四节点四面体单元仅在四个顶点处取为节点，其编号为i,j,m,p 。每个单元的计算简图如图3-7所示。

在位移法中，取节点位移为基本未知量，四节点四面体单元共有十二个自由度(位移分量)，其节点位移列阵为

{}[

]

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 =??????????????=δδδδδ

其子矩阵 {}[]i i

i

i w v u =δ (i,j,m)

相应的节点力列阵为

{}[]

T

p m

j i

e F F F F F -

图3-7 空间四面体单元

其子矩阵 {}[]T

i i i i W V U F =

一、单元法位移函数

结构中各点的位移是坐标x 、y 、z 的函数。当单元足够小时，单元内各点的位移可用简单的线性多项式来近似描述，即

??

?

??

+++=+++=+++=z y x w z y x v z y x u 121110087654321αααααααααααα (3-49) 式中1α，2α，…，12α是十二个待定系数，它们可由单元的节点位移和坐标确定。假定节点i,j,m,p 的坐标分别为(i x i y i z )、(j x j y j z )、(m x m y m z )、 (p x p y p z )，将它们代入(3-49)式的第一式可得各个节点在x 方向的位移

???

?

???

??????

?+++=+++=+++=+++=p p p p m m m m j j j j i i i i z y x u z y x u z y x u z y x u 4321432143214321αααααααααααααααα (3-50)

解上述线性方程组，可得到1α，2α，3α，4α，再代入(3-50)式，得

]

)()()()[(61

p p p p p m m m m m j

j j j j i i i i i u z d y c x b a u z d y c x b a u z d y c x b a u z d y c x b a V

u +++-+++++++-+++=

(3-51) 其中V 为四面体ijmp 的体积，a i ,b i ，…，c p ,d p 为系数。

p

p

p

m m m j j j i i i z y x z y x z y x z y x V 1111=

(3-52)

p p p m m m

j

j j

i z y x z y x z y x a = 111j

j

i m m p p y z b y z y z = p

p

m m j j i z x z x z x c 111= 1

11p

p

m m j j i y x y x y x d = (i,j,m,p) (3-53)

为了使四面体的体积V 不致为负值，单元四个节点的标号i,j,m,p 必须按照一定的顺序：在右手坐标系中，要使得右手螺旋在按照i →j →m 的转向转动时，向p 的方向前进，象图3-1中单元那样。

用同样方法，可以得出其余二个位移分量：

]

)()()()[(61

p p p p p m m m m m j

j j j j i i i i i v z d y c x b a v z d y c x b a v z d y c x b a v z d y c x b a V

v +++-+++++++-+++=

(3-54) ])()()()[(61

p p p p p m m m m m j

j j j j i i i i i w z d y c x b a w z d y c x b a w z d y c x b a w z d y c x b a V

w +++-+++++++-+++=

(3-55) 综合表达式(3-51)、(3-54)及(3-55)，可以将位移分量表示成为

{}[]{}

[

]

{}

e

p m

j i

e

T

IN IN IN IN N w v u

f δδ===][ (3-56)

其中I 是三阶的单位矩阵，[N]为形函数矩阵，而各个形函数为

?

??

+++-=+++=),(6/)(),(6/)(p j V

z d y c x b a N m i V z d y c x b a N i i i i j i i i i i (3-57) 和平面问题相似，(3-49)式中的系数1α，5α，6α代表刚性移动0u ，0v ，0w ；系数2α，7α，12α代表常量的正应变；其余6个系数反映了刚性转动x w ，y w ，z w 和常量剪应变。这就是说，12个系数充分反映了单元的刚体位移和常量应变。同时，可以证明：由于位移模式是线性的，两个相邻单元的共同边界在变形过程中 ，始终是相互贴合的，使

得离散的模型变形中保持为连续体。这样，选用的位移函数满足收敛的充分必要条件，保证了有限单元法解答收敛于精确解。

二、载荷移置

空间问题的单元载荷移置和平面问题一样，也是根据静力等效原则，将不作用在节点上的集中力、体力、面力移置成作用在节点上的等效节点载荷。其通用公式的形式和平面问题也是一样的，只不过多出一维空间分量。

1. 集中力

设单元上某点(x,y,z)作用有集中力{}T

x

y

z P P P P ??=??

则仍然得到等效节点载荷

{}{}P N R T ][= (3-58)

这里 {}T p p

p m m m j j j i i i

e

Z Y X Z Y X Z Y X Z Y X R ][=

2. 分布体力

单元上作用有分布体力{}T Z Y

X

P ][=，则

{}{}dV P N R T e ?=][ (3-59)

其中dV 是单元中的微分体积，对于直角坐标系上式为

{}{}dxdydz

p N R e

T e ???=][ (3-60)

3. 分布面力

单元的某一边界面 S 上作用有一分布面力{}[

]T

Z Y

X

P =

则 {}{}

dA P N R T

e

?

=][

其中dA 是边界面S 上的微分面积。

4. 常见载荷的移置

上列公式是空间问题载荷移置的通用公式。对于四节点四面体单元，由于其采用线性位移模式，采用直接计算虚功的方

法求出节点载荷比较简单。下面介绍常见的二种载荷的移置。

(1) 重力

四面体单元的自重为W ，作用在质心C 处(如图3-8)。为求得节点载荷X i ,Y i ,Z i ,可分别假想发生1*=i u ,1*=i v 或

1*=i w 的虚位移。

在1*=i u 或1*=i v 时，整个单元上各点的均没有z 方向上的虚位移，重力W 不做功，所以X i =Y i =0。

当1*=i w 时，jmp 面上各点的虚位移为零，即0*

=b w ，又因bi bc 4

1

=

,所以有 41*=

c w ， 4

W

Z i -= 对于其余三个节点可得同样结论，于是有

{}

T

e

i W R ?????

?

-=400(i,j,m,p) (3-61)

即，对于四节点四面体单元承受的重力载荷，只需要把共

4

1

移置到每个节点上即可。 (2) 界面压力

设四面体的一个边界面ijm 上受有一线性分布的压力P ，共在三个节点上的强度分别为q i ,0,0。很容易看出，该力向p 点移置的等效节点力为零。由水力学知，总压力ijm q P i ?=

3

1

,作用于ijm 面上的d 点，d 点到ij 边和im 边的距离分别为m 到ij 及j 到im 边的距离的1/4。于是可得

{}

T

ijm i T

e

i q P P P R ??

?

????=??????=0212

1

161044

2

(3-62) 所得各节点载荷的方向和分布力的方向相同，要求各节点载荷分量还需乘上相应的方向余弦。

图3-8 重力移置

由上述面力移置结果，可求出任意线性分布面的等效节点载荷。如在ijm 面受有线性分布面力在各点强度分别为q i ，q j ，q m ，时，在i 节点的等效载荷为

ijm m j i i q q q P ?++=

)2

1

21(61 (i ，j ，m) (3-63)

三、应力应变矩阵

空间问题几何方程为

{}T

z y x z y x z u x w y

w z u x v y u z w y v

x u ??

??

????+????+????+????????=?????

??????

???????????=γγγεεεε 将四面体单元之位移表达式(3-52)、(3-54)和(3-55)代入几何方程，即得单元应变。用节

点位移可表示为

{}{}[]

{}e

p m

j i

e B B B B B δδε--==][ (3-64)

式中应变矩阵子矩阵为6×3矩阵：

??????

???

??????

?????=i i

i i i i

i i i i b d c d b c d c b V B 0

000

000061][ (i,j,m,p) (3-65) 由上式可以看出，每一个单元的应变矩阵是一个常量矩阵；因此，采用线性位移模式

的四面体单元是常应变单元。这与平面问题中的三角形单元是一样的。而与平面问题的不同之处仅在于应变矩阵的阶数不同。

将表达式(3-16)代入空间问题的物理方程，即可得出用单元节点位移表示的单元应力：

{}{}[]{}{}e e S B D D δδεσ][][][=== (3-66)

式中弹性矩阵[]D 为

???

???????????

????????

???????

?---------=)1(22100

0)1(2210000)1(221000111111][μμμμμμ

μμ

μμμμ称

对D 应力矩阵 [

]

p m j i

S S S S ][=S (3-67)

令 μμ

-=11A ， )

1(2212μμ--=

A

则

???

???

???

??????

?????-+-=i i

i i i i b A d A c A d A b A V E S 22222i

2i i

1i 1i 1i i 1i 1i 1i 0

00c A d c A b A d A c b A d A c A b )21)(1(6)

1(][μμμ (i,j,m,p) (3-68) 显然，式(3-68)中各元素均为常量，应力矩阵[Ｓ]是常量矩阵，所以，四面体单元是

常应力单元。

四、单元刚度矩阵

空间问题的单元刚度由虚功方程导出。假设该单元发生某虚位移，相应节点虚位移为

{}*e

δ。此时相应的虚应变为

{}{}e

B **

][δε=

将上式及式(3-66)代入虚功方程，有

{}

{}{}

{}dxdydz B D B F e

T e

v

e

T e

δδδ]][[)]([)(**???=

通过与平面问题一样的处理，并注意到矩阵[B]中的元素为常量，可以得到

{}{}{}{}e

e e

T e

v

T e K V B D B dxdydz B D B F δδδ][]][[][]][[][===??? (3-69)

式中，e K ][为单元刚度矩阵：

???==e

T T e V B D B dxdydz B D B K ]][[][]][[][][ (3-70)

将式(3-64)和(3-68)式代入，可以得出

???

?

??

?????

???-------=pp pm

pj

pi

mp mm mj

mi jp jm jj ji

ip im ij ii e K K K K K K K K K K K K

K K K K K ][ (3-71)

其中，e rs K ][为3×3阶方阵：

212121*********()

(1)[]()36(1)(1)()r s r s r s r s r s r s r s e rs r s r s r s r s r s r s r s r s r s r s r s r s r s r s b b A c c d d A b c A c c A b d A d b E K A c b A b c c c A b b d d A c d A d c V

A d b A b d A d c A c d d d A b b c c μμμ++++??-??=++++??+-??++++??

(r,s=i,j,m,p) (3-72)

有了单元节点力和节点位移之间的关系之后，通过分析每个节点的平衡条件可得到

{}{}∑∑∑==e

r

e p

m j i s s

e rs

R K ,,,][δ

这个矩阵形式的方程实际上代表了关于r 节点三个坐标轴方向的力平衡方程式。将关

于结构物所有节点的线性方程式集合起来，可以得到

{}{}R K =δ][

式中{}δ代表整个结构的节点的位移，是所求之基本未知量；{}R 代表整个结构的节点载荷；

][K 为整体刚度矩阵，其是由每个单元刚度矩阵升阶后组集得到，即

∑==NE

e e K K 1

][][

其为3NP 阶方阵。显然，对每一个子矩阵，应有

∑==NE

e e rs rs K K 1

][][

和平面问题一样，][K 是对称、带状、稀疏矩阵，在消除刚体位移之后，它是正定的。

由平衡方程组可以解出节点位移，随后即可求得所需节点和单元应力。

五、形成四面体的对角线划分方法

在实际计算中，用一系列的四面体来组合成一个空间物体，这个形象是很难想象的。但是如果先用一系列较为直观的六面体(图3-9)来划分弹性体，然后由计算机来将这些六面体及三棱柱划分为若干个四面体，则要方便得多。同时减少许多准备及输入工作，也为将来结果分析带来方便。

现在介绍一种适合计算机进行自动划分四面体的方法——对角线划分法。 1. 将六面体划分为四面体的方法

通过连接六面体上一些四边形的对角线，可以把一个六面体划分为五个或六个四面体。为叙述方便，先将六面体的八个角点进行局部编号，编号原则是先顶面后底面，对于顶面或底面的节点来说，则是先前后后，从左到右(见图3-9)排列。

图3-9 六面体和三棱柱

(1) 将一个六面体划分为五个四面体

这种方法是先过六面体的一些四边形的对角线，从六面体的四个角上切下四个四面体，最后剩下中心的一个四面体，共得五个六面体单元。选择被切下的角点不同，有二种不同的划分结果，如图3-10(a )和(b )所示。我们分别称之为A5型划分和B5型划分。A5型划分所得五个四面体为1246，1347，1467，1567，4678；而B5型划分则得到1235，2348，2358，2568，3578五个四面体。

以上二种划分方法的共同特点是，六面体二对面四边形的剖分对角线是交叉的。这就使得如果一个六面体按A5型划分，那么与之相邻的各个六面体必定要按B5型划分。

(a) (b)

图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为两个“三棱柱”的分界面

图3-13 A6型划分，以折面3456为两个“三棱柱”的分界面

图3-14 B6型划分拆成两个“三棱柱”

我们看到，A6型和B6型划分，由于其相对四边形的对角线“平行”，而剖分大对角线35和45并不在六面体表面上，其相邻的六面体可以全部采用A6型划分或B6型划分，两种划分也可以交替使用。一个六面体划分为六个四面体，各四面体的体积大小一般较为均匀；但是在相等的六面体数目下，A6型和B6型划分所产生的四面体单元的总数，要比A5型和B5型产生的多六分之一。

此外，在A6型与B6型划分中，如果离散体的节点整体编号是按本节开头所述，从上到下，从左到右连续进行的；同时每个六面体八个节点的整体编号的大小次序与其局部编号的大小次序相一致(由小到大)的话，那么在划分中，对底面上任一节点，与它构成四面体的三个节点中的最小号码，不会比其正上方那个节点的号码更小。这是由于在A6型及B6型划分中，注意到在连各个四边形对角线时，不使节点编号之差较大的三个节点出现在同一个三角形中的结果。例如图3-11中，对于1357四边形，我们连接了35对角线而使1、7两节点分别属于两个三角形中。这样的划分能获得一个带宽较窄的刚度矩阵。

(3) 编号推算

如果将六面体的八个顶点的节点整个编号置于数组D[１：８]中，而将前述图中的局

部编号１～８理解为数组D[１：８]的下标时，于是上述问题就转化为：要在有八个元素的数组中按一定规律，每次取四个元素构成一个四面体单元的节点编号问题。对于A6型和A5型划分所得四面体顶点编号的规律性进行一些分析之后，可以导出下列公式，分别表示按预定规律划分成的四面体的各节点编号：

)

1,0;2,1,0;2,1()]4)(1(2/)3(5([)]2/)5(2)(1()()1(2([)]

2/)1(1)(1()13([)]

)(1()1(31([===++-+-++-++-+++-++++-+-++-+-+m J I J I m J J I m D J J I m J I J J I m D J J I m J I m D J I m J I m D (3-73) 通过直接代入数字检验，知道 m=0对应着A6型划分所得的六个四面体，m=1则对应着A5型划分所得的五个四面体，(此时I=J=2形成的四面体应舍去)。例如m=0的情况，当I=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]三个四面体，而I=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[１：８]中两个添加点对应元素为零，并规定，所得四面体的四个节点编号中若有一个号码是零，就表示为空单元，不进行单元编号，从而只留下三个四面体单元。具体地说，如图3-9所示的三棱柱，可以认为是图3-12中六面体以折面2376为分界面拆开的两个三棱柱中的一个——三棱柱234678。于是将数组D[1：8]中的第一个和第五个元素置零，其他元素则存入相应的节点整体编号。

实际计算中还可能遇到五节点的五面体，同样可以通过增添三个节点，使之构成一个六面体。此时应在D[1：8]的适当位置补上三个零，使用公式(3-73)后，便得两个非空的四面体单元节点号数组。

最后，应当指出，在我们的划分中，都没有要求六面体某个四边形是平面的，即没有要求四个顶点在同一个平面内。因为在划分中出现的四边形，都是用连对角线的办法一分为二的，因此，任一四边形都可以是由二个折面构成的，也就是说，八个角点的位置可以是空间任意安置的。这对于将弹性体划分为一系列六面体、三棱柱来说是便利的。