1.1 什么是有限元方法

有限元方法(Finite Element Method, FEM)是一種求解由偏微分方程描述或可表示為泛函極小化問題的數值方法。感興趣的域被表示為有限單元(finite elements)的集合。有限元中的逼近函數是根據所求物理場的節點值確定的。FEM將一個連續的物理問題轉化為節點值未知的離散化有限元問題，并得到一個線性方程組，求解該方程組就可以獲得待求的物理量。有限元內部的值可以使用節點值恢復。

值得一提的是，FEM的兩個特點：

在有限元上的物理場的分段近似提供了很好的精度，即使用簡單的近似函數(增加單元的數量，我們可以達到任意精度)。

局部逼近將導致離散問題的方程組稀疏，這有助于解決具有大量節點未知數的問題。

1.2 FEM是如何工作的

為了說明有限元方法是如何工作的，我們列出如下幾個有限元的主要求解步驟

離散化。第一步是將求解域劃分為有限單元。有限元網格通常由預處理程序生成。網格的描述由幾個數組組成，其中主要是節點坐標和單元的連通性。

選擇插值函數。插值函數用于對單元上的場變量進行插值，這里通常選擇多項式作為插值函數。多項式的次數取決于單元上的節點數。

尋找單元屬性。建立待求函數的結點值與其他參數聯系起來的矩陣方程，此步驟有多種方法可以使用，但最方便的方法是伽遼金方法(Galerkin method)。

組裝單元方程。要找到整個求解區域的全局方程組，必須將所有的單元方程組合在一起。在組裝過程中需要用到單元之間的連通性。同時，在求解之前還需要施加邊界條件（邊界條件不在單元方程中考慮）。

求解全局方程組。有限元整體方程組具有典型的稀疏、對稱和正定性質，可采用直接法和迭代法來求解該方程組。待求函數在結點處的值被當做最終的解。

計算其他量。對于很多情況，我們需要計算其他的量，例如在力學問題中，除了位移外，還需要考慮應變和應力，這些位移是在求解整體方程組后得到的。

1.3 有限元方程的表述

有很多方法可以將問題的物理形式轉換成有限元離散后的形式。如果物理問題是由偏微分方程描述的，則經常使用的是伽遼金方法(Galerkin method)。如果文理問題是由泛函最小化描述的，則經常使用變分公式(variational formulation)將其轉換為有限元方程。

1.3.1 伽遼金方法(Galerkin method)

讓我們先從簡單的一維問題開始，介紹使用伽遼金方法的有限元公式。假設我們需要求解下面方程的數值解
$$a\frac{d^{2}u}{d{x}^2}+b=0,0\le x\le 1\text{\hspace{1em}(1)}$$
邊界條件為
$$\begin{array}{l}{u|}_{x=0}=0\\ {a\frac{du}{dx}|}_{x=2L}=R\end{array}$$
其中$u$是待求解量。我們將使用下圖所示的兩個線性一維有限元來解決這個問題。

首先我們考慮圖1.1右半部分的有限元，該單元擁有兩個節點，則近似函數$u(x)$可以表示為
$$\begin{array}{l}u=N_1{u}_1+N_2{u}_2=[N]\{u\}\\ [N]=\left[\begin{array}{cc}{N}_1& N_2\end{array}\right]\\ \{u\}=\left\{\begin{array}{cc}{u}_1& u_2\end{array}\right\}\end{array}$$
這里的$N_i$也稱為形函數(shape functions)
$$\begin{gathered} N_1=1?\frac{x?x_1}{x_2?x_1} \\ {N}_2=\frac{x?x_1}{x_2?x_1} \end{gathered}$$
使用這些形函數在結點出的值來插值$u(x)$。我們需要從全局離散方程中求解出未知的結點值$u_1$和$u_2$。

