溫馨提示×
您好，登錄后才能下訂單哦！
點擊 登錄注冊 即表示同意《億速云用戶服務條款》
您好，登錄后才能下訂單哦！
有限元法最初的概念源自用結構力學的方法解決彈性力學的問題，后來擴展到各種用微分方程描述的學科。最初的有限單元法多用于工程計算，隨著計算機性能的不斷增強以及對圖形表現真實性的進一步需求，有限單元法也逐漸應用到了圖形的領域。
在圖形學中，可變形的物體通常建模為連續介質的三維物體，其中最重要的三個量是：位移，應力，應變。
對于一維的情況，如下圖：
由Hooke's law 很容易得到：
其中E是楊氏模量，描述的是材料的剛度。引入應力
下面來定義一些符號。
一個三維可變形物體主要由一個未變行狀態的體（undeformed shape, rest shape,initial shape）和一些材料屬性來定義.Rest shape為一個連續的三維空間子集，記為
當有力施加到rest shape 上的時候，x的位置移動到 p(x),因為 p(x) 定義為所有材質上的點，那么 p(x)為
為了在三維空間中使用Hooke 定律，需要用新的方法來表示。在大部分情況下，變形體內的應變都不是一樣的。比如一個一端固定的懸臂梁，上側是受拉，下側則是受壓。則應變其實是材料坐標的函數，記為
應變可以表示為一個 3*3 空間上的張量（tensor）：
有兩種從位移域計算應變張量的方法：
前者稱為格林非線性應變張量，后者是科希線性應變張量。科希張量少了后面的二次部分，無法獲取旋轉量。位移域的梯度矩陣如下：
逗號后面的字母表示對其你偏導，求出梯度矩陣后很容易得到應變張量。
應力 Stress
應力指的是單位面積的受力，和應變一樣，應力也可用一個3階張量表示：
求得應力之后，想要求得某個面上的手里，則還需要知道這個面的法向，然后通過下面的公式求得：
本構關系
本構關系表示的應力和應變之間的關系。在三維空間中，Hooke 定律表示為：
在這里應力和應變都是3階的張量，對與各向同性的材料（isotropic materials）,Hooke 定律有：
這里的 v 是楊氏模量，取值范圍是[0...1/2), 描述固體材料抵抗形變能力的物理量。
將牛頓第二定律運用到彈性材料中的質點，有
由于我們不知道整個物體的質量，對與一個單位體積的立方體dxdydz，化為單位體積質量（也就是密度）的形式：
這里的 f(x) 指的是作用在x的 內力+外力合。
現在要做的是計算內力。對于一個變形體內的一個單位立方體dxdydz，垂直X方向的應力如下
應力是單位面積的力，則應力乘以面積就可以得到內力，再除以體積則得到x方向內力合為：
逗號表示 空間導數。
則由內部應力產生的內力合可表示為：
則完整的偏微分方程可以化為：
通過這個方程，我們就可以通過求出內力外力和，繼而求出加速度，然后更新位置。
整理一下這個求解方法的pipeline：
1.密度和外力已知；
2.定義材料坐標 X 和世界坐標 P；
3.計算位移域 u =p-x;
4.計算位移域梯度，然后應變張量就可以獲得（格林應變張量或者柯希應變張量）；
5.根據Hooke's law 計算應力域；
6.根據牛頓第二定律計算出加速度；
7.利用顯示或者隱式積分，更新P(x)。
這里還有兩個概念需要說明一下。
材料線性：材料的應力和應變關系符合簡單的Hooken law.(線性關系)
幾何線性：應變張量使用的是柯希張量。
材料線性+集合線性得到的結果是動力學偏微分方程線性，偏微分方程線性的好處就是好求解。在deformate object 模擬中，大部分都是假設成線性的，這樣能夠使計算簡單，獲得更快的計算速度，但是對于大變形或者旋轉，集合上的線性簡化會引起一些很怪異的現象，這在后面會介紹一些解決的方法。
在上面得到的偏微分方程中，p是一個連續向量空間，怎么描述物體中點位移會非常困難，因為物體有無限的質點，而我們需要求解解析解。
有一個方法可以無限近似的方式逼近真實解，這個方法的核心思想就是Finite Element Method（FEM），在力學中稱為有限單元法。通過一些算法，將deformate object 的體mesh化分為網格單元mesh，網格單元之間不會有重疊，對于每一個單元，向量空間都可以根據單元頂點的位置寫出解析式。
網格剖分的相關的文章：專注網格剖分 - TetGen,NETGEN,Steller
一個球體剖分的例子：
剖分后：
在有限單元法中，最簡單的方法就是用四面體作為有限單元，這樣的網格mesh稱為四面體mesh。
在單元中，如果形變域為一個常量的話會更加簡單，但是這樣會導致應力為0，然后就無法模擬變形物體了。
對于如上圖的一個四面體單元，x0,x1,x2,x3是在完全不受外力情況下的位置（rest shape），在形變狀態下的位置為p0,p1,p2,p3，單純的平移變換并不會產生內力，所以假設x0 = p0 = 0；
對與四面體中的某一點x,可以用重心坐標插值來表示：
變形之后，這一點的重心坐標是不變的，可以記為：
兩式連立，消去b，得到
這是一個線性關系，[x1,x2.x3]矩陣的逆可以事先求出。
由于p(x)是線性的，所以有
又位移域 u = p-x，則
使用格林應變張量
再利用Hooken law，求出應力，然后乘以面積就可以得到內力。
例如，對于面(0,1,2)面上方的力如下：
三角形的兩個向量叉乘得到三角形的面積。
最后，我們只要吧面力離散到三角形的各個頂點就可以了。
就如我們看到的，整個過程都非常簡單，只是根據之前的一個狀態不斷地去計算各種變量，同時我們可以很輕易的將Hooken law 替換成其他的非線性模型。
下面給出整個模擬過程的偽代碼.
一行行解釋。
1-4：對有限元mesh中的所有單元進行初始化，初始位移和rest shape重合。
5-7：預處理每個單元的Xi；
8-11：這里開始主模擬循環。每次循環開始都要重新設置頂點上的力，初始設置為外力和；
13-14：計算出位移梯度；
15-16：用格林應變張量計算出應力；
17-22：計算單元中每個面上的力，然后將力離散到面上的每一個頂點；
24-26：到這一步，單元上的上的每一個頂點上的力都已經搞定了，接下來是更新每個點的位置；
28：更新顯示。
算法的流程很像Mass - Spring - System. 雖然力的計算很耗時，但算法并沒有因此更難理解或實現。
對于更穩定的模擬，一定要使用隱式積分。
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
免責聲明：本站發布的內容（圖片、視頻和文字）以原創、轉載和分享為主，文章觀點不代表本網站立場，如果涉及侵權請聯系站長郵箱：is@yisu.com進行舉報，并提供相關證據，一經查實，將立刻刪除涉嫌侵權內容。
