有限元理論基礎及Abaqus內部實現方式研究系列28：幾何非線性的T.L.和U.L.描述方法

SnowWave02

瀏覽：5970 評論：22 收藏：7

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

1 概述

有限元理論基礎及Abaqus內部實現方式研究系列28：幾何非線性的T.L.和U.L.描述方法的圖3 本系列文章研究成熟的有限元理論基礎及在商用有限元軟件的實現方式，通過

（1）基礎理論

（2）商軟操作

（3）自編程序

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

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

有限元理論基礎及Abaqus內部實現方式研究系列28：幾何非線性的T.L.和U.L.描述方法的圖4

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

通用結構有限元求解器iSolver介紹視頻：

有限元理論基礎及Abaqus內部實現方式研究系列28：幾何非線性的T.L.和U.L.描述方法的圖6

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

==第28篇：幾何非線性的T.L.和U.L.描述方法==

物理空間的很多量和方程都是相對某個參照系或者時刻點的，幾何非線性的方程也不例外。K.J.Bathe教授1979年首次提出了應用于工程問題的幾何非線性理論，他將幾何非線性理論公式分為完全拉格朗格式(T.L. Total Lagrangrian格式)和更新拉格朗日格式(U.L. Updated Lagrangian格式)兩種公式描述方式，物理方程中的所有的物理量統一在一種描述方式下表示，最終建立幾何非線性的求解方程。此后，這兩種描述方式稱為商用結構有限元最廣泛采用的幾何非線性描述方法方法，此文將簡單介紹一下物理量、網格的Lagrangina和Euler描述，虛功原理的T.L.和U.L.描述方法，然后通過和自主結構求解器iSolver的比對，驗證Abaqus的幾何非線性的描述方法。

1.1 初始構型和當前構型

我們先從構型的含義說起，幾何非線性理論由三維連續體的虛功方程，提出了適用大位移、大轉動幾何非線性的解法。若把物體視為由無數質點所構成，并把這些質點稱為材料點，一個物體從初始狀態由于受到外部載荷運動，運動到另一個狀態的過程中，每一個時刻，物體中所有材料點的位置的集合，定義為該物體的一個構型（configuration），且材料點不會憑空消失。

物體的材料點在0時刻占有的區域稱為初始構型C0（initial configuration），材料點的位置在在當前t時刻占有的區域稱為當前構型Ct（current configuration）。譬如下面的Abaqus和iSolver計算的幾何大變形下的懸臂梁，初始時刻，懸臂梁在一個水平平面上，由未變形的體單元組成初始構型。而下壓的任意時刻t，懸臂梁將彎曲，由變形的體單元組成當前構型。

1.2 物理量的Lagrangian和Euler描述方式

為描述各個物理量，可以為初始構型和當前構型選取不同的參考坐標系。當然，一般情況下，無論是初始構型還是當前構型，都以空間中固定不變的全局直角坐標系O來參考，此時，C0下任意一個材料點a在該全局坐標系O下記為X，因為材料點不會憑空消失，在t時刻C1構型下，必然會移動到某一位置A，該點在同一個全局坐標系下的坐標為x。

在數學上，X和x就相當于兩種自變量，對于參考坐標系O中任意的變量y，就有兩種表示方法。

1.2.1 Lagrangian描述

第一種：表示為X的函數：y=f(X)。此時我們稱為Lagrangian（Lagrangian英文后面統一簡稱為L）坐標描述方法。典型物理量：Green應變。

注：Green應變也稱為Green-Lagrange 應變，猜測可能就是因為是Lagrange描述的原因。

1.2.2 Eulerian描述

第二種：表示為x的函數：y=f(x)。此時我們稱為Eulerian（Eulerian英文名后面統一簡稱為E）坐標描述方法。典型物理量：真實應變。

1.2.3 兩種描述分析

這兩者表述的是同一個模型，舉個簡單的例子，譬如材料點隨時間t的運動為x=2*X*t。那么在坐標系O下，在t時刻的速度的Lagrangian和Eulerian描述方法分別為：

v(X)=2*X和v(x)=2*x/2t=x/t

這兩種描述方法都是相對同一坐標系O的物理量的表示方法，那么得到的數值結果完全一致，只不過寫成函數的時候自變量不一致一樣。

1.3 網格的Lagrangian和Eulerian描述

在講T.L.和U.L.之前，我們先解釋一下L.代表的Lagrangian網格的含義。前面所述物理量的坐標描述也有L和E之分，網格也有L和E之分。這兩者是不同的概念。在L和E網格定義中，表示的是兩種網格隨時間變化的位置定位方式，而物理量的坐標描述僅僅是變量的表示方式，這兩者沒有直接關系。為了說明L和E網格的區別，舉一個簡單的例子，譬如二維平面的運動如下：x=X+tY, y=Y，那么一個長方形運動后的變形為一個平行四邊形，具體的材料點a(在全局坐標系O中坐標為X，Y)變形后為A（在O中坐標為x，y）。

網格主要由節點組成，而節點在全局空間O中的坐標決定了網格的形狀。節點的坐標有兩種取法：

1.3.1 網格的Lagrangian描述

無論在哪個時刻，取同一個材料點變形后的xi，yi坐標，此時如果一開始網格點在a點，那么變形后必然在A點，而且，如果用L描述，那么可以得到，Xi，Yi=常數。那么得到的網格將類似下方，隨構型變形而自動更新變化。這就是Lagrangian網格。

1.3.2 網格的Eulerian描述

無論在哪個時刻，取同一個材料點變形前的Xi，Yi坐標，此時如果一開始網格點在a點，那么變形后還是在a點。那么得到的網格將類似下方，固定在空間不變。這就是Eulerian網格。

