ABSTRACT

The locking phenomenon will appear when the commonly used finite elements are applied to the solution of nearly incompressible problems in three dimensions. It is necessary to use some special methods. The penalty function conforming finite element method is an effective method to overcome this locking phenomenon since it is simple for the realization of the resulting program and easy to determine the penalty number and it also does not change the functional stationary value properties. In this paper, the computing format of penalty function finite element method is carefully derived, the conditions for success of the resulting method is analyzed and the effectiveness and robustness of this method are finally verified by some numerical experiments for nearly incompressible elasticity problems. The quality of the mesh used in three-dimensional finite element analysis has a great effect on the accuracy and computational efficiency. If the isotropic grids can be used in the practical calculations, the method will have better convergence.

Keywords:Nearly Incompressible Problems, Locking Phenomenon, Penalty Function Finite Element, Mesh Quality, Reduced Integration

三维近不可压缩弹性问题的罚函数有限元分析

肖映雄^*，周  磊

湘潭大学土木工程与力学学院，湘潭

Email: ^*xyx610xyx@xtu.edu.cn, 455611862@qq.com

收稿日期：2014年10月13日；修回日期：2014年11月7日；录用日期：2014年11月11日

摘  要

对三维近不可压缩弹性问题，利用常规有限元进行求解时会出现体积闭锁现象，需要采用某些特殊的方法。罚函数协调有限元法具有程序实现简单、罚数易于确定以及不改变泛函驻值性质等特点，是克服体积闭锁现象的一种有效方法。本文，针对混合边界条件的三维近不可压缩问题，详细推导了罚函数有限元法的计算格式，分析该方法实施成功的条件，并通过数值实验验证了该方法对解决体积闭锁现象的有效性和鲁棒性。在三维有限元分析中，剖分网格的质量将对计算精度和求解效率产生很大影响，实际计算时若能采用各向同性网格，则对问题的分析将具有更好的收敛性。

关键词

近不可压缩问题，闭锁现象，罚函数有限元，网格质量，减缩积分

1. 引言

对弹性近不可压缩问题，利用常规有限元进行分析时会出现体积闭锁现象，求解所得位移值偏小，甚至趋向于0。解决体积闭锁现象的有效方法很多，它们按照变分原理可以大致分为两类，一类是基于H-R变分原理的混合有限元方法[1] -[3] ，另一类是基于能量极小化原理的标准有限元方法，如非协调元法[4] -[7] ，高阶协调元法[8] [9] ，罚函数法[10] ，减缩积分法[11] [12] 等。这两类方法中，基于能量泛函极小的离散变分问题更容易被解决，并且有限元法形成的代数方程组的系数矩阵是正定的，将给数值求解带来很大方便。解决体积闭锁现象一般采用非协调元法，但在实际应用时，非协调有限元尽管具有自由度少的优点，但该方法对网格十分敏感。若采用高阶单元逼近原问题，则需要更多的计算机存储单元，具有更高的计算复杂性，对三维情形尤为突出。减缩积分会使单元刚度阵秩降低，常引起多余的零能模式。罚函数有限元法是解决体积闭锁现象的一种有效方法，它具有程序实现简单、罚数易于确定及不改变泛函驻值性质等特点，但对三维情形还未见相关系统的研究工作。

本文，针对混合边界条件的三维近不可压缩弹性问题，详细推导了罚函数有限元法的计算格式，分析了该方法实施成功的条件，并通过数值实验验证了罚函数有限元法是解决体积闭锁现象的有效方法。同时，我们还系统研究了有限元网格质量对计算精度的影响，数值结果表明：当采用基于曲面CVT网格的Delaunay四面体网格时，罚函数有限元解具有更好的收敛性，求解所得位移值随单元数的增加而稳定地收敛于理论解，具有更好的鲁棒性。

2. 模型问题及体积闭锁现象

考虑如下混合边界条件的三维弹性问题：

(1)

其中，，且及，为位移向量，为外力，为边界上的面力，为边界的单位外法线向量，拉梅常数和可用杨氏模量和泊松比表示为

(2)

这里，泊松比，当时，称相应的弹性力学问题为近不可压缩问题。

问题(1)对应的变分问题可描述为：

求，使得对任意的，下式总成立：

, (3)

其中，

这里，为线性微分算子，为弹性矩阵。

设是区域上的四面体网格剖分，其中为上所有剖分单元的最大直径。引入如下次拉格朗日有限元空间：

