有限元素法

有限元素法（finite element method）或有限單元法、有限元法，是一種將求解域或連續體劃分為眾多個有限大小的單元體的集合（離散化）為若干個有限大小的單元體的集合，以解答求解域或連續體問題的數值方法。在物理學上，即為使用有限元素分析物理現象，有限元法用於求解微分方程組或積分方程組數值解。有限元素（finite element）或有限元，就是被研究的結構域或介質域所分解成的有限數量小几何子區域。

在解偏微分方程的過程中，主要難點是如何構造一個方程來逼近原本研究的方程，並且該過程還需要保持數值穩定性。目前有許多處理的方法，他們各有利弊。當區域改變時（就像一個邊界可變的固體），當需要的精確度在整個區域上變化，或者當解缺少光滑性時，有限元方法是在複雜區域（像汽車、船體結構、輸油管道）上解偏微分方程的一個很好的選擇。

為了解決問題，有限元素法將大型物理系統細分為更小、更簡單的部分，稱為有限元。這是通過在空間維度上進行特定的空間離散化來實現的，該離散化是通過構建對象的網格實現的：解決方案的數值域具有有限數量的點。邊值問題的有限元素法公式化最終形成了一個代數方程組。該方法在域上近似未知函數。然後，將對這些有限元建模的簡單方程式組合成一個對整個問題進行建模的較大方程式系統。然後，有限元素法通過最小化關聯的誤差函數，使用來自變異演算的變異方法來近似求解。

由研究人員建立的有限元素網格，然後使用套裝軟體找到解決磁性問題的方法。顏色表示分析人員已為每個區域設置了材料屬性，在此圖中，淺藍色的鐵磁成分（可能是鐵），淺紫色的是空氣。將整個物理系統細分為更簡單的部分具有以下優點：

- 精確表示複雜的幾何形狀。
- 可以描述多樣的材料特性。
- 輕鬆表示整體解決方案。
- 精確描述局部現象。

該方法的工作流程包括

（1）將問題的域劃分為子域的集合，每個子域由一組元素方程表示為原始問題，然後（2）系統地將所有元素方程組重組為用於最終計算的全域方程組。

在上面的第一步中，元素方程是簡化過的方程，可以局部地近似要研究的原始復雜方程組，其中原始方程通常是偏微分方程。為了求此方程式的近似解，通常將有限元素法作為伽遼金法的特例來處理。用數學語言來說，該過程是將殘差和加權函數取內積，並將該積分設為零。簡而言之，它是通過將試驗函數擬合到偏微分方程中來最小化近似誤差的過程。殘差是由試驗函數引起的誤差，權重函數是投影殘差的多項式逼近函數。該過程消除了偏微分方程中的所有空間導數，從而使偏微分方程局部近似為一組穩態問題的代數方程，或是一組用於瞬態問題的常微分方程。如果基礎偏微分方程是線性的，則元素方程也是線性的，反之亦然。穩態問題中出現的代數方程組，便利用數值線性代數方法求解，而瞬態問題中出現的常微分方程組則使用其他數值方法（例如歐拉方法或Runge-Kutta法）通過數值積分來求解。

## 歷史

有限元法最初起源於土木工程和航空工程中的彈性和結構分析問題的研究。它的發展可以追溯到Alexander Hrennikoff（1941）和Richard Courant (1942)的工作。這些先驅者使用的方法具有很大的差異，但是他們具有共同的本質特徵：利用網格離散化將一個連續區域轉化為一族離散的子區域，通常叫做元.Hrennikoff的工作離散用類似于格子的網格離散區域; Courant的方法將區域分解為有限個三角形的子區域，用於求解來源於圓柱體轉矩問題的二階橢圓偏微分方程. Courant的貢獻推動了有限元的發展，繪製了早期偏微分方程的研究結果。

有限元方法的發展開始於五十年代中後期使用在機身框架和結構分析上，並於六十年代通過斯圖加特大學的和柏克萊加州大學的在土木工程中的應用工作中積累經驗。

