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

==概述==

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

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

iSolver介紹視頻：

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

==第18篇：幾何非線性的應變 ==

上一章我們介紹了幾何非線性的物理含義，這章接著介紹幾何非線性中重要的物理量：應變。物體在受到外力作用下會產生一定的變形，變形的程度稱應變。理解應變的物理含義對有限元結果的判別和試驗的對應有重要意義。線性情況下的應變非常容易計算，但如果存在幾何非線性，就不是所有人都知道Abaqus等商軟結果中的應變E代表的是什么了。本章首先從位移、變形和應變的區別說起，然后再具體介紹一下幾何非線性下的應變是如何度量的。

1.1 位移 Displacement

一個物體從初始狀態C0受到外部載荷運動到另一個狀態C，由于材料不能憑空消失，那么C0下的任意的一個物質點a（設坐標為X）必然在C下有個對應點A（設坐標為x）。位移表示兩個位置坐標相減，為u=x-X。

1.2 變形 Deformation

位移無法表示出變形情況，譬如大黃蜂變形為汽車時，身體360度轉動，位移很大，但沒有任何的變形。

變形必須要由兩個點之間材料纖維的拉伸才能決定。對初始構想C0任意的兩點a、b，變形后設為為AB，那么AB矢量和ab之間的關系可以由

                                                                         dx=F*dX

來表示。

其中F定義為deformation gradient matrix。

1.3 應變 Strain

顯然F與X和x的坐標系相關，dX和dx都是在特定三維空間笛卡爾坐標系下的矢量值。它是個張量，有9個量。F在理論上可以表示物體受力后的變形情況，但實際中試驗沒法測出這么多分量，因為拉伸試驗等測試的幾何量只能是位移，長度等有限的幾個標量，這些基本量是無法得到F的全部分量的。所以只能尋找新的量來表示物體變形情況。

一種很自然的想法就是取模，也就是取F的絕對值，稱為拉伸率strech ratio：

                                                               r=||dx||/||dX||

橡膠等材料，變形比較大，r可以達到2以上，那么在試驗上很容易測量，但對金屬等變形比較小的情況，一倍的力和兩倍的力得到的r只是1.0001和1.0002，差異非常小，依然很難標定。

所以，為了測量和計算方便，人為假定了另一個表示物體變形情況的變量，就是應變，因為是人為設定的，所以應變的取法有很多種，主要考慮兩個條件的便捷性：

（條件a）在試驗上可以簡單的表示為位移，長度等基本量，且工況不同時，差異較大，能容易反應出物體加載情況。

（條件b）在理論上可以方便地和位移、應力一起形成最終正確的求解方程。

1.4 應變的度量選擇

為了簡單起見，我們只針對一維情況來說明應變度量選擇，應變的一般公式可以表示為r的一個函數：

                                                                       Strain=f(r)

由于應變是反應物體變形的量，那么當不變形時，也就是lambda=1時，顯然Strain=0。

1.4.1 工程應變

當變形特別小時，也就是r~=1時，假定應變和位移是線性的，那么取伸長量和原始長度的比值。

也就是：

                                                           Strain=u/X=(x-X)/X =r-1

x和X分別為最終和最初的長度。一倍的力和兩倍的力得到的應變變為了0.0001和0.0002，差兩倍關系，很容易區分出來。

此時，顯然條件a很容易滿足，同時，在理論上將應變和位移看成是線性關系去求剛度陣，受力后得到的位移也容易計算，且結果和試驗一致，因此，條件b也滿足。

1.4.2 真實應變

當存在幾何大變形的情況，變形比較大時，此時x-X已經和X比擬了。

此時工程應變中的應變和位移是線性關系這種假定求出的位移已經和實際偏差較大了，因此，需要采用新的應變度量方式。一種選擇是應變取從初始時刻到最終時刻點這段時間的累加，顯然，這種選擇就與材料點在整個時間段內的變形路徑有關。

問題：那么從上圖的X到最終的x，我們取哪個路徑呢？

答：取真實的路徑。

實際情況是什么路徑，就應該取什么路徑，此時在真實路徑下得到的應變就是真實應變。當然真實路徑不一定總是直接從X直線到x點的，譬如下方

(1)  如果受力是FA，那么的確是X直線到x點A路徑；

(2)  但如果是受力FB，那么可能是路徑B，只不過試驗上總是采用FA來做拉伸實驗罷了。

在路徑A下，可以得到

也就是真實應變和拉伸率的對數相等，這也是為什么真實應變也稱為對數或者自然應變的原因。

1.4.3 Green應變

真實應變在一維下很容易計算，但是在三維情況下，Strain和F相關，上面的積分就沒法直接求了，只能得到一個類似

這種表達式求積分，時間的度量由增量來表示，真實應變由應變增量累加得到。因此，在計算機中每次求真實應變必須花費大量時間。

• 注：實際中有限元的幾何非線性求應變增量需要計算局部模態，而模態分析的求解時間成本很大。

一般情況只能老老實實積分計算，但對大位移大轉動小應變的特殊情況，發現應變如果取為Green應變，計算精度不受影響，同時，三維的應變也可以簡單的表示出來。所以對小應變情況，可以采用Green應變，一維形式如下：

                                                                 Strain=0.5（r^2-1）

1.5 殼的應變度量

如果是殼單元，那么用的是局部坐標系，假設為L坐標系，初始構型下存在坐標系變換dX=F2*dL。類似的，幾何大變形中任意一個構形也有類似的變換。此時變形量F或者真實應變就和L坐標系相關了。

具體的也可看一下第十四篇：殼的應力方向或者下面配合iSolver求解器的視頻解說

http://www.yqgqt.org.cn/college/video/c12884 20.14 理論系列文章14：殼的應力方向

==總結==

本章首先從位移、變形和應變的區別說起，然后通過一維的簡單例子具體介紹了幾何非線性下的應變的度量方式。

至于線性的應變在商軟或者自主軟件中在單元中實現，可以參考下面的視頻：

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

深入淺出有限元：基礎理論->Abaqus操作->matlab編程

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

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

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

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