结构有限元计算的奥秘：Locking与F-Bar单元解析

在超弹材料的结构有限元仿真中，单元的锁死现象（Locking）时常困扰着工程与研发人员，是什么样的机理导致了这种数值失真，我们如何来通过科学的方法解决，从而得到正确的计算结果，是我们一直在探索的。而F-Bar数值方法似乎就可以很好的解决这个问题，今天我们来了解一下什么是F-Bar单元，它是如何完美解决超弹大变形的锁死问题的。

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 integration, Incompatible，Assumed strain，Mixed formulation，MITC等等。对于近似不可压缩材料的体积锁死，最常见的对应方法是采用B-bar单元，这种单元的新式变种是F-bar单元，这是本文下面要讲的内容。

2. B-bar单元[3,4]

一般的，单元内任意一点的应变ε表为单元节点i的位移u的函数。在B-bar单元中，将应变分解，其中为单元的平均体积应变或单元中心的体积应变值为偏差应变值。相应地，单元内任意一点的应变用下式计算

这里将体积应变做平均计算的方法缓和了体积锁死。这种方法实现起来非常简单，计算费用不高，效果也不错，因此得到非常广泛的应用。但是这种方法视乎不能满足LBB条件[5]。

3. F-bar单元[6,7]

相应于B-bar单元对应变的加法分解，单元中的任意一点的变形梯度（deformation gradient）可以分解为体积变形梯度和偏差变形梯度两个部分 F=Fs Fv。 理论上来说这种分解方式是严密的，而B-bar单元中的应变加法分解是一种近似，只是在应变很小时与变形梯度的乘法分解相近。因此把B-bar单元推广到 F-bar单元是非常自然的，在变形较大时采用F-bar单元精度更好。

在F-bar单元中，单元中的任意一点的变形梯度表述为，J=detF. 其中

可以是单元中Jocobian J的平均值 ，也可以指单元中心的Jocobian值。

3.1  Total Lagrange算式展开

在Total Lagrange计算中，单元中变形能由Green应变E和第2 Piola-Kirchhoff 应力S的乘积在初始构型V中的积分计算得出，其中中的采用上述变形梯度表达式。如此，变形能变分为，单元的consistency stiff matrix可以由变形能变分的微分展开得出。

3.2 Updated Lagrange算式展开

在Updated Lagrange计算中，单元中变形能由Almansi应变A和Cauchy应力T的乘积在现在构型v中的积分计算得出。

, 其中。把这些变换带入3.1的相应算式即可得出Updated Lagrange法的内力，单元的consistency stiff matrix计算公式。

3.3 微小变形计算的算式展开

如前所述，F-bar单元用于大变形计算。在微小变形计算中应用B-bar单元就可以了。不过由于应变的加法分解可以理解为变形梯度乘法分解的近似。那么可以从F-bar单元公式推出B-bar公式来。虽然从实用的角度来看这样做似乎意义不大。参见文献[9]。

3.4 F-Bar-patch 单元[8]

F-Bar单元不能应用于三角形一次和四面体一次单元，因为这些单元内的应变本身就是均匀分布的。为了将这一方法扩展到三角形一次和四面体一次单元，文献[8]提出了将单元分解为相邻的patch的方法，称为F-Bar-patch单元。

4. 算例: Cook's membrane [10,11]

Cook's membrane benchmark见上图，锲形版材质设定为近似不可压缩的Mooney Rivlin超弹性材: 其·弹性能函数, 给定参数C10=1.0，C01=0.925，D1=0.000245。 采用单元为六面体F-Bar单元，B-Bar单元和一般的等参六面体单元。计算结果如下：

（a）等参单元：受力方向最大位移值=0.2706

（b）B-bar单元：受力方向最大位移值=3.898

（c）F-Bar单元：受力方向最大位移值=4.918

通过数值的实例计算，可以明显的看出，B-Bar单元在超弹材料的计算上要明显优于等参数单元，而F-Bar单元，在变形较大时是优于B-Bar单元的。值得一提的，由中国人主导研发的大型通用有限元软件 WELSIM 已经支持了F-Bar单元，这为我们在工程分析中得到更加精确的结果提供了有力工具。

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删