溫故而知新。
--《論語》
這里只討論基於位移(注1)的有限元方法,它最終建立的是關於位移為未知量的方程組。也就是說,在推導過程中涉及到的其他未知物理量都要轉化為位移的表達。
推導過程中用到的關鍵物理定律,虛功原理或最小勢能原理。由於這兩者在某些情形下是等價的,這里只以虛功原理為例進行說明。所謂的虛功原理,是說對於靜態平衡的系統,力(注2)在虛位移上所做的功等於應力在虛應變上所作的功[1]。
本段來具體解釋虛功原理中涉及到的概念和方程。先來看功。在標量的情形下,功等於力和其作用距離的乘積,而在向量的情形下,功等於力與物體位移(這兩者都是向量)的內積,標量的結果可以看作它的一種特殊情形。值得一提的是,由於體力,面力以及外力和物體的具體形狀都給定,所以虛功原理中的力可以看作是已知的。再來看應力和應變。直觀地理解,應力可以看作物體的內力(注3)(鮑老師語),應變是在應力作用下引起的變形的一種描述,它包括物體內一點鄰域內變形前后長度和角度的變化量[2]。在彈性力學中給出了應力和應變之間的關系,稱為本構關系,它是有限元推導過程中一個必不可少的關系(注4);由於討論的是基於位移的有限元,如果能夠給出應變與位移之間的函數關系(稱為幾何關系,推導過程中另一個關系),加上剛提到的本構關系,那么應力在虛應變上所作的功就可以表述為位移的函數。最后來看“虛”這個概念,虛位移指的是符合約束條件的任意的無窮小位移。
根據虛功原理列出相應的方程組后,經過移項等處理后,得到關於位移的函數乘以虛位移等於零的結果。注意到虛位移的任意性,得到關於位移的函數為零(這一點可以給出證明[3],它也是數學分析中一個經典的證明,有興趣的同學可以自己證明下)。這就是有限元要建立的最終結果。
細心的童鞋可能會發現,上面的敘述只是一個理論上的框架,其中某個未知變量關於位移的函數只是一筆帶過。那么,具體到編程實現上,該如何給出相應的函數呢?首先,給出整個物體上關於位移(它是一個連續量)的函數在絕大部分情況下是不可行的。退而居其次,在離散的空間中進行定義,如果離散的程度很小,就可以近似連續情形下的結果。這是數值求解實際問題的一個通用的步驟,稱之為連續問題離散化。具體到有限元中,就是將一個連續區域划分成彼此不相交的區域(二維情形下是三角形或者四邊形,三維情形下是四面體或者六面體),稱之為單元或者網格,單元之間通過節點相連,有限元最終建立的就是關於單元節點處位移的方程組。而單元中除了節點外其他的點處的位移,可以基於節點處的位移值進行插值得到。其次,即便將問題域離散化,那么采用什么樣的函數(稱之為形函數)來擬合單元內某點處的位移呢?最容易想到的,自然是多項式函數[4],主要是它容易理解,計算簡單,逼近的精度也能達到要求的優點。當然,函數系的選取要結合具體的問題。如果邊界條件是周期邊界條件,經常使用三角函數來進行逼近。現在有了單元內任一點的位移表達式,那么相應的應變和應力也就可以表達出來,自然地,應力在虛位移上所作的功也可以表示出來。注意的是,這里要使用積分的表達式,而且最終的結果中是關於位移的二次項表達。
再來看等式的另一端,目的是將之表達為節點處位移的函數,再根據虛位移的任意性,就可以將問題轉化為關於位移的方程組。僅以體力為例進行說明,給定單元內任一點處的位移表達式,根據等功原理,即體力在單元上所作的功將其等價為集中力力在單元節點處所作的功,就將體力的功轉化到單元節點處的表達。
主體框架搭建完畢,簡言之,最終建立了一個關於單元節點處位移的方程組。剩下的一些關鍵問題,比如從單元剛度矩陣到總體剛度矩陣的構建(涉及到坐標系的旋轉變換,計算機圖形學的基礎內容),線性方程組的求解,積分的數值運算,力的疊加等等,其中第一個問題是編程時必須考慮的問題,總剛矩陣最后的結果使用稀疏矩陣來存儲,后兩個問題是數值計算研究的主要內容。至於說力的疊加,對每個單元使用虛功原理,求出單元上每個節點相應的內力,添加到節點力的矩陣中。也就是說,具體的實現就成了一個數學問題。 上述問題如果不加任何限制的話,在數學上是有無數多組解的(如果是線性方程組,對應的系數矩陣奇異)。為了得到唯一的解,需要增加約束條件,通常在有限元法中對邊界施加,稱為邊界條件。以傳熱問題為例,邊界條件又分為邊界節點的溫度固定(第一邊界條件),邊界節點處溫度的梯度固定(第二邊界條件),或者兩者的組合(第三邊界條件)。不較真地可以這樣認為,上述問題在數學上都是可解的。 盡管從理論上問題是可解的,但是從程序實現上仍然有很多地方需要注意。由於方程組的系數矩陣是一個大規模的對稱稀疏矩陣,在存儲時采用三個向量來存儲,一個用來存儲某一行的非零元素的個數,一個用來存儲該行非零元素所在的列,一個用來存儲對應的某行非零元素的具體數值[5];在將邊界條件考慮進去以后,對應的方程組也要相應地進行修改,對於第一邊界條件,直觀的想法自然是直接對其賦值,但是在程序的具體實現過程中,使用矩陣乘大數法;至於大規模稀疏線性方程組的求解,以前常常以迭代法為主,隨着並行和多核的發展,目前直接法也得到了越來越多的關注;而針對大規模的非線性問題,則使用Newton-Rapshon方法將其轉化為一系列線性方程組的求解。 感謝耐心地看到這里,最近一個關於有限元的邊界條件問題仍舊困擾着我,一般的教材講到有限元邊界條件時,都是在問題求解區域的邊界上施加限制, 那么能否在一部分邊界上施加條件,一部分內部節點上施加條件呢?從程序實現的角度來分析是可行的,另一個棘手的問題是這樣施加的條件,有唯一解嗎?更進一步,施加什么樣的條件,數學上能夠保證問題有唯一解?希望大家不吝賜教。 總結:有限元法用到的知識點,虛功原理(剛度矩陣的構造和整體方程的構造),轉換矩陣(局部坐標系和整體坐標系的轉換),疊加原理(力的疊加,單元剛度矩陣構造整體剛度矩陣),邊界條件(實際問題及唯一解的需求)。
附:線性問題中三角形單元涉及到各個量的轉換關系圖[6]:

