XFEM数值实现 | 一维强不连续问题

今日主要内容：以一维强不连续问题为引子，介绍XFEM理论如何数值实现？

写在前面的话

接触扩展有限元（XFEM）也有一段时间了，大部分时间停留在商业软件XFEM的操作学习，但是究其根本，仍对理论涉及不深，在现成的理论层面进行Copy-Copy，论文中也是对已有公式进行Copy...

公式中的一些具体概念到底什么意思？浑然不知！比如，XFEM经典位移模式：    ，国内的教材或者已有的文献中对于该公式的解读停留在增加了自由度，仔细一想，在哪增加了自由度？增加了多少自由度？怎么将自由度添加进入位移模式中？这些细节问题往往是我们常常忽略的！

上述的问题，说大不大，说小不小。写论文时，将公式复 制即可，不影响发表，但若是对XFEM这一理论想要进一步学习，了解它数值实现过程的话，就需要对它的很多细节自己去“过一遍”，用简单的模型去验证，通过自己编程或者手算的方式来真正掌握其实现过程，最好可以“把玩它”~

在本篇推文中，木木将结合最简单的一维杆单元强不连续问题，带着大家慢慢认识XFEM数值实现过程。由于问题模型较简单，故不用程序编制，手算即可，与常规有限元结果做对比，过程中也会对XFEM原始的位移模式进行修正。精彩内容，请继续往下浏览~

问题模型

如下图所示，1D杆模型（每个节点1个自由度），左端固定，右端施加    位移，试用XFEM和FEM求得其位移。

FEM方法：划分4个单元，6个节点，整体6个自由度；

XFEM方法：划分3个单元，4个节点，整体6个自由度（常规4个自由度，2号单元为Heaviside富集单元，每个节点附加1个自由度）；

 

FEM求解

In the FEM analysis, the stiffness matrix is obtained for each bar element as

and the total stiffness matrix of FEM model can be obtained as

由已知的边界条件：    ，可求得其余节点位移    

XFEM求解

In X-FEM analysis, the discontinuity is modeled using the Heaviside step function by enrichment of the nodal points 2 and 3, in which  . The standard and enriched parts of the shape function    , the vector of Cartesian shape function derivatives    , and the stiffness matrix     for these three elements are as follows; for element (1)

Heaviside step function具体的形式已在模型图中标明，单元中心为0，以单元中心为原点进行阶跃。

For element (2):

      For element (3):

  for element (1)；     for element (3)。对于这两个式子等于0，我是这样想的，在1号单元和3号单元中，只对节点2和节点3进行加强，在单元上是没有施加富集操作的，所以对于节点2，那就是    ，对于节点3，那就是    。

The total stiffness matrix of X-FEM model can be obtained as

由已知的边界条件：    ，可求得其余节点位移      。

最终位移结果：

目前存在的问题

在求解单元刚度矩阵时，对于富集节点处的应变矩阵积分运算我还没有搞懂，具体是怎样计算的？如：    。

XFEM位移模式的修正

在以上形函数和应变矩阵公式中，我都用到了    ，对节点处的位移值进行了修正，即减掉节点处富集函数在节点处的取值。

结论

1. XFEM在满足精度要求的前提下，可用较少单元模拟强不连续问题；

2. 由于加入了富集自由度的影响，XFEM单元刚度矩阵规模较大，当模型较大时，形成的整体刚度矩阵较FEM小的多；

3. XFEM在面对裂纹这种强不连续问题时，不需要对裂纹处进行特殊单元设置，即分析过程中，网格独立于裂纹域，可大幅提高程序分析时间。

以上内容主要参考来自伊朗谢里夫大学Amir R. Khoei大佬编写的《EXTENDED FINITE ELEMENT METHOD THEORY AND APPLICATIONS》，这本书对于XFEM的理解很有启发性，带入新概念时，会一点一点带入进来，让人读了还想读。