廣義有限元方法

廣義有限元方法（Generalized Finite Element Method，GFEM），是基於傳統有限元方法（FEM）的一種數值分析方法，能夠更加優雅、準確地解決材料交界面和斷裂力學等非連續問題。這類問題在劃分單元網格時，應用傳統有限元方法往往需要使網格邊界與非連續界面儘可能重合，帶來極大的網格依賴（mesh-dependence）。對於複雜邊界問題及形狀變化問題（例如拓撲優化和裂紋擴展等），傳統有限元需要花費大量的時間建立單元網格以及重建網格，同時還需要引入計算量處理裂尖端（crack tips）來避免奇異剛性矩陣。廣義有限元正是在這樣的背景下誕生，它對傳統有限元的形函數（shape functions）進行擴展，在包含非連續的單元網格內使用強化方程（enrichment functions），來準確模擬非連續的物理場，這種方法又可稱為擴展有限元（extended finite element method，XFEM）。

## 背景

### 非連續性

非連續性，在有限元方法中的定義，指物理場量或其梯度在極小的區間內發生了極大的改變，其中「極小的區間」是相對於觀察定義域而言，「物理場量「包括位移場（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是採用了廣義有限元演算法的商業有限元軟體。