有限元分析（FEA）是個什麼東東

一、有限元能幹什麼

通常能通過第二牛頓定律能得出這樣的方程：

我們把這種方程稱之為 常微分方程 ，即 自變數只和時間 有關係，和空間位置沒關係 ，換句話說，常微分方程描述的是單質點的變化規律，如：某個物體在重力作用下做自由落體運動，下落距離隨時間變化的規律；火箭在發動機推動下在空間飛行，飛行的軌道等等。常微分方程一般是把研究物件當成一個質點或者剛體，研究整體的運動規律。

值得高興的是，線性常微分方程相對比較好解決，我們可以通過 傅立葉變換 （ J Pan：傅立葉變換後面的到底有什麼小祕密 ）和 拉普拉斯變換 （ J Pan：從另一個角度看拉普拉斯變換 ）將微分方程變成代數方程，進而得出準確的 解析解 。

事實證明，對於繁複紛雜的大自然，只用常微分方程是不夠的，因為有些研究物件不能簡化成質點，一個典型的例子就是琴絃：

琴絃是一個柔性體，在彈撥的時候每個點振動都是不一樣的，一根均勻的弦，其自由振動的數學方程是什麼呢？假定 表示 點在 時刻的位移，如果取琴絃中一個微元進行力學分析，就可以得到琴絃應遵守的方程為：

可以發現，這個方程的 自變數不僅僅是時間，還有空間座標 ，我們把這種方程稱之為偏微分方程。也就是說，偏微分方程能描述的連續體的各個點隨時間連續變化的情況，本質上就是一種“場”的描述了。

嚴格來說， 自然界的各種現象，都需要用偏微分方程來描述 ，只不過有些情況下為簡單起見，我們忽略了位置項帶來的影響，如火箭發射的時候，在強烈的振動下殼體會發生變形的，理論上要完整描述火箭的狀態，是需要用偏微分的，但是實際上我們還是把它簡化成一個質點來描述，因為這樣可以簡化很多工作。

當然，這個琴絃除了滿足上面的方程，還有其它約束，比如說琴絃的兩端要固定在支架上，琴絃內部與外界是通過這個邊界聯絡起來的，邊界的變化是會影響琴絃的波動的，我們把這些稱之為 邊界條件 ，比對於琴絃，邊界條件一般為：

； ；

除了邊界，琴絃的波動顯然還和它的“歷史”有關。兩根同樣的弦，一根在薄刀背的敲擊下發生的聲音就比較刺耳，另一根在手指的彈撥下就比較和諧，兩根弦由於初始時刻的振動情況不一樣，後來振動下發出來的聲音就不一樣，我們把初始時刻的狀態稱之為 初始條件 ： ； ；

通常邊界條件（含邊界條件和初始條件）有三類，一類是變數本身的約束，如 ，稱之為狄利克雷邊界條件；一類是變數的導數有約束，如 ，稱之為諾依曼邊界條件；還有第三類就是前兩者的混合，成為羅賓邊界條件。

遺憾的是，對於偏微分方程，除了少數極簡單的情況下能計算出解析解（形式還很複雜），絕大多數情況下我們還無法很好的求解——有限元就是在這種情況下應運而生，簡要來說， 有限元法就是為了求解偏微分方程！

二、偏微分方程如何求解（有限元發明以前）

前面我們說了，對於大多數偏微分方程而言，我們是沒有辦法求得解析解的，只能退而求其次，尋求個近似解，即便是這個近似解的過程也是相當不容易的，方法也很簡單、粗暴—— 試湊法 ！比如我們現在我們假定一個位移函式，其中包含若干個待定係數：

當然這些函式也不是隨意假定的，它需要滿足位移邊界條件，這個靠 來滿足， 主要用來模擬非邊界條件的點，它們在邊界處取值應該為零，否則和 疊加後會破壞位移邊界條件。 為待定係數。。

在給定的外力作用下，在滿足位移邊界條件的各組位移中，實際存在的一組位移應使總勢能為最小，因此該式又稱為 最小勢能原理 。

