第8章 空间实体单元

1、第第8 8章章 空间实体单元空间实体单元 8.1 概述概述 许多工程实际问题属于空间问题。用有限元法分析空间许多工程实际问题属于空间问题。用有限元法分析空间 问题和分析平面问题在原理、思路和解题方法完全相同，问题和分析平面问题在原理、思路和解题方法完全相同， 基本未知量仍然是节点位移。不同的是单元具有三维特点。基本未知量仍然是节点位移。不同的是单元具有三维特点。 节点位移在节点位移在x、y、z三个坐标轴方向都有分量：三个坐标轴方向都有分量：u、v、w。 它的基本方程比平面问题要多，有它的基本方程比平面问题要多，有3个平衡方程，个平衡方程，6个几何个几何 方程，方程，6个物理方程。个物理方程。分

2、析方法仍然是先进行单元分析，再分析方法仍然是先进行单元分析，再 进行系统分析，最后求解系统的节点平衡方程，解算内力进行系统分析，最后求解系统的节点平衡方程，解算内力 或应力。或应力。 空间离散化后的单元模型主要有：四面体单元、长方空间离散化后的单元模型主要有：四面体单元、长方 体单元、直边六面体单元、曲边六面体单元，如图体单元、直边六面体单元、曲边六面体单元，如图8-1所示。所示。 (a)(b)(c)(d) 图图8-1 空间实体单元模型空间实体单元模型 (a) 图为图为4节点四面体单元，是空间问题最简单的单元，也节点四面体单元，是空间问题最简单的单元，也 是是常应变、常应力单元常应变、常应力单

3、元。类似平面问题三节点三角形单元。类似平面问题三节点三角形单元 进行分析。进行分析。 ？ (b)图图w为长方体单元，可以类似平面四节点矩形单元进行为长方体单元，可以类似平面四节点矩形单元进行 分析。分析。 (c)图为任意八节点六面体单元，可以类似平面四节点任意图为任意八节点六面体单元，可以类似平面四节点任意 四边形等参元进行分析。四边形等参元进行分析。 (d)图为图为20节点曲边六面体单元，可以类似平面八节点曲边节点曲边六面体单元，可以类似平面八节点曲边 四边形等参元进行分析四边形等参元进行分析 。 8.2 4节点四面体常应变单元节点四面体常应变单元 1、位移模式、位移模式 如图如图8-2所示

4、，取四面体的所示，取四面体的4个顶点个顶点 i,j,m,n 为节点。每为节点。每 一个节点有一个节点有3个位移分量，即个位移分量，即 ),(nmjiwvu T iiii （8-1） 单元节点位移向量为单元节点位移向量为 T nnnmmmjjjiii TT n T m T j T i e wvuwvuwvuwvu （8-2） 与平面问题式（与平面问题式（2-12）类似，假定单元内一点的位移）类似，假定单元内一点的位移 分量为坐标的线性函数分量为坐标的线性函数 zayaxaaw zayaxaav zayaxaau 1211109 8765 4321 （8-3） 将式（将式（8-3）的第）的第1式应

5、用于式应用于4个结点，则个结点，则 i j m n x y z 图图 8-2 nnnn mmmm jjjj iiii zayaxaau zayaxaau zayaxaau zayaxaau 4321 4321 4321 4321 （8-4） 由此可解出由此可解出a1a4，再代回到式（，再代回到式（8-3）的第）的第1式，与式（式，与式（2- 19）的第）的第1式类似，有式类似，有 nnmmjjii uNuNuNuNu（8-5） 式中形函数具有与式（式中形函数具有与式（2-18）类似的形式：）类似的形式： ),()( 6 1 nmjizdycxba V N iiiii （8-6） 其中其中 ),

