有限單元法(The Finite Element Method)

簡介

       有限元法最初的概念源自用結構力學的方法解決彈性力學的問題，后來擴展到各種用微分方程描述的學科。最初的有限單元法多用于工程計算，隨著計算機性能的不斷增強以及對圖形表現(xiàn)真實性的進一步需求，有限單元法也逐漸應用到了圖形的領域。

連續(xù)介質力學 Continuum Mechanics

       在圖形學中，可變形的物體通常建模為連續(xù)介質的三維物體，其中最重要的三個量是：位移，應力，應變。

         對于一維的情況，如下圖：

      由Hooke's law  很容易得到：

其中E是楊氏模量，描述的是材料的剛度。引入應力和應變，則式子可簡化為：

下面來定義一些符號。

一個三維可變形物體主要由一個未變行狀態(tài)的體（undeformed shape, rest shape,initial shape）和一些材料屬性來定義.Rest shape為一個連續(xù)的三維空間子集，記為，其中的點稱為點的材質坐標（material coordinates）。

當有力施加到rest shape 上的時候，x的位置移動到 p(x),因為 p(x) 定義為所有材質上的點，那么 p(x)為上的一個向量空間，記向量空間 u(x) = p(x) - x 為位移域（displacement field）.

應變 Strain

為了在三維空間中使用Hooke 定律，需要用新的方法來表示。在大部分情況下，變形體內的應變都不是一樣的。比如一個一端固定的懸臂梁，上側是受拉，下側則是受壓。則應變其實是材料坐標的函數(shù)，記為，如果物體未發(fā)生形變，應變就是0，單純的空間位移無法引起形變，所以應變也為0，因此應力其實是由位移域決定的。在三維空間中，位移域有三個部分：，每一個分量都可以對x,y,z方向微分。應變就不是一個單值了，這個非常重要，因為在一個三維物體中的一個點，可以在一個方向受壓的同時，在另一個方向受拉。

應變可以表示為一個 3*3 空間上的張量（tensor）：

有兩種從位移域計算應變張量的方法：

前者稱為格林非線性應變張量，后者是科希線性應變張量?？葡埩可倭撕竺娴亩尾糠?，無法獲取旋轉量。位移域的梯度矩陣如下：

逗號后面的字母表示對其你偏導，求出梯度矩陣后很容易得到應變張量。

應力 Stress

應力指的是單位面積的受力，和應變一樣，應力也可用一個3階張量表示：

求得應力之后，想要求得某個面上的手里，則還需要知道這個面的法向，然后通過下面的公式求得：

本構關系

本構關系表示的應力和應變之間的關系。在三維空間中，Hooke 定律表示為：

在這里應力和應變都是3階的張量，對與各向同性的材料（isotropic materials）,Hooke  定律有：

這里的 v 是楊氏模量，取值范圍是[0...1/2), 描述固體材料抵抗形變能力的物理量。

動力方程 Equation of Motion

將牛頓第二定律運用到彈性材料中的質點，有

由于我們不知道整個物體的質量，對與一個單位體積的立方體dxdydz，化為單位體積質量（也就是密度）的形式：

這里的 f(x) 指的是作用在x的 內力+外力合。

現(xiàn)在要做的是計算內力。對于一個變形體內的一個單位立方體dxdydz，垂直X方向的應力如下

應力是單位面積的力，則應力乘以面積就可以得到內力，再除以體積則得到x方向內力合為：

逗號表示 空間導數(shù)。

則由內部應力產生的內力合可表示為：

則完整的偏微分方程可以化為：

通過這個方程，我們就可以通過求出內力外力和，繼而求出加速度，然后更新位置。

整理一下這個求解方法的pipeline：

1.密度和外力已知；

2.定義材料坐標 X 和世界坐標 P；

3.計算位移域 u =p-x;

4.計算位移域梯度，然后應變張量就可以獲得（格林應變張量或者柯希應變張量）；

