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

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

§ 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 V

V

其子矩阵

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

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 z

则仍然得到等效节点载荷

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 7， 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 -

(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]

E3i B j B m B p e(3-64)

式中应变矩阵子矩阵为6X 3矩阵:

b i00

0C i0

1 00d i

[B i] c (i，j，m,p)(3-65)

6V c i b i0

0d i q

d i0b i

由上式可以看出，每一个单元的应变矩阵是一个常量矩阵；因此，采用线性位移模式的四面体单元是常应变单元。这与平面问题中的三角形单元是一样的。而与平面问题的不同之处仅在于应变矩阵的阶数不同。

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

e e

[D] D[B] [S] (3-66)

式中弹性矩阵［D］为

应力矩阵[D]

11

00

12

2(1)

000

000

称

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) (3-68)

A2b i0

A2d i A2G

0A2b 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阶方阵：

b r b s A 2(GC s d r d s )

有了单元节点力和节点位移之间的关系之后，通过分析每

个节点的平衡条件可得到

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

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

[B]

K H K i

. ij

K.

K-

II [K]e

J 1

K mi

JJ

K mj

K pi K pj

K im K ip

K jm ?

K

mm K

mp

(3-71)

K

pm

K pp

Ab 「C s AWs Ab r d s A z d r b s :s A 2(b r b s d r d s ) AC r d s A z d r C s

Adi s AC r d s d r d s A 2(bR s C r C s )

[K ]e

E(1

)

AC r b s Ad r b s

AbrC s AbA

rs

36(1 )(1 )V

[K rs ]e s

e s i,j ,m, p

R r

e

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

[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) 将一个六面体划分成六个四面体

六面体和三棱柱

(a)

(b)

将六面体划分成六个四面体有很多种划分方法。这里介绍两种，如图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(l1)J(1 m)(l J)]

D[m(3l J1)(1 m)(1l 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(3J) /2 (1m)(4 l J)]

m 0,1)

(l 1,2; J0,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) 后，便得两个非空的四面体单元节点号数组。

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