其中，为单元上次数不超过的多项式的全体。图1给出了剖分单元中四面体线性、二次和三次元插值节点的分布情况。

变分问题(3)的次元提法为：求，使得对任意的，下式总成立：

(4)

设变分问题(4)的次有限元离散化线性系统的矩阵形式为

(5)

由于，利用Korn不等式，易于证明有限元方程(5)的解存在且唯一。但对近不可压缩问题，通常的低阶协调元(如线性元、二次元)解不再收敛到原问题的解或达不到最优收敛阶，即出现所谓的体积闭锁现象。下面，举几个例子来进行验证。

算例1(纯位移边界弹性力学问题)该数值例子取自文献[13] 。设求解域，，选取精确解为，其中：

图1. 四面体线性、二次和三次单元插值节点的分布情况

对区域进行一致四面体网格剖分，取及不同的进行数值实验，表1列出了线性元解和精确解之间的绝对误差和相对误差，表2列出了二次元和精确解之间的绝对误差、相对误差以及相应收敛阶的比值。表格中的“TCR”表示理论收敛阶的比值，对二次元情形，它可写为：，

“NCR”表示数值收敛阶的比值。从表中结果可知，当时，对线性元情形，我们有或

，即使当很小时，其相对误差总是接近于1，无法得到任何收敛性结果；对二次元情形，有限元解虽然收敛，但达不到最优收敛阶。这就是近不可压缩材料的所谓体积闭锁现象。

算例2(悬臂梁问题)考虑如图2所示细长梁，设长，高度和厚度相等即。弹性模量，在端部固支，在悬臂梁上表面全长施加向下的分布载荷，不考虑自重。采用四面体一致网格进行剖分，其中图3给出了两种类型的网格剖分图。由于，采用材料力学的解作为近似理论解，如对轴线上的点和点，它们在方向上的位移解分别为和。分别采用4节点线性元和10节点二次元进行计算，数值结果总结如表3所示。

图2. 悬臂梁几何结构示意图

图3. 悬臂梁网格剖分

表3. 悬臂梁轴线上点和点挠度(单位：m)

算例3(Cook膜问题)考虑如图4所示的三维Cook膜问题，其左侧面固定，右侧面上作用有大小为且垂直向上的切应力，弹性模量，不计自重。采用两种类型四面体网格进行剖分，分别采用4节点线性元和10节点二次元进行计算，相应的数值结果总结如下。

Case 1：一致网格

对区域进行一致四面体网格剖分。表4给出了特征点随问题规模及不同下方向上的线性元和二次元解。

Case 2：基于CVT网格的Delaunay四面体网格

对作基于曲面CVT网格的Delaunay四面体网格剖分[14] ，图5给出了含不同单元数的两种网格剖分图。表5给出了特征点随问题规模及不同泊松比下方向上的线性元和二次元解。

从上述算例2和3的数值结果可知，对于近不可压缩问题，当采用线性元时，即使采用各向同性网格或Delaunay四面体网格，所求得的位移值均偏小，出现闭锁现象；当采用二次元时，特征点处的位移能得到一定程度的改善，但对有些情形，如悬臂梁问题，计算结果还不是很令人满意，需要研究新的数值方法。另外，实际计算时若能采用各向同性网格，对问题将具有更好的收敛性。

本文，我们将利用罚函数有限元法来解决三维弹性材料的体积闭锁现象，该方法因其具有程序实现简单、罚数易于确定及不改变泛函驻值性质等特点，是克服体积闭锁现象的一种有效方法。下面，我们针对混合边界条件的三维近不可压缩弹性问题，详细推导罚函数有限元法的计算格式，分析该方法实施成功的条件，并通过数值实验验证方法的有效性。

图4. 三维Cook膜问题

表4. 特征点在一致网格下的有限元解 (单位：mm)

图5. 两种基于曲面CVT网格的Delaunay四面体网格

表5. 特征点在Delaunay四面体网格下的有限元解

3. 罚函数有限元分析

3.1. 罚函数有限元格式的建立

设三维近不可压缩弹性问题的能量泛函为

引入

和

其中，是平均应力(也称体积应力)，是体积应变，在小变形下它们存在下述本构关系：

或

这里，，是体积模量，它可用和表示为。于是，可引入偏斜应力向量和偏斜应变向量，具体如下：

而和之间又存在如下关系：

利用以上各式，弹性应变能可写为：

(6)

上式右端第一项和第二项分别表示偏斜应变能和体积应变能。将式(6)代入，可得：

