有限元方法及ANSYS平台的力學模擬（一）

今年寒假，寫寫關於有限元理論的數學力學基礎以及基於ANSYS平台的力學模擬。（本文內容冗長，多半是學習的筆記和心得，後面如果自己在這個領域繼續研究，就再寫一個關於自己觀點小短篇）話不多說，開始了~

材料力學和彈性理論的知識在專欄的其他文章中寫到了，這裡就不贅述啦，直接進入主題，如果力學還有問題，那抓緊時間補補吧~

先引入一個小問題，要求計算它的位移場：

一般的材料力學理論不要求計算位移場，一般只用計算某點位移就可以了，這個位移場的確定，首先應該取微元體，建立微元滿足的微分方程，容易得到上圖右側的控制方程。1D問題，2階常微分方程，應該有兩個邊界條件，定解條件。這個問題的解析解十分容易獲得，只需要移項，積分兩次，就可以得到原函數，再根據B.C.定待定係數就可以了。但是大多B.C.複雜的工程問題，沒有解析解或者很難求解，需要我們採用數值方法，這裡先引入兩個數值方法，差分法與加權殘差法，和解析解對比可知，精確度都比較高。

方法一：差分法

第一步，用二階差分代替二階微分，將差值帶成L/5得到上圖中的四個等式。原邊界條件化為U0=0與U4=U5。將方程轉化為矩陣形式並求解，可以得到：

conclusion：差分的解法將微分方程轉化為線性方程組，便於求解。

如果提高差分的個數，將提高精確度：

方法二：加權殘差法

首先找到一個試函數，試函數必須首先滿足邊界條件，然後試函數帶入控制方程，（一般情況下，不能滿足控制方程，如果恰好滿足控制方程，則恰好找到解析解），得到差值，由於函數中存在待定常數，確定待定常數的方法就是使得差值最小（為0）時，所確定的待定常數。

現在用試函數法來求解這個問題：

上面的算例是一個試函數時的計算結果，接下來利用兩個試函數進行計算：

計算結果中明顯的發現，兩個函數的結算結果非常接近解析解。

接下來對各個方法進行比較：

在逼近的過程中，可以採用全域逼近，或者子域逼近，全域逼近有良好的解析性質，但是難度大，甚至無法實現，子域逼近操作簡單，逼近程度高，但是存在不可導點。

下面對2D，3D以及簡單，複雜形狀進行逼近：

得出的結論是，3D的複雜的結構，很難進行全域逼近，而可以實現子域逼近。有限單元法的核心思想，也就是將整體複雜的幾何構型，離散化為簡單的單元，事先解決單元的問題，就等於解決了各種類型複雜的工程結構的問題。

根據幾何構型與力學條件，創造出不同的單元，分別有如下單元可供選擇：

從1D的單元開始介紹，首先引入彈簧單元：

對於彈簧單元，取彈簧的兩端為兩個節點，再構造兩個列陣，分別稱之為節點位移列陣和節點力列陣。根據力學的平衡方程和彈簧的本構方程，聯立，得到彈簧的剛度方程。聯想到結構力學中兩種超靜定問題的求解思路，力法，即柔度法，根據變形協調關係求解，位移法，即剛度法，根據力的平衡方程求解，這裡的剛度矩陣類似於位移法中的剛度矩陣，位移法方程是平衡方程，此處也如此。

接下來考慮兩個彈簧串聯的情況，當兩個彈簧串聯時，有三個節點，在連接點處施加外力F2，根據兩個彈簧的平衡方程和本構方程分別分析兩個彈簧。

將兩個彈簧的方程分別用矩陣表示，將矩陣擴充，多餘位置設置為0，再將擴充後的矩陣相加，得到系統總體的剛度方程，根據系統總體的平衡方程和邊界條件，求解位移於反力。

有了彈簧單元的剛度方程，接下來討論桿單元的剛度方程，由於桿在受力後的彈性階段，擁有和彈簧類似的「正比」本構關係，因此，在求解桿的剛度方程時，可以對彈簧進行剛度替換即可。

將k替換成EA/L，便得到了單個桿的剛度方程，類似的單元剛度矩陣，單元節點力，單元節點位移，就依次導出了。接下來對剛才的串聯彈簧，通過剛度代換的方式，嘗試求解串聯桿的剛度方程，通過這些，就大致了解了桿單元的單元剛度矩陣。

那麼現在還有一些問題，比如每一個一維的桿件的剛度方程已經建立了，但是，一個複雜的結構並不是一維的，例如結構力學問題，它由許許多多的一維桿構成，那麼每個桿用一個坐標系顯然不合適，我們要通過坐標變換得到統一的坐標公式。

本身的桿是一個1D的建模，現在取新的坐標系，稱之為整體坐標系。現在就產生了兩個列陣，一個稱之為局部坐標系中的節點位移，一個稱之為整體坐標系中的節點位移。現在尋求這兩者之間的關係，如上圖。

局部坐標的節點位移列陣與整體坐標的節點位移列陣之間的關係就是坐標變換矩陣。用Te表示。這裡有順序問題，順序是：整體坐標系中的節點位移列陣左乘坐標轉換矩陣，可以得到局部坐標系中的節點位移列陣。接下來，為了得到整體坐標系中的剛度矩陣，將局部坐標中的節點位移用整體坐標的節點位移用坐標變化表示，再等式兩邊左乘一個坐標轉換矩陣的轉至，就得到了新的剛度方程，即整體坐標系下的剛度方程。整體坐標系下新的節點位移列陣，節點力列陣，剛度矩陣，也依次得到。

若是3D情況，則依次類推，通過坐標轉換矩陣的方式實現。