有限元學習必須知道：隱式與顯式有限元演算法

導語：對於

有限元

的學習，最最關鍵的其實不是本構方程，也不是屈服準則，而在於對求解演算法的理解。本人根據長期學習經驗，在美國作訪問學者學習經歷，純屬興趣，在這裡做一點膚淺的總結，完全原創。

隱式

與顯式有限元最大的區別在於是否迭代，是否所有的物理量在同一時刻獲得。採用隱式迭代求解平衡方程（位移、速度和加速度）、而不管是否用隱式與顯式的方法（前向或者後向尤拉求解方法）求解本構方程（應力和應變）叫做

隱式有限元

；用顯式時間積分的方法求解叫做顯式有限元。

首先，對於本構方程的求解，通常分為前向和後向尤拉演算法

。對於後向

尤拉演算法

求解彈塑性問題，所有的物理量（包括等效塑性應變增量、N+1迭代步的應變和應力以及相關依賴於solution的狀態變數）均是同時求解獲得，因為涉及到多個物理量，而通常情況下他們是相互依賴、相互成為函式，所以必須透過牛頓迭代同時求解幾個方程（如採用試應力方程、

屈服函式徑向返回演算法

（對於各向異性，也叫回映演算法，最近點的投射演算法）聯合求解等效塑性應變增量）。對於前向尤拉，直接由N時刻的應力和應變求出N+1時刻的應力和應變，無需迭代。

其次，對於平衡方程的求解，通常分為隱式和顯式有限元演算法

。對於隱式有限元演算法，由應力平衡方程+邊界條件變分之後獲得的剛度方程KU=F，隱式求解必須引入雅可比矩陣（二次收斂、隻影響計算速率、不影響數值精度；K又稱為雅可比），其是實時更新的，是N+1時刻的應力、應變以及狀態變數（如損傷內變數）的函式，隱式求解是很robust的，確保了計算精度，但是不足之處在於計算非常expensive，每次迭代都要計算K的逆矩陣，也容易產生數值收斂性問題，目前解決的方法有弧長法、粘性阻尼法等，個人認為粘性阻尼法效果最好。

對於顯示演算法，採用時間積分，用t+1時刻的積分點應力、應變，獲得t+1時刻的節點位移，無需迭代求解，也不需要雅可比矩陣（應力對應變偏導數）；如果硬是要有，連續雅可比，基於本構模型而不是剛度方程推導近似的連續雅可比。對於顯示演算法，單元高斯積分點應力、應變的求解可用前向或者後向尤拉方法，然後透過時間積分求取節點位移。本質上，平衡方程中位移的迭代求解與本構方程中的應力、應變求解沒有關聯，這點很容易造成誤解，很多時候將前、後尤拉演算法視為顯式和隱式的區別，大大錯誤。通常應用較廣的顯示演算法紐馬克法、威爾遜-sita法，其中改變紐馬克法中的兩個引數，可以實現隱式與顯式求解，其中alpha=0。5和beta=0是中心差分法（二階精度）。目前一個大的誤區認為只有顯示演算法可以求解動力學問題，隱式只能求解準靜態問題（如低速衝擊），alpha=0。5和beta=0。25就是隱式，所有的物理量在t+1時刻同時求解，通常ABAQUS軟體中所說的隱式動力學求解採用了

斯坦福大學Hilbe

r、HUGHES院士（現在

德克薩斯大學奧斯丁分校

）和加州大學伯克利分校Taylor院士提出的無條件穩定隱式差分演算法，可以求解低速動力學問題，缺點是不適合含阻尼的求解、計算效率不高；alpha=0。5和beta=0時的紐馬克法更適合求解

動力學

問題，主要原因在於比隱式求解計算效率更高，不足之處在於其是條件穩定，時間增量過大位移解容易震盪，根本原因是差分演算法的條件穩定導致的，時間增量必須非常小（其值越大，一方面不穩定、另一方面計算誤差也更大），其依賴於波速、

彈性模量

和最小單元網格尺寸，這是顯式演算法計算最耗時的地方。相對於隱式演算法，顯式演算法的功能更強大，通常計算依賴於率的變形和應力，也可以求解穩態問題，如alpha=0。5和beta=0時，對於剛度方程中引入阻尼矩陣後，叫做動態鬆弛法，可以解決靜力學問題的一些穩態問題（如重力、

預應力

引起的初始應力）。此外，一些準靜態的剪下自鎖問題，本質上有解，但是用牛頓法求解失效，中心差分引入質量矩陣後，可以獲得正常的解。需要注意的是，時間積分演算法通常採用Lumped集中對角質量矩陣而不是一致質量矩陣，以提高計算效率。

總體來說，由於計算效率的問題，隱式時間積分演算法ABAQUS-Standard特別適合於低速衝擊問題；對於高速衝擊問題，由於存在不連續非線性接觸的動響應過程，隱式演算法解決不好，使用顯式時間演算法ABAQUS-Explicit更好。此外，對於

瞬態

和穩態熱傳導問題，半離散的拋物線方程，中心差分法可以較好獲得溫度分佈。

對於依賴於率的粘塑性問題（對於本質上的粘性材料），與彈塑性材料的根本區別在於，一般來說是一致性條件不滿足（排除彈塑性材料在高溫下的軟化問題，對於這種問題，屈服條件也可以滿足），即屈服條件不滿足，N+1時刻的物理量不用回映到N+1時刻的屈服面上，粘塑性模型成為過應力模型，顯示和隱式演算法都可以求解。對於依賴於率的本構模型，其可解決模擬高速衝擊、爆炸、彈道射擊問題時存在的動態應變區域性化問題（對於動態問題，平衡方程喪失雙曲線特性；對於靜態問題，平衡方程失去橢圓性），解決網格尺寸效應，其實質上是引入了適當的阻尼遲滯效應。需要注意的是，

對於大變形(又稱為有限變形)問題

，Cauchy應力率和速度梯度（包括客觀和對稱的扭曲張量率D、不客觀和反對稱的spin旋轉張量W兩個部分）均是不客觀的，為解釋剛體旋轉（如純剪下變形就包含剛體旋轉），在共旋座標系下面求解真實應力和應變，應力和應變積分求解的時候應首先求解客觀性的Jaumann應力率（相對於真實應力，空間座標系），相對於Second Piola-Kirchhoff應力是Truesdell率（材料座標系）。

ABAQUS軟體對於大變形問題已經做了旋轉

。

對於一些耦合場問題，由於計算量非常大，同時要求解太多物理量，如熱-流-固耦合，要求解位移、壓力、溫度，採用純隱式演算法或顯式演算法基本不太獲得收斂或準確的解，這時候可採用混合的隱式與顯式有限元格式mixed implicit-explicit partitoning方法，將剛度矩陣和阻尼矩陣分成兩個部分，在同一區域採取不同演算法，提高計算效率和精度、穩定性和收斂性。