有限元理論基礎及Abaqus內部實現方式研究系列27： Abaqus內部計算和顯示的應變

SnowWave02

瀏覽：5729 評論：31 收藏：6

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

==概述==

有限元理論基礎及Abaqus內部實現方式研究系列27： Abaqus內部計算和顯示的應變的圖3 本系列文章研究成熟的有限元理論基礎及在商用有限元軟件的實現方式，通過

（1）基礎理論

（2）商軟操作

（3）自編程序

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

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

有限元理論基礎及Abaqus內部實現方式研究系列27： Abaqus內部計算和顯示的應變的圖4

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

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

有限元理論基礎及Abaqus內部實現方式研究系列27： Abaqus內部計算和顯示的應變的圖5

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

==第27篇：Abaqus內部計算和顯示的應變 ==

應變度量在線性情況下沒有差異，但在幾何非線性情況下有不同的度量方式。系列文章18：幾何非線性的應變中，我們提到了幾何非線性常用的三種應變：工程應變、Green應變、對數應變，其中工程應變和Green應變分別用于線性和小應變情況，對數應變用于幾何大應變情況，那么在Abaqus的幾何大應變中，是否真的就簡單的采用了對數應變呢？Abaqus后處理顯示的應變和計算的應變是一樣的嗎？我們將舉殼單元和體單元兩個例子通過iSolver和Abaqus的結果對比來驗證幾何大應變Abauqs內部計算用的應變和后處理顯示的應變。前面有朋友也問幾何非線性時iSolver計算用的應變是什么？趁此文章我們也將證明iSolver計算時不同的單元用不同的應變，每種應變的選取方式和Abaqus一致。

1.1 對數應變的計算方法

當存在幾何大變形的情況，變形比較大時，此時終止時刻位移x-初始時刻X已經和X比擬了。

對數應變取從初始時刻到最終時刻點這段時間的累加，與材料點在整個時間段內的變形路徑，是一種積分行為。

（1）在一維情況下，設拉伸率為r：

這個等式其實隱含了應變是沿x方向的，如果在三維空間可以認為有一個主方向n=（1,0,0）’，因為應變是二維張量，那么上式就是：

（2）在三維情況下，存在三個主方向n1,n2,n3，拉伸率分別為r1，r2，r3：

剩下的問題就是怎么求這三個主方向，這個主要是求一個與變形F相關的3X3矩陣的特征值和特征向量，具體可參考其它論文書，其實和Abaqus后處理中顯示應變時有三個Principal的值是一樣的求法。

1.2 變形率積分的計算方法

由上面計算方法發現每次計算對數應變都需要求一個特征值和特征向量，在數值計算中，特征問題的求解耗時相對較多，且計算相對復雜（一般人都是認為計算復雜才采用別的應變，個人不太認可），而實際許多非線性材料中，都有這樣一個規律，就是彈性應變都相對較小，譬如典型的鋼材料，楊氏模量為2.1e11Pa，屈服應力為235Mpa，那么達到屈服時的應變為235e6/2.1e11=1e-3，同時，典型的應力應變本構曲線如下圖，在塑性段譬如C點的彈性應變和屈服應變差異并不大。

因此，Abaqus中假定內置的所有材料都滿足彈性應變相對較小，此時，理論可以證明，對數應變可以簡單的取為變形率D的積分：

上式無法直接得到數學表達式，但在有限元中，可以通過增量形式累加。

1.3 調試Abaqus內部計算應變的方法

由于對數應變和試驗最接近，因此Abaqus后處理中的E都是用對數應變來顯示的，Abaqus為了進一步提示對數應變，直接在后處理中如果選了NLGeom=On，應變的顯示從E變味了LE，但Abaqus幾何非線性實際計算的應變并不完全一致。

