FEM10-位移/压力混合插值单元

简介

在“FEM08-广义坐标有限元”中，介绍了将一次函数作为试函数，采用标准形式的基于纯位移变量的单元（standard pure displacement-based element）进行有限元计算的方法。对于不可压缩材料，存在体积自锁问题(volumetric locking)。位移/压力混合插值单元，增加了压力插值，可以有效解决体积自锁问题。

体积自锁

体积应变ε_v(volumetric strain)是指材料在外界载荷作用下产生的体积改变量与原体积的比值。

ε_v = ΔV/V; where, ΔV = Change in volume, and V = Original volume.

如图所示的3D立方体微元，在三个坐标轴方向的正应变分别为ε_xx、ε_yy和ε_zz。根据体积应变的定义，将ΔV/V展开，并略去二阶以上的无穷小量，可以得到：

ε_v = ε_xx+ε_yy+ε_zz (体积应变等于三个正交方向的正应变之和)。

下面通过一个例子，来说明基于纯位移变量的单元在处理不可压缩材料时，表现出来的体积自锁现象。

如图所示的不可压缩平面应变单元，2、3、4三个节点固定，给1号节点施加强制位移。为了保持体积不变，施加的位移应满足：v1=-u1

根据“FEM08-广义坐标有限元”中介绍的方法，写出单元的位移函数。由于2、3、4三个节点固定，这些节点的位移为0，所以位移函数只与1号节点的位移(u1,v1)有关。再根据位移函数，推导出正应变，计算体积应变。

对于不可压缩材料，受到载荷变形后体积保持不变，体积应变等于零。

ε_v = 0； 得到：u1=v1=0

以上分析采用的单元是标准形式的基于纯位移变量的单元，最终得到的结论表明：只有当施加的位移为零时，才能维持体积不变，保持不可压缩性。换句话说：如果要移动1号节点，需要施加无穷大的力。单元刚度非常大，表现出锁定状态，这种现象就是体积自锁。

对于各向同性材料，如果泊松比等于0.5，就是完全不可压缩状态。泊松比越接近0.5，材料越接近不可压缩状态。针对前面的平面应变单元的例子，在Adina中采用基于纯位移变量的单元进行分析。给1号节点施加位移v1=-u1=0.1，杨氏模量E=100；泊松比分别为0.49和0.4999999，对应的支反力分别为22.9和1.96E6. 由此可见：材料趋近于不可压缩状态时，基于纯位移变量的单元的刚度趋近于无穷大，出现体积自锁，验证了前面的结论。

位移/压力混合插值单元

想象一下，一个不可压缩的球体，表面受到均匀的压力P。由于不可压缩性，不管P多大，都不会有任何变形，球的任何位置的位移都是零。基于纯位移变量的单元，只把位移作为变量，应变、应力和压力都是位移的导出量。因此，对于不可压缩材料，基于纯位移变量的单元，无法区分位移都是零而压力不同的状态，也就无法准确计算压力。

对于近似不可压缩材料，基于纯位移变量的单元很难准确计算压力，其困难程度取决于材料接近不可压缩的程度。细化网格可以得到比较精确的结果，但所需的网格数量，通常要比处理同类问题的可压缩材料所需的网格数量大得多。

位移/压力混合插值单元(Displacement/Pressure mixed interpolation)，是在基于纯位移变量的单元的基础上，增加了压力P作为变量。在引入压力P作为变量之前，首先需要了解这个压力(Pressure)到底是指什么？

Pressure is simply the negative of hydrostatic stress. 压力其实就是静水应力的负值。我们知道，拉伸状态的应力为正，压缩状态的应力为负，而在描述压力/压强的时候，虽然是压缩的状态，但工程应用中一般不写负号。比如，标准大气压是101kPa，而不写成-101kPa。所以，为了既符合工程应用的习惯，又不破坏数学的严谨性，这里规定压力是静水应力的负值。

解决了正负号的问题，再来解释什么是"静水应力"。Any state of stress can be decomposed into a hydrostatic (or mean)stress σ_m and a deviatoric stress s. 任何应力状态可以分解为静水应力（即平均应力）σ_m和偏应力s。

压力 P 是静水应力 σ_m 的负值，即 σ_m = -P；同时，压力 P 和偏应力 S，可以分别用体积应变 ε_v 和偏应变 ε' 来表示。

其中，K是体积模量，G是剪切模量。

回顾一下，“FEM05-虚位移原理”这节内容中，应力τ没有分解，内力的虚功只有一项。现在我们把应力τ拆分为压力 P 和偏应力 S，对应的虚功也相应地拆分为两项。这样就得到了同时带有位移和压力插值的虚位移原理的表达式。

再根据位移插值函数和压力插值函数，推导出单元的刚度矩阵。这类单元，称为"位移/压力混合插值单元"。

现在，重新采用位移/压力混合插值单元，对前面那个不可压缩平面应变单元的例子进行求解。

推导出应变-位移矩阵。

假设单元内部压力是一个常数，即：H_p=1, P=p0. 当然，也可以采用其他假设给出内部压力变化的插值函数。不同的压力插值函数，产生不同的单元类型。

杨氏模量E=100, 泊松比ν=0.49, 剪切模量G=E/[2(1+ν)]=33.56, 体积模量K=1666.67.

根据前面推导的应变-位移矩阵，并在单元内对dx和dy双重积分，计算单元刚度矩阵。

最后，根据KU=R，将各个节点的位移代入，就可以计算出每个节点受到的外力。节点位移：v1=-u1=0.1，其他都是0.

在Excel中就可以完成以上矩阵的乘法。计算结果显示：1号节点在u和v两个方向受到的外力分别是-2.3306和2.3306，则该节点的合力是：3.296. 比前面Adina中采用基于纯位移变量的单元得到结果22.94小了不少。说明采用位移/压力混合插值方法，确实解决了体积自锁造成刚度偏大的问题。现在，将Adina仿真模型中的单元改为位移/压力混合插值单元，得到的节点受力情况如下图。泊松比为0.49时，1号节点外力是3.296，与前面理论分析Excel计算结果一致；泊松比为0.4999999时，1号节点外力是3.274. 说明泊松比趋近于0.5时，计算结果趋于稳定，没有出现刚度无穷大的情况。因此，位移/压力混合插值单元，对于不可压缩或近似不可压缩的情况，是适用的。

Reference:

• Finite Element Procedures by Klaus-Jürgen Bathe 轩建平 译
• A First Course in Finite Elements by Jacob Fish & Ted Belytschko