1. 什么是体积锁死现象
图1 平面应变模型
在图1所示的又两个三角形单元组成的平面应变问题中，如果变形体材料是不可压缩的，那么我们不管在加载点施加多大的力，从有限元法计算得到的所有节点的所有位移都为零。因为任意位移都会导致两个三角形中的一个的体积变化。这种有限元单元精度不足以表现所需体积应变的问题就是体积锁死现象。
图2 受内压的线弹性厚壁圆筒；理论解（点线）和有限元计算结果的比较（1）
图2给出了另外一个例子。这里显示了受内压的厚壁圆筒在弹性圆筒的泊松比为0.3时，有限元计算结果逼近理论解；而当泊松比逼近0.5时，有限元计算结果与理论解相差极大（实际上当泊松比逼近0.5时，有限元计算位移逼近于零）。
1.1 为什么会发生锁死现象
这里给出一个非常粗糙的解释: 我们通常可以把变形体的变形功W分解为几个部分，比如说体积变形部分和非体积变形部分
W = W1 + W2
由于数值计算的误差，当不同部分变形功的差相差很大时，将会带来很大的舍入误差，使计算结果变得很不可靠。具体来说，由于
1） 单元尺寸带来的误差: 壳，梁单元是其典型。由于壳单元面内的尺寸和壳厚相差很大，使得面外应变与面内应变的计算精度相差很大，带来剪切锁死现象。由于细分后网格的单元边长和壳厚将会减小, 剪切锁死现象可以通过细分网格来减轻。
2）由于材质带来的误差: 不可或近似不可压缩材料如泊松比趋于0.5的弹性材料，塑性材料，这时的体积应变对应于巨大的变形能，因而带来体积锁死。体积锁死现象不能靠细分网格来回避。
3）由于物理现象不同带来的误差。比如说在做热-变形耦合分析时，温度（在三角形单元中线性分布）和应变（在三角形单元中均匀分布）计算精度带来的误差。
对于该问题的进一步解释建议参看文献2），在哪里给出了通俗和准确的解释。
1.2 如何减轻锁死
可以通过单元技术来减轻锁死的影响。如Selective reduce integeration, ·Incompatible，Assumed strain，Mixed formulation，MITC等等。对于近似不可压缩材料的体积锁死，最常见的对应方法是采用B-bar单元，这种单元的新式变种是F-bar单元，这是本文下面要讲的内容。
2. B-bar单元3,4)
一般的，单元内任意一点的应变$\varepsilon$表为单元节点i的位移·ui·的函数$\varepsilon=\Sigma_iB_iu_i$. 在B-bar单元中将应变分解$\varepsilon=\varepsilon_s+\bar{\varepsilon}_v$。其中$\bar{\varepsilon}_v$为单元的平均体积应变$\bar{\varepsilon}_v=\sum_i\varepsilon_{vi}$或单元中心的体积应变值$\bar{\varepsilon}_v=\varepsilon_{vo}$. $\varepsilon_s$为偏差应变值。相应地，单元内任意一点的应变用下式计算
$\varepsilon=\sum_i{(B_{si}+\bar{B}_{vi})u_i}$
这里将体积应变做平均计算的方法缓和了
体积锁死。这种方法实现起来非常简单，计算费用不高，效果也不错，因此得到非常广泛的应用。但是这种方法视乎不能满足LBB条件5)。
3. F-bar单元6,7)
相应于B-bar单元对应变的加法分解，单元中的任意一点的变形梯度（deformation gradient）可以分解为体积变形梯度和偏差
变形梯度两个部分: $F=F_sF_v$。 理论上来说这种分解方式是严密的，而B-bar单元中的应变加法分解是一种近似，只是在应变很小时与变形梯度的乘法分解相近。因此把B-bar单元推广到 F-bar单元是非常自然的，在变形较大时采用F-bar单元精度更好。
在F-bar单元中，单元中的任意一点的变形梯度表述为$\bar{F}=(\bar{J}/J)^{1/3}F, J=detF$. 其中$\bar{J}$可以是单元中Jocobian J的平均值$\bar{J}=\frac{1}{V}\int_{V}JdV$，也可以指单元中心的
Jocobian值。
3.1 Total Lagrange算式展开
在Total Lagrange计算中，单元中变形能由Green应变E和第2Piola-Kirchhoff应力S的乘积在初始构型V中的积分计算得出$W=\int_{V}ESdV$, 其中$E=\frac{1}{2}(\bar{F}^T\bar{F}-I)$中的$\bar{F}$采用上述变形梯度表达式。如此，变形能变分为$\int_{V}\delta ESdV$, 单元的consistency stiff matrix可以由
变形能变分的微分$d(\int_{V}\delta ES)dV$展开得出。
3.2 Updated Lagrange算式展开
在Updated Lagrange计算中，单元中变形能由Almansi应变A和Cauchy应力T的乘积在现在构型v中的积分计算得出。$W=\int_{v}ATdv$, 其中$\bar{F}^TA\bar{F}=E$, $T=\frac{1}{J}\bar{F}S\bar{F}$。把这些变换带入3.1的相应算式即可得出Updated Lagrange法的内力，单元的consistency stiff matrix计算公式。
3.3 微小变形计算的
算式展开
如前所述，F-bar单元用于大变形计算。在微小变形计算中应用B-bar单元就可以了。不过由于应变的加法分解可以理解为变形梯度乘法分解的近似。那么可以从F-bar单元公式推出B-bar公式来。虽然从实用的角度来看这样做似乎意义不大。参见文献（9）
3.4 Fbar-patch 单元（8）
Fbar单元不能应用于三角形一次和四面体一次单元，因为这些单元内的应变本身就是均匀分布的。为了将这一方法扩展到三角形一次和四面体一次单元，文献8）提出了将单元分解为相邻的patch的方法，称为Fbar-patch单元。
4. 算例: Cook's membrane （10,11)
Cook's membrane benchmark见上图，锲形版材质设定为近似不可压缩的Mooney Rivlin超弹性材: 其·弹性能函数$W=C_{10}(\bar{I}_1-3)+C_{01}(\bar{I}_2+3)+\frac{1}{D_1}(J-1)^2$, with $C_{10}=1.0, C_{01}=0.925, {D_1}=0.000245$. 采用单元为六面体F-bar单元,B-bar单元 和一般的等参六面体单元。计算结果如下
（a）等参单元:受力方向最大位移值=0.2706
（b）B-bar单元:受力方向最大位移值=3.898
（b）F-bar单元:受力方向最大位移值=4.918
参考文献
1）Allan F. Bower: Applied Mechanics of Solids, CRC, 2010, p535
2） Gangan Prathap: FINITE ELEMENT ANALYSIS AS COMPUTATION， http://www.cmmacs.ernet.in/cmmacs/pdf/ch08.pdf
3）TJR Hughes: Ceneralization of selective integration procedure to anisotropic and nonlinear media, Intl. J. Numer. Methiods Engng., 15(1980), 1413-1418
4) JC Simo, KS Pister, RA Taylor: Variational and projection methods for the volume constraint in finite deformation elatoplasticity, Comp. Mthods Appl. Mech. Eng., 51(1985), 177-208
5) F Brezzi, M Fortin: Mixed and Hybrid Finite Element Method, Springer, 1991
6) B Moran, M Oritz, CF Shin: Formulation of implicit finite-element methods for multiplicative finite deformation platicity, Intl. J. Numer. Methiods Engng., 29(1990), 483-514
7) EA de Souza Neto, D Peric, M Dutko, DRJ Owen: F-bar based linear triangles and tetraheral for finite element strain analysis of nearly incompressible solids. Part I: Formulation and benchmarking Intl. J. Numer. Methiods Engng., 33(1996), 3277-3296
8) EA de Souza Neto, D Peric, M Dutko, DRJ Owen: Design of simple low order finite element for large strain analysis of nearly incompressible solids, Intl. J. Numer. Methiods Engng., 62(2005), 353-383
9）CD Foster, T Mohammad Nejad,;
Trilinear Hexahedra with integral-averaged volumes for nearly incompressible nonlinear deformation, Engineering, 7(2015), 765-788
10) Cook, R. D.; Improved two-dimensional finite element; ASCE J. Struct. Div., ST9, 1851-1863 (1974).
11) Brink, U., and E. Stein,
“On some Mixed Finite Element Methods for Incompressible and Nearly Incompressible Finite Elasticity,
” Computational Mechanics, vol. 19, pp. 105–119, 1996.
