基于T5单元的体积不可压缩问题光滑有限元法

王思照，张仪萍

（1.浙江大学建筑工程学院，浙江杭州310058；2.浙江省水利水电勘测设计院，浙江杭州310002）

基于T5单元的体积不可压缩问题光滑有限元法

王思照1，2，张仪萍1

（1.浙江大学建筑工程学院，浙江杭州310058；2.浙江省水利水电勘测设计院，浙江杭州310002）

提出5节点的四面体单元（T5），将无网格法的面积权重应变光滑法和光滑有限元法应用于该5节点四面体单元，提出用于解决三维体积不可压缩线弹性体的算法：基于节点光滑域的选择性体积权重应变光滑模型（T5-p NVW/NVW）.数值算例显示，四节点四面体单元采用基于节点的光滑有限元法（T4-NS）无法完美解决体积锁定，相比于T4-NS法，利用提出的T5-p NVW/NVW模型能够较精确地解决体积锁定问题，并完美解决应力的棋盘式波动.

体积锁定；体积权重应变光滑；5节点的四面体单元（T5）；αFEM；压力波动

体积锁定是有限元分析的一个难点，它是指当处理材料不可压缩问题时（即当材料泊松比趋于0.5或体积模量与切变模量比趋于无穷时），有限元方程的位移解出现病态和棋盘式的压力解[1-3].自有限元发展以来，诸多学者对其进行研究，其中颇具代表性的解决方案有：一致降阶积分法[4-7]、B-bar法[2]、u/p混合元法[2，8]、平均节点压力法[9-12]等.

为了提高有限元法三角单元的计算效率，Liu等[13-17]结合了无网格法的应变光滑法[18-20]，提出一系列基于不同积分子域的光滑有限元法，包括：基于单元子域的CS-FEM[13]、基于边域的ES-FEM[21]、基于节点光滑域的NS-FEM[22]、基于三维面域的FS-FEM[23]（三维四面体单元）等.其中，采用一个光滑子单元域的CS-FEM具有免于锁定的特性，该特性类似于采用一点高斯积分的FEM，虽然能够得到较精确的位移解，但是存在额外零能模态和不均匀的压力分布等缺陷.基于光滑子单元域的选择性积分法（selective cell-based SFEM）[24-25]虽然能够解决体积锁定和额外的零能模态，但其解锁本质是基于一个单元子域的光滑有限元法（CS-FEM），无法解决压力解的波动问题；另一个天然的免于体积锁定的光滑算法是NS-FEM[26-27]，不存在额外的零能模态，无法解决压力的棋盘式分布，当采用较规则网格时，采用该算法能够给出问题的应变能上限.还有一些光滑算法包括：αFEM[28-29]和基于不同光滑域的选择性算法（a domain-based selective scheme）[21]均基于NS法的解锁本质，无法解决压力的棋盘式波动.

Nguyen-Xuan等[30]研究指出NS-FEM的解锁特性不完美，通过二维模态分析[1-2，31-33]发现，NSFEM法采用三角单元并无法解决沙漏模态，在解决一些对约束条件较苛刻的体积不可压缩问题时，NS-FEM法的解锁特性不彻底，既无法解决压力的波动，而且会出现位移解的不精确.Nguyen-Xuan 等[30]通过在FEM法的三角单元内增加一个节点并采用气泡型形函数[34]（bubble function），结合ESFEM法提出bES-FEM[30]法用于解决二维体积锁定问题，能够较有效地解决体积锁定和压力的棋盘式振荡.Wu等[35-37]提出一种新型的ME-FEM的无网格四面体单元，用于解决体积锁定问题，能够较好地解决压力波动，但该单元使用的隐式形函数对于简单问题的计算效率没有有限元法高.Zhang等[38]提出基于5节点的四边形单元和9节点的六面体单元，用于解决二维和三维体积不可压缩问题，均较好地解决了体积锁定，获得了较精确的位移解和压力解.

本文通过在四面体单元（T4）内部增加一个节点，采用气泡型形函数（bubble function）构造出一种5节点四面体单元（T5）.基于该T5单元，结合无网格法的面积权重应变光滑技术和NS-FEM法，提出一种解决体积锁定的算法，用于解决普遍算法中无法解决的应力波动问题.本文提出的T5-p NVW/NVW模型具有αFEM的特点，能够调节偏斜应变能以达到提高解的准确性的目的，通过一系列数值算例检验了该单元在解决体积不可压缩问题时的效能.

