等参单元

5．等参单元
本章包括以下内容： 5.1等参单元的基本概念 5.2四边形八节点等参单元 5.3等参单元的单元分析 5.4六面体等参单元
5.1等参单元的基本概念
在进行有限元分析时，单元离散化会带来计算误差，主要采用两种方法来降低单元离散
化产生的误差：1）提高单元划分的密度，被称为h 方法（h-method ）；2）提高单元位移函数多项式的阶次，被称为p 方法（p-method ）。

在平面问题的有限单元中，我们可以选择四结点的矩形单元，如图5-1所示，该矩形单元在x 及y 方向的边长分别为2a 和2b 。

图5-1 四结点矩形单元
同第三章的方法类似，将单元的位移模式选为，
xy a y a x a a u 4321+++= xy a y a x a a v 8765+++=
（5-1）
可得到，
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 +++=
（5-2）
形态函数为， )1)(1(41b y a x N i --=
)1)(1(41b y a x N j -
+
=
)1)(1(41b y a x N m ++= )1)(1(4
1b
y a
x N p +
-
=
（5-3）
上述单元位移模式满足位移模式选择的基本要求： 1）反映了单元的刚体位移和常应变， 2）单元在公共边界上位移连续。

在矩形单元的边界上，坐标x 和y 的其中一个取常量，因此在边界上位移是线性分布的，由两个结点上的位移确定。

与三结点三角形单元相比，四结点矩形单元的位移模式是坐标的二次函数，能够提高计
算精度，但也有显著的缺点，两种单元的比较如下。

表5-1 三结点三角形单元与四结点矩形单元比较
如果任意形状的四边形四结点单元采用矩形单元的位移模式，则在公共边界上不满足位移连续性条件。

为了既能得到较高的计算精度，又能适应复杂的边界形状，可以采用坐标变换。

图5-2任意四结点四边形单元
图5-3四结点正方形单元
在图5-2所示的任意四边形单元上，用等分四条边的两族直线分割四边形，以两族直线
的中心为原点，建立局部坐标系),(ηξ，沿ξ及η增大的方向作为ξ轴和η轴，并令四条边上的ξ及η值分别为1±。

为了求出位移模式，以及局部坐标与整体坐标之间的变换式，在局部坐标系中定义一个四结点正方形单元，如图5-3所示。

参照矩形单元，四结点正方形单元的位移模式为，
44332211u N u N u N u N u +++=
44332211v N v N v N v N v +++=
（5-4）
其中， )1)(1(411ηξ--=N )1)(1(412ηξ-+=N )1)(1(413ηξ+-=N )1)(1(4
14ηξ++=
N
（5-5）
四个结点的坐标为),(i i ηξ，定义新的变量， ξξξi =0，ηηηi =0 (i=1，2，3，4) (5-6)
形态函数表示为，
)1)(1(4
100ηξ++=
i N
（i=1，2，3，4） （5-7）
把ξ及η作为任意四边形单元的局部坐标，把（5-4）的位移模式和（5-7）的形态函数
用于任意形状的四边单元，可得：
1）在四个结点处可以得到结点的位移；
2）在单元的四条边上，位移线性变化，保证了单元公共边界上位移的连续性。

因此给出任意四边形单元的结点位移就能得到整个单元上的位移，（5-4）的位移模式就是所要找的正确的位移模式。

把局部坐标与整体坐标的变换式也取为，
44332211x N x N x N x N x +++=
44332211y N y N y N y N y +++=
（5-8）
将坐标变换式用于任意四边形单元，可得： 1）在四个结点处给出结点的整体坐标，
2）在四条边上的整体坐标是线性变化的。

只要给出任意四边形单元四个结点的整体坐标，用（5-8）式就可以建立局部坐标系中的正方形单元和整体坐标系中的任意四边形单元之间的坐标变换关系。

把图5-3中的局部坐标系中的正方形单元称为基本单元。

把图5-2中的在整体坐标系中的任意四边形单元看作由基本单元通过坐标变换得来的，
称为实际单元。

单元几何形状和单元内的未知量采用相同数目的结点参数以及相同的插值函数进行变换，称为等参变换。

采用等参变换的单元，称为等参单元。

由于形态函数i N ，正好反映了单元形状的变化，也称为形函数（Shape function ）。

采用等参单元，使我们可以在局部坐标系中的规则单元上进行单元分析，然后在映射到实际单元上。

等参单元同时具有计算精度高和适用性好的特点，是有限元程序中主要采用的单元形式。

5.2四边形八节点等参单元
为了更好地反映物体内的应力变化，适应曲线边界，在弹性力学平面问题的分析中经常使用四边形八节点等参单元。

如图5-4所示，由于每条边上增加了一个结点，单元的边是一条二次曲线，可以更好地适应曲线边界，
图5-4四边形八结点单元
图5-5 八结点基本单元
对于等参单元，先在图5-5所示的八结点基本单元上进行分析。

八结点单元一共有16
个已知的结点位移分量，基本单元中取如下的位移模式： 2
8272652
4321ξηηξηξηξηξa a a a a a a a u +++++++=
2
82
72
6524321ξηηξηξηξ
ηξb b b b b b b b v +++++++=
（5-9）
该位移模式实际上是一个双二次函数，待定系数由结点位移分量确定。

在单元的每条边
上，局部坐标1±=ξ或1±=η，位移是局部坐标ξ或η的二次函数，完全由边上的三个结点的位移值确定，所以这个位移模式满足位移连续性条件。

实际单元内的位移用形函数表示为，
i i i u N u ),(8
1ηξ∑
==
i i i v N v ),(8
1
ηξ∑
==
（5-10）
其中的形函数为： )1)(1)(1(411-----=ηξηξN
)1)(1)(1(4
13---+=
ηξηξN
)1)(1)(1(415-+++=ηξηξN )1)(1)(1(417-+-+-=ηξηξN )1)(1(212
2ηξ--=N )1)(1(2126ηξ+-=N )1)(1(2124ξη+-=N
)1)(1(2
128ξη--=
N
将形函数归纳为，
⎪⎪⎪
⎩
⎪⎪⎪⎨⎧=+-=+-=+++=)
8,4()1)(1(21)
6,2()1)(1(21)7,5,3,1()
)(1)(1(41
),(22
i i i N i i i i i i i ξξηηηξηηξξηηξξηξ （5-11）
形函数),(ηξi N 在单元的i 结点上的值为1，在其它结点上的值均为0。

坐标变换式采用如下相似的公式，
i
i i i
i i y N y x N x ),(),(8
1
8
1ηξηξ∑
∑
===
=
（5-12）
将1=ξ代入公式（5-12），可以得到单元345边在整体坐标下的参数方程： f
e d y c b a x ++=++=ηηηη22
（5-13）
可见在整体坐标系中，单元的边是一条抛物线或退化为一条直线。

图5-6 ANSYS提供的Plane82单元
如图5-6所示，ANSYS提供的PLANE82单元是一个四边形八结点等参单元，局部坐标定义为s和t，如图5-7所示。

PLANE82单元可以退化为三角形六结点单元。

图5-7 Plane82的基本单元
ANSYS理论手册中给出的PLANE82单元的位移模式如图5-8所示，位移模式与公式（5-10）展开后是一样的。

图5-8 Plane82单元位移模式
5.3等参单元的单元分析
在本节，以平面问题的四边形八结点等参单元为例，介绍构造等参单元的单元刚度矩阵的基本过程。

弹性力学平面问题的单元刚度矩阵为，
[]⎰⎰=tdxdy B D B K T
e
]][[][
单元的应变为，
e
B }]{[}{δε=
单元的结点位移，
{}T
e
v u v u v u ]...
[882
211=δ
将形函数代入后，可以得到应变的矩阵表达式，
{}⎪⎪⎪⎪⎭
⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⎥⎥⎥⎥⎥⎥
⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂=8811882211
8
2
1821
.....
...0
...0v u v u x
N y
N x
N y
N x
N y
N
y N y
N y N x
N x
N x N
ε （5-14）
可得应变矩阵的分块矩阵，
⎥⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂=x
N y
N
y
N x
N
B i i
i i
i 0
0][ （5-15）
由于等参单元的形函数是局部坐标),(ηξ的函数，因此应变矩阵[B]也是局部坐标),(ηξ的函数。

形成等参单元的单元刚度矩阵需要在整体坐标系中对局部坐标的函数进行积分，包括以下三个基本步骤：1）计算用局部坐标表示的形函数),(ηξi N 对整体坐标x 、y 的偏导数；2）将整体坐标系中的面积积分转换为在局部坐标系中的面积积分；3）用数值积分计算出单元刚度矩阵中的元素。

（一）计算形函数对整体坐标x ，y 的偏导数
由于局部坐标与整体坐标之间存在坐标转换关系，因此形函数N i
是局部坐标的函数，同时也可以看作是整体坐标的函数。

由复合函数求导法则可得：
η
η
η
ξξξ∂∂∂∂+
∂∂∂∂=
∂∂∂∂∂∂+∂∂∂∂=∂∂y y N x x N N y y N x x N N i i i i i i
（5-16）
⎪⎪⎭
⎪⎪⎬⎫
⎪⎪⎩⎪⎪⎨⎧∂∂∂∂⎥⎥⎥⎥⎦⎤⎢⎢⎢
⎢⎣⎡∂∂∂∂∂∂∂∂=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂y
N x N y x y x
N N i
i
i i ηηξξηξ 定义，
⎥⎥⎥⎥⎦
⎤
⎢⎢⎢
⎢⎣⎡∂∂∂∂∂∂∂∂=ηηξξy x y x
J ][ （5-17）
⎪⎪⎭
⎪⎪⎬⎫
⎪⎪
⎩⎪⎪⎨⎧∂∂∂∂=⎪⎪⎭
⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂y N
x N J N N i
i
i i
][η
ξ （5-18）
矩阵[J]称为雅可比矩阵（Jacobian Matrix ），单元的整体坐标可以形函数来表示，因此用
坐标变换公式可以计算雅可比矩阵。

⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢
⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂=88
2
2
11
821821......
...
...][y x y x y x N N N N N N J ηη
η
ξξ
ξ （5-19）
由（5-18）可得， ⎪⎪⎭
⎪⎪⎬⎫
⎪⎪⎩⎪⎪⎨
⎧∂∂∂∂=⎪⎪⎭⎪⎪⎬⎫
⎪⎪⎩⎪⎪⎨⎧∂∂∂∂-η
ξ
i
i i
i N N J y
N
x N 1][ （5-20）
⎥⎥⎥⎥⎦
⎤
⎢⎢⎢
⎢⎣⎡∂∂∂∂-∂∂-
∂∂==-ξη
ξηx x
y y
J J J J 1][]
[*1
（5-21）
将（5-21）代入（5-20）得到， ⎪⎪⎭
⎪⎪⎬⎫
⎪⎪⎩⎪⎪⎨
⎧∂∂∂∂+∂∂∂∂-∂∂∂∂-
∂∂∂∂=⎪⎪⎭⎪⎪⎬⎫
⎪⎪⎩⎪⎪⎨⎧∂∂∂∂ηξξ
ηηξξηi i
i i i
i
N x N x N y N y J y
N
x N 1 （5-22）
（二）将整体坐标系中的面积积分转换为在局部坐标系中的面积积分 []⎰⎰=tdxdy B D B K T
e
]][[][
（5-23）
在整体坐标系中，面积微元为x 方向和y 方向微矢量的叉乘的模量，
y d x d dA
⨯= （5-24）
ηηξξ
d x d x x d ∂∂+∂∂=
ηηξξ d y d y y d ∂∂+∂∂=
η
ξξ
ηη
ξηη
ξξηηξξd d y x y x d x d y d x d x dA )()
()(
∂∂∂∂-
∂∂∂∂=
∂∂+∂∂⨯∂∂+∂∂=
ηξd d J dA = （5-25）
代入（5-23），得到单元刚度矩阵在局部坐标系中的积分公式：
[]
ηξd d J t B B B D B B B K T
e
]...
][[]...
[82
1
82
11
1
11⎰
⎰--=
（5-26）
单元刚度矩阵中的任意一个分块矩阵的积分公式为，
ηξd d J t B D B K s T
r rs ]][[][][1
1
1
1
⎰
⎰--=
（5-27）
（三）用数值积分计算出单元刚度矩阵中的元素
等参单元刚度矩阵的每个元素都是局部坐标的函数，在有限元程序中不用解析的办法来计算局部坐标系中的积分，而采用数值积分方法。

通常采用高斯积分方法计算单元刚度矩阵中的元素。

高斯积分方法预先定义了积分点和相应的加权系数，求出被积分的函数在指定积分点上的数值，加权后求和，就得到了该函数的积分。

这种方法具有比较高的计算精度。

已经证明，采用n 个积分点的高斯积分可以达到2n-1阶的精度，也就是说，如果被积分的函数是2n-1次多项式，用高斯积分可以得到精确的积分结果。

例如，ANSYS 提供的PLANE82单元采用的积分点数为22⨯。

关于高斯积分方法的更多内容参见，王勖成、邵敏编著“有限元方法基本原理和数值方法”的4.5节数值积分方法。

将作用在单元上的外载荷同样表示为局部坐标的函数，就可以在局部坐标下完成单元的
载荷移置。

体力移置的公式为， ηξd d J t p N R T
e
}{][}{1
1
1
1
⎰
⎰--=
（5-28）
面力移置的公式也类似，例如在1=ξ的边上受到面力作用，
ηξξd J
t P N R T
e
1
11
1
}{][}{==-⎰
=
（5-29）
在点),(00ηξ集中力移置的公式为， }{][}{),(00P N R T
e
ηξ=
（5-30）
5.4六面体等参单元
多数弹性力学问题需要按照三维空间问题来求解。

三维弹性力学问题的有限元法的基本步骤与平面问题的步骤一样，包括单元离散化、选择单元位移模式、单元分析、整体分析和方程求解。

在分析三维问题时，所选择的单元主要为四面体单元和六面体单元。

每个单元节点上定义有三个位移分量u 、v 、w 。

三维问题有限元法有以下两个主要难点：
1）单元划分比较复杂
无法采用人工方法完成复杂三维实体的单元划分，需要有功能强大的单元划分程序，从CAD 模型直接生成离散的单元网格。

现在的有限元软件可以读入IGES 、STL 等格式的图形交换文件。

六面体单元的计算精度比较高，但是对于复杂三维实体无法实现六面体单元的自动划分。

采用四面体单元能够实现单元自动划分，但是四面体单元的计算精度比较低。

2）计算规模大
三维问题的单元数目大，节点自由度多，导致计算规模大，对计算机硬件的要求很高。

为缩短计算时间，有许多问题需要采用巨型计算机，如CRAY ，或并行计算机。

常用的三维等参单元有六面体八节点等参单元和六面体二十结点等参单元。

等参单元的位移模式和坐标变化式采用相同的形函数，
i
i n
i i i n
i i i n
i w N w v N v u N u ),,(),,(),,(111ζηξζηξζηξ∑
∑
∑
====
==
（5-31）
i
i n
i i i n
i i i n
i z N z y N y x N x ),,(),,(),,(1
11ζηξζηξζηξ∑
∑
∑
====
==
（5-32）
ANSYS 提供的Solid45单元就是六面体八节点等参单元，每个节点有代表x 、y 、z 三个方向位移的三个自由度（DOF ，Degree of Freedom ），可以退化为五面体棱柱和四面体单元，如图5-9所示。

单元局部坐标为r 、s 、t ，如图5-10所示，六面体八结点等参单元的基本单元如图5-11所示。

图5-9 ANSYS 提供的Solid45单元
图5-10 Solid45的基本单元 图5-11 八结点基本单元 六面体八节点等参单元的基本单元如图5-11所示，其形函数为，
)8,...,1()
1)(1)(1(8
1=+++=
i N i i i i ζζηηξξ
（5-33）
其中，i i i ζηξ,,为结点的局部坐标。

如图5-12所示，ANSYS 提供的Solid95单元是六面体二十节点等参单元，每个节点有
代表x 、y 、z 三个方向位移的三个自由度，可以退化为五面体棱柱、五面体金字塔形和四面体单元。

Solid95单元的基本单元如图5-13所示。

图5-12 ANSYS 提供的Solid95单元
图5-13 Solid95的基本单元
图5-14 二十结点基本单元
与六面体八结点等参单元相比，六面体二十结点等参单元能更好地适应不规则的形状，计算误差比较小，基本单元如图5-14所示，其形函数为，
)
8,...,2,1()2)(1)(1)(1(8
1=-+++++=
i N i i i i i i i ζζηηξξζζηηξξ （5-34）
)19,17,11,9()1)(1)(1(412
=++-=i N i i i ζζηηξ )20,18,12,10()1)(1)(1(412
=++-=i N i i i ξξζζη
)16,15,14,13()
1)(1)(1(4
12=++-=
i N i i i ηηξξζ
其中ζηξ,
,i i 为单元结点在局部坐标系中的坐标。

单元刚度矩阵为，
[]dxdydz B D B K T
e
]][[][⎰⎰⎰
=
（5-35）
按照上节介绍的等参单元分析的基本步骤可以得到三维单元的单元刚度矩阵。

雅可比矩阵为， ⎥
⎥⎥⎥⎥
⎥⎥⎦⎤⎢⎢⎢⎢⎢
⎢⎢⎣
⎡∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂=ζζ
ζηηηξξ
ξz y x z y x
z y x J ][ （5-36）
形函数对整体坐标的偏微分可以用雅可比矩阵表示为形函数对局部坐标的偏微分， ⎪⎪⎪⎭
⎪⎪⎪⎬⎫
⎪⎪⎪⎩
⎪⎪⎪⎨⎧∂∂∂∂∂∂=⎪⎪⎪⎭⎪⎪
⎪⎬⎫
⎪⎪⎪⎩⎪⎪⎪⎨⎧∂∂∂∂∂∂-ζηξi i
i i i i
N N N J z
N y N x N 1][ （5-37）
将（5-32）式代入（5-36）式可以计算出雅可比矩阵，
⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡⎥⎥⎥
⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂=⎥⎥⎥⎥⎥⎥⎥
⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂=∑
∑
∑∑
∑
∑∑
∑
∑=========
n n
n
n n n i n
i i
i n
i i
i n i i i
n
i i i n
i i i n
i i i n
i i i n
i i i n i z y x z y x z y x N N N
N N N N N N z N
y N x N z N y N x N z N y N x N J ......
......][2
22111
21
21211
1
1
111
111ζζ
ζ
ηηηξ
ξξζζ
ζ
ζ
ηηξξξ
（5-38）
利用雅可比矩阵的行列式，将整体坐标系下的积分转换为在局部坐标系下的积分。

在整体坐标系中的体积微元为，
dxdydz z d y d x d dV =⨯⋅=)(
微矢量在局部坐标系中表示为，
321 ζζηηξξd x d x d x
x d ∂∂+∂∂+∂∂=
321 ζζ
ηηξξd y d y d y y d ∂∂+∂∂+∂∂=
321
ζζ
ηηξξd z d z d z z d ∂∂+∂∂+∂∂=
其中3,
2,1
为局部坐标系中ςηξ,,
方向上的单位向量。

121
323
ζηηζηξξζζηζ
ηηξξηζξζξηξηξd d z y d d z y d d z y d d z y d d z y d d z y z d y d ∂∂∂∂-∂∂∂∂+∂∂∂∂+∂∂∂∂-∂∂∂∂-∂∂∂∂=⨯
ζηξζηξζ
ζ
ζ
ηηηξξξd d d J d d d z y x z y x z y x dxdydz =∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂=
ζηξd d d J B D B K T
e
]][[][][1
1
1
1
11
⎰
⎰⎰---=
（5-39）
最后，用高斯积分计算出单元刚度矩阵。

同样，用上节中类似的公式就可以在局部坐标下完成单元的载荷移置。

体力移置的公式为， ζηξd d d J t p N R T
e
}{][}{1
1
1
1
1
1
⎰
⎰⎰---=
（5-40）
在1=ξ的面上受到面力作用，面力移置的公式为， ζηξξd d J
t P N R T
e
1
11
1
11
}{][}{==--⎰
⎰=
（5-41）
在点),,(000ζηξ集中力移置的公式为， }{][}{),,(000P N R T
e
ζηξ=
（5-42）。