广义有限元方法

广义有限元方法（Generalized Finite Element Method，GFEM），是基于传统有限元方法（FEM）的一种数值分析方法，能够更加优雅、准确地解决材料交界面和断裂力学等非连续问题。这类问题在划分单元网格时，应用传统有限元方法往往需要使网格边界与非连续界面尽可能重合，带来极大的网格依赖（mesh-dependence）。对于复杂边界问题及形状变化问题（例如拓扑优化和裂纹扩展等），传统有限元需要花费大量的时间建立单元网格以及重建网格，同时还需要引入计算量处理裂尖端（crack tips）来避免奇异刚性矩阵。广义有限元正是在这样的背景下诞生，它对传统有限元的形函数（shape functions）进行扩展，在包含非连续的单元网格内使用强化方程（enrichment functions），来准确模拟非连续的物理场，这种方法又可称为扩展有限元（extended finite element method，XFEM）。

## 背景

<image>

### 非连续性

非连续性，在有限元方法中的定义，指物理场量或其梯度在极小的区间内发生了极大的改变，其中“极小的区间”是相对于观察定义域而言，“物理场量“包括位移场（displacement field）、温度场（temperature field）、势场（potential field）等。物理场量的非连续往往被称为强非连续（strong discontinuities），其梯度的非连续被称为弱非连续（weak discontinuities），前者的表现方式包括断裂力学中裂缝两边的位移差，后者可举例为复合材料受力时在材料交界面的位移梯度差或是位移转折（kink）。

### 建模方法

非连续问题的建模方法，可从基本原理上分为两类。第一类方法依赖于多项式近似空间（polynomaial approximation spaces）和单元网格与非连续界面的重合程度。其中最高效的多项式近似空间构造方法即为传统有限元方法，同时还有一些无网格法给多项式近似空间的构造带来不同的可行性。然而该类方法在处理复杂边界问题时，需要花费大量的时间建立单元网格，对于形状变化问题（例如拓扑优化和裂纹扩展等），求解过程中重建网格不可避免，对于断裂力学还需要引入计算量处理裂尖端来避免奇异刚性矩阵，给求解过程带来极大不便，求解问题被严重限制。即便如此，在广义有限元未被推出时代，传统有限元被广泛用在分析非连续问题上，尤其在土木工程中对河床、片岩、叶理、节理、裂纹对管道及地基等基础设施建造的影响分析中起到重要作用。

随着计算机硬件的提升，工程学上对力学建模的精确程度要求逐渐提高，又由于传统有限元分析非连续问题带来的种种复杂性，广义有限元应运而生。这就属于非连续问题的第二种建模方法——强化有限元多项式近似空间，使得非连续的建模可以独立于单元网格。这种方法在传统有限元的公式中，加入一项特殊的函数，称为强化项（enrichment），该强化项包含解中非连续性的信息，即物理场或其梯度的在非连续界面两端的差值，详细原理参见下一章节。该概念基于有限元单位分解法（partition of unity method），最初于1995年J.M. Melenk的博士论文中被提出，随后Melenk与他的博士导师、有限元领域的地基人物Ivo Babuška于1996年整理发表。同一时期，另一所学校的博士生Duarte也在毕业论文中发表了相似理论，称为hp云方法（hp cloud method）。随后，广义有限元在断裂力学中的应用立即被深入研究，由Belytschko及其同事延伸成扩展有限元方法，由于最终被证明广义有限元与扩展有限元其实是同样的一种方法，该两个名词可互换使用。广义有限元的重要特性即在于解决了传统有限元中求解过程与单元网格的依赖性，网格的划分可以更加自由，同时，诸如裂纹扩展、拓扑优化、凝固过程（solidification）等涉及到非连续界面不断发生改变的复杂的问题，该方法都能更有效率、更准确地计算。

## 基本公式

在广义有限元的基本公式中，第一项与传统有限元完全一致，第二项是其对传统有限元近似空间的扩展，即强化项，此处为了简化，只使用了一项强化项作为说明。

$$u(x)\approx u^h(x) = \underbrace{\sum_{i\in I}N_i(x)u_i}_\text{std.FEM} + \underbrace{\sum_{i\in I}N^*_i(x)\psi(x)a_i}_\text{enrichment}$$

其中$I$为定义域内所有节点的集合；$N_i$与$N_i^*$是传统有限元中的形函数，他们大部分情况下会相等（这里区别表示是由于，在后期一些高阶广义有限元的发展中，一些研究者曾尝试在这两项中使用不同阶数的形函数）；系数$u_i$是标准有限元的节点所对应的自由度（degrees of freedom）;$a_i$是强化项引入的额外节点自由度，对应强化项作用的强化节点（enriched nodes），即是包含非连续的单元网格的节点。

$\psi(x)$作为该公式中最重要的部分，是结合对解当中非连续性的预了解而建立的强化函数（enrichment function）。$\psi(x)$与$N^*_i(x)$的乘积，即运用了传统有限元拉格朗日形函数的单位分解性质（partition of unity property）：

$$\sum_{i\in I}N^*_i(x)=1$$

使$\psi(x)$仅仅作用在包含非连续界面的网格单元中。

下面解释最初的广义有限元概念中，使用的强化函数$\psi(x)$：

$$\psi(x)=H(\phi(x))= \begin{cases} <pre> 0\quad \phi(x)\leqslant 0 \\ 1\quad \phi(x)> 0\\ </pre> \end{cases}$$

$$\psi(x)=abs(\phi(x))=|\phi(x)|$$

其中，$\phi(x)$被称为Level set函数，用于描述非连续界面的位置：

$$\phi(x)= \begin{cases} <pre> -\phi_d(x) \quad \forall x \in \Omega^- \\ +\phi_d(x)\quad \forall x \in \Omega^+\\ </pre> \end{cases}$$

$\Omega^-$和$\Omega^+$是定义域$\Omega$被非连续界面分割出的子域。当中$\phi_d(x)$是每个点$x$到非连续界面$\Gamma$的最短距离：

$$\phi_d(x) = \min_{x_{\Gamma}\in \Gamma}\|x-x_{\Gamma} \|,\quad x\in \Omega$$

对于一个一纬杆的弱非连续问题，如果使用1个单元，并给每个节点添加两个强化函数，用上述基本公式构建的广义有限元近似空间，其解$u^h$可用矩阵形式表示为：

$$u^h(x) = [ \, N_1 \quad N_1\psi_{11} \quad N_1\psi_{12}\quad| N_2 \quad N_2\psi_{21} \quad N_2\psi_{22} ]\, \begin{bmatrix} <pre> u_1 \\ a_{11} \\ a_{12} \\ u_2 \\ a_{21}\\ a_{22} \end{bmatrix} </pre> = \textbf{NU}$$

## 自身问题

虽然广义有限元的网格独立性很大程度上提高了非连续的建模，可它基本公示中的强化项却引入了新的问题。

- Blending elements的计算误差难以克服。我们称包含非连续界面的单元为reproducing element，而与其相邻的单元则为blending element。在后者中，由于强化项的影响导致其单元的单位分解性质丢失或不完全，因而无法准确计算出节点场值。
- 强化项给狄利克雷边界条件施加带来困难。
- 一些强化项会导致奇异刚性矩阵。

## 发展

- 高阶有限元：提高了其计算精确度和收敛率
- Interface-Enriched Generalized Finite Element Method
- Discontinuity-Enriched Finite Element Method

## 软件

Abaqus是采用了广义有限元算法的商业有限元软件。