將節點處的值和形函數帶入到式(1)的微分方程中，得到如下的近似形式
$$a\frac{d^2}{d{x}^2}[N]\{u\}+b=\psi \text{\hspace{1em}(2)}$$
由于式(2)是有限元中的近似形式，所以殘差$\psi$是一個非零的數。伽遼金方法通過將式(2)的項乘以形函數，在單元上積分并令其等于0最小化殘差：
$${\int }_{x_1}^{x_2}[N{]}^{T}a\frac{d^2}{d{x}^2}[N]\{u\}dx+{\int }_{x_1}^{x_2}[N{]}^{T}bdx=0$$
對其進行分部積分，得到該單元上微分方程的一種離散形式
$${\int }_{x_1}^{x_2}{\left[\frac{dN}{dx}\right]}^{T}a\left[\frac{dN}{dx}\right]dx\{u\}?{\int }_{x_1}^{x_2}[N{]}^{T}bdx?\left\{\begin{array}{c}0\\ 1\end{array}\right\}a{\frac{du}{dx}|}_{x=x_2}+\left\{\begin{array}{c}1\\ 0\end{array}\right\}a{\frac{du}{dx}|}_{x=x_1}=0\text{\hspace{1em}(3)}$$
為了表示方便，我們使用如下幾個記號

$$\begin{array}{rl}& [k]={\int }_{x_1}^{x_2}{\left[\frac{dN}{dx}\right]}^{T}a\left[\frac{dN}{dx}\right]dx\\ & \{f\}={\int }_{x_1}^{x_2}[N{]}^{T}bdx+{\left\{\begin{array}{l}0\\ 1\end{array}\right\}a\frac{du}{dx}|}_{x=x_2}?{\left\{\begin{array}{l}1\\ 0\end{array}\right\}a\frac{du}{dx}|}_{x=x_1}\end{array}$$

則式子(3)可以寫成
$$\boxed{[k]\{u\}=\{f\}}$$
在固體力學中，$[k]$稱為剛度矩陣(stiffness matrix)；$\{f\}$稱為載荷向量(load vector)。在上述我們考慮的偏微分問題中，我們分別計算兩個長度為$L$單元的剛度矩陣和載荷向量。
$$[k_1]=[k_2]=\frac{a}{L}\left[\begin{array}{cc}1& ?1\\ ?1& 1\end{array}\right]$$

$$\{f_1\}=\frac{bL}{2}\left\{\begin{array}{c}1\\ 1\end{array}\right\},\{f_2\}=\frac{bL}{2}\left\{\begin{array}{c}1\\ 1\end{array}\right\}+\left\{\begin{array}{c}0\\ R\end{array}\right\}$$

對于整個求解域上2個單元，3個結點的全局方程組可以通過組裝單元方程組獲得。在此實例中，第二個結點時單元共用結點。組裝全局方程組為
$$\frac{a}{L}?\begin{array}{ccc}1& ?1& 0\\ ?1& 2& ?1\\ 0& ?1& 1\end{array}??\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}?=\frac{bL}{2}?\begin{array}{c}1\\ 2\\ 1\end{array}?+?\begin{array}{c}0\\ 0\\ R\end{array}?$$
施加邊界條件$u(x=0)=0$后，最終的全局方程為
$$\frac{a}{L}?\begin{array}{ccc}1& 0& 0\\ 0& 2& ?1\\ 0& ?1& 1\end{array}??\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}?=\frac{bL}{2}?\begin{array}{c}0\\ 2\\ 1\end{array}?+?\begin{array}{c}0\\ 0\\ R\end{array}?$$
結點處的值$u_i$可以通過求解上面的線性方程組得到。$u$在有限元內任一點處的值可以通過形函數計算。當$a=1,b=1,L=1,R=1$時有限元的解如下圖所示。

在此例中，精確解是二次函數，而有限元方法中在每個單元上使用的是簡單的線性函數。通過增加單元的個數或使用更加復雜的形函數可以得到更加精確的數值解。值得注意的是，在結點處有限元方法得到的解是精確的（僅針對此例中）。如果精確解是二次函數，則有限元方法使用線性形函數可在結點處獲得精確值；而如果精確解是三次函數，則有限元方法使用二次形函數可在結點處獲得精確值，以此類推。

總結

以上是生活随笔為你收集整理的有限元分析简介及伽辽金法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。