1 5节点四面体单元（T5）

考虑一个四面体单元[39-41]，中心加一节点，第5节点形函数为

其中，xe（i）、ye（i）、ze（i）为四面体单元各角点的坐标值.T5单元其他4个节点的形函数为

2 基于节点光滑域的体积权重应变光滑法（T5-NVW）

问题域ΩNVW的划分是基于节点的，通过将第5节点与各个边中点和面中心点相连，将T5单元划分为4个部分（见图1），每个部分占用一个角节点，因此各个节点光滑域为与相同节点k相邻的所有部分的合集.基于节点的应变光滑域有3类，包括内部节点光滑域、边节点光滑域和角节点光滑域.典型的内部节点光滑域k的光滑域体积为

基于节点k的应变光滑算子为

式中：dI为节点I的位移列向量，I为与节点k的光滑域相关的所有节点，uh（x）为位移向量表示光滑后的应变位移算子）为光滑的应变位移矩阵，

图1 基于节点的体积权重应变光滑域（T5-Nv W）Fig.1 Illustration of node-based volume-weighted strain smoothing domain for T5-NVW

其中

总体刚度矩阵为

其中，nNVW为节点单元域总数.

3 基于不同积分域的选择性积分法（T5-p NVW/NVW）

将三维线弹性问题材料矩阵D[21-22]分解为与形变和体积相关两部分D1和D2.其中，μ=E/[2×（1+ν）]为切变模量，与D1相关；λ=2νμ/（1-2ν）为体积模量，与D2相关.对于三维问题，材料矩阵分解为

T5-p NVW/NVW模型的刚度矩阵为

4 算例分析

4.1 三维悬臂梁

考虑一个一端固支顶部受均布荷载的三维悬臂梁例子，用于比较分析本文所提方法在解决体积锁定问题的精确性.如图2所示，悬臂梁尺寸为5×1×1，顶部受均布荷载的合力f=1，泊松比ν=0.499 999 9，弹性模量E=1000.

图2 顶部受均布荷载的三维悬臂梁几何尺寸及5节点四面体单元网格划分Fig.2 Domain discretization using five-node tetrahedral element for three-dimensional cantilever beam subjected to uniform pressure at top

为了分析解的准确性，本文加入4节点四面体单元（T4）采用NS-FEM法进行比较分析.对于模型T5-p NVW/NVW，参数p的确定如图3所示.图中，Een为应变能.采用3套网格（具有相同的划分比例5∶1∶1）计算问题应变能，根据αFEM法的计算原理可知，3条曲线的交点即为问题的估计最佳应变能Een≈0.038 5，对应最优参数p=0.288.从确定参数p的3套曲线图可以看出，参数p在定义域（0，1]范围变化时，只是进行应变能的微调搜索过程，未出现由“体积锁定”引起的应变能无界状态（突变巨大且无界），这证明参数p对是否发生“体积锁定”没有任何影响.如表1所示为当网格划分精细时，三维悬臂梁应变能的收敛趋势.可见，相比于T4-NS法，T5-p NVW/NVW（p=0.288）模型能够较快地找到精确应变能解.如表2所示为三维悬臂梁顶部端节点A（见图2）的z向位移uz（八节点六面体单元采用选择性降阶积分H8-SRI，计算得出参考值为uz=-0.19，其中负号表示方向为z轴负方向，采用网格为100×20×20），T5-p NVW/NVW（p=0.288）相对于T4-NS法[22]得到较精确的位移解.如图4所示为当ν=0.499 999 9，位移放大倍数为5时，三维悬臂梁的位移图及压力云图.图4显示T5-p NVW/NVW模型和T4-NS法在解决体积锁定问题时，均能够获得较准确的位移解.从压力云图可以发现，T4-NS在端部的压力偏高，T5-p NVW/NVW模型的压力解比T4-NS法更光滑.

图3 顶部受均布荷载的三维悬臂梁的参数p的确定（T5-p Nv W/Nv W）Fig.3 Parameter determination of three-dimensional cantilever beam subjected to uniform pressure at top for T5-p NVW/NVW

4.2 三维冲压问题

考虑一三维冲压（punch）问题，尺寸及网格划分如图5所示，E=1000.顶部1/4区域受z向位移uz=-0.03.

表1 顶部受均布荷载的三维悬臂梁应变能（网格i≡i×10×2×2×6）Tab.1 Strain energy obtained using different methods for three-dimensional cantilever beam subjected to uniform pressure at top（mesh i represents i×10×2×2×6）

