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

==概述==

本系列文章研究成熟的有限元理論基礎及在商用有限元軟件的實現方式，通過

（1）   基礎理論

（2）   商軟操作

（3）   自編程序

三者結合的方式將復雜繁瑣的結構有限元理論通過簡單直觀的方式展現出來，同時深層次的學習有限元理論和商業軟件的內部實現原理。

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

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

自主結構有限元求解器iSolver介紹視頻：

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

==第22篇：幾何非線性的剛度矩陣求解==

幾何非線性在界面上是很容易設置的，但商軟內部的處理相當復雜，我們從最基本的剛度矩陣的求解出發，看看在幾何非線性設置后，剛度矩陣具體是怎么實現的。本文首先介紹幾何非線性下的剛度矩陣的理論推導和計算機求解方法，說明理想的求解方式的困難點和猜測Abaqus內部的解決方法。最后利用一個簡單的算例通過對比iSolver和Abaqus的結果，部分驗證我們對Abaqus幾何非線性的剛度矩陣的實現方式的猜測。

1.1 幾何非線性的剛度矩陣推導理論

在前面17章：幾何非線性的物理含義中，我們提到如果是非線性系統，應變能W隨t的變化就是個非線性過程。每個時刻點可以求出一個斜率，這個斜率最終會形成當前時刻點的剛度矩陣。

求導后得到的剛度K：

也就是剛度矩陣將分為兩塊：

（1）   上式的前面一部分稱為材料剛度陣，依然是以前的BDB形式，只不過B換成了當前時刻的應變位移矩陣

（2）   后面新增項一般稱為幾何剛度陣，在Abaqus中稱為初始應力矩陣（initial stress stiffness）。

1.2 幾何非線性的剛度矩陣計算機求解

1.2.1 理想的求解方式

理論上受力曲線是一條光滑曲線，計算機沒法求解曲線上每個時刻點的結果，只能求解部分有限間隔點的結果。非線性問題不是一條直線，所以需要多次迭代才能實現。

非線性問題求解有多種方法主要分為以下幾類：增量法、迭代法、增量迭代法和弧長法。

 

在增量迭代法，此時時刻t和t+Δt分別表示為當前增量步的開始和結束時間。此時W求導表示為增量的形式：

上面的應力應變理論上將都應該取當前增量步的結束時間的結果。先將應力增量轉換為應變和本構關系，得到：

1.2.2 解題困難點

1.2.3 困難點的解決方法

解決上述困難點的基本原則：對一個非線性問題來說，在增量步中的斜率K即使計算不精確，只要不是偏差太大，依然收斂，而且滿足精度范圍。K只影響迭代的次數，不影響結果的精度。關于這個原則，具體的解釋和算例驗證可以看我們下面的視頻：

http://www.yqgqt.org.cn/college/video/c13034 Abaqus用戶子程序UMat詳解與開發工具（未完，待續）：2.4 整體理論02-Jacobian和應力矩陣的含義

當然這個K的近似也不能偏差太大，要不然迭代次數太多或者根本無法收斂。本著這個基本原則，兩個困難點的解決方法分別如下：

1.2.4 Abaqus中的內部實現方式猜測

當然，我們不知道Abaqus的內部計算方式，但可以通過兩個外部接口可以猜測Abaqus的內部實現方式。

（1）困難點1：應變增量與增量步結束時刻的位移增量的關系，Abaqus可以通過單元函數實現，內部單元函數接口和Abaqus的子程序UEL的接口是一致的，如下表示：

其中，Abaqus的UEL的U和dU表示的是當前增量步最后時刻的全量和增量的緣故。剛度矩陣只與增量步的量有關，迭代步的量只是對增量值的修正，不直接反應到剛度矩陣中，固UEL輸入參數沒有迭代相關的任何量。Abaqus無論是內部單元還是外部單元都是通過輸入增量步最后時刻的位移全量和增量來求應變增量。具體的UEL說明可參看下方視頻：

http://www.yqgqt.org.cn/college/video/c14948 深入淺出有限元：基礎理論->Abaqus操作->matlab編程

 （2）困難點2：應力增量的求法通過材料函數實現，內部材料函數接口和Abaqus的另一個UMAT接口是一致的，UMAT接口如下：