6、(nmji 1 1 1 1 1 1 1 1 1 nn mm jj i nn mm jj i nn mm jj i nnn mmm jjj i yx yx yx d zx zx zx c zy zy zy b zyx zyx zyx a （8-7） nnn mmm jjj iii zyx zyx zyx zyx V 1 1 1 1 （8-8） 在式（在式（8-8）中，）中，V为四面体的体积。为使其计算值不为负，为四面体的体积。为使其计算值不为负， 单元的节点（单元的节点（i,j,m,n）编号次序应遵循右手法则。）编号次序应遵循右手法则。(p4) 采用同样的方法，可得采用同样的方法，可得 nnmm

7、jjii vNvNvNvNv（8-9） nnmmjjii wNwNwNwNw （8-10） 将式（将式（8-5）、（）、（8-9）（）（8-10）统一用矩阵式表示，可得与）统一用矩阵式表示，可得与 平面问题式（平面问题式（2-20）类似的公式）类似的公式 （8-11） e N w v u f 式中式中N为单元形函数矩阵，其维数为为单元形函数矩阵，其维数为312。进一步可写。进一步可写 为与平面问题式（为与平面问题式（2-21）、（）、（2-22）类似的子块形式）类似的子块形式 （8-12） nmji NNNNN 其中，子矩阵其中，子矩阵 （8-13）),( 00 00 00 nmjiIN N

8、N N N i i i i i 式中，式中，I为为3阶单位矩阵。阶单位矩阵。 2、应变矩阵、应变矩阵 在空间问题中，每点有在空间问题中，每点有6个应变分量。几何方程为：个应变分量。几何方程为： T T zxyzxyzyx z u x w y w z v x v y u z w y v x u （8-14） 将式（将式（8-11）（）（8-13）和（）和（8-6）代入上式，得）代入上式，得 （8-15） e nmji e BBBBB 式中式中 （8-16） ),( 0 0 0 00 00 00 6 1 nmji bd cd bc d c b V B ii ii ii i i i i 上式（上式（

9、8-15）、（）、（8-16）与平面问题式（）与平面问题式（2-24）（）（2-26） 类似。与平面问题三节点三角形单元相同，在四节点四面类似。与平面问题三节点三角形单元相同，在四节点四面 体单元中，体单元中，B的元素都是常量，因此是常应变单元。的元素都是常量，因此是常应变单元。 2、应力矩阵、应力矩阵 三维问题的应力应变关系也可写为式（三维问题的应力应变关系也可写为式（2-8）的矩阵）的矩阵 形式形式 D（8-17） 与平面问题不同，这里与平面问题不同，这里和和分别由分别由6个分量组成，个分量组成， 弹性矩阵弹性矩阵D是一个是一个66的矩阵：的矩阵： T zxyzxyzyx （8-18） （

10、8-19） T zxyzxyzyx （8-20） )1 (2 21 00000 )1 (2 21 0000 )1 (2 21 000 1 11 1 1 1 )21)(1 ( )1 ( 称 对 E D 将式（将式（8-20）所表示的）所表示的D和式（和式（8-15）、（）、（8-16）所表）所表 示的示的B代入式（代入式（8-22），并将），并将S定成分块矩阵的形式，有定成分块矩阵的形式，有 将式（将式（8-15）代入式（）代入式（8-17），得），得 e S（8-21） 应力矩阵应力矩阵S为为 （8-22）BDS 由于由于D、B都是常数矩阵，因此应力矩阵都是常数矩阵，因此应力矩阵S也是常也是常

11、 数矩阵。也就是说，单元中的应力分量也是常数。数矩阵。也就是说，单元中的应力分量也是常数。 （8-23） nmji SSSSS 式中式中 （8-24）),( 0 0 06 22 22 22 11 11 11 3 nmji dAdA cAdA bAcA dcAbA dAcbA dAcAb V A BDS ii ii ii iii iii iii ii 其中其中 （8-25） )21)(1 ( )1 ( )1 (2 21 1 321 E AAA 3、单元刚度矩阵、单元刚度矩阵 仿照平面问题中的推导，可得单元平衡方程仿照平面问题中的推导，可得单元平衡方程 ee Fk（8-26） 单元刚度矩阵具有与式

12、（单元刚度矩阵具有与式（2-33）类似的形式）类似的形式 （8-27） v T dVBDBk 式中，式中，k是一个是一个1212的矩阵。由于的矩阵。由于B、D都是常数矩都是常数矩 阵，所以阵，所以k也是一个常量矩阵。并且也是一个常量矩阵。并且 （8-28）VBDBk T 写成分块矩阵的形式，有写成分块矩阵的形式，有 （8-29） nnnmnjni mnmmmjmi jnjmjjji inimijii kkkk kkkk kkkk kkkk k 式中子矩阵式中子矩阵krs为为33的矩阵的矩阵 VBDBk s T rrs ),(nmjisr （8-30） )( )( )( 36 22121 212

13、21 21212 3 srsrsrsrsrsrsr srsrsrsrsrsrsr srsrsrsrsrsrsr ccbbAdddcAcdAdbAbdA cdAdcAbbddAcccbAbcA ddAdbAbcAcbAddccAbb V A 4、等价节点力向量、等价节点力向量 式（式（8-26）中的单元等价节点力也包括体积力、表面力、）中的单元等价节点力也包括体积力、表面力、 集中力几部分。体积力与表面力的计算公式与平面三角形集中力几部分。体积力与表面力的计算公式与平面三角形 单元公式（单元公式（2-36）、（）、（2-37）类似：）类似： V V Te V dVpNF（8-31） （8-32）

14、 S S Te S dSpNF 对于简单情形，也可采用静力等效原则简化计算。对于简单情形，也可采用静力等效原则简化计算。 进一步的整体平衡方程的建立（即结构刚度矩阵、结进一步的整体平衡方程的建立（即结构刚度矩阵、结 构等价节点力列阵的组集）、位移约束条件的引入、线性构等价节点力列阵的组集）、位移约束条件的引入、线性 方程组的求解等，和平面问题有限元法一样，不再赘述。方程组的求解等，和平面问题有限元法一样，不再赘述。 8.3 二十结点六面体等参数单元二十结点六面体等参数单元 由于精度高，容易适应不同边界，在平面问题中常选由于精度高，容易适应不同边界，在平面问题中常选 用了八节点四边形等参数单元。

15、与此类似，在三维问题中用了八节点四边形等参数单元。与此类似，在三维问题中, 常选用二十节点六面体等参数单元。常选用二十节点六面体等参数单元。 如图如图8-3所示，在整体坐标系的二十节点六面体的实际单元与中所示，在整体坐标系的二十节点六面体的实际单元与中 心在局部坐标的原点、边长为心在局部坐标的原点、边长为2的立方体基本单元相对应。的立方体基本单元相对应。 12 3 4 5 6 7 8 9 10 11 12 13 14 15 1617 18 19 20 x y z 6 12 18 1 2 3 4 5 78 9 10 11 13 14 15 16 17 19 20 = -1 = 1 = 1 = -

16、1 = 1 图图 8-3 实际单元实际单元基本单元基本单元 1、位移模式、形函数和坐标变换式、位移模式、形函数和坐标变换式 = -1 20 1 20 1 20 1ii iiiii i i wNwvNvuNu （8-33） （8-34） 20 1 20 1 20 1i i i i i iiii zNzyNyxNx 形函数的表达式如下：形函数的表达式如下： （8-35） )20,19,18,17(4/ )1)(1)(1 ( )16,15,14,13(4/ )1)(1)(1 ( )12,11,10, 9(4/ )1)(1)(1 ( )8, 2 , 1(8/ )2)(1)(1)(1 ( 00 2 00

17、 2 00 2 000000 iN iN iN iN i i i i 式中式中 iii 000 , 位移模式和坐标变换式可写为如下形式：位移模式和坐标变换式可写为如下形式： 根据几何方程，单元中的应变为根据几何方程，单元中的应变为 2、应变矩阵、应变矩阵 ee BBBB 2021 （8-36） 其中其中 （8-37） Te wvuwvuwvu 202020222111 应变矩阵的子矩阵应变矩阵的子矩阵Bi为：为： （8-38）)20, 2 , 1( 0 0 0 00 00 00 i x N z N y N z N x N y N z N y N x N B ii ii ii i i i i 根

18、据复合函数求导规则，有根据复合函数求导规则，有 i i i i i i N N N J z N y N x N 1 （8-39） 式中式中J-1为雅可比矩阵为雅可比矩阵J的逆矩阵。的逆矩阵。 J的表达式为的表达式为 202020 222 111 2021 2021 2021 zyx zyx zyx NNN NNN NNN zyx zyx zyx J （8-40） 3、应力矩阵、应力矩阵 单元中的应力为单元中的应力为 e BDD 单元应力矩阵单元应力矩阵S为为 （8-41） BDS 4、单元刚度矩阵、单元刚度矩阵 单元刚度矩阵可写成单元刚度矩阵可写成 （8-42） 1 1 1 1 1 1 | d

19、ddJBDBk T k是一个是一个6060的矩阵，式（的矩阵，式（8-42）通常采用高斯法进行）通常采用高斯法进行 积分。积分。 5、等价节点力矩阵、等价节点力矩阵 等价节点力计算公式如下：等价节点力计算公式如下： （1）体积力）体积力 设单位体积力是设单位体积力是 ， 则等价节点力为则等价节点力为 T vzVyVxV pppp dddJ p p p N F F F F Vz Vy Vx i Vzi Vyi Vxi Vi | 1 1 1 1 1 1 （8-43） （2）表面力）表面力 设某边界面上作用表面力设某边界面上作用表面力 ，则等价节点，则等价节点 力为力为 T SzSySxS pppp

20、 dS P P P N F F F F sz Sy Sx i Szi Syi Sxi Si （8-44） 设该边界面对应基本单元设该边界面对应基本单元 =1 的面。由数学公式，结构的面。由数学公式，结构 坐标下曲面微元坐标下曲面微元dS对应于单元坐标下的微元面积式为对应于单元坐标下的微元面积式为 （8-45） dd yxyxxzxzzyzy dS 222 )()()( 则式（则式（8-44）可写为单元坐标系下的积分公式）可写为单元坐标系下的积分公式 （8-46） 1 1 1 1 sx sy sx isi p p p NF dd yxyxxzxzzyzy 222 )()() ( 以上为以上为 =

21、1的表面力计算公式。对于其他表面力，可类的表面力计算公式。对于其他表面力，可类 似处理。式（似处理。式（8-46）通常也采用高斯法进行积分计算。）通常也采用高斯法进行积分计算。 8.4 空间轴对称单元空间轴对称单元 许多工程构件，其几何形状、约束条件及所受的荷载许多工程构件，其几何形状、约束条件及所受的荷载 都对称于某一轴，因而所有的位移、应变和应力分量也都都对称于某一轴，因而所有的位移、应变和应力分量也都 对称于该轴。这类问题称对称于该轴。这类问题称空间轴对称问题空间轴对称问题。 对空间轴对称问题，采用圆柱坐标系，对空间轴对称问题，采用圆柱坐标系，r表示径向坐表示径向坐 标，标，z表示轴向坐

22、标，任一对称面为表示轴向坐标，任一对称面为rz面。在有限元分析时，面。在有限元分析时， 采用绕对称轴旋转一周的采用绕对称轴旋转一周的轴对称环形单元轴对称环形单元。把轴对称的工。把轴对称的工 程构件与程构件与rz坐标平面正交的截面划分为一些三角形，例如坐标平面正交的截面划分为一些三角形，例如 其中一个是其中一个是ijm（图（图8-4）。由这些三角形旋转一周就构成）。由这些三角形旋转一周就构成 上述环形单元。上述环形单元。 图图8-4 轴对称构件及三角形环形单元轴对称构件及三角形环形单元 r (u) z (w) i j 当然，也可以划分为当然，也可以划分为矩形矩形，构成矩形截面的环形单元。总，构成

23、矩形截面的环形单元。总 之，单元截面形状，可采用之，单元截面形状，可采用三角形、四边形等平面有限元三角形、四边形等平面有限元 法所用单元形状法所用单元形状。本节只讨论三角形截面的环形单元，其。本节只讨论三角形截面的环形单元，其 它截面形状环形单元可参照本节方法进行分析。它截面形状环形单元可参照本节方法进行分析。 各个环形单元（后简称单元）间用位于三角形顶点的各个环形单元（后简称单元）间用位于三角形顶点的 环形铰联系起来，就成了离散结构模型。我们把环形铰在环形铰联系起来，就成了离散结构模型。我们把环形铰在 rz平面的交点（即三角形顶点）称为单元的节点（如平面的交点（即三角形顶点）称为单元的节点（

24、如 i,j,m ）。因此）。因此对轴对称问题进行有限元分析只需在对轴对称问题进行有限元分析只需在rz平面平面 划分网格划分网格，就像平面问题在，就像平面问题在xy平面上划分网格一样。照此平面上划分网格一样。照此 思路，轴对称空间问题被大大简化。思路，轴对称空间问题被大大简化。 1、位移模式、位移模式 对于如图对于如图8-4所示轴对称三角形环形单元，是由如图所示轴对称三角形环形单元，是由如图8- 5所示平面上的三角形绕对称轴轴回旋一周得到的。考虑到所示平面上的三角形绕对称轴轴回旋一周得到的。考虑到 问题的轴对称特点，可以借助于平面问题三节点三角形单问题的轴对称特点，可以借助于平面问题三节点三角形

25、单 元位移模式进行单元分析。元位移模式进行单元分析。 图图8-5 r z i j m ri rj rm 由于轴对称，由于轴对称，环向位移恒等于零环向位移恒等于零。 只有径（只有径（r）向位移和轴（）向位移和轴（z）向位移，）向位移， 它们恰在它们恰在rz坐标平面上。设径向位移为坐标平面上。设径向位移为 u，轴向位移为，轴向位移为w。对于图。对于图8-5所示情形，所示情形， 借助平面问题的三角形单元，取位移模借助平面问题的三角形单元，取位移模 式为式为 ？ zrw zru 654 321 （8-47） 代入节点位移后，可解出代入节点位移后，可解出a1a6，再代入上式，得，再代入上式，得 （8-4

26、8） mmjjii mmjjii wNwNwNw uNuNuNu 其中，形函数其中，形函数 （8-49）),()( 2 1 mjizcrba A N iiii 式中式中A, ai, bi, ci，与平面问题三角形单元的对应公式（，与平面问题三角形单元的对应公式（2- 15）、（）、（2-17）一致，区别仅仅是将那里的）一致，区别仅仅是将那里的x, y换成这里的换成这里的 r, z。 式（式（8-48）也可写为式（）也可写为式（2-20）同样的形式）同样的形式 e mji e NNNNf（8-50） 式中，式中， Ni= NiI(i, j, m)，其中，其中I为为2阶单位矩阵。阶单位矩阵。 2、

27、应变矩阵、应变矩阵 根据弹性力学理论，空间轴对称问题的几何方程为根据弹性力学理论，空间轴对称问题的几何方程为 r w z u z w r u r u rz z r （8-51） 与平面问题中与平面问题中只含有只含有3个应变分量不同，这里个应变分量不同，这里含有含有4 个应变分量个应变分量。将。将u, w的表达代入式（的表达代入式（8-51），得），得 e mji e BBBB（8-52） 式中式中 ),( 0 0 0 2 1 mji bc c f b A B ii i i i i （8-53） 其中其中 ),(mji r zc b r a f i i i i （8-54） 由式（由式（8-52

28、）（）（8-54）可见，矩阵中含有变量）可见，矩阵中含有变量r, z，因此，因此 它不是常数矩阵。即它不是常数矩阵。即轴对称问题的三角形环形单元不是常轴对称问题的三角形环形单元不是常 应变单元应变单元。 3、应力矩阵、应力矩阵 根据弹性力学理论，空间轴对称问题的应力根据弹性力学理论，空间轴对称问题的应力-应变关系为应变关系为 D rz z r （8-55） 式中式中D是轴对称问题的弹性矩阵是轴对称问题的弹性矩阵 1 21 000 1 11 1 1 1 )21)(1 ( )1 ( 称 对 E D（8-56） 将式（将式（8-52）代入式（）代入式（8-55），得），得 ee SBD（8-57）

29、应力矩阵应力矩阵 S=DB 显然，除显然，除 rz外，单元中其他应力分量不是常数。外，单元中其他应力分量不是常数。 4、单元刚度矩阵、单元刚度矩阵 轴对称问题的单元刚度矩阵可由式（轴对称问题的单元刚度矩阵可由式（8-27）计算）计算 V T V T dzrdrdBDBdVBDBk 由于被积函数由于被积函数 与无关，故在三角形截面的环单元的积分可与无关，故在三角形截面的环单元的积分可 简化为在三角形截面上的积分：简化为在三角形截面上的积分： A T rdrdzBDBk2（8-58） 式（式（8-58）中，）中，B不是常数，而由式（不是常数，而由式（8-53）、（）、（8- 54）确定。式（）确定

30、。式（8-53）、（）、（8-54）中，）中，ai, bi, ci是（是（P32） （1）一般公式）一般公式 （2）近似计算）近似计算 可见，单元刚度矩阵计算比平面三角形单元麻烦。须对其进行数值可见，单元刚度矩阵计算比平面三角形单元麻烦。须对其进行数值 计算或近似计算。现在，讨论面积分计算或近似计算。现在，讨论面积分A g(r, z)drdz的计算问题。的计算问题。 虽为常数，由式（虽为常数，由式（2-17）计算，但）计算，但 fi 是是r, z的函数。式（的函数。式（8- 58）被积函数）被积函数BTDBr是是r, z的函数，暂简写为的函数，暂简写为g(r,z)。式。式 （8-58）可写为）

31、可写为 图图8-6 r z i j m 对于一个典型三角形截面对于一个典型三角形截面ijm（图（图8- 6），过），过j点作垂直于点作垂直于r轴的直线，将轴的直线，将 三角形三角形ijm，分成两个三角形：，分成两个三角形：A1， A2。于是，有。于是，有 A2 A1 A drdzzrgk),(2 12 ),(),(),( AAA drdzzrgdrdzzrgdrdzzrg d e rz rz e f rz rz dzzrgdrdzzrgdr )( )( )( )( 2 1 3 1 ),(),(（8-59） 完成上述积分可采用数值方法，有专完成上述积分可采用数值方法，有专 门的二重数值积分程序可

32、供引用。门的二重数值积分程序可供引用。 也可采用近似方法完成：也可采用近似方法完成： 当单元较小时，常把各个单元中的当单元较小时，常把各个单元中的r, z近似看作常数，并且近似看作常数，并且 分别等于各单元形心的坐标，即分别等于各单元形心的坐标，即 )( 3 1 ),( 3 1 mjimji zzzzzrrrrr i j m A2 A1 r z d e f Z1(r) Z2(r) Z3(r) 则式（则式（8-54）成为常数：）成为常数： r zc b r a ff i i i ii （8-60） 这样，就可把各个单元近似地当做常应变单元，式（这样，就可把各个单元近似地当做常应变单元，式（8-5

33、8） 变成变成 2BDBArk T （8-61） 式中式中B是从式（是从式（8-52） （8-54）将）将B中的中的r, z用用r, z代替后得出的。代替后得出的。 若将单元刚度矩阵若将单元刚度矩阵k写为子块形式写为子块形式 mmmjmi jmjjji imijii kkk kkk kkk k 其中，子矩阵可近似为其中，子矩阵可近似为 ArBDBk srrs 2（8-62） 写成显式写成显式 srsrsrssr srrrssrsrsrsrsr rs bbAcccbAfbcA bcAfbcAccAbffbAffbb A Ar k 221 2121 3 )( )()( 2 ),(mjisr （8-63） 式中式中A1, A2, A3 由式（由式（8-25）确定。）确定。 5、等价节点力、等价节点力 轴对称问题的等价节点力可用类似平面问题的式（轴对称问题的等价节点力可用类似平面问题的式（2- 36）、（）、（2-37）写出。对于作用于三角形环单元上的体积）写出。对于作用于三角形环单元上的体积 力、表面力的等效结点力为：力、表面力的等效结点力为： rdrdzpNF A V Te V 2（8-64） （1）体积力的等价节点力）体积力的等价节