假设已被划分为个单元，在每一个单元上，位移函数可用单元节点位移列阵表示，即，而，这里，为形函数矩阵，为应变矩阵。于是，我们可得到离散形式的总位能表达式为：

设，其中，为单元节点自由度和结构节点自由度的转换矩阵，为结构节点位移列阵。令和，则有

(7)

设，，我们有

(8)

由可得有限元求解方程：

(9)

或

(10)

这里，为罚函数法中的罚数，且当时，；子块矩阵和分别为

当时，为了使有限元方程(9)或(10)有非零解，子块矩阵必须是奇异的，则有：秩，其中，是结构的单元总数，是计算子块矩阵所采用的积分点数，是中被积函数中独立行(或列)数，易知。另一方面，为保证有限元方程(9)有解，系数矩阵

必须是非奇异的，即是满秩，则有：，其中，是计算的积分点数，它可以不同于，是中被积函数中独立的行(或列)数，对于三维问题，。这样，罚函数有限元方法实施成功的条件是

(11)

和

(12)

其中，是系统的总自由度数。有关罚函数有限元更为详尽的描述可参见文献[15] 。下面，以四面体二次单元为例，在计算和时为保证条件(11)和(12)成立，讨论相应积分点数的选取方案。

为方便起见，考虑0狄利克雷边界条件的三维近不可压问题，求解区域，对进行一致四面体网格剖分，并设每个方向的剖分段数分别为，和，取二次单元(如图6所示)进行分析，计算时Hammer积分点数为，计算时Hammer积分点数为。不妨设，整个结

构的单元总数为，的阶次为，取不同的，验算条件(11)(或

)和条件(12)(或)是否成立，具体如表6所示。一般地，若及，则条件(12)显然成立；但若，则条件(11)已不再成立，非奇异。

3.2. 数值实验及结果分析

将上面建立的罚函数二次元应用于第二节中三个算例的分析，以验证其有效性和鲁棒性。对算例1，当时，在不同网格规模下罚函数二次元解和精确解之间的误差及收敛阶比较总结如表7所示。对算例2，当时，在不同网格规模下悬臂梁轴线上点和点的罚函数二次元解分别如表8所示。

对算例3，当时，在基于曲面CVT网格的Delaunay四面体网格剖分下罚函数二次元解及与其它几种方法的比较如图7所示，相应的位移场如图8所示。

图6.单元插值节点积分点积分点

图7. 特征点在方向上的几种有限元解的比较

图8. Cook膜问题位移场，其中：左图为线性元，右图为罚函数二次元

从上述数值结果可知，对近不可压缩问题，罚函数二次元具有更好的收敛性。例如，对纯位移边界

弹性力学问题，当和时，协调二次元解的相对误差为0.144264，数值收敛阶为3.716，

而罚函数二次元的相对误差为0.024652，数值收敛阶为5.135；对悬臂梁问题，当时，在网格剖分下，协调二次元解，而罚函数二次元解；对Cook膜问题，当采用基于曲面CVT网格的Delaunay四面体网格时，罚函数二次元解具有更好的收敛性，所选特征点的位移值均随单元数的增加而稳定地收敛于理论解。

4. 结语

有限元方法是求解三维近不可压弹性问题的一类重要的离散化方法，若利用通常的低阶协调元方法会出现体积闭锁现象。罚函数有限元法是解决体积闭锁现象的一种行之有效的数值方法。该方法由于将一般罚函数方法中的罚数与体积模量相关联，解决了罚数难以确定的问题，若计算时能确保与罚数相关的子块刚度矩阵为奇异矩阵，则相应的有限元解一致收敛于精确解。另外，实际计算时，为提高精度，应尽可能地采用各向同性网格，特别是基于曲面CVT网格的Delaunay四面体网格。

本文仅对二次元进行了罚函数有限元实验及结果分析，若采用更高阶单元，则计算精度将会得到大大改善，但其计算复杂性将增大。另外，对近不可压缩问题，体积闭锁现象解决之后，下一步将是如何快速求解相应的有限元离散化代数系统，这是提高其有限元分析整体效率必须解决的问题之一。针对罚函数有限元分析中形成的大型的、稀疏的和高度病态的正定方程组，如何设计高效求解算法是一个非常困难的问题，这将是我们今后进一步研究的问题。

致  谢

感谢国家自然科学基金资助项目(10972191)和湖南省自然科学基金资助项目(14JJ2063)对本论文工作的资助。

参考文献 (References)

NOTES

^*通讯作者。