（原創，轉載請注明出處）

==概述==

本系列文章研究成熟的有限元理論基礎及在商用有限元軟件的實現方式。有限元的理論發展了幾十年已經相當成熟，商用有限元軟件同樣也是采用這些成熟的有限元理論，只是在實際應用過程中，商用CAE軟件在傳統的理論基礎上會做相應的修正以解決工程中遇到的不同問題，且各家軟件的修正方法都不一樣，每個主流商用軟件手冊中都會注明各個單元的理論采用了哪種理論公式，但都只是提一下用什么方法修正，很多沒有具體的實現公式。商用軟件對外就是一個黑盒子，除了開發人員，使用人員只能在黑盒子外猜測內部實現方式。

                                             

一方面我們查閱各個主流商用軟件的理論手冊并通過進行大量的資料查閱猜測內部修正方法，另一方面我們自己編程實現結構有限元求解器，通過自研求解器和商軟的結果比較來驗證我們的猜測，如同管中窺豹一般來研究的修正方法，從而猜測商用有限元軟件的內部計算方法。我們關注CAE中的結構有限元，所以主要選擇了商用結構有限元軟件中文檔相對較完備的Abaqus來研究內部實現方式，同時對某些問題也會涉及其它的Nastran/Ansys等商軟。為了理解方便有很多問題在數學上其實并不嚴謹，同時由于水平有限可能有許多的理論錯誤，歡迎交流討論，也期待有更多的合作機會。

iSolver介紹：

http://www.yqgqt.org.cn/college/video/c12884

==第13篇：顯式和隱式的區別==

CAE求解方法一般有兩種，分別為顯式（Explicit）和隱式（Implicit）。顯式求解算法基于動力學方程，當前時刻的位移只與前一時刻的速度和位移相關，求解過程中無需迭代；而隱式求解基于虛功原理，一般需要進行迭代計算。

在Abaqus中，顯式求解和隱式求解一般都會采用增量求解，即將分析步分割為若干個增量步，在當前增量步達到平衡時計算下一個增量步。

1. 顯式(Explicit)

在顯式求解過程中，每個增量步內不需要進行迭代求解，且無需形成切線剛度矩陣，故每個增量步內計算量相對于隱式求解方法消耗較小，一般與單元規模成正比。但增量步長也不能過大，一般不超過模型最小自由振蕩周期的1/10，否則容易導致計算結果發散。

2. 隱式(Implicit)

在隱式求解過程中，每個增量步都需要進行平衡迭代，需要形成切線剛度矩陣，計算量相對較大，一般與單元規模和迭代收斂速度相關。隱式求解的收斂速度和穩定性根據選擇迭代方法的不同而不同。因此，需要針對模型特性選擇合適的增量步長，保證計算結果的收斂。

綜上，無論是顯式求解還是隱式求解，都需要根據模型和求解問題合理設置分析步的增量步長和求解方法，保證分析的精度和質量。雖然這兩種求解方法已經是有限元的基本動力學求解方法，但由于有限元本身的復雜性，往往很多人都難以理解兩者的區別和顯式為何發射，本文將以一個簡單的算例配合代碼實現來直觀的解釋一下隱式和顯式的區別。

1.1 前向歐拉（Forward Euler）和后向歐拉（Backward Euler）

前向歐拉和后向歐拉分別是顯式和隱式的一個典型計算方法，本文將引用這兩個方法來盡量直觀地展現顯式求解和隱式求解的區別。

前向歐拉：

fn+1 (x) = fn (x) + h*f'n (x)

后向歐拉：

fn+1 (x) = fn (x) + h*f'n+1 (x)

1.2 算例

現在以一個微分方程算例為例：

y'(x) = -y+x+1

且有y(0) = 1。     

1、前向歐拉法

使用前向歐拉法，可得：

yn+1 = yn + h*(-yn +xn +1)

顯然，這個式子不需要做任何額外運算就能從yn出yn+1，因此每步迭代計算量小。但h在一個最大限制，具體見下面的推導過程：

yn+1 = (1-h)yn + h*(xn +1)

可得

yn+1 = (1-h)2 yn-1 + h*(xn +1) + h*(xn-1 +1)

則得

yn+1 = (1-h)n +O(s2 )

 

當h<2時，y收斂，但h>=2時，y將發散。

2、后向歐拉法