基於五十年代至六十年代大型水壩計算研究的實踐經驗，1965年，中國計算數學專家馮康發表了《基於變分原理的差分格式》一文，奠定了有限元計算方法的嚴格數學理論，為後世有限元計算方法的實際應用提供了理論保證。且馮康教授的「有限元法」嚴密理論體系是先於西方的，是國際公認的當代計算數學的一項重大成就，不同的是馮康教授只是從數學方面提出有限元法的。

## 有限元概念

### 單元

單元（Element）是由節點組成的幾何體，如三角形單元，四面體單元等。

### 節點

節點（Node）是單元幾何體的端點、頂點或特定點，單元的各物理量變化均體現在節點上，例如在彈性力學問題中，一個有兩個節點的線單元的質量集中在兩個節點上，受力也只能作用在節點上，變形也用節點的位移表示。

### 自由度

節點自由度（Degree of Freedom，簡寫 DoF），是節點上變量的個數，例如用位移法解結構問題時節點自由度為3，表示單個節點上三個坐標方向上的位移，又例如熱分析時節點自由度為1，表示某個節點處的溫度值。

### 網格

網格（Mesh）是由多個單元通過共用節點組成的單元網路，用以表示待解問題域。

## 分析方法

以下用有限元分析解決兩個簡單問題，更一般的問題可以類似的推導出來。

P1是一個較簡單的一維問題

$$\mbox{ P1 }:\begin{cases} u''(x)=f (x) \mbox{ in }(0,1), \\ u (0)=u (1)=0, \end{cases}$$

其中$f$是已知函數, $u$是關於$x$的未知函數, $u''$是$u$對$x$的二階導數。

二維比較簡單的問題是狄利克雷問題

$$\mbox{P2 }:\begin{cases} u_{xx}(x,y)+u_{yy}(x,y)=f(x,y)& \mbox{ in } \Omega, \\ u=0 & \mbox{ on } \partial \Omega, \end{cases}$$

其中$\Omega$是$(x,y)$平面上的連通開區域，它的邊界$\partial \Omega$是良好的（例如，光滑流形或多邊形）, $u_{xx}$和$u_{yy}$分別表示$x$和$y$的二階導數。問題P1能夠通過計算不定積分而直接解決。然而，解決邊值問題的這一方法只有在空間維數為1時才可用，並且不能推廣到高維問題以及形如$u+u''=f$的問題。出於這種考慮，我們將用有限元方法解決P1並將其推廣至問題P2.

我們的描述分為兩步，每步都反映了用有限元解決邊值問題的本質。

- 將原問題描述為它的弱形式，或變分形式。這一步很少或不需要計算。
- 離散化，將弱形式在有限維空間離散化。

這兩步之後，我們可以構造一個大型有限維線性方程，線性方程的解就是原邊值問題的逼近解。然後，這一有限維問題由電腦求解。

### 弱解形式

#### P1的弱解形式

第一步是將問題P1和P2轉化為他的等價變分形式，或弱解形式。

如果$u$是問題P1的解，那麼對任何滿足邊界條件的光滑函數$v$ ,有

（1）$\int_0^1 f (x)v (x) \, \mathrm{d}x = \int_0^1 u''(x)v (x) \, \mathrm{d}x$

相反如果 $u$對任何光滑函數$v (x)$滿足$u(0)=u (1)=0$和（1）。

那麼$u$是P1的解。對於二次可導函數$u$證明這一點是非常容易的（利用中值定理）。

通過對（1）的右側使用分部積分，可以得到

（2）$\begin{align} <pre> \int_0^1 f (x)v (x) \, \mathrm{d}x & = \int_0^1 u''(x)v (x) \, \mathrm{d}x \\ & = u'(x)v (x)|_0^1-\int_0^1 u'(x)v'(x)\, \mathrm{d}x \\ & = -\int_0^1 u'(x)v'(x)\, \mathrm{d}x = -\phi(u,v) </pre> \end{align}$

其中假設$v (0)=v (1)=0$。

#### P2的弱解形式