5.根據(jù)Hooke's  law 計算應力域；

6.根據(jù)牛頓第二定律計算出加速度；

7.利用顯示或者隱式積分，更新P(x)。

這里還有兩個概念需要說明一下。

材料線性：材料的應力和應變關系符合簡單的Hooken law.(線性關系)

幾何線性：應變張量使用的是柯希張量。

材料線性+集合線性得到的結果是動力學偏微分方程線性，偏微分方程線性的好處就是好求解。在deformate object 模擬中，大部分都是假設成線性的，這樣能夠使計算簡單，獲得更快的計算速度，但是對于大變形或者旋轉，集合上的線性簡化會引起一些很怪異的現(xiàn)象，這在后面會介紹一些解決的方法。

有限單元離散

        在上面得到的偏微分方程中，p是一個連續(xù)向量空間，怎么描述物體中點位移會非常困難，因為物體有無限的質點，而我們需要求解解析解。

        有一個方法可以無限近似的方式逼近真實解，這個方法的核心思想就是Finite Element Method（FEM），在力學中稱為有限單元法。通過一些算法，將deformate object 的體mesh化分為網格單元mesh，網格單元之間不會有重疊，對于每一個單元，向量空間都可以根據(jù)單元頂點的位置寫出解析式。

       網格剖分的相關的文章：專注網格剖分 - TetGen,NETGEN,Steller

      一個球體剖分的例子：

剖分后：

常應變四面體網格 Constant Strain Tetrahedral Meshes

        在有限單元法中，最簡單的方法就是用四面體作為有限單元，這樣的網格mesh稱為四面體mesh。

        在單元中，如果形變域為一個常量的話會更加簡單，但是這樣會導致應力為0，然后就無法模擬變形物體了。

對于如上圖的一個四面體單元，x0,x1,x2,x3是在完全不受外力情況下的位置（rest shape），在形變狀態(tài)下的位置為p0,p1,p2,p3，單純的平移變換并不會產生內力，所以假設x0 = p0 = 0；

對與四面體中的某一點x,可以用重心坐標插值來表示：

變形之后，這一點的重心坐標是不變的，可以記為：

兩式連立，消去b，得到

這是一個線性關系，[x1,x2.x3]矩陣的逆可以事先求出。

由于p(x)是線性的，所以有

又位移域 u = p-x，則

使用格林應變張量

再利用Hooken law，求出應力，然后乘以面積就可以得到內力。

例如，對于面(0,1,2)面上方的力如下：

三角形的兩個向量叉乘得到三角形的面積。

最后，我們只要吧面力離散到三角形的各個頂點就可以了。

就如我們看到的，整個過程都非常簡單，只是根據(jù)之前的一個狀態(tài)不斷地去計算各種變量，同時我們可以很輕易的將Hooken law 替換成其他的非線性模型。

下面給出整個模擬過程的偽代碼.

一行行解釋。

1-4：對有限元mesh中的所有單元進行初始化，初始位移和rest shape重合。

5-7：預處理每個單元的Xi；

8-11：這里開始主模擬循環(huán)。每次循環(huán)開始都要重新設置頂點上的力，初始設置為外力和；

13-14：計算出位移梯度；

15-16：用格林應變張量計算出應力；

17-22：計算單元中每個面上的力，然后將力離散到面上的每一個頂點；

24-26：到這一步，單元上的上的每一個頂點上的力都已經搞定了，接下來是更新每個點的位置；

28：更新顯示。

算法的流程很像Mass - Spring - System. 雖然力的計算很耗時，但算法并沒有因此更難理解或實現(xiàn)。

對于更穩(wěn)定的模擬，一定要使用隱式積分。

參考

SIGGRAPH 2008 Course - Real Time Physics Class Notes

SIGGRAPH 2012 Course - FEM Simulation of 3D Deformable Solids: A practitioner’s guide to theory, discretization and model reduction