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

==概述==

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

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

iSolver介紹視頻：

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

==第17篇：幾何非線性的物理含義 ==

幾何非線性，即考慮大變形對于結構平衡位置的影響。爆炸沖擊、沖壓成型、大型結構件的彎曲等都含有幾何非線性問題，幾何非線性也是現代結構有限元商業軟件的必備發展方向。在Abaqus中只要簡單的在Step中勾選NL Geom這個開關就行。

K.J.Bathe教授1979年提出的幾何非線性理論也是目前應用于有限元分析最廣泛的幾何非線性力學。但要想在自主結構有限元程序中編程加入幾何非線性的理論，遠遠不是Abaqus或者Ansys表面上的看起來只要加一個NLGeom=on/off這個開發那么容易。同時，要讓自研程序的幾何非線性做到和商軟結果接近遠比線性或者材料非線性難，我們在iSolver編寫幾何非線性的過程中也發現，除了剛度矩陣的修改，增量迭代法的自動步長選取，收斂準則等都有極大的影響。從本章開始，將介紹幾何非線性的一些理論和Abaqus的實現方式，同時通過iSolver的程序驗證Abaqus的實現方式。本章將介紹幾何非線性的簡單的物理含義，并通過幾何非線性的懸臂梁Abaqus和iSolver的小應變情況的結果，從直觀上理解幾何非線性和線性的差異。配合本章的視頻解說和操作演示可看下方：

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

2.1 幾何非線性的物理含義

2.1.1 從線性到幾何非線性

一個物體從初始狀態A由于受到外部載荷運動，如果現在已知了另一種狀態B的位移，那么其它的任意狀態C的位移怎么求？

如果能直接從B的位移乘以一個常量就得到C，那么這個系統就是線性系統。譬如下面的800mm的懸臂梁問題，在Abaqus中用線性計算，載荷F和位移u是直線關系。

載荷1N的時候Abaqus計算得到最大位移時1.177mm。

那么載荷是1000N時是多少？顯然，不用計算也知道就是1177mm。 但1177這個值明顯有問題，已經超過了梁的長度，按生活經驗判斷這個梁估計都斷了或者極端扭曲了，所以這種情況需要用幾何非線性來計算。

2.1.2 幾何非線性簡單物理含義

虛功原理如下：

在物理上可解釋能量守恒原理，即在某一個時刻點，假定在外力作用下有個虛擬的位移，那么外力在虛擬位移下做的虛功=內部應變能的變化相同。

為了更好的理解上面的物理解釋，如果我們把當成真實的位移，那么外層加上對時間的積分，可以理解為外力在虛擬位移下做的虛功=內部應變能的一段小時間內對應變能的積分：

那么就和我們高中時學過的小球受重力作用后勢能、動能相互轉換是一樣的原理。

舉個簡單的有限元例子，譬如下方一個減縮積分S4R單元的右端受X方向兩個力F：

得到的位移為U（相當于虛位移），那么力F做的功是W=2F*U。

增加的應變能為S*E*V在對所有時刻點t的積分，S、E、V都是當前時刻的值。

應變E的取法有很多種，采用真實應變，那么E取為位移U(t)和長度L(t)的比值，按虛位移的定義，虛位移必然相對原始長度比較小，也就是L(t)=L0+U(t)可以用L0代替，E=U(t)/L0，如果是線性系統，U(t)=U*t/總時間，積分很容易計算出來，得到應變能V=S*U*b*h，因為內力和外力平衡，減縮積分S4R面積內的所有點的應力和中心點一樣，所以S=F/截面積=2F/(b*h)，此時V=2F*U=W。

如果是非線性系統，那么應變就沒法簡單的用E=U(t)/L0，W隨t的變化就是個非線性過程。每個時刻點可以求出一個斜率，這個斜率最終會形成當前時刻點的剛度矩陣。

如果是對當前時刻的體積積分，那么對W求導就很困難，因為V也是與時間有關的，可以選擇一個不變的初始構型V0，此時應力和應變也需要做相應的變化，我們假定分別變為了S和E。

也就是剛度矩陣將分為兩塊，上式的前面一部分依然是以前的BDB形式，只不過B換成了當前時刻的應變位移矩陣，而后面新增項一般稱為幾何剛度陣，在Abaqus中稱為初始應力矩陣（initial stress stiffness）。

2.1.3 幾何非線性的計算機求解方式

理論上受力曲線是一條光滑曲線，計算機沒法求解曲線上每個時刻點的結果，只能求解部分有限間隔點的結果。非線性問題不是一條直線，所以需要多次迭代才能實現，而不再考慮u=F/K這種一次性就能計算的簡單問題。非線性問題求解有多種方法主要分為以下幾類：增量法、迭代法、增量迭代法和弧長法。

具體的理論和Abaqus實現過程可參考我們以前系列文章：第四篇：非線性問題的求解。介紹Abaqus在非線性分析中采用的數值計算的求解方法。

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

2.2 懸臂梁的幾何非線性算例

此次驗證，依然使用自研求解器iSolver與Abaqus計算結果對比的方式。

2.2.1 算例說明

我們采用上面提到的懸臂梁的例子，我們采用殼單元來模擬整個梁。

參數如下：

尺寸：X方向長度L=800，Y方向寬度b=20，Z方向厚度h=20。

材料：Young’s Modulus E=11000。

點載荷：F=1000。

網格：在X方向劃分20份，Y方向劃分8份（為了避免Abaqus和iSolver中都會出現的沙漏現象，厚度方向盡量多劃分網格）。

Step中設置NLGeom打開，并且步長=0.25。

• 我們發現就算不設置這個步長，采用自動步長，Abaqus也依然會自動調整為0.25的步長，為了和iSolver對比，我們都設置0.25，避免因為收斂判據的不同造成的結果差異。

如果線性情況下，上面已經得到了就是1177，現在幾何非線性情況，我們采用小應變的幾何非線性。

• 注：這個問題雖然位移很大，但依然還是個小應變問題，小應變單元的計算速度遠遠高于大應變，所以我們只用小應變單元來模擬。

2.2.2 Abaqus結果

Abaqus中設置單元為S4R5。

得到最后一個增量步的位移為628.4：

同時可以發現應變為0.03，和1相比是小量，所以可以用小應變來模擬。

2.2.3 iSolver結果

在運行結束后，得到的位移為627.9，和Abaqus相差0.08%。

2.3 遺留問題

Abaqus應力判據為：

其中，

為最大應力殘差；

為誤差因子，默認為5e-3。

q為平均時間應力，即.msg中的Average Force值，在Abaqus的幫助文檔中有具體的計算公式，但我們按照公式第一個增量步的結果和Abaqus完全一致，后面的后面的增量步差異很大，這是導致我們的迭代次數和Abaqus沒法完全一致的原因，有知道q第二增量步怎么計算的大神還請不吝賜教。

==總結==

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

視頻解說及算例操作演示如下：

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

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

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。介紹基于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

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

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

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