yn+1 = yn + h*(-yn+1 +xn+1 +1)

這個式子不能直接得出yn+1, 必須做進一步計算得到

yn+1 = (yn + h*(xn+1 +1)) / (1+h)

當h>0時，y無條件收斂。

 

3、結果比較

假設n=30

（1）當h = 1.9時，計算結果如下圖所示,顯式和隱式都和理論接近：  

        

（2）當h = 2.1時，計算結果如下圖所示，顯式求解發散了：

 ==總結==

本文簡單介紹了顯式和隱式的區別，本質在于采用不同的物理學平衡方程，因此在不同的物理學問題也有不同的表現。文中給出一個數學算例，并分別利用前向歐拉和后向歐拉求解，以求直觀表現顯式和隱式在求解過程中的差異，以及增量步長對求解結果的影響。

如果有任何其它疑問或者項目合作意向，也歡迎聯系我們：

snowwave02 From www.yqgqt.org.cn

email: snowwave02@qq.com

以往的系列文章：

第一篇：S4殼單元剛度矩陣研究。介紹Abaqus的S4剛度矩陣在普通厚殼理論上的修正。

http://www.yqgqt.org.cn/content/post/338859

第二篇：S4殼單元質量矩陣研究。介紹Abaqus的S4和Nastran的Quad4單元的質量矩陣。

http://www.yqgqt.org.cn/content/post/343905

第三篇：S4殼單元的剪切自鎖和沙漏控制。介紹Abaqus的S4單元如何來消除剪切自鎖以及S4R如何來抑制沙漏的。

http://www.yqgqt.org.cn/content/post/350865

第四篇：非線性問題的求解。介紹Abaqus在非線性分析中采用的數值計算的求解方法。

http://www.yqgqt.org.cn/content/post/360565

第五篇：單元正確性驗證。介紹有限元單元正確性的驗證方法，通過多個實例比較自研結構求解器程序iSolver與Abaqus的分析結果，從而說明整個正確性驗證的過程和iSolver結果的正確性。

http://www.yqgqt.org.cn/content/post/373743

第六篇：General梁單元的剛度矩陣。介紹梁單元的基礎理論和Abaqus中General梁單元的剛度矩陣的修正方式，采用這些修正方式可以得到和Abaqus梁單元完全一致的剛度矩陣。

http://www.yqgqt.org.cn/content/post/403932

第七篇：C3D8六面體單元的剛度矩陣。介紹六面體單元的基礎理論和Abaqus中C3D8R六面體單元的剛度矩陣的修正方式，采用這些修正方式可以得到和Abaqus六面體單元完全一致的剛度矩陣。

http://www.yqgqt.org.cn/content/post/430177

第八篇：UMAT用戶子程序開發步驟。介紹基于Fortran和Matlab兩種方式的Abaqus的UMAT的開發步驟，對比發現開發步驟基本相同，同時采用Matlab更加高效和靈活。

http://www.yqgqt.org.cn/content/post/432848

第九篇：編寫線性UMAT Step By Step。介紹了線性UMAT的接口功能和關鍵接口變量的含義，并通過簡單立方體靜力分析的算例詳細說明了基于Matlab線性UMAT的開發步驟。

http://www.yqgqt.org.cn/content/post/440874

第十篇：耦合約束（Coupling constraints）的研究。介紹了耦合約束的定義和用途，具體闡述了Abaqus中運動耦合約束和分布耦合約束的原理。

http://www.yqgqt.org.cn/content/post/531029

第十一篇：自主CAE開發實戰經驗第一階段總結。結合自研有限元求解器iSolver第一階段開發的實戰經驗，從整體角度上介紹自主CAE的開發難度、時間預估、框架設計、編程語言選擇、測試、未來發展方向等。

http://www.yqgqt.org.cn/content/post/532475

第十二篇：幾何梁單元的剛度矩陣。研究了Abaqus中幾何梁的B31單元的剛度矩陣的求解方式，以L梁為例，介紹General梁用到的面積、慣性矩、扭轉常數等參數在幾何梁中是如何通過幾何形狀求得的，根據這些參數，可以得到和Abaqus完全一致的剛度矩陣，從而對只有幾何梁組成的任意模型一般都能得到Abaqus完全一致的分析結果，并用一個簡單的算例驗證了該想法。

http://www.yqgqt.org.cn/content/post/534362