FEM的數值精度與高階單元

21 人贊了文章

之前聽過清華FEM的Mooc前九章，知道了FEM的最基本的概念，介紹了FEM的流程，簡單的單元研究，但實際上開始用ABAQUS操作以後，發現在ABAQUS中，常常採用的單元，比C3D20R這種三維二階減縮積分單元，分析結果是摸不清楚的，看了庄茁老師講的實體單元，有了一些頭緒，大致懂了什麼叫剪切自鎖，什麼叫沙漏，但也並沒有給出一些更具體的數學表達，於是回頭看看曾老師講的第十和十一章，有關數值問題的精度和高階單元的介紹，補完這一系列的課程，再加一些相關的筆記，順便推薦一本書《ABAQUS常見問題解答》。

在正式寫筆記之前，先補一下參數單元的知識：

參數單元就是長得不規則的單元，因為真實物理性狀複雜，不可能每個都長得四四方方。對於不同的參數單元，希望從規則單元通過變換來得到參數單元，以2D問題為例，首先引入兩個坐標系，物理系和基準系，基準系是建立在單元上的（可變），而物理系是建立在空間上的（不變），將建立在不同單元上的基準系得到的變數映射到建立在空間上的系中，可以增加不同形狀單元得到變數的通用性，下面分析需要哪些變換：

由此可見，需要變換的有：

<1>坐標；<2>偏導數；<3>積分變數（面積微元）；

下面首先來說<1>坐標的映射：

學過斜角直線坐標的張量分析，就應該對坐標變換非常熟悉了，現在空間中有一個位置P，可以用一個向量P來表示，這個向量在不同的坐標系中有不同的坐標，即：

並且分量之間一定存在唯一確定的轉換關係：

首先，這些都是理論上行的通的，然而，在這裡我們不用精確的表達式演算法，而用插值的方式來獲得每一個基準繫上的點在物理系中的坐標（以2D4Node為例）：

回憶一下形函數的構造，根據自由度數確定待定係數，插值函數從Pascal三角形中獲得，這裡是完全相同的道理，所以選取的坐標插值函數為：

接下來，通過四個節點處的八個坐標值來確定坐標插值函數中的待定係數：

其實仔細考察，便會發現，這和2D4Node的形函數是完全一致的，只不過一個是通過節點位移插值出位移場，一個是在坐標變換中通過節點坐標插值出整個場的坐標。

最後關於坐標變換，給出整個場的變換表達式：

關於<2>的偏導數變換，相信學過高等數學的人都會，引入雅各比矩陣，直接上圖不費話了：

最後是<3>積分變數的變換，或者說面積微元的變換，這裡需要用到一些行列式和外積的幾何性質，以及微分幾何的一點內容：

第一個表達式，可以從幾何意義入手，叉乘的大小代表平行四邊形的面積，這個不是絕對值，也不是行列式，是給一個向量取模。那麼具體應該如何實施呢，通過將兩個向量表達成以i和j為基地的空間中，再進行外積，我寫一下我的理解，從幾何角度分析的（如果不對請指正）：

到這裡，三個代換就做完了，這裡在第一個坐標變換中引入了坐標的插值，根據坐標的插值和位移的插值的相對大小，參數單元分為：

下次的補充內容再說高斯積分，下面是正式的課程記錄筆記。