其中的上述四個關鍵參數中，輸入的STRESS是增量步開始時刻的應力全量，STRAN和DSTRAN分別是增量步開始時刻和結束時刻的應變全量和應變增量，需要通過這三個量來求出結束時刻的應力增量，累加到STRESS上，從而輸出到程序中，而求解的過程猜測無論是Abaqus內部材料還是UMAT都采用了開始時刻的材料本構關系。譬如下方為塑性材料的UMAT的應力增量求解的matlab代碼。

1.3 算例驗證

從UEL和UMAT兩個接口只能猜測Abaqus內部單元幾何非線性的剛度矩陣的實現方式，為了證明我們的猜想，我們創建了一個幾何非線性的懸臂梁的簡單算例，采用iSolver和Abaqus分別計算，而iSolver的兩個困難點按我們前面所述實現，最終的位移結果如下：

兩者結果相差0.08%，基本一致，從而部分驗證了本文對Abaqus幾何非線性下剛度矩陣實現方式的猜測，不能說完全驗證，畢竟幾何非線性除了上面提到的兩個困難點，想要和商軟完全一致，涉及的問題很多，但至少在如果上面兩個困難點的解決方法和Abaqus不一致，那么得到的結果也不會一致。該算例的視頻解說及iSolver操作演示如下：

http://www.yqgqt.org.cn/college/video/c12884 20理論系列文章17-幾何非線性的物理含義。

1.4 聯系方式

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

SnowWave02 From www.yqgqt.org.cn

email: snowwave02@qq.com

1.5 以往的系列文章

以往的系列文章：

1.5.1 ========第一階段========

第一篇：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。介紹基于Matlab線性零基礎，從零開始Step by Step的UMAT的編寫和調試方法，幫助初學者UMAT入門。

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

第十篇：耦合約束（Coupling constraints）的研究。介紹Abaqus中耦合約束的原理，并使用兩個簡單算例加以驗證。

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

1.5.2 ========第二階段========

第十一篇：自主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

第十三篇：顯式和隱式的區別。介紹了顯式和隱式的特點，并給出一個數學算例，分別利用前向歐拉和后向歐拉求解，以求直觀表現顯式和隱式在求解過程中的差異，以及增量步長對求解結果的影響。

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

第十四篇：殼的應力方向。簡單介紹了一下數學上張量和Abaqus中殼的應力方向，并說明Abaqus這么選取的意義，最后通過自編程序iSolver來驗證殼的應力方向的正確性。

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

第十五篇：殼的剪切應力。介紹了殼單元中實際的和板殼近似理論中的剪切應力，也簡單猜測了一下Abaqus的內部實現流程，最后通過一個算例來驗算Abaqus中的真實的剪切應力。

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

第十六篇：Part、Instance與Assembly。介紹了Part、Instance與Assembly三者之間的關系，分析了Instance的網格形成原理，并猜測Abaqus的內部組裝實現流程，隨后針對某手機整機多part算例，通過自編程序iSolver的結果比對驗證我們的猜想。

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

第十七篇：幾何非線性的物理含義。介紹了幾何非線性的簡單的物理含義，并通過幾何非線性的懸臂梁Abaqus和iSolver的小應變情況的結果，從直觀上理解幾何非線性和線性的差異。

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

第十八篇：幾何非線性的應變。首先從位移、變形和應變的區別說起，然后通過一維的簡單例子具體介紹了幾何非線性下的應變的度量方式，并給出了工程應變、 真實應變、Green應變三者一維情況下在數學上的表達方式。

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

第十九篇：Abaqus幾何非線性的設置和后臺。首先介紹了幾何非線性一般的分類，然后詳細說明了Abaqus中幾何非線性的設置方式和常用單元的分類，最后以一個殼單元的簡單算例為對象，可以發現應變理論、Abaqus和iSolver三者在線性、小應變幾何非線性和大應變幾何非線性三種情況下都完全一致，從而驗證Abaqus幾何非線性后臺采用的應變和我們的預想一致。

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

第二十篇：UEL用戶子程序開發步驟。本文首先簡單的討論了UEL的一般含義，并詳細的介紹了基于Fortran和Matlab兩種方式的UEL的開發步驟，對比發現開發步驟基本相同，但Matlab更加高效和靈活。

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

1.5.3 ========第三階段========

第二十一篇：自主CAE開發實戰經驗第二階段總結。從實戰角度介紹自主CAE在推廣和工程化應用的過程中的體會，同時說明一個CAE平臺最重要的兩個特點：可擴展和易維護。

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