表2 顶部受均布荷载的三维悬臂梁节点A处的z向位移uz（网格i≡i×10×2×2×6）Tab.2 Tip deflection at point A obtained using different methods for three-dimensional cantilever beam subjected to uniform pressure at top（mesh i represents i×10×2×2×6）

图6采用均匀10×10×5×6网格，比较分析T4-CS/NS法和T5-CS/NVW模型（为了方便对比分析，T5-p NVW/NVW模型参数直接采用p=1，因此退化为T5-CS/NVW法；由于T4-NS法的计算结果波动较大，采用T4-CS/NS法进行比较，T4-CS/NS法为采用基于T4单元的CS-FEM计算偏斜应变能，NS-FEM计算体积应变能）在材料近似不可压缩问题下的位移变形图和压力云图的光滑度.由图6（a）、（b）可知，当泊松比为0.49时，采用2种方法均能够得到位移解，不存在压力波动.当泊松比为0.499 999 9时（趋于0.5），T4-CS/NS法在处理体积锁定问题时无法得到精确的位移解，顶部位移变形图（见图6（c））不如T5-CS/NVW法光滑（见图6（d）），还出现了应力的不均匀波动现象；相比于T4-CS/NS法，T5-CS/NVW模型得到了较准确的位移解和光滑的压力解.

5 结 论

（1）本文基于无网格ME法和基于节点的NSFEM法，提出应用于有限元法的T5单元，给出相应的显式形函数.

（2）结合无网格法的面积权重应变光滑法与光滑有限元法的CS法和NS法，提出基于节点光滑域的三维选择性体积权重应变光滑模型（T5-p NVW/NVW）.数值算例显示，采用该模型能够很好地解决体积锁定问题引起的位移解不准确和棋盘式压力波动.

（3）T5-p NVW/NVW具有类似αFEM的功能，通过参数p来调整系统的刚度使其更加接近系统的真实刚度，从而提高T5-p NVW/NVW法的解的准确性.参数p的确定至少需要2套网格，在一定程度上降低了计算效率，高效、简便地确定参数p是下一步研究的重点.采用高效准确的FS-FEM法、三维ES-FEM法或其他算法代替p NVW计算偏斜应变能是一种选择.p的调节与体积锁定没有任何关系，p在定义域（0，1]内变化均对是否发生“体积锁定”没有任何影响.实际工程问题时，可以直接采用T5-CS/NAW计算，不需要确定p，也可以满足工程要求.如果需要提高精度，在网格、载荷、结构的几何形式确定之后，可以采用本文方法确定参数p.

图4 顶部受均布荷载的三维悬臂梁位移变形图及压力云图Fig.4 Deformation and distribution of pressure for three-dimensional cantilever beam subjected to uniform pressure at top

图5 三维冲压问题的几何尺寸及T5单元网格划分Fig.5 Domain discretization using five-node tetrahedral element for three-dimensional punch problem

图6 三维冲压问题位移变形图（放大倍数为50）及压力云图Fig.6 Deformation and distribution of pressure for three-dimensional punch problem in compressible and nearly incompressible case

（4）本文提出的T5单元采用三维四面体网格，适用于进一步研究大变形及网格极不规则的问题.

（

）：

[1]HUGHES T J R.The finite element method[M].New York：Dover Publications，2000.

[2]BATHE K J.Finite element procedures[M].New Jersey：Prentice-Hall，1996.

[3]李忠学.梁板壳有限单元中的各种闭锁现象及解决方法[J].浙江大学学报：工学版，2007，41（7）：1119-1125.

LI Zhong-xue.Strategies for overcoming locking phenomena in beam and shell finite element formulations [J].Journal of Zhejiang University：Engineering Science，2007，41（7）：1119-1125.

[4]HUGHES T J R.Reduced and selective integration techniques in the finite element analysis of plates[J].Nuclear Engineering and Design，1978，46（1）：203-222.

[5]FLANAGAN D P，BELYTSCHKO T.A uniform strain hexahedron and quadrilateral with orthogonal hourglass control[J].International Journal for Numerical Methods in Engineering，1981，17（5）：679-706.

[6]BELYTSCHKO T，BACHRACH W E.Efficient implementation of quadrilaterals with high coarse-mesh accuracy[J].Computer Methods in Applied Mechanics and Engineering，1986，54（2）：279-301.