（1）先說第十章的前三節，和文章題目關係不大，簡單做一下筆記。分別是講了帶寬（見下圖），形函數與剛度矩陣的性質，邊界條件的處理，研究總體剛度矩陣的帶寬是為了節省計算機的存儲空間的，我目前並不怎麼關係怎麼存儲；邊界條件的處理，是講除了我們在矩陣位移法裡面對於零邊界劃掉行列以外，為了保持原有剛度矩陣各行各列的編號，採用了一些利於計算機編程的處理思路，比如置一法，乘大數法，罰函數法，並且他們可以處理非零邊界，這個在有限元軟體上是已經實現過的，我目前不自己編程，也不怎麼需要了解；另外的就是剛度矩陣的性質，其實考試前早都背過了，再寫一寫，先說單元剛度矩陣，性質有一下幾條：Kii代表在i處產生單位一的位移時，在i節點所施加的力的大小；Kij（非對角線元素）表示在j處產生單位一的位移，在i處所施加的力的大小；害怕第二條記混淆吧，沒關係，第三條告訴我們們Kij=Kji，即剛度矩陣是對稱陣；然而這個矩陣也是奇異的，從求解剛度矩陣前需要帶入邊界條件就可以得知，剛度矩陣這個平衡方程是允許產生剛體位移的；既然允許產生剛體位移，這個矩陣就不能是正定的，因為產生的如果是剛體位移，其應變能應該依然為零，可見，剛度矩陣具有半正定的特徵；最後回到剛度矩陣的由來，推導它通過的是最小勢能原理，最小是能原理是平衡方程在變分條件下轉化而來的，剛度矩陣因此具有平衡的特性，在每一行對剛度矩陣的元素求和，結果均為零，也不難通過這個性質得到它的行列式為零，又反推出奇異性。另外還有形函數的性質：Ni的含義是i處產生單位一的位移，其他位置不產生位移時整個的位移場；對Ni進行求和，結果為常數一，感覺像是權函數對吧，實際上是通過帶入剛體位移得到的這個性質。

（2）第四節重點來了，關於收斂性性討論。收斂的意思是指「單元的尺寸趨向於0時，FEM所給出的解趨向於真實解，不僅要討論單元內部，還要討論單元與單元之間」，另一方面，我們應該盡量使得結果單調收斂。

思路很明確，我們從收斂這個結果，一步一步的推出它的必要條件，第一個是單元尺寸趨向於零時，勢能的表達式依然存在，首先對勢能的表達式進行展開：

位移函數必須連續，當體積趨向於零時，省略的高階部分都趨向於零，為了使得勢能存在，應該存在常應變項與常位移項，另外，位移函數在單元之間保持連續，才能保證應變的存在。

然後就有兩個收斂準則要給出來：

<1>完備性要求：如果在勢能泛函中所出現位移函數的最高階導數為m階，則單元內部位移場函數應該為m階完備的多項式；

<2>協調性要求：如果在勢能泛函中所出現位移函數的最高階導數為m階，則單元邊界位移場應該具有m-1階連續導數；

其實，之前在我們確定位移模式的時候已經默認了這些要求，只不過沒有具體的提出來，現在給出位移模式中函數的選取規則：

<1>唯一確定性原則：節點的自由度個數等於待定係數個數，從低階到高階進行選取，選取完備的多項式提高精度（參考帕斯卡三角形）；

<2>單元內部的完備性：包含常數項（常位移）和完備的一次項（常應變）；

<3>單元之間的協調性：對於 連續性單元（存在一階連續導數）較難保證；

（3）第五節， 第六節，關於 單元與 單元，說一下我的理解，這裡說的 型單元，是針對協調性要求來講的，由模型的勢能泛函決定。對於平面問題和空間問題，勢能泛函中出現了應變，而應變為位移的一階導數，所以只要求在邊界上是零階導數連續，也就是函數本身連續，所以平面和空間問題都是 型單元，一階單元和二階單元都滿足 連續性。對於梁單元，板和殼單元，勢能泛函中出現了位移的二階導數，則在邊界上要求一階導數連續，梁，板，殼屬於 型單元。單元的完備性往往容易滿足，但是對於 型問題的協調性不太容易滿足。滿足完備性和協調性要求的單元稱之為協調單元，協調單元在尺度趨向於零時，FEM分析的結果趨向於真實的結果。

（4）第十章後面幾節課得到的結論：

<1>FEM算出的位移整體上偏小；

<2>控制誤差的h方法和p方法，h方法即增大網格密度，p方法即增加網格的多項式階次；

後面講的h和p方法我們可能每天都在使用，只是不確定它的名字。

總結一下，今天寫了參數單元的變換原理，以及第十章的內容，關於有限元收斂性，精度和誤差控制，想一下學懂真的很難，日後再仔細體會吧。下一次補上高斯積分和十一章的高階單元。大神門多多指點我們初學者吧。