由于個人的原因最近一直在回顧有限元的一些內容，再加上之前一直覺得自己寫作能力有很大的問題，也算是自己給自己一個總結，就想非常認真的寫一篇有限元的科普文章，為廣大結構工作者能在使用有限元程序的時候能更加好的理解軟件到底在算些什么，順便看看自己的寫作能力還有沒有救，也希望廣泛地接受各位大神給我提點建議或者說一起討論討論。

在我學習有限元的過程中，也是參考了各行各業或者不同網站的資源，首先我覺得何曉明教授上傳到B站上的《有限元基礎編程》是非常不錯的一個系列課程，但是需要你有一點有限元和Matlab以及編程的基礎。

知乎上的Dr. Stein大神的文章也是非常不錯的一個有限元資料，但是他講的非常數學也非常專業，需要多讀讀。

再就是如果能翻墻的話Allan Bower 的網站也是非常好的資料并且給了不同版本的代碼，Brown University畢業的學生跟ABAQUS軟件有千絲萬縷的聯系（當然也是聽說），去年上課的時候教授是Brown University畢業的大神，他說比較早期的ABAQUS是Brown University的學生參與了早期的設計，因此編程思路上跟Allan Bower教授網站上給的代碼非常相似，包括形成總剛度矩陣的方式等。

如果真的想好好理解有限單元法，最好再去看看數值分析或者微分方程相關的基礎知識。后面如果有什么想起來比較好的資源在補充。

首先，從本質上來說，有限元是一種求近似解的思路，最簡單的明了的了解就是對于定義域在[0,π]上，f(x)=sin(x)是一條曲線，如果想求f(x)在定義域上的面積，最簡單的辦法就是積分。

但是如果我不知道具體的函數表達式，僅僅知道若干個點在定義域上所對應的函數值（f(0)=0, …,f(π/2) =1, …, f(π)=0）, 只要知道的點的個數足夠多，那么我也可以把不同的點連成直線然后求面積，也可以得到一個近似的積分值。

那么，有限元到底是怎么應用到結構分析上的，在我看來首先要理解一個應力（stress），應變（strain），位移（displacement用u來表示）以及位移的導數

之間的關系。

用最簡單的例子來解釋，假設一個軸向受力單元（1D-bar element），材料的彈性模量是E，截面面積是A，受作用力P之前的位置是x1，x2，受力以后的位置是x1’，x2’，如圖1：

△??圖1

待解的未知量u1=x1’-x1，u2=x2’-x2，由于作用了一個力在單元上，單元會有一個拉長Δu = (x2’-x2)- (x1’-x1)，再帶入材料力學（我隱約記得是第二章）給出的應變的表達式：

假設這個單元非常非常小，也就是對等式的右邊取極限，就可以將應變轉換為導數的形式：

在一維線彈性單元，應力應變之間遵循胡克定律：

再將上式稍微調整一下位置，再把P換成一般常用的字母小f來表示，Govern equation（實在不知道應該怎么翻譯）可以表達為：

當然這是最最最簡單的一種情況，復雜一點的單元沒準會在后續中介紹（如果有的話，手動狗頭）。

上面這個方程確實很好解，只需要把f移到右邊，兩邊同時積分，就可以得到線彈性情況下的精確解，但是比如涉及到穩定性問題（stability），那么Govern equation就會變成：

我隱約記得之前某個結構群里討論桿端彎矩對桿件失穩（buckling）有沒有影響，如果是按照我的理解那肯定是有影響，但是似乎最后有人說了一句沒有影響，然后問問題的大哥就說好的，那就是不影響。有限單元法王勖成

現在回想一下，似乎生活節奏變快以后大家會自動過濾跟自己想法不同的聲音，然后選擇性相信自己想相信的內容，行吧，扯遠了，如果對穩定性問題有疑慮，我強烈（夾帶私貨）的推薦你們看一下W.F.Chen和E.M.Lui的《Structural stability: theory and implementation》。

回到上面的Govern equation，上面的等式實質上是一個x, y(x), y''(x)的微分方程，如果用待定系數法也可以得到精確解，但是解起來很麻煩解一道題差不多就用了好幾頁A4白紙。