[7]FREDRIKSSON M，OTTOSEN N S.Fast and accurate 4-node quadrilateral[J].International Journal for Numerical Methods in Engineering，2004，61（11）：1809-1834.

[8]BREZZI F，FORTIN M.Mixed and hybrid finite element methods[M].New York：Springer，1991.

[9]BONET J，BURTON A J.A simple average nodal pressure tetrahedral element for incompressible and nearly incompressible dynamic explicit applications[J].Communications in Numerical Methods in Engineering，1998，14（5）：437-449.

[10]DOHRMANN C R，HEINSTEIN M W，JUNG J，et al.Node-based uniform strain elements for three-node triangular and four-node tetrahedral meshes[J].International Journal for Numerical Methods in Engineering，2000，47（9）：1549-1568.

[11]KRYSL P，ZHU B.Locking-free continuum displacement finite elements with nodal integration[J].International Journal for Numerical Methods in Engineering，2008，76（7）：1020-1043.

[12]LAMICHHANE B P.Inf-sup stable finite element pairs based on dual meshes and bases for nearly incompressible elasticity[J].IMA Journal of Numerical Analysis，2009，29（2）：404-420.

[13]LIU G R，DAI K Y，NGUYEN-THOI T.A smoothed finite element for mechanics problems[J].Computer Methods in Applied Mechanics and Engineering，2007，39（6）：859-877.

[14]LIU G R，NGUYEN-XUAN H，NGUYEN-THOI T.A theoretical study of S-FEM models：properties，accuracy and convergence rates[J].International Journal for Numerical Methods in Engineering，2010，84（10）：1222-1256.

[15]LIU G R.A G space theory and weakened weak（W2）form for a unified formulation of compatible and incompatible methods.Part I：theory and Part II：application to solid mechanics problems[J].International Journal for Numerical Methods in Engineering，2010，81（9）：1093-1156.

[16]NGUYEN-XUAN H，BORDAS S，NGUYEN-DANG H.Smooth finite element methods：convergence，accuracy and properties[J].International Journal for Numerical Methods in Engineering，2008，74（2）：175-208.

[17]BORDAS S，RABCZUK T，NGUYEN-XUAN H，et al.Strain smoothing in FEM and XFEM[J].Computers and Structures，2010，80（23）：1419-1443.

[18]HUERTA A，FERNANDEZ-MENDEZ S.Locking in the incompressible limit for the element free-Galerkin method[J].International Journal for Numerical Methods in Engineering，2001，51（11）：1361-1383.

[19]CHEN J S，WU C T，YOON S，et al.A stabilized conforming nodal integration for Galerkin mesh-free methods[J].International Journal for Numerical Methods in Engineering，2001，50（2）：435-466.

[20]YOO J W，MORAN B，CHEN J S.Stabilized conforming nodal integration in the natural-element method [J].International Journal for Numerical Methods in Engineering，2004，60（5）：861-890.

[21]LIU G R，NGUYEN-THOI T，LAM K Y.An edgebased smoothed finite element method（ES-FEM）for static，free and forced vibration analyses of solids[J].Journal of Sound and vibration，2009，320（4）：1100-1130.

[22]LIU G R，NGUYEN-THOI T，NGUYEN-XUAN H，et al.A node-based smoothed finite element method for upper bound solution to solid problems（NS-FEM）[J].Computers and Structures，2008，87（1）：14-26.

[23].NGUYEN-THOI T，LIU G R，LAM K Y，et al.A face-based smoothed finite element method（FS-FEM）for 3D linear and nonlinear solid mechanics problems using 4-node tetrahedral elements[J].International Journal for Numerical Methods in Engineering，2009，78（3）：324-353.

[24]NGUYEN-THOI T，LIU G R，DAI K Y.Selective smoothed finite element method[J].Tsinghua Science and Technology，2007，12（5）：497-508.

[25]NGUYEN-XUAN H，BORDAS S，NGUYEN-DANG H.Addressing volumetric locking and instabilities by selective integration in smoothed finite elements[J].Communications in Numerical Methods in Engineering，2009，25（1）：19-34.

[26]王建明，张刚，戚放，等.基于光滑有限元法的体积锁定研究[J].山东大学学报：工学版，2012，42（3）：94-99.WANG Jian-ming，ZHANG Gang，QI Fang，et al.Study of volumetric locking based on the smoothed finite element method[J].Journal of Shandong University：Engineering Science，2012，42（3）：94-99.