結構有限元中一般都是Lagrangian網格，這也是為何T.L.和U.L.方法稱為Lagrangian方法的原因。

1.4 虛功原理的T.L.和U.L.描述方式

有限元中，T.L.和U.L.用于非線性問題中的虛功原理的描述方式，是整個有限元求解方程的描述方式，而前面物理量的描述方式是虛功原理中各個量的描述方式。

非線性問題在有限元中都是以增量法來求解，增量法的基本思想就是將施加外載荷的過程分割為若干個時間增量步，假定這個時間增量步是固定的，dt,2dt,3dt,….，每一步只施加一個比較小的載荷增量，總載荷引起的位移就是所有增量步位移增量的累加。

具體的理論和Abaqus實現過程可參考我們以前系列文章4：非線性問題的求解。

由虛功原理在(n+1)dt時刻應變能為：

上面每個量都是都是參考(n+1)dt時刻的坐標的，此時將相對參考時刻的坐標寫在左下方，將要求的時刻點寫在左上方，那么就是：

1.4.1 T.L.描述

注意，有些讀者認為T.L.只能用于小應變，但這種說法個人認為并不正確，這里的推導沒有涉及是否是大應變還是小應變，同時，后面文章我們將證明U.L.和T.L描述是可以等價轉換的，既然U.L.可用于大應變，那么T.L.也適用于大應變、大轉動、大位移情況。

1.4.2 U.L.描述

以當前時刻ndt構型為參考構型，用Eulerian坐標x描述物理量，那么就是U.L.方法（Updated Lagrangian方法）

上式的轉換數學上并不嚴格，還差一個應變對時間的導數，嚴格推導和系列文章22：幾何非線性的剛度矩陣求解類似，需要將上述表達式先做線性化。

1.5 Abaqus的內部描述方式驗證方法

在Abaqus所有的幫助手冊中，找不到任何T.L.還是U.L.的信息，由上節所述，我們只能反過來，通過查看Abaqus的變量的度量方式，來猜測Abaqus用了T.L.還是U.L.方式，Abaqus文檔中多處提到了各種變量的度量，但并不一定真實，而且前后明顯矛盾，最好的方式還是直接查Abaqus程序的結果。

注意，此處不一定能采用后處理模塊中顯示的應變量。為了和試驗對應，Abaqus幾何非線性后處理顯示是自然應變和Cauchy應力，它將程序中的相關度量在顯式時做了轉換，與程序中的真正用于計算的度量并不一致。詳細說明可以查看本系列文章27： Abaqus內部計算和顯示的應變。

我們通過UMAT等子程序來查看Abaqus子程序接口中的應變應力度量，同時和iSolver采用同樣度量的結果比對，大體猜測Abaqus幾何非線性的描述方式。

1.5.1 Abaqus中T.L.的描述

Abaqus中S4R5、STRI3、STR65、S4RS、S8R5、S3RS、B33，B23等單元都采用T.L.的描述方式。采用一個簡單的算例，證明如下：

1.5.1.1 算例介紹

參數如下：

尺寸：5X1，厚度0.1。

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

左側兩個節點固支。

右側兩個節點每個加集中力1e5，x方向。

劃分為一個單元。

1.5.1.2 Abaqus應變結果

Abaqus中采用S4R5單元，幾何非線性單元NLGeom打開，得到一個迭代步后結果：

1.5.1.3 iSolver應變結果

iSolver中也采用S4R5的度量方式，同時采用T.L.方式，得到的應變度量如下，可發現，和Abaqus完全一致。

1.5.2 Abaqus中U.L.的描述

Abaqus中C3D8/C3D8R、S4/S4R、S3/S3R、B31，B21等非線性單元都采用T.L.的描述方式。采用一個上述同樣的算例，證明如下：

1.5.2.1 Abaqus應變結果

幾何非線性開關NLGeom=On，同時單元類型改為S4R。得到應變：

1.5.2.2 iSolver應變結果

iSolver中也采用S4R的度量方式，同時采用U.L.方式，得到的應變度量如下，可發現，和Abaqus完全一致。

1.6 T.L.和U.L.描述方式的本質

既然Bathe教授在開辟幾何非線性時就創立了T.L.和U.L.兩種虛功原理的描述方式，而且有限元商軟也是按這兩種方式實現的，那Abaqus中為何從頭到尾沒提T.L.還是U.L.的描述方式呢？

我們猜測這是Abaqus有意為之。這兩種虛功原理的描述本來就是描述的同一個物理對象，所以理論上你采取任何一種都可以得到最終的正確解，但無論哪種描述都需要它用到物理量的描述方式結合起來滿足一定的要求，譬如功的共軛、能量守恒、動量守恒等條件。同時，因為有限元中物理量的描述方式可以相互轉換的，為了有限元的結果更能和實際試驗的度量對應，Abaqus也會在T.L.中采用Euler描述的物理量的度量方式或者U.L.中采用Lagrangian描述的物理量的度量方式，譬如它本身提供的UMAT是真實應變和真實應力，但你完全可以不用管UMAT輸入的DSTRAN等應變值，而用它提供的另一個參數變形梯度DFGRD1按Lagrangian描述方式實現Green應變（當然，這樣得到的Jacobian矩陣DDSDDE和應力STRESS都沒問題，但你改不了Abaqus內部利用STRESS來計算初始應力剛度和本身的材料剛度矩陣的算法，所以還需要Abaqus后面提供更多的接口才行），或者在UEL采用Green應變和2nd PK力來實現剛度陣AMATRX和內力RHS，也就是說已經突破了傳統的T.L.和U.L.的界限，所以干脆Abaqus就不再提T.L.還是U.L.方式了。