雪城第一98k，
美國雪城大學(xué)PhD在讀；高級(jí)摸魚運(yùn)動(dòng)員；有限元科研工作者；PUBG端游愛好者。
由于個(gè)人的原因最近一直在回顧有限元的一些內(nèi)容，再加上之前一直覺得自己寫作能力有很大的問題，也算是自己給自己一個(gè)總結(jié)，就想非常認(rèn)真的寫一篇有限元的科普文章，為廣大結(jié)構(gòu)工作者能在使用有限元程序的時(shí)候能更加好的理解軟件到底在算些什么，順便看看自己的寫作能力還有沒有救，也希望廣泛地接受各位大神給我提點(diǎn)建議或者說一起討論討論。
在我學(xué)習(xí)有限元的過程中，也是參考了各行各業(yè)或者不同網(wǎng)站的資源，首先我覺得何曉明教授上傳到B站上的《有限元基礎(chǔ)編程》是非常不錯(cuò)的一個(gè)系列課程，但是需要你有一點(diǎn)有限元和Matlab以及編程的基礎(chǔ)。
知乎上的Dr. Stein大神的文章也是非常不錯(cuò)的一個(gè)有限元資料，但是他講的非常數(shù)學(xué)也非常專業(yè)，需要多讀讀。
再就是如果能翻墻的話Allan Bower 的網(wǎng)站也是非常好的資料并且給了不同版本的代碼，Brown University畢業(yè)的學(xué)生跟ABAQUS軟件有千絲萬縷的聯(lián)系（當(dāng)然也是聽說），去年上課的時(shí)候教授是Brown University畢業(yè)的大神，他說比較早期的ABAQUS是Brown University的學(xué)生參與了早期的設(shè)計(jì)，因此編程思路上跟Allan Bower教授網(wǎng)站上給的代碼非常相似，包括形成總剛度矩陣的方式等。
如果真的想好好理解有限單元法，最好再去看看數(shù)值分析或者微分方程相關(guān)的基礎(chǔ)知識(shí)。后面如果有什么想起來比較好的資源在補(bǔ)充。
首先，從本質(zhì)上來說，有限元是一種求近似解的思路，最簡單的明了的了解就是對于定義域在[0,π]上，f(x)=sin(x)是一條曲線，如果想求f(x)在定義域上的面積，最簡單的辦法就是積分。
但是如果我不知道具體的函數(shù)表達(dá)式，僅僅知道若干個(gè)點(diǎn)在定義域上所對應(yīng)的函數(shù)值（f(0)=0, …,f(π/2) =1, …, f(π)=0）, 只要知道的點(diǎn)的個(gè)數(shù)足夠多，那么我也可以把不同的點(diǎn)連成直線然后求面積，也可以得到一個(gè)近似的積分值。
那么，有限元到底是怎么應(yīng)用到結(jié)構(gòu)分析上的，在我看來首先要理解一個(gè)應(yīng)力（stress），應(yīng)變（strain），位移（displacement用u來表示）以及位移的導(dǎo)數(shù)
之間的關(guān)系。
用最簡單的例子來解釋，假設(shè)一個(gè)軸向受力單元（1D-bar element），材料的彈性模量是E，截面面積是A，受作用力P之前的位置是x1，x2，受力以后的位置是x1’，x2’，如圖1：
待解的未知量u1=x1’-x1，u2=x2’-x2，由于作用了一個(gè)力在單元上，單元會(huì)有一個(gè)拉長Δu = (x2’-x2)- (x1’-x1)，再帶入材料力學(xué)（我隱約記得是第二章）給出的應(yīng)變的表達(dá)式：
假設(shè)這個(gè)單元非常非常小，也就是對等式的右邊取極限，就可以將應(yīng)變轉(zhuǎn)換為導(dǎo)數(shù)的形式：
在一維線彈性單元，應(yīng)力應(yīng)變之間遵循胡克定律：
再將上式稍微調(diào)整一下位置，再把P換成一般常用的字母小f來表示，Govern equation（實(shí)在不知道應(yīng)該怎么翻譯）可以表達(dá)為：
當(dāng)然這是最最最簡單的一種情況，復(fù)雜一點(diǎn)的單元沒準(zhǔn)會(huì)在后續(xù)中介紹（如果有的話，手動(dòng)狗頭）。
上面這個(gè)方程確實(shí)很好解，只需要把f移到右邊，兩邊同時(shí)積分，就可以得到線彈性情況下的精確解，但是比如涉及到穩(wěn)定性問題（stability），那么Govern equation就會(huì)變成：
我隱約記得之前某個(gè)結(jié)構(gòu)群里討論桿端彎矩對桿件失穩(wěn)（buckling）有沒有影響，如果是按照我的理解那肯定是有影響，但是似乎最后有人說了一句沒有影響，然后問問題的大哥就說好的，那就是不影響。
現(xiàn)在回想一下，似乎生活節(jié)奏變快以后大家會(huì)自動(dòng)過濾跟自己想法不同的聲音，然后選擇性相信自己想相信的內(nèi)容，行吧，扯遠(yuǎn)了，如果對穩(wěn)定性問題有疑慮，我強(qiáng)烈（夾帶私貨）的推薦你們看一下W.F.Chen和E.M.Lui的《Structural stability: theory and implementation》。
回到上面的Govern equation，上面的等式實(shí)質(zhì)上是一個(gè)x, y(x), y''(x)的微分方程，如果用待定系數(shù)法也可以得到精確解，但是解起來很麻煩解一道題差不多就用了好幾頁A4白紙。
因此，有些時(shí)候就需要提高計(jì)算速度，但是又可以得到相對精確的解，而且一定程度上最好能概念化結(jié)構(gòu)化的交給計(jì)算機(jī)去做，有限單元法（Finite element analysis）就出現(xiàn)了。
首先，假設(shè)我們對于微分方程（Partial differential equation）：
的近似解解由一組基函數(shù)(basis function)的疊加組成：
為基函數(shù)，
為待解的未知系數(shù)，除非我們假設(shè)的基函數(shù)跟實(shí)際解的形式相同（或包含關(guān)系），那么我們根據(jù)假設(shè)的解與精確解之間一定會(huì)有誤差。
若將近似解帶入原方程，原微分方程就會(huì)一些誤差，用R來表示：
因此，加權(quán)余量法(Method of Weighted Residual) 要求我們?nèi)舨荒苁刮⒎址匠痰挠嗔縍=0，最起碼要保證在定義域上的余量R盡可能地小。
為此，引入一個(gè)試探函數(shù)（test function）v，對微分方程的余量R與試探函數(shù)v的積在定義域上求積分，再令積分為0。
具體為什么這么做可以參考一下知乎上的大神Dr.Stein的系列文章：
變分法簡介Part 1.（Calculus of Variations） - 知乎 (zhihu.com)
https://zhuanlan.zhihu.com/p/20718489
對于上式，我的理解是經(jīng)過上面的處理，相對減弱了對于微分方程解的限制。
假設(shè)一個(gè)小球從1點(diǎn)滾到2點(diǎn)，我們想要求小球最快下降的時(shí)間，那么時(shí)間既是坐標(biāo)（x,y）的函數(shù)，同時(shí)也是小球運(yùn)動(dòng)路徑的函數(shù)，那對于這個(gè)問題，精確解就是圖2左邊，兩點(diǎn)之間直線最短。
但是，假設(shè)這個(gè)路徑很難直接計(jì)算出來，那么我們稍微放松一下限制條件，小球反正1點(diǎn)2點(diǎn)定了，你要是不愿意沿著直線走那也行吧，但是只要走的別太離譜最后到達(dá)2點(diǎn)的時(shí)間大差不差就行，這樣就可以更輕松的得到解，但是精度上可能大于第一種直接求得的解。
如果放到結(jié)構(gòu)分析里來可以這么理解，一根桿件在P作用下的第一失穩(wěn)模態(tài)（buckling mode）的精確解是y=Asin(πx/L)，但是如果我們假設(shè)他的解的形式是y=a0+a1*x+a2*x^2，也能得到一個(gè)大差不差的解，但是可能會(huì)略大于精確解，這種情況叫做stiffened effect（加強(qiáng)效應(yīng)？我也不知道怎么翻譯了，但是知道是啥事就行了）。
我個(gè)人的理解就是桿件一定會(huì)沿著它本身最不利的情況破壞，我們預(yù)設(shè)的情況如果能滿足邊界條件，那么就會(huì)求得一個(gè)略大的近似解。
OK，啰嗦半天就是希望大家能看懂我想表達(dá)的意思，再回到之前的Govern equation上，表達(dá)式中u(x)（trial function）是精確解，假設(shè)近似解是由一系列的基函數(shù)（Basis function）組成，即：
戴帽子的u表示的就是近似解，那么之前的Govern equation就可以重新表達(dá)為：
上式中對我來講c一般是一個(gè)常數(shù)，因此材料的模量和截面性質(zhì)一般是常數(shù)，就算涉及到材料非線性或者屈服，我們依然可以別的處理方法完成計(jì)算。
如果你好奇為什么這兩項(xiàng)等于零，可以去認(rèn)真的讀一讀我發(fā)在上面鏈接的文章，我們在這里討論的是定端變分的問題，因此在定義域的端點(diǎn)上的變分一定會(huì)等于零。這里如果實(shí)在不好理解可以先放一放，回頭再看可能就看懂了。
這樣處理的好處就是讓Trial function(u)和Test function(v)的derivative order（導(dǎo)數(shù)階？）趨于一致。
在之前提過Test function(v)是由我們定，那好，我直接令他等于近似解
，這就是傳說中的伽遼金法 (Garlekin method) 。那么弱形式的終形態(tài)就變?yōu)椋?
當(dāng)然，除了弱形式，還有強(qiáng)形式(Strong form)。由于弱形式用的比較多，這里就主要介紹弱形式。
到這里，如果帶入我們之前的最早例子里面的參數(shù)c = EA，f = P，那么方程的左邊就是應(yīng)變能，右邊就是外力做功。可以假設(shè)每一個(gè)單元都是一個(gè)小彈簧，結(jié)合回顧初中學(xué)胡克定律時(shí)候的內(nèi)容，是不是感覺一切都聯(lián)系起來了。
再將近似解帶入我們之前假設(shè)的基函數(shù)的形式：
這種寫法看著還是有點(diǎn)迷惑，但是作為僅僅是一篇科普文，那我換一種更加直接了當(dāng)?shù)恼f法：
上式中
是跟
相關(guān)，
是常數(shù)，那么我們的近似解就可以分離為一組常數(shù)（節(jié)點(diǎn)位移）乘以
的形式，這里的
就是傳說中的型函數(shù)（Shape function）。
如果把
提出來，左右兩邊同時(shí)約去，再把求和號(hào)拿出來，最后的形式就可以表達(dá)為：
這樣就把問題轉(zhuǎn)移到求
上，注意一下，
是跟材料截面相關(guān)的常數(shù)，
是我們的未知量。重寫成一種更為熟知的形式就是：
在此再說一下型函數(shù)，型函數(shù)是我們根據(jù)已知推未知的一種方式，假設(shè)已知我們節(jié)點(diǎn)值，但是節(jié)點(diǎn)中間的值我們并不知道，因?yàn)槲覀円呀?jīng)將我們的問題離散化了，也就是把一個(gè)整體劃分成很多小的單元(mesh)。
你可以理解成，已知兩個(gè)人身高，第三個(gè)人的身高介于兩人之間，你就可以根據(jù)Shape function大差不差的猜一下第三個(gè)人的身高。這聽上去很有道理，如果你已知某個(gè)人身高介于姚明易建聯(lián)之間，求出來的誤差可能在幾厘米，但是如果告訴你某個(gè)人身高介于姚明和郭敬明之間，你是不是想弄死我，因?yàn)槲艺f了一句廢話，而且可能誤差也非常大，手動(dòng)狗頭。
故事從另一個(gè)角度告訴我們網(wǎng)格密度的正確性：如果兩個(gè)節(jié)點(diǎn)之間比較接近，那么節(jié)點(diǎn)值之間的差值較小，那么你使用插值的誤差就較小。
關(guān)鍵問題來了，怎么猜？你可以假設(shè)姚明和郭敬明之間的第三個(gè)人遵循一次函數(shù)之間的關(guān)系，也就是說從姚明頭頂連一條線到郭敬明的頭頂，那就可以根據(jù)離郭敬明近還是離姚明近來預(yù)測一下身高。那如果你假設(shè)姚明和郭敬明之間的第三個(gè)人遵循拋物線的關(guān)系，那你就需要謹(jǐn)慎了，特定情況下計(jì)算出來的身高可能比我郭還要矮。
故事從另一個(gè)角度告訴我們兩個(gè)道理：
1、當(dāng)你的計(jì)算結(jié)果看著不是那么對，需要考慮是不是單元選的有問題。
2、所有的有限元計(jì)算結(jié)果需要結(jié)合實(shí)際的物理或者結(jié)構(gòu)常識(shí)來對比。
假設(shè)一個(gè)軸向受力構(gòu)件，受一個(gè)大小為P的外力，我將其離散為3個(gè)單元4個(gè)節(jié)點(diǎn)，圖3所示。我們的未知量就是在每個(gè)節(jié)點(diǎn)處的位移(u1, u2, u3, u4) ，邊界條件 (Boundary Condition) ：假設(shè)最右端是固定端 (u1 = 0) ，求解其余三個(gè)未知量。
在這里我假設(shè)每個(gè)相鄰的點(diǎn)的位移之間相互有聯(lián)系，而且是線性的關(guān)系，線性關(guān)系非常好理解，在單元1中的位移u(x)由其左右兩邊的節(jié)點(diǎn)位移u1和u2線性插值得到，你可以理解成一個(gè)在節(jié)點(diǎn)1處u = u1+0*u2，但是隨著向右邊移動(dòng)，直到節(jié)點(diǎn)2 u = 0*u1+u2，如圖4所示：
我們假設(shè)每個(gè)單元的長度是h，在此h=L/3，在此處i是型函數(shù)的編號(hào)，j是節(jié)點(diǎn)的編號(hào)，對于第一個(gè)單元來說，非零項(xiàng)僅存在于i,j=1,2。
將型函數(shù)帶入單元?jiǎng)偠染仃嚨谋磉_(dá)式，就可以化簡為：
對于單元2，i,j=2,3，帶入單元?jiǎng)偠染仃?
關(guān)于矩陣組合(matrix assemble)，如圖5所示，如果實(shí)在是記不清了，麻煩回去看一下矩陣位移法。
終于快到最后了，對于等式右邊在節(jié)點(diǎn)123處，沒有外力，所以：
此時(shí)將計(jì)算結(jié)果帶入全局剛度矩陣 (Global Stiffness Matrix)和外力矩陣：
如果這個(gè)時(shí)候嘗試去求解矩陣，那么一定會(huì)提示你剛度矩陣不可逆，因?yàn)闆]有施加邊界條件。
大家可以想象一下，一根軸向受拉構(gòu)件，沒有約束，受P大小的作用拉力，那豈不是會(huì)一直做剛體運(yùn)動(dòng)(rigid body movement)。因此，我們需要處理一下邊界條件。
矩陣的本質(zhì)就是由一系列的方程組成，如果想給節(jié)點(diǎn)1賦值，可以令u1的系數(shù)等于1，u2, u3, u4的系數(shù)等于0，然后令結(jié)果等于1，那么最終的矩陣就會(huì)變?yōu)椋?
后續(xù)剩下的內(nèi)容就是非常簡單的線性代數(shù)運(yùn)算了~