系統在外力的作用下，發生位移，產生變形。位移可以是各種各樣的，但必須滿足位移的邊界條件。但是，滿足位移邊界條件的位移位移有無窮多組，我們應該選哪一組呢？——那就是當 總勢能取極小值 時的位移，這種現象，我們稱之為 最小作用量原理 ，這是自然界中一個普遍成立的原理。

實際我們是怎麼做的呢？我們知道，不同的 組合就會獲得不同的位移，這樣我們就可以列出總勢能的表示式，其中只包含待定係數 ，由前面我們說的最小勢能原理， 只有一組是真實的解，那就是當總勢能變分 為零的時候，怎麼做到的呢？很簡單，將 分別對 求導取零就行了，這個時候偏微分方程組就變成線性方程組了，這種解偏微分的方法稱之為利茲法。

如果我們所取的位移不僅滿足位移邊界條件，而且根據它們求得的應力還滿足應力邊界條件(不要求滿足平衡方程），這種方法稱之為伽遼金方法，這種方法對位移函式的要求較高，但計算量小一些。

上面說的都太抽象，我們舉個具體的例子，注意下面的例子不用看懂，只是為了演示傳統求解偏微分方程是多麼費勁！

我們舉一個利茲法的例子，可能大家會更容易理解一些：兩端簡支的等截面樑，受均勻分佈載荷 作用如圖所示，試求解樑的撓度 。

首先構造位移試函式

顯然這個試函式是滿足位移邊界條件

； ；

則總的勢能為：

整理得到

根據 ，得到

所以

再次強調，這個例子看思路即可！

三、有限元法的基本思路

根據第二節的描述，我們知道傳統解決偏微方程的方法（如利茲法），採用的試函式 非常複雜，而且是在全域定義的，工程上的物體（零元件或者系統）及其邊界條件都太複雜了，我們幾乎不能夠構建一個很好的基於整個彈性體（全域）的位移函式。

那怎麼辦呢？——還記得圓周率是怎麼計算的嗎？其中有一種幾何方法是這樣的，就是用正多邊形等效圓形，就像切西瓜一樣，將圓切成有限個等腰三角形，示意圖如下：

每個等腰三角形的面積

由所有的等腰三角形組成的正多邊形的面積為：

我們知道， ，所以當 時， 正對變形就趨向於圓形了，此時可獲得圓的面積為 。也就是說我們可以採用“ 化整為零 ”的思想，將複雜不易求解的東西切成簡單的容易計算的幾何形狀來進行等效。

也就是說，當整體曲面或者曲線較為複雜時，我們就可以把它切開，切開的每一個規則的小塊稱之為單元，單元與單元之間通過節點聯絡起來，整個彈性體就被劃分成了有限個單元，簡稱有限元。對一個規則的單元（子域）再假設位移函式就簡單多了，比如可以採用線性的試函式 ，試函式的要求也不高（連續性階次較低），缺點就是計算量變大了，不過我們現在有計算機了，這都不是事。

簡而言之， 有限元分析就是“化質為量”，就是將解微分方程這個難題（質）轉化為大量的數值計算（量）。

四、有限元法的數學基礎——降維

很多童鞋都讀過劉慈欣的《三體》，裡面提到了一個很火的概念，叫做“ 降維打擊 ”，什麼意思呢？就是將攻擊目標本身所處的空間維度降低，致使目標無法在低維度的空間中生存從而毀滅目標。其實“降維”這種做法我們早就使用了，比如地圖就是典型的降維的結果，因為我們關心的是經緯度、不關心每塊地區的海拔高度，這樣就可以將三維的地球資訊用二維的地圖表示。

有限元分析其實也就是降維分析，將“無限維”降低到“有限維”進行計算。以上都是抽象的描述，接下來我們看看在數學上怎麼描述。舉個簡單的例子，在每個節點，我們選擇線性分段函式作為試函式 ，它的特點是在節點 處取值為1，在 ( )處取0，其他位置線性變化，形狀曲線如下：

數學表示式如下：

因為 長得像個帽子，所以一般稱之為hat function。這些函式有個特點就是“作用範圍”很小，它們只在節點周圍的若干個單元上有值，其它地方一概為0。接下來要引入一個非常重要的概念： 作為基函式，可以張成空間 ， 中任一函式都可以表示成：

什麼意思，大家估計都玩過積木，基礎模組只有幾種，通過各種不同的空間組合，卻能拼出各種房子、動物等等。

前面我們說的 就可以看成是積木中的基礎模組，我們稱之為基函式；積木搭建的各種房子、動物等，我們稱之為基函式張成的空間 ； 的線性組合就是 中一個函式向量。這個點我們要多說一下，因為對於沒有接觸過泛函分析的人，確實不太容易理解。我們都很熟悉泰勒分析：

這個方程啥意思呢？——多項式 為基函式組成一個空間，空間裡的任意函式（足夠光滑）都可以由這些基函式線性疊加而來。再比如傅立葉分解也是這個意思：任意周期函式（滿足狄利克雷條件），都可由正餘弦函式疊加而來，而係數就是被分解函式與基函式的內積，也就是被分解函式在基函式上的投影！

說到投影大家都很熟悉了，比如X光透射，就是光透過三維物體留下的二維影像！

那函式 是不是也可以通過計算其在空間 的投影來計算投影后係數 呢？——當然是可以的，比如我們可以這麼計算：

其中， 為 在空間 的投影， 為空間 任意函式向量。既然 與 中所有的函式 正交，則說明 是 在空間 中的投影，具體請看下圖：

由於 的任意性，為方便起見，我們可以選擇基函式 ，於是可以得到：

由於 在空間 中，我們可以假設

所以

這個方程看著太複雜，我們假設：

於是方程

就變成了

注意，這是一個 線性方程組， 為 個未知量，我們通過解這個線性方程組，就可以獲得基函式的係數 進而獲得函式 的近似。式中 稱為負載向量，由於歷史原因， 一般稱之為質量矩陣。

接下來我們就詳細看一下這個線性方程組具體長什麼樣。我們先畫出相鄰的兩個hat functions：

由前面的分析，採用辛普森積分法則，可以得到：

其中 ， ， 和 是half hat，因此需要單獨計算 ， ,同樣可以計算 如下：

同樣可以得到： 。

我們對質量矩陣進行組裝，可以得到：

MATLAB程式碼如下：

function M = MassMat1D(x)
n = length(x)-1; % number of subintervals
M = zeros(n+1,n+1); % allocate mass matrix
for i = 1:n % loop over subintervals
h = x(i+1) - x(i); % interval length
M(i,i) = M(i,i) + h/3; % add h/3 to M(i,i)
M(i,i+1) = M(i,i+1) + h/6;
M(i+1,i) = M(i+1,i) + h/6;
M(i+1,i+1) = M(i+1,i+1) + h/3;
end

採用梯形積分法則，可以可以計算負載向量 如下：

完整矩陣形式為：

MATLAB程式碼如下：

function b = LoadVec1D(x,f)
n = length(x)-1;
b = zeros(n+1,1);
for i = 1:n
h = x(i+1) - x(i);
b(i) = b(i) + f(x(i))*h/2;
b(i+1) = b(i+1) + f(x(i+1))*h/2;
end

現在我們來驗證一下演算法是否正確，假定 ,程式碼如下：

function y = Foo(x)
y=x.*sin(x)

我們來計算一下投影函式和原函式分別長什麼樣，程式碼如下：

n = 5 ;% number of subintervals
h = 1/n; % mesh size
x = 0:h:1; % mesh
y = Foo(x);% calc f(x)
M = MassMat1D(x); % assemble mass
b = LoadVec1D(x,@Foo) % assemble load
Pf = M\b; % solve linear system
plot(x,Pf,x,y); % plotprojection

計算結果如下，可見投影函式基本和原函式一致，但是誤差比較大，這是因為我們節點去太少了，減小 後重合度會更高。

五、如何獲得“弱形式”的解

前面我們計算了函式 的基函式空間的投影情況，我們說有限元是為了解決偏微分方程求解的，那 和偏微分方程又有什麼關係呢？這就牽扯到微分方程的“弱形式”了！ 獲得微分方程的“弱形式”解，可以說是有限元法最重要的概念之一 ，我們看一下是怎麼一回事！假設我們微分方程具有下列形式：

邊界條件為：

其中 和 為給定的函式， ， ， 和 是給定常數。具有給定熱源的一維熱傳導方程就具備這種形式。微分方程描述的是微觀狀態，但是我們一般更容易掌握的是巨集觀的現象，所以呢，方程我們一般喜歡把微分方程轉化成積分方程，比如對上式兩邊同時積分：

乍一看，皆大歡喜，實際上有問題，問題在哪呢？——以上方程表示在整個計算域 中， 方程的積分值相等，也就是 平均值相等 ，與原來要求處處相等相比，顯得“太弱了”！那怎麼辦呢？——我們可以找一個試函式 （任意的），將原方程改寫成如下方式：

當然， 不能表現太糟糕，至少要保證兩邊積分都是收斂的，如果上述方程對於所有滿足條件的 都成立，這個要求太強了，顯然當 時這個方程是成立的，那是不是隻有這種情況下才成立呢？——實際上就是的，這可以通過變分法證明。

下一個問題就來了，既然 有千千萬萬，那我們到底選擇哪一個？總不能都選擇吧！——前面我們說了，有限元就是為了解決偏微分方程的求解問題，思路就是將複雜的空間切成很多小區域，每個區域上假設一個試函式test function（或者叫形函式shape function） ，然後通過求取試函式的係數 來獲得函式 的近似。試函式有一個特點，它只在一個很小的區間上有數值，其他地方都為零。如果我們在每個點上都用試函式測試一遍，這樣至少能保證所單元上都成立的！——這種將試函式和形函式選擇為同一個函式的方法，是一個俄國人最先使用的，我們稱之為伽遼金（Galerkin）法。

這樣我們就得到了如下的方程：

這就是所謂的原微分方程的 弱形式 ！之所以“弱”，是因為這個方程的解只能保證在每個單元上方程左右兩邊的平均值相等，而不是原微分方程要求的每一點都相等。注意此時的 已經與原微分方程的解 有了差異，因為我們現在是用有限維空間去近似了原無限維空間。

接下來我們要繼續變形了，要稍微用點數學知識——不知道大家還記不得分部積分法：假設有函式 和 ，則很容易得到 乘積的微分為：

為二階無窮小，可略去，所以可以得到：

兩邊積分得到：

現在我們令 ， ，則有

所以我們可以得到：

通過這種變化，我們對 的要求由2階降到了1階，將形函式 從0階增加到了1階，由於形函式是我們自己設計的，這就降低了對原微分方程解的光滑性要求。將邊界條件代入上式，可得展開的微分方程弱形式：

我們可以假設

代入展開的微分方程弱形式就可以得到：

其中 和 為 矩陣， 和 為 向量，具體表達式如下：

這樣我們就將微分方程轉化成了線性矩陣方程：

其中 為 矩陣，稱之為剛度矩陣， 為負載矩陣，和前面的定義一樣。和前面的計算相似，可以獲得 的完整矩陣形式如下：

為邊界條件帶來的等效矩陣，顯然只有 或者 時 才有取值，很容易求得 ， ，即

所以可以得到：

MATLAB程式碼如下：

function A = StiffMat1D(x,kappa)
n = length(x)-1;
A = zeros(n+1,n+1);
for i = 1:n
h = x(i+1) - x(i);
A(i,i) = A(i,i) + 1/h; 
A(i,i+1) = A(i,i+1) - 1/h;
A(i+1,i) = A(i+1,i) - 1/h;
A(i+1,i+1) = A(i+1,i+1) + 1/h;
end 
A(1,1) = A(1,1) + kappa(1);
A(n+1,n+1) = A(n+1,n+1) + kappa(2);

很容易得到：

MATLAB程式碼如下：

function b = LoadVec1D(x,f,kappa,g)
n = length(x)-1;
b = zeros(n+1,1);
for i = 1:n
h = x(i+1) - x(i);
b(i) = b(i) + f(x(i))*h/2;
b(i+1) = b(i+1) + f(x(i+1))*h/2;
end
b(1) = b(1) + kappa(1)*g(1);
b(n+1) = b(n+1) + kappa(2)*g(2);

舉個例子：假設我們現在有一維導熱方程如下：

邊界條件為：

也就是左端是恆定溫度邊界條件，右端是絕熱邊界條件。對照之前我們理論推導，設定的邊界條件應該具有下列形式：

兩者貌似不一致，怎麼回事？——我們之所以用諾依曼邊界條件進行推導，是因為微分方程弱形式進行分部積分的時候會用到諾依曼條件，那假如是別的邊界條件怎麼辦呢？比如本例中的就有狄利克雷條件： ，這就需要一點數學技巧了！

我們知道，物理模型的量一般都是有界的，不會漫天飛，因此，我們可以認為 也是一個有限值，現在假如 非常大，比如說達到 ，那 要非常接近0才行，於是 ，所以只要我們設定 即可；

第二個邊界條件本身就是諾依曼條件，就簡單了，直接設定 即可；

下面我們通過MATLAB來看一下，熱源的程式碼為：

function y = Source(x)
y = 0.03*(x-6)ˆ4; % heat source

求解結果如下：

%%solve Poission equations 
h = 0.1; % mesh size
x = 2:h:8; % mesh
kappa = [1.e+6 0];
g = [-1 0];
A = StiffMat1D(x, kappa);
b = LoadVec1D(x, @Source, kappa, g);
U = A\b;
plot(x,U)

顯然， 在 是取值為-1，目視在 是導數為0，是符合我們邊界條件設定的。

六、二維、三維有限元計算

前面說的都是一維情況，其節點形函式一般為：

插值函式為：

對於二維幾何域，其形函式基本類似，如下圖所示：

選用了 和 兩個方向的線性函式，與一維的hat function類似，每個函式在點 上的值為 1，但在其他非相鄰點 上的值為零，此時插值函式變成了：

情況會複雜一點，但是基本思想和處理方式是一致的。當然，我們最常用的還是二維和三維的情況，對於二維和三維的線性函式，最常見的單元所示：

前面說的都是，一階插值，就是節點之間都是線性的，特點是簡單，缺點也有，那就是精度偏低，我們可以通過提高階數來提升精度，比如採用二階，此時單元的邊和麵都可以是彎曲的，如下圖：

因此，在有限元中，有兩種提高精度的方法，一種就是細化網格，此時會數值解更逼近真解，我們稱之為h-refinement，也就是h法；另一種就會保持網格不變，選擇更高階數的形函式，稱之為p-refinement，也就是p法。

最後說一點，人們對於有限元有兩種對立的看法，一種認為有限元無所不能，上算天文，下算地理，給個電腦就能工作，這種人一般沒有親自做過數值計算；另一種認為有限元就是個花架子，除了花花綠綠的圖片能裝點門面外，沒啥太大的參考價值，因為可能經歷了太多別人給的模擬結果還不如自己的經驗來的準確；兩邊都有道理，因為隨著商業CAE軟體的普及，大多人不需要怎麼培訓，任何模型給他，點選幾下滑鼠就能出結果，這會給人一種錯覺，認為模擬無所不能，這其實忽略了一個重要的問題，那就是 給出結果 和 給出正確的結果 是兩回事——而正確的結果是什麼，可能壓根就沒人知道！

真的想做好模擬，首先要知道碰到的問題數學模型是什麼，邊界條件時什麼，其中邊界條件往往更復雜，難以獲得，這就需要評估其中誤差有多大；其次，要知道數值計算的基本原理和過程，知道自己的模型輸入到電腦之後，資料是怎麼處理的，要有基本的預判。 模擬的作用不是復現現實，而是幫助工程師理解現實，模擬只能代替人計算，不能代替人思考，而人的而價值，就在於他會思考，有想象力和創造力。