注1:假設通過求解線性方程組得到位移的結果,那么根據幾何關系,就可以得到應變的分布,再根據本構關系,能夠得到應力的分布。
注2:這里的力包括體力和面力。體力,指的是分布在物體體積上的力,如重力等,它在有限元中對應於體積載荷,推導過程中要將其轉化為對應的單元節點力,文中有述;面力,即作用在物體表面上的力,在有限元中對應於邊界分布的情形,和體力一樣,推導時也要將其轉化為對應的單元節點力。
注3:關於應力和內力的轉換,依舊使用虛功原理來實現,只不過這個時候將其應用於某一個單元內,而不是整個求解區域。這樣就可以求解出該單元每個節點的內力向量,它表達為一個矩陣(線性變形下為常數)乘以應力。節點內力和位移之間的關系通過單元剛度矩陣來表達[6]。
注4:在一維情形下,彈簧的受力(應力)拉伸或者壓縮(應變,變化前后的長度差和原始長度的比值)滿足本構關系,即f=kx。在二維情況下,應力和應變之間的轉換,Omega=M*Delta,只不過這里應力和應變是向量形式,對應的M是矩陣形式。如果M是常數矩陣的話,那么相應的力學問題成為線性問題,如果M中的元素是關於位移的函數時,則該問題是非線性問題。
參考文獻:
[1] http://zh.wikipedia.org/wiki/%E8%99%9B%E5%8A%9F#.E8.99.9B.E5.8A.9F.E5.8E.9F.E7.90.86
[2] 《物理學與偏微分方程第二版(上冊)》,第五章:彈性力學,李大潛等著,北京:高等教育出版社,2006年
[3] 《微分方程數值解法(第四版)》,李榮華等著,北京:高等教育出版社,2009年
[4] http://zh.wikipedia.org/wiki/%E9%80%BC%E8%BF%91%E7%90%86%E8%AE%BA
[5] Helfenstein R, Koko J. Parallel preconditioned conjugate gradient algorithm on GPU[J]. Fuel & Energy Abstracts, 2012, 236(15):3584-3590. 4.1 Matrix storage
[6] http://wenku.baidu.com/link?url=wNatQMjgOjH-viTSBwPu-xNT1a46mPl9q0bf6logQQL0iXVl6urOlzEoQ0PwlwnB4LujELtOcR1PixcZY__TB2lMBFwVDVNagX7kKfHvcgu (同濟大學孫遠韜)第7講_三角形單元 P24 以及 Step3 單元剛度矩陣