f當我們使用格林恆等式來表示式（2）， P2可以 $u$的積分型式表示，在此定義 $\phi(u,v)$

$$\int_\Omega fv\,ds = -\int_\Omega \nabla u \cdot \nabla v \, ds \equiv -\phi(u,v),$$

此處 $\nabla$代表梯度，即為二維平面上的內積。另外 $\,\!\phi$可以轉為內積空間 $H_0^1(\Omega)$，且 $\Omega$的一次微分函數 $\partial \Omega$為零。我們也可以假設 $v \in H_0^1(\Omega)$ （詳見索伯列夫空間）也可以顯示解的存在性和唯一性。

## 離散化

P1 和 P2 通過上述過程被離散化，並簡化為 子問題 (3)。基本思路是將無限維線性問題替換掉：

找到 $u \in H_0^1$ 使

$$\forall v \in H_0^1, \; -\phi(u,v)=\int fv$$

表示唯有限維度的形式：

子問題(3) Find $u \in V$ such that

$$\forall v \in V, \; -\phi(u,v)=\int fv$$

此處 $V$ 是$H_0^1$的一個有限維度線性子空間。

$V$有許多可能的形式，但對於有限元素而言，在此將假定$V$存在於分段多項式函數的空間中。

### 對於 P1

在此在區間 $(0,1)$之中選擇 $n$個$x$的可能值 $0=x_0<x_1<\cdots<x_n<x_{n+1}=1$ 接著定義$V$ 為：

$$V=\{v:[0,1] \rightarrow \mathbb R\;: v\mbox{ is continuous, }v|_{[x_k,x_{k+1}]} \mbox{ is linear for } k=0,\dots,n \mbox{, and } v(0)=v(1)=0 \}$$

令e $x_0=0$ 和 $x_{n+1}=1$。觀察到在 $V$ 之中的函數根據微積分的基本定義是不可微分的。當然，當 $v \in V$，則通常不定義 $x=x_k$，$k=1,\ldots,n$的導數，但導數事實上存在於每一個$x$的位置，並可以利用這些導數來進行部分積分運算。

### 對於 P2

$V$ 是屬於 $\Omega$的一系列函數。在右圖中，圖片下半部是一個15邊形的平面$\Omega$的三角分割，以及該多邊形的分段線性函數（圖片上半部彩色部份），即 在三角分割所形成的每個三角形上呈線性；空間$V$則由在所在的三角分割的每個三角形上的函數線性組合而成。

我們希望當下面的三角形網格變得越來越精細，離散子問題（3）的解在某種意義上將收斂到原始邊界值問題P2的解。為了測量此網格的細度，三角分割由一很小的實數$h > 0$所表示。此參數將與三角分割中最大或平均三角形的大小有關。當我們提高三角分割的精度時（分割出更多三角形），分段線性函數的空間$V$應會隨$h$變動。因此，在某些文獻中會以$V_h$來代表。

## 外部連結

- MIT Video Lecture on the Finite Element Method
- Multiphysics Glossary （Glossary of Multiphysics and Finite Element Modeling terms by COMSOL）
- NAFEMS—The International Association for the Engineering Analysis Community
- IFER —Internet Finite Element Resources - an annotated list of FEA links and programs
- Workshop "The Finite Element Method in Biomedical Engineering, Biomechanics and Related Fields"
- Finite Element Analysis Resources- Finite Element news, articles and tips
- COMSOL Multiphysics Finite Element Analysis Software - Official site
- CAD, Finite Element Analysis（Abaqus,Ansys）, CAE, Programming - FEM, CAD, Programming, discussion forums
- Finite Element Books- books bibliography
- Mathematics of the Finite Element Method
- Finite Element Methods for Partial Differential Equations
- FEM AVI-gallery at CompMechLab site, St.Petersburg State Polytechnical University, Russia
- Intro to FEA
- Introduction to FEA for EM modeling（includes list of currently available software）
- Finite Element modeling of light feapower
- propagation
- World Association of Fatigue, Durability and Fracture Mechanics - Fatigue for Finite Element Models