弹性力学空间问题与体单元

1、6 弹性力学空间问题与体单元6. 1 空间问题的基本描述（1）基本变量 空间问题的有限元法与平面问题的有限元法的原理和求解思路类似，只是离散化是在三维空间进行，基本变量均是三维变量。 空间物体上的一点的位移矢量在空间直角坐标系中记为工程应变有应力列阵有（2） 几何方程 在小变形条件下，描述应变场与位移场关系的几何方程如下简记为L 三维微分算子矩阵。（3）物理方程 对于无初始应力、无初应变的各向同性线弹性体，应力和应变间存在比例关系，即广义胡克定律，记为或若用应变的线性组合表示应力，有引入弹性矩阵D，则其中（4）平衡方程记作用在物体内的体力矢量为对于空间问题，同1. 3节的分析，可导出描述应力与

2、体力之间关系的静力平衡方程简记为上述式整理得到位移为未知函数的平衡方程6. 2 有限元公式一般空间问题有限元离散化形成的单元为三维体单元。将空间问题有限元离散的第 i 个结点位移记为若单元有m个结点，则单元结点位移列阵为注意上式中单元内点号r（1m）与总点号i 不同。单元内任一点的位移由下式定义则单元内应变为单元内应力为 为了建立有限元公式，设单元在结点上受的单元结点力为 空间问题，每个结点有三个分量，其中第i 个结点的结点力为为了得到单元刚度矩阵公式，将势能原理用与一个单元，此时设单元上只受单元结点力，根据式（6.7）、（6.11）、（6.12），单元的总势能为（Ve为单元体积）由最小势能原

3、理可知，真解使取最小值。即 当单元受体力 f ，面力 作用时，由虚功原理得体力、面力的等效结点力分别为由 ，得此时单元结点载荷列阵为6. 3 常应变四面体单元四结点四面体单元是最简单的空间单元，如图 6. 1所示。四个角点为结点，分别以 i、j、k 和 m 表示，各结点的坐标分别结点的编号规定如下： 在右手坐标系中，当右手螺旋按照ij k转向时，拇指指向 m 结点。空间单元中每个结点有3个方向位移分量，记为单元结点位移列阵为6. 3. 1 位移函数 单元变形时，单元内各点有三个位移分量，各个位移分量一般为坐标的函数。 对于四结点四面体单元，其单元内部位移假定为坐标的线性函数 位移函数式（6.

4、22）中的12个待定常数可用单元的结点位移来表示。把各结点的坐标分别代入式（6. 22）得第一式中上式改写为引入式（6.26）中的系数 ai 、bi 等由 A 的伴随矩阵确定。上式中V 是四面体的体积。定义：设是一个 n 阶矩阵称 为的伴随矩阵，其中 Aij 是行列式中元素 aij的代数余子式。第 i 行第 j 列系数ai 、bi 、ci、di代入式（6.22）的第一式即式中 Ni、Nj、Nk、Nm为各结点的形函数。形函数统一写为同理用结点位移和形函数描述的y、z方向的位移分量为把式（6.28a）、（6.30）、（6.31）统一写为矩阵形式简写为I 三阶单位矩阵。由于位移函数是线性的，在相邻单

5、元的接触面上，位移显然是连续的，因此四面体单元是协调元。6. 3. 2 单元应变对于空间问题，一点的六个分量和位移场的关系由下式描述。当给定位移场后，把 代入式（6.1），可得单元内任一点的应变。其中 Be 为单元应变转换矩阵，其子块为63的矩阵 由此可见 Be 中的所有元素均是与坐标值无关的常数，说明四面体单元是常应变单元。 由于单元位移假定为现行函数，且应变仅为与位移的一阶导数，因此单元内任一点的应变是常数。6. 3. 3 单元应力把用结点位移表示的应变代入空间问题的物理方程（6. 2b）中，可得单元内人一点的应力其中Se为单元应力转换矩阵，其子块为63的矩阵 同样，分析Sie可以看出，一

6、旦单元确定， Sie也就确定了，此时单元内的应力仅依赖于结点位移，应力也是常数。因此四结点四面体单元是常应力单元。6. 3. 4 单元刚度矩阵当得到单元应变矩阵Bie 和单元应力矩阵Sie 后，将其代入式（6.17），可计算出单元刚度矩阵Ke，为1212的矩阵由于矩阵Bie和Sie为常量矩阵，有单元刚度矩阵的结点分块形式为令 和 ，上式中各分块为对于具体的载荷，可由 计算出单元结点等效结点载荷，而单元刚度矩阵的组装过程同平面问题中是类似的。6. 4 常用的三维单元常应变四面体单元对边界的适应能力较强，但计算精度低。空间单元还有五面体和六面体单元。6. 4. 1 体积坐标 为了四面体单元形函数构