[27]王建明，樊现行，裴信超，等.光滑节点域有限元法[J].山东大学学报：工学版，2013，43（2）：54-61.WANG Jian-ming，FAN Xian-xing，PEI Xin-chao，et al.Node-based smoothed cells based on finite element method[J].Journal of Shandong University：Engineering Science，2013，43（2）：54-61.

[28]LIU G R，NGUYEN-THOI T，LAM K Y.A novel FEM by scaling the gradient of strains with scaling factorα（αFEM）[J].Computational Mechanics，2009，43（3）：369-391.

[29]LIU G R，NGUYEN-THOI T.A novel alpha finite element method（αFEM）for exact solution to mechanics problems using triangular and tetrahedral elements [J].Computer Methods in Applied Mechanics and Engineering，2009，197（45）：3883-3897.

[30]NGUYEN-XUAN H，LIU G R.An edge-based smoothed finite element softened with a bubble function（bES-FEM）for solid mechanics problems[J].Computers and Structures，2013，128（1）：14-30.

[31]HUECK U，SCHREYER H，WRIGGERS P.On the incompressible constraint of the 4-node quadrilateral element[J].International Journal for Numerical Methods in Engineering，1995，38（18）：3039-3053.

[32]HACKER W L，SCHREYER H L.Eigenvalue analysis of compatible and incompatible rectangular fournode quadrilateral elements[J].International Journal for Numerical Methods in Engineering，1989，28（3）：687-704.

[33]KIDGER D J，SMITH I M.Eigenvalues and eigenmodes of 8-node brick elements[J].Communications in Applied Numerical Methods，1992，8（3）：193-205.

[34]ARNOLD D，BREZZI F，FORTIN M.A stable finite element for the stokes equations[J].Calcolo，1984，21（4）：337-344.

[35]WU C T，HU W.Meshfree-enriched simplex elements with strain smoothing for the finite element analysis of compressible and nearly incompressible solids[J].Computer Methods in Applied Mechanics and Engineering，2011，200（45）：2991-3010.

[36]WU C T，HU W，CHEN J S.A meshfree-enriched finite element method for compressible and near-incompressible elasticity[J].International Journal for Numerical Methods in Engineering，2012，90（7）：805-938.

[37]WU C T，HU W.A two-level mesh repartitioning scheme for the displacement-based lower-order finite element methods in volumetric locking-free analyses [J].Computational Mechanics，2012，50（1）：1-18.

[38]ZHANG Y P，WANG S Z，CHAN D.A new fivenode locking-free quadrilateral element based on smoothed FEM for near-incompressible linear elasticity [J].International Journal for Numerical Methods in Engineering，2014，100（9）：633-668.

[39]TIMOSHENKO S P，GOODIER J N.Theory of elasticity[M].3rd ed.New York：McGraw-Hill，1970.

[40]罗伯特D.库克.有限元分析的概念与应用[M].关正西，译.4版.西安：西安交通大学出版社，2007.

[41]徐芝纶.弹性力学简明教程[M].北京：高等教育出版社，2008.

Nearly incompressible linear elasticity using five-node tetrahedral element based on smoothed finite element method

WANG Si-zhao1，2，ZHANG Yi-ping1

（1.College of Civil Engineering and Architecture，Zhejiang University，Hangzhou 310058，China；2.Zhejiang Design Institute of Water Conservancy and Hydroelectric Power，Hangzhou 310002，China）

A new five-node tetrahedral element（T5）was proposed.The area-weighted strain smoothing technique and the smoothed finite element method were introduced into T5 element.A volumetric lockingfree scheme for three-dimensional tetrahedral meshes was proposed，which is the node-based selective domain-based strain smoothing scheme（T5-p NVW/NVW）.The benchmark numerical examples show that the proposed method can solve the volumetric locking and the pressure oscillation compared to the node-based smoothed FEM using the four-node tetrahedral element（T4-NS）.

volumetric locking；volume-weighted strain smoothing；five-node tetrahedral element（T5）；αFEM；pressure oscillation

TB115；O 343

A

1008-973X（2015）10-1967-07

2014-08-21. 浙江大学学报（工学版）网址：www.journals.zju.edu.cn/eng

浙江省科技创新团队资助项目（2010R50037）.

王思照（1989—），男，硕士生，从事软土地基处理及有限元数值计算等的研究.ORCID：0000-0003-3091-026X.

E-mail：21212076@zju.edu.cn

张仪萍，男，教授.ORCID：0000-0001-8537-8181.E-mail：zhangyiping@zju.edu.cn