Abaqus真正計算的變量度量可以通過它的子程序的輸入參數獲取。在Abaqus中，增量步即代表時刻點，可以查看增量點時刻的子程序輸入來猜測Abaqus的內部量描述方式。UMAT子程序中，在材料本構函數中要利用應變增量和當前應力等物理量更新應力，查看UMAT等子程序的接口：

可知其中STRAN和DSTRAN分別表示當前增量步最后時刻的應變全量和增量。具體的介紹也可參考下面視頻講解：

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

變形率D在一維上代表對數應變的導數，但三維上并不是對數應變的導數，這是有很大區別的，同時，可以利用iSolver分別采用上述兩種應變度量和Abaqus子程序接口的結果比對來確認Abaqus計算的應變是哪種度量。所以下面我們將找一個體單元和一個殼單元的例子來驗證到底Abaqus計算和顯示的應變是什么。

1.4 體單元的例子

1.4.1 算例介紹

體單元算例參數如下：

尺寸：5X1X0.1。

材料：Young’s Modulus 1e8， Poisson Ratio 0.3。

左側四個節點固支。

右側四個節點約束位移為5,1,1。

劃分為一個殼單元C3D8R。

幾何非線性開關NLGeom=On，且控制只迭代一次。

1.4.2 Abaqus的應變

Abaqus中采用殼的UMAT子程序進行計算。

（1）顯示應變：Abaqus計算完畢后得到導入結果，在后處理中查看，應變E11=0.6819,E22=5.621e-3如下：

（2）計算應變：Abaqus中采用UMAT子程序，利用我們的子程序調試插件DUS調試UMAT，在Visual Studio中查看dStran的值，發現在計算完應變后，進入UMAT時，顯示如下，可得E11=0.6667， E22=0：

1.4.3 iSolver的應變

猜測體單元中，abaqus顯示的是對數應變，而計算還是用的變形率。為了證明這點，在iSolver以同樣方式建模，只是iSolver中采用自帶材料進行計算，材料參數和UMAT的輸入完全一致。

計算中，iSolver也采用變形率計算方式，得到的應變顯示E11=0.6667， E22=0，如下，可發現和Abaqus完全一致。

1.5 殼單元的例子

1.5.1 算例介紹

殼單元算例參數如下：

尺寸：5X1，厚度0.1。

材料：Young’s Modulus 1e8， Poisson Ratio 0.3。

左側兩個節點固支。

右側兩個節點每個加集中力6.73128E+006，x方向。

劃分為一個殼單元S4R。

幾何非線性開關NLGeom=On，且控制只迭代一次。

1.5.2 Abaqus的應變

Abaqus中采用殼的UMAT子程序進行計算。

（1）顯示應變：Abaqus計算完畢后得到導入結果，在后處理中查看，應變E11=8.528e-1,E22=-5.173e-1如下：

（2）計算應變：Abaqus中采用UMAT子程序，利用我們的子程序調試插件DUS調試UMAT，在Visual Studio中查看dStran的值，發現在計算完應變后，進入UMAT時，E11=8.528e-1,E22=-5.173e-1，調試如下：

可以發現殼單元Abaqus的計算應變和顯示應變一樣，猜測都是對數應變。

1.5.3 iSolver的應變

iSolver中采用自帶材料進行計算，材料參數和UMAT的輸入完全一致。

為了計算和Abaqus完全一致，iSolver也采用對數應變計算方式，得到的應變顯示如下，可發現和Abaqus完全一致。

==總結==

由上可以看到，在實際計算中，對體單元，Abaqus和iSolver都采用變形率積分方式來計算應變，對殼單元，Abaqus和iSolver都采用對數應變。一般理論書都認為Abaqus是因為對數應變計算復雜才采用別的應變，但個人認為應該不是這個原因，因為Abaqus對體單元為了顯示對數應變，依然重新計算了一遍，說明Abaqus體單元采用變形率是有其它原因的，具體什么原因我也沒研究清楚，歡迎探討。

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

snowwave02 From www.yqgqt.org.cn

email: snowwave02@qq.com

以往的系列文章：