7、造方便，介绍定义在单元上的自然坐标体积坐标。体积坐标的引入将简化高次四面体单元的分析过程。如图6-2所示，一常应四面体变单元，P(x , y , z)为单元内任意一点。P点与四个角点的连线把原来的四面体分割成四个小四面体。V1四面体Pjkm 的体积；V2四面体Pkmi 的体积；V3四面体Pmij的体积；V4四面体Pijk的体积。四个四面体体积分别记为令式中 V e四面体ijkm 的体积。考虑到以及V2、V3、V4等。对比式（6.29）可以得到则任意一点 P 在单元内的位置可以由 L1L4 来确定，称 L1L4 为体积坐标，由 V1+V2+V3+V4=V e ，有 可见体积坐标仅有3个独立，当P

8、点在四面体某个表面三角形时，体积坐标退化为面坐标。由式（6. 28）、（6. 26）有结合式（6. 25）并引入式（6. 42），得到6. 4. 2 高次四面体单元 常应变四面体中的各点应力为常量，但实际工程中的各点应力是随着坐标变化的。为了反映真实情况，可以假定高次位移函数，相应的单元称为高次单元。虽然高次单元的形函数比较复杂，但当引入体积坐标以后，高次四面体单元的形函数可仿照平面问题中高次单元的构造方法。（1） 10结点二次四面体单元二次四面体单元，如图6. 3所示，除了四个角点为结点，各条棱边的中点也均为结点，共10各结点。其位移函数为 其中 Ni 为各结点的形函数，用体积坐标来描述，有

9、角结点棱边结点这样定义的形函数在结点上仍然满足以及，对于任意一点 二次位移函数对坐标求一次导后，得到线性函数。 因此10结点二次四面体单元的应变随坐标线性变化，称为线形应变四面体单元。 此单元计算精度大大高于4结点常应变四面体单元。（2） 20结点四面体单元 对于三次四面体单元，如图6. 4所示，共有20各结点，包括4个角点、12个棱边3等分点以及4个表面形心。各结点的形函数用体坐标表示为角结点：表面形心结点：棱边 3 等分点： 由于位移函数为三次多项式，20结点四面体单元的应变和应力是坐标的二次函数。6. 4. 3 三棱柱单元 三维求解域为界面形状复杂的柱体，可以选用三棱柱五面体单元来解决界

10、面形状不规则的问题。图6.5为六结点线性三棱柱单元，利用其几何特性，结合面坐标和方向的拉格朗日插值函数，可得各结点的形函数（1） 6结点线性三棱柱单元式中 Li 三角形面积坐标。式中 Li三角形面积坐标。（2）15结点二次三棱柱单元 15结点五面体单元除了角点外，每边的中点也为结点，如图6.6所示，为自然坐标。矩形边中点该单元的三棱柱形函数为角结点三角形边中点6. 4. 4 长方体单元（1） 8结点长方体单元 如图6. 7所示，8结点长方体单元的特点是结点坐标有如下规律单元中心坐标 不同的单元长、宽、高是不同的。 为了将单元规范，引入单元自然坐标在这种坐标下结点坐标值为 单元各结点的形函数由三

11、个方向的两点拉格朗日插值函数构成 类似于二维问题，可以有三个方向的等距三点插值函数构成更复杂的拉格朗日单元，这样的长方体有27个结点。由于结点过多不适用，常见的有20结点与21结点的单元。（2） 20结点长方体单元 如图6. 8所示，20结点长方体单元较之8结点单元增加了12条棱边的中点，即 仍采用式（6.53）的单元自然坐标系。用自然坐标描述的单元结点形函数为i =18结点 长方体单元精度明显高于同样结点数的四面体单元，但长方体形状难以适应工程结构的复杂外形。 通过几何变换，将立方体单元映射为可变形的六面体单元就成为必然采用的方法。6. 4. 5 六面体等参元四结点单元形状灵活，但精度很低。

12、长方体单元精度高，但形状难以全部离散为长方体单元。解决的方法是采用等参元将单元自然坐标系中的立方体单元映射为实际空间中的任意六面体单元。 形函数是在自然坐标中的母单元中定义的，如同上节一样。 将立方体的母单元映射为实际问题的六面体单元是通过如下坐标变换实现的。其中形函数矩阵变换式（6.56）实现了 的一一对应。同时，单元位移场设定为 其中u、v、w表示为、的显函数，应变分量需有位移对x、y、z 求偏导得出：其中涉及形函数的偏导 共24个。形函数Ni 对直角坐标和自然坐标的偏导数间关系为其中雅可比矩阵 对任一点（ ， ），J均可有式（6.60）算出来。由式（6.60），有至此，由式（6.11）所定义的Be 元素均可算出来。采用数值积分，如高斯积分，可计算出单元刚度矩阵其中积分在母单元上进行，且K=2424矩阵。 当单元表面上作用有外载荷，在计算单元等效结点在何时，积分式中会涉及面极微元的变换。如在=常数的面上有