因此，有些時候就需要提高計算速度，但是又可以得到相對精確的解，而且一定程度上最好能概念化結構化的交給計算機去做，有限單元法（Finite element analysis）就出現了。

首先，假設我們對于微分方程（Partial differential equation）：

且有邊界條件：

的近似解解由一組基函數(basis function)的疊加組成：

為基函數，為待解的未知系數，除非我們假設的基函數跟實際解的形式相同（或包含關系），那么我們根據假設的解與精確解之間一定會有誤差。

若將近似解帶入原方程，原微分方程就會一些誤差，用R來表示：

因此，加權余量法(Method of Weighted Residual) 要求我們若不能使微分方程的余量R=0，最起碼要保證在定義域上的余量R盡可能地小。

為此，引入一個試探函數（test function）v，對微分方程的余量R與試探函數v的積在定義域上求積分，再令積分為0。

具體為什么這么做可以參考一下知乎上的大神Dr.Stein的系列文章：

變分法簡介Part 1.（Calculus of Variations） - 知乎 (zhihu.com)

https://zhuanlan.zhihu.com/p/20718489

對于上式，我的理解是經過上面的處理，相對減弱了對于微分方程解的限制。

引用經典變分法的例子來說就是：

△??圖2

假設一個小球從1點滾到2點，我們想要求小球最快下降的時間，那么時間既是坐標（x,y）的函數，同時也是小球運動路徑的函數，那對于這個問題，精確解就是圖2左邊，兩點之間直線最短。

但是，假設這個路徑很難直接計算出來，那么我們稍微放松一下限制條件，小球反正1點2點定了，你要是不愿意沿著直線走那也行吧，但是只要走的別太離譜最后到達2點的時間大差不差就行，這樣就可以更輕松的得到解，但是精度上可能大于第一種直接求得的解。

如果放到結構分析里來可以這么理解，一根桿件在P作用下的第一失穩模態（buckling mode）的精確解是y=Asin(πx/L)，但是如果我們假設他的解的形式是y=a0+a1*x+a2*x^2，也能得到一個大差不差的解，但是可能會略大于精確解，這種情況叫做stiffened effect（加強效應？我也不知道怎么翻譯了，但是知道是啥事就行了）。

我個人的理解就是桿件一定會沿著它本身最不利的情況破壞，我們預設的情況如果能滿足邊界條件，那么就會求得一個略大的近似解。

OK，啰嗦半天就是希望大家能看懂我想表達的意思，再回到之前的Govern equation上，表達式中u(x)（trial function）是精確解，假設近似解是由一系列的基函數（Basis function）組成，即：

戴帽子的u表示的就是近似解，那么之前的Govern equation就可以重新表達為：

上式中對我來講c一般是一個常數，因此材料的模量和截面性質一般是常數，就算涉及到材料非線性或者屈服，我們依然可以別的處理方法完成計算。

因此，上式的左邊就可以重新寫為：

再對其進行分部積分處理：

同時帶入邊界條件：

如果你好奇為什么這兩項等于零，可以去認真的讀一讀我發在上面鏈接的文章，我們在這里討論的是定端變分的問題，因此在定義域的端點上的變分一定會等于零。這里如果實在不好理解可以先放一放，回頭再看可能就看懂了。

上面處理過最終的形式就化簡為：

這樣處理的好處就是讓Trial function(u)和Test function(v)的derivative order（導數階？）趨于一致。

在之前提過Test function(v)是由我們定，那好，我直接令他等于近似解

，這就是傳說中的伽遼金法 (Garlekin method) 。那么弱形式的終形態就變為：

當然，除了弱形式，還有強形式(Strong form)。由于弱形式用的比較多，這里就主要介紹弱形式。

到這里，如果帶入我們之前的最早例子里面的參數c = EA，f = P，那么方程的左邊就是應變能，右邊就是外力做功。可以假設每一個單元都是一個小彈簧，結合回顧初中學胡克定律時候的內容，是不是感覺一切都聯系起來了。

再將近似解帶入我們之前假設的基函數的形式：

這種寫法看著還是有點迷惑，但是作為僅僅是一篇科普文，那我換一種更加直接了當的說法：

