（２）在進入非線性的世界前，先來複習一下線性有限元素法
第二篇想先回來談一下「線性有限元素法」。有許讀者會有明明這是「非線性」的連載，為什麼還跑出來線性的內容？其實這有兩個理由。第一個理由是就算已經實際上跑過模擬，但還是可能對於「有限元素法」的理論背景不清楚。另一個理由是，雖然談的是「非線性模擬」，其實基礎還是建立在「線性模擬」之上；因此在正式說明非線性模擬之前，還是要複習一下縣性的問題，以加深理解的印象。
有限元素法（ＦＥＭ：Finite Element Method）
到底有限元素法是甚麼呢？這裡想先從這個介紹起。現在很多物理行為藉由電腦來求得解答的狀況是越來越常見，這也是大家都很熟悉、被稱為「模擬（Simulation）」的方法。而使用電腦來求解物理現象的手法，現在這個世界中也有非常多種。其中在結構行為的模擬上最廣泛被使用的求解方法，正是這個有限元素法。
那麼這個有限元素法是甚麼呢？簡單來說，就是「將想要分析模擬的對象分割成有限個數的元素，利用微分方程式找出『近似解』的方法」。說的更具體一點，就是利用「柯西（Cauchy）第一運動定律」、或者是以近似的方式解出平衡方程式，就是有限元素法。想要更加明白柯西第一運動定律的人，可以去找有關連續體力學的書籍來研讀。
方程式１
順便補充一下，不同的書籍會將方程式１寫成不同的表示方式。
那麼被分析模擬的對象，是怎樣透過有限元素法來分析模擬其行為的呢？接下來就想來請大家確認一下其相關流程了。基本的流程大略會是分成以下的各個階段：
１．離散化（建立網格＝Mesh）
２．建立元素的剛性矩陣（Stiffness Matrix）
３．建立全體的剛性矩陣
４．套入邊界條件（Boundary Condition）
５．矩陣運算求解（解出聯立一次方程式）
６．計算應力等結果
這個流程，在各位讀者有使用的有限元素法模擬軟體中，不管是哪一種軟體，基本上都是一樣的。也是這個緣故，接下來就要一一來看看這些程序的內容。
１．離散化（建立網格）
為了使用有限元素法進行模擬，不管怎樣都需要執行「切割連續體」的步驟。順便補充一下，在連續體分析學中，所謂的連續體是指「位移連續變化」的意思。而將連續體分割成各個元素（Element），就是一種「離散化」。
例如像圖１這樣，當出現了一個任意形狀的分析模擬對象時，就是要將其變成圖１右側那樣的分割元素的集合體。
圖１ 將被分析模擬對象離散化
換句話說，就是將原本的連續體領域稱為Ｖ，而經過元素分割之後、個個元素所組合起來的領域稱為Ｖｅ；而離散化之後的模型全體領域，就等於元素領域的集合體。
方程式２
換句話說，這種離散化，就是一種近似表現原始物體模型的方法。
順帶一提，在進行離散化之際所建立的元素，其形狀／種類並不會只有一種。如果是平面元素的話，通常會是三角形或四邊形；而立體元素的話，通常會是四面體或六面體。關於元素本身的詳細內容，這裡就割愛了；唯一要說的是在立體元素切割上，六面體的一階元素在精度與計算資源花費的關係上會比較均衡，而通常為模擬分析者所喜愛，但是六面體的自動切割網格會相當困難，一般都是選擇二階的四面體（而且幾乎都是唯一選擇）。
圖２ 離散化時建立的元素種類
那麼在元素表面上出現的變形，基本上就是用構成元素的整體變形來表現，那麼在元素內部的連續性（變形連續性）又該怎麼樣來表示呢？這裡就要使用被稱為「形狀函數」的函數來近似表示。雖說，單看底下這個方程式，可能會搞不懂其意義；總之就是請理解為怎麼使用形狀函數的一個範例。
圖３ 四節點之四邊形元素
方程式３ （形狀函數的一部分）
順便補充一下，一般在後處理中表示應力分布等等結果的輪廓圖（等高線圖），也都是直接沿用這個形狀參數。此外，上述這些元素的處理，都會將正交（Normal）座標系統轉換成整體座標系統的座標來進行。
２．建立元素的剛性矩陣（Stiffness Matrix）
所謂的切割網格，就是在重新定義被分析模擬對象的形狀。不過這也僅僅是在定義被分析模擬對象的形狀而已；接下來則是必須賦予這個定義好的形狀怎樣的材料性質（硬或軟等等性質）。這就和元素剛性有關係，接下來就稍微解釋這個部分。
本篇一開頭說明過使用柯西第一運動定律來求得近似解就是有限元素法。說的再具體一點，根據這個運動定律，給予應力／應變的關係式、形狀、位移的邊界條件、應力的邊界條件，就能知道物體的位移／變形。而且通常我們要處理的實際物件形狀都相當複雜，要直接手算是相當困難，因此都必須使用電腦、執行有限元素法的方式來計算。
那麼我們就來看看要賦予這個柯西第一運動定律的各個條件吧。
應力應變方程式
其中應力應變的關係式，正如方程式４所示：
方程式４
或許會覺得這個方程式好像在哪裡看過對吧？不過這個方程式，如果用張量（Tensor）的形式來表示的話，其內部看起來就很複雜了。以立體的線性彈性體而言，其內部會長的像底下那樣。這個用Ｄ來表示的應力應變關係矩陣，通常也會被稱為Ｄ矩陣。
方程式５
假設這個分析模擬對象是由等方向性線性彈性材料所構成的話，其物性就能以楊氏係數（Ｅ）與蒲松比（ν）來表示，也就有可能把方程式５寫成方程式６那樣的形式。從方程式６當中，就能理解為什麼在模擬軟體的材料物性中需要輸入楊氏係數與蒲松比了。
方程式６
從以上各個方程式來看，或許可以讓大家對應力應變的關係有些概念。再來就是要看應變與位移的關係了。
應變變形關係式
應變與變形（位移）的關係，則如以下所示：
方程式７
順便補充一下，如果是線性彈性分析的問題，應變可以用以下微分方程式表示：
方程式８
此外，變形可以利用節點變形與形狀函數來表示成以下的方程式：
方程式９
以這些方程式為基礎，舉例來說，就能將四節點四邊形元素的變形與應變表示成以下的形式：
方程式１０
形狀函數跟元素的形狀與階數有關，所以方程式７中的Ｂ矩陣的內部也都跟元素的形狀與階數相關。
接下來就要將有限元素的公式建立起來了。
使用前述的柯西第一運動定律，可以將物件全體的受力關係表示成以下的方程式：
方程式１１
然後套用虛功原理，就可以推導出方程式１２：
方程式１２
然後利用分部積分法與高斯發散原理，就可以將方程式改寫成以下的形式：
方程式１３
再假設這是靜態問題，將速度ai等於０。省略中間詳細的算式簡化推導、再加以離散化，就能得到方程式１４：
方程式１４
再將前述的應力應變關係、變形應變關係與形狀函數等等方程式帶入方程式１４中加以整理，並且將虛位移套入這個式子也成立，就會再得到方程式１５：
方程式１５
方程式的左邊有剛性矩陣與變形向量，右邊則是外力向量，就可以看出這個方程式實際是被整理成[K]{u}＝{f}的樣子。這個方程式本身是聯立一次方程式，就可以根據這個聯立一次方程式來求出某個外力負載下造成的變形結果。這些式子的詳細積分過程也一併割愛，總之就是利用高斯積分等數值階分方法來計算。元素的應力或應變，為了高斯數值積分的計算方便，通常是計算被稱為「高斯積分點」的取樣位置上的結果。
至於方程式１５中左邊的部分：
方程式１６
就是各元素的剛性矩陣，如果假設其為平面的四節點四邊形元素（一階元素），這個矩陣舊會長得像方程式１７的樣子
方程式１７
此時這個矩陣的大小，就由節點數ｘ節點自由度來決定。對於一般的連續體元素而言，節點的自由度只有平移自由度而已；以平面元素而言，就是每個節點只有Ｘ與Ｙ的兩個自由度，可以得知剛性矩陣的大小就是８ｘ８。
至此就建立好了元素的剛性矩陣，再來就要建構全體的剛性矩陣了。
３．建立全體的剛性矩陣
全體的剛性矩陣，簡單來說就是「將各元素的剛性疊加組合起來」，就能建立。
舉例來說，如同前面提及的平面四邊形一階元素，整體模型只有兩個元素的話，就有以下的兩個剛性矩陣。
圖４ 建立全體剛性矩陣
方程式１８
如方程式１８所示，將兩個矩陣合成起來，就能建立全體的剛性矩陣，而將變形向量、負載向量表現出來的話，會像是底下這個方程式的樣子：
方程式１９
矩陣當中０的位置，就是兩個要素之間彼此不相關的部分，而相加起來的部分，則是兩個要素之間彼此相關的部分。
４．套入邊界條件
到這裡幾乎就可以準備開始求解了，但還稍微少了一點資訊。這個缺少的資訊，就是所謂的「邊界條件」。而邊界條件在很多模擬軟體的使用者介面上，會以「限制（拘束）條件」或「負載條件」來表示。
如圖５所示，就是將左邊的兩個節點限制拘束住，而右邊的兩個節點施以任意的負載ｆ，此時就會建立出方程式２０。此時在沒有負載（未施加外力）的節點自由度上，帶入０的數值。
圖５ 套入邊界條件
方程式２０
之後，再將變形（位移）為零處的行列消除，就會讓未知數的數目與方程式的數目相符起來。只不過，一般常用的商用模擬軟體，並不會使用如上述的教科書等級理論求解方式，往往會用其他的形式來求解就是了（詳細的原理此處也就割愛）。
５．矩陣運算求解（解出聯立一次方程式）
建立好全體剛性矩陣與負載向量後，就能求解聯立一次方程式，並且計算出未知數的變形u。當然，實際上要求解的模擬模型，其自由度數量非常非常多，不可能用手來計算。此時就得依靠電腦程式來求解了。像這樣用來求解聯立一次方程式的程式，一般會稱為「矩陣求解器」。
矩陣求解器其實有各式各樣的理論，很多的商用模擬軟體則是用好幾種求解器組合起來計算。而矩陣求解器基本上大致可以分成被稱為「直接法」與「反覆疊代法」兩種類型。
直接法的求解器，是直接消去聯立方程式的形式，只要矩陣不是非正定（non-positive-definite），就一定能解出結果（但條件不好的話，是有可能造成求解精度惡化），而且能獲得穩定的解，所以很多商用模擬軟體預設都是使用直接法來求解。相反地，直接法求解器的缺點，就是很吃記憶體、因而計算時間也會變長，遇到大規模的問題時，就需要小心注意。
至於反覆疊代法，則是會先推測出一個初期值，然後反覆疊代計算直到解答收斂的程度合理為止，是一種求出近似解的方法。反覆疊代法與直接法相比，常常會有需要的記憶體與計算時間都比較少的好處，但是遇到矩陣的狀態不佳，也會讓反覆疊代的次數變多而讓計算時間變長。換句話說，反覆疊代法的求解效率要看分析模擬的問題內容而定。如果無法很好地掌握求解器的特徵的話，使用直接法會更加確實一點。
６．計算應力等結果
矩陣運算求解可以得到的結果，就是元素的節點變形（位移）。那麼要如何從節點的變形中，獲得大家常常想看到的應力、應變的結果呢？那就是根據前面介紹過的應力應變關係式（{σ}＝[D]{ε}）與變形應變關係式（{ε}＝[B]{u}）了。
以底下這個元素為例：
圖６ 計算應力等結果
首先是將求解出來的節點變形與Ｂ矩陣，求出元素中積分點的應變。元素積分點的應變值，就是從高斯數值積分的結果求出。獲得應變值後，就可以再利用Ｄ矩陣得到積分點的應力值。
以上雖然有點冗長，但這就是線性有限元素法程式的計算流程。
但這也只不過是說明大概的流程而已，對於專門寫作模擬程式的人而言，這是再也理所當然不過的內容。對於想要更詳細知道其內容的讀者，可以參考有限元素法的專門書籍。
＝＝＝
沒有留言:
張貼留言