2.3　扩展有限单元法的基本原理

XFEM是基于单位分解的思想在常规FEM位移模式中加入加强函数，即阶跃函数和裂尖渐进位移场函数，从而能够反映裂纹面不连续性和裂尖奇异性。与常规FEM相比，XFEM有如下优势：①网格与结构的几何尺寸和内部结构相互独立，因此划分网格时无需把裂纹作为单元边界，裂尖无需作为单元结点；②无需为反映裂尖应力奇异性而进行高密度的网格划分；③分析裂纹扩展时，无需重新划分网格。

2.3.1　单位分解法

单位分解法（Partition of Unity Method-PUM）是利用分片的思想逼近局部函数，然后再将各片进行“黏合”，进而形成对函数的全局逼近。单位分解法是由Duarte[130]和Melenk[131]在1996年先后提出来的。在求解区域Ω内，对Ω分为相互交叉的子区域ΩI，引入函数φI（x），使得φI（x）与子区域ΩI相对应，函数φI（x）只在ΩI内非零，并满足单位分解条件，即在域内任一点x，函数φI（x）满足

基于单位分解的思想，可以在常规有限元的位移模式中添加一些加强函数，从而模拟局部不连续特性。XFEM就是基于单位分解的思想而发展起来的，从而保证了计算的收敛性；由于形函数在单元内部仍然具有单位分解特性，从而使得XFEM继承了常规有限元的优点，是一种能够有效分析不连续问题的数值方法。

2.3.2　水平集法

水平集法（level set method-LSM）是Osher和Sethian首先提出的一种确定界面位置和追踪界面移动的数值技术。该方法是将曲线或曲面嵌入到比高一维空间的水平集函数Ψ[X（t），t]。在一定的速度场驱动下，通过求解水平集方程实现曲线或曲面的边界运动分析与跟踪。

移动界面γ（t）⊂R2可表示为函数Ψ[X（t），t]的水平集曲线：R2×R→R，其中γ（t）=｛X∈R2：Ψ（X，t）=0｝。

γ（t）的移动可由Ψ的演化方程得到[132]：

式中：F为界面上点X∈γ（t）在界面外法线方向的速度。

初始条件Ψ（X，0）一般选用符号距离函数表示，在界面一侧其值大于零，在另一侧其值小于零，界面处等于零。

若裂纹不是将区域分成两个离散的子区域，则需采用两个相互垂直的水平集函数描述裂纹，即法向水平集函数Ψ[X（t），t]和切向水平集函数φ[X（t），t]。Ψ[X（t），t]和φ[X（t），t]均为符号距离函数。水平集函数ΨX（t），t可表示为：

式中：X为考察点；X*为裂纹面上的点；XΓ为裂纹面上距离X最近的点；n为裂纹面上的单位外法线矢量。如图2.4（b）所示，若X位于γ（t）所定义的裂纹的上方，上式符号取正，否则取负。同理，可以定义φ[X（t），t]。

由图2.4（a）可知，Ψ=0且φ＜0为裂纹面；Ψ=0且φ=0为裂纹波前。对于含有两个裂尖的裂纹，需定义两个切向水平集函数。为了计算方便，通常取其最大值，即φ=max（φi）（i=1，2），定义单一的切向水平集函数。

水平集法是追踪界面移动的一种十分方便有效的数值方法，它有如下优点：

（1）可以在笛卡儿网格上对演化中的曲线曲面进行数值计算而不必对曲线曲面参数化；

（2）可以方便地追踪物体的拓扑结构改变。

2.3.3　位移模式

如图2.5所示为含一任意裂纹的扩展有限元计算网格，网格的划分独立于裂纹的位置。XFEM的位移模式是在FEM的位移模式上加入加强函数来反映不连续性。对于裂纹问题，分别采用广义的Heaviside函数（阶跃函数）和裂尖渐进位移场函数加强裂纹贯穿单元和裂尖单元。

如图2.5所示问题的位移逼近可以表示为

式中：ui为结点位移；aj和bαk为加强结点的附加变量；Ni（x）、Nj（x）和Nk（x）为结点形函数；I为离散域内所有结点集；J为形函数支撑域完全被裂纹切割的结点集，为H（x）函数加强的结点集；K为K*与其邻近结点（与K*中结点共享单元的结点）的结点集，为裂尖函数φα（x）加强的结点集；R（x）为线增函数，该函数的引入是为了消除裂尖函数加强引起的混合单元（见图2.5中灰色区域单元）[116]。

对于裂纹贯穿单元的加强函数H（x），取广义的Heaviside函数（阶跃函数）。在裂纹的一侧为+1，另一侧取-1，即

式中：x为考察点；x*为裂纹上距离考察点x最近的点；n为x*点处裂纹的单位外法线向量。

线增函数R（x）的定义为

式中：φα（x）为裂尖加强函数，一般根据裂尖渐进位移场函数确定，反映裂尖应力的奇异性和不连续性。对于各向同性弹性体，φα（x）一般采用下列形式的裂尖分支函数：r和θ为裂尖处的局部极坐标。

在式（2.34）中的四个分支函数中，反映了沿裂纹表面的间断性，其他三个函数是为了改善裂尖附近解的奇异性[133]。

2.3.4　支配方程

位移模式构造后，就可以和常规有限元方法一样，由虚功原理推导其支配方程。

在外力作用下，裂纹面可能处于张开或闭合状态。当裂纹面处于闭合状态时，需要考虑裂纹面间的接触条件，否则裂纹面间会发生嵌入现象。下面推导分析接触问题的XFEM支配方程。

根据XFEM的位移逼近，式（2.31），裂纹面的相对位移uc可以表示为

式中：d为单元结点参数列阵，结点位移存在该列阵前面，结点加强变量放在后面；T为整体坐标到裂纹面局部坐标的转换矩阵；d与Nc和裂纹单元的结点加强方式有关。

其中

对于H（x）函数加强的结点

对于φα（x）函数加强的结点

式（2.15）中的增广型的Lagrange乘子场可以写成如下的矩阵形式：

式中：M为Lagrange乘子场的插值形函数。

其中，

对于H（x）函数加强的结点

对于φα（x）函数加强的结点

另外，增广型的Lagrange乘子矩阵为

离散虚功方程式（2.27），并结合式（2.31）中的XFEM位移模式，则分析摩擦接触问题的XFEM支配方程为

其中，d=[u a b]T为结点未知量向量。

式中：u为常规结点未知量，a和b为结点附加未知量；K为劲度矩阵，来自式（2.27）的第1项；Kc为裂纹面接触劲度矩阵，来自式（2.27）的第2项；K和Kc的叠加为系统的整体劲度矩阵；F为整体荷载列阵。

K通过单元劲度矩阵组装而成，单元劲度矩阵为

式中：分别为常规应变矩阵、广义的Heaviside函数加强和裂尖函数加强对应的加强应变矩阵。

Kc由裂纹面单元劲度矩阵组装而成

F由单元荷载列阵组装而成

其中

2.3.5　加强函数的偏导数计算

在应变矩阵的计算过程中，需计算加强函数的偏导数。广义的Heaviside函数的偏导数等于零，下面给出裂尖加强函数对整体坐标的偏导数计算过程。

φα对局部坐标（图2.6）的导数为

由式（2.34）中φα的具体形式，可以得到：

由图2.6中的局部坐标可知：

将式（2.54）和式（2.53）代入式（2.52），可得和的具体形式为

结合局部坐标与整体坐标的关系，就可以得到裂尖加强函数对整体坐标的偏导数如下：

将式（2.55）中的各分量代入式（2.56）即可得到裂尖加强函数对整体坐标的偏导数。