上式中

是跟相關，是常數，那么我們的近似解就可以分離為一組常數（節點位移）乘以的形式，這里的就是傳說中的型函數（Shape function）。

如果把

提出來，左右兩邊同時約去，再把求和號拿出來，最后的形式就可以表達為：

這樣就把問題轉移到求

上，注意一下，是跟材料截面相關的常數，是我們的未知量。重寫成一種更為熟知的形式就是：

這樣是不是親切多了（手動狗頭）。

在此再說一下型函數，型函數是我們根據已知推未知的一種方式，假設已知我們節點值，但是節點中間的值我們并不知道，因為我們已經將我們的問題離散化了，也就是把一個整體劃分成很多小的單元(mesh)。

你可以理解成，已知兩個人身高，第三個人的身高介于兩人之間，你就可以根據Shape function大差不差的猜一下第三個人的身高。這聽上去很有道理，如果你已知某個人身高介于姚明易建聯之間，求出來的誤差可能在幾厘米，但是如果告訴你某個人身高介于姚明和郭敬明之間，你是不是想弄死我，因為我說了一句廢話，而且可能誤差也非常大，手動狗頭。

故事從另一個角度告訴我們網格密度的正確性：如果兩個節點之間比較接近，那么節點值之間的差值較小，那么你使用插值的誤差就較小。

關鍵問題來了，怎么猜？你可以假設姚明和郭敬明之間的第三個人遵循一次函數之間的關系，也就是說從姚明頭頂連一條線到郭敬明的頭頂，那就可以根據離郭敬明近還是離姚明近來預測一下身高。那如果你假設姚明和郭敬明之間的第三個人遵循拋物線的關系，那你就需要謹慎了，特定情況下計算出來的身高可能比我郭還要矮。

故事從另一個角度告訴我們兩個道理：

1、當你的計算結果看著不是那么對，需要考慮是不是單元選的有問題。

2、所有的有限元計算結果需要結合實際的物理或者結構常識來對比。

讓我們回到一維軸向受力單元問題上：

△??圖3

假設一個軸向受力構件，受一個大小為P的外力，我將其離散為3個單元4個節點，圖3所示。我們的未知量就是在每個節點處的位移(u1, u2, u3, u4) ，邊界條件 (Boundary Condition) ：假設最右端是固定端 (u1 = 0) ，求解其余三個未知量。

在這里我假設每個相鄰的點的位移之間相互有聯系，而且是線性的關系，線性關系非常好理解，在單元1中的位移u(x)由其左右兩邊的節點位移u1和u2線性插值得到，你可以理解成一個在節點1處u = u1+0*u2，但是隨著向右邊移動，直到節點2 u = 0*u1+u2，如圖4所示：

△??圖4

我們假設每個單元的長度是h，在此h=L/3，在此處i是型函數的編號，j是節點的編號，對于第一個單元來說，非零項僅存在于i,j=1,2。

將型函數帶入單元剛度矩陣的表達式，就可以化簡為：

對于單元2，i,j=2,3，帶入單元剛度矩陣

同樣，單元3，i,j =2,3

△??圖5

關于矩陣組合(matrix assemble)，如圖5所示，如果實在是記不清了，麻煩回去看一下矩陣位移法。

終于快到最后了，對于等式右邊在節點123處，沒有外力，所以：

此時將計算結果帶入全局剛度矩陣 (Global Stiffness Matrix)和外力矩陣：

如果這個時候嘗試去求解矩陣，那么一定會提示你剛度矩陣不可逆，因為沒有施加邊界條件。

大家可以想象一下，一根軸向受拉構件，沒有約束，受P大小的作用拉力，那豈不是會一直做剛體運動(rigid body movement)。因此，我們需要處理一下邊界條件。

矩陣的本質就是由一系列的方程組成，如果想給節點1賦值，可以令u1的系數等于1，u2, u3, u4的系數等于0，然后令結果等于1，那么最終的矩陣就會變為：

后續剩下的內容就是非常簡單的線性代數運算了~

總結

以上是生活随笔為你收集整理的到底什么是有限单元法？的全部內容，希望文章能夠幫你解決所遇到的問題。

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