各向異性屈服準則的UMAT子程序二次開發研究

喬順成，吳建軍，展學鵬

（1.中航工業西安飛機工業集團有限責任公司，陜西 西安 710072；2.西北工業大學 機電學院，陜西 西安 710072）

ABAQUS以其強大的非線性迭代計算和前后處理功能而廣泛應用于材料彈塑性有限元分析，它是現階段在各個工程領域廣泛使用的大型通用有限元軟件之一。ABAQUS[1]有大量的單元庫和求解模型供用戶使用，而且用戶也能夠通過這些模型求解絕大多數問題。但是實際問題是相當復雜的，ABAQUS不可能直接處理所有可能的問題[2-4]，所以，有必要給用戶提供二次開發的接口來解決實際中出現的新問題。在用戶自定義材料（UMAT即 User-defined Material Mechanical Behavior）研究中，很多學者通過UMAT子程序二次開發，將新的本構模型用于實際問題的有限元分析。例如：李平[5]等通過數值分析模擬了普通碳鋼連續熱軋的過程，將率相關的各向同性硬化熱變形本構模型通過UMAT子程序用于有限元仿真，分析了軋制過程中的溫度和應力、應變之間的耦合關系。楊曼娟[6]提出了五參數的柳玉起屈服準則二次開發開發程序，并將Rankine準則的Mohr-Coulomb模型通過UMAT子程序嵌入ABAQUS，最后通過單元測試和實際工程算例驗證了二次開發的正確性和適用性，惟一不足的就是作者并沒有附上二次開發的Fortran程序代碼。唐榕蔚[7]將雙剪統一強度理論通過UMAT子程序嵌入到 ABAQUS軟件中，用于巖土材料的彈塑性有限元分析，作者最后還附上了二次開發程序代碼，為以后新本構模型的二次開發研究者提供了參考和指導。除此之外，還有許多學者通過有限元二次開發來解決一些實際的工程問題[8，9]。

ABAQUS嵌入了Mises各向同性屈服準則、Tresca各向同性屈服準則、Hill48各向異性屈服準則，前兩個各向同性屈服準則只能用于各向同性有限元分析，Hill48各向異性屈服準則可以用于各向異性有限元分析。但是，Hill48屈服準則形式簡單，對屈服面描述粗略，不涉及材料的晶體結構理論和晶體塑性理論，僅對r＞1（r指的是材料厚向異性指數）的各向異性材料屈服狀態和塑性變形預測較好[10]，所以符合Hill48屈服準則的各向異性有限元分析必然存在局限和不足，不能精準地預測復雜屈服狀態，不能精確地描述復雜的塑性變形，許多新特征、新規律不能被發現。例如，倪向貴等[11]將Yld89、Yld91屈服準則本構模型通過二次開發應用于有限元分析，研究了符合不同各向異性屈服準則（Hill48、Yld89、Yld91）本構關系的方盒件拉深成形，并將三種屈服準則的模擬結果與實驗進行了比較，研究證明符合Hill48屈服準則的方盒件拉深成形仿真與實驗結果相差最大，成形精度最差。

為了彌補之前屈服準則的缺陷和不足，Barlat等[12]提出了Yld2004-18p屈服準則，該準則需要RDTD平面上沿軋制方向每15°的單向屈服應力σY、厚向異性指數r，還有雙向拉伸屈服應力σb和厚向異性指數rb，可以精確地預測拉深成形制件不同方位6或8個“制耳”，也可以精確地預測面內流動應力和厚向異性指數的變化，更加全面地反映了塑性流動各向異性和應力-應變響應，是目前描述FCC（面心立方）和BCC（體心立方）材料各向異性性能最準確的屈服準則之一。目前，大多數有限元軟件都沒有嵌入這個屈服準則本構模型。

為了使得材料加工工藝和塑性成形有限元分析更加精確，更加精準地預測屈服狀態，更加合理地描述塑性變形行為，拓展ABAQUS在各向異性有限元分析領域的應用，很有必要利用ABAQUS提供的UMAT子程序二次開發接口，建立編譯運行環境，采用Fortran語言編程，將高級各向異性屈服準則本構模型嵌入ABAQUS，并應用于彈塑性有限元分析。

本文構建了Yld2004-18p各向異性屈服準則本構模型，計算了Yld2004-18p屈服準則的各向異性系數，通過UMAT子程序二次開發，將Yld2004-18p各向異性屈服準則本構模型嵌入ABAQUS，結合有限元模型進行實例分析，驗證UMAT子程序二次開發的正確性，同時提出一個通用的、柔性的二次開發結構模式。

1 Yld2004-18p各向異性屈服準則及其系數

為彌補之前的Yld89、Yld91、Yld96屈服準則的不足，F.Barlat等[12]提出了更高級的屈服準則——Yld2004-18p屈服準則，該準則給應力偏張量s又加了兩個線性變換，這兩個線性變換共包含有18個各向異性系數。該屈服準則如下：

式中：指數m與材料的晶體結構類型有關，對于BCC（體心立方）材料，m=6，對于FCC（面心立方）材料；S?′，S?″是與應力偏張量s?′，s?″有關的主值；s?′，s?″由應力偏張量s通過線性變換得到，如下式：

式中：C′和C″是應力張量線性變換矩陣，如下式：

T也是線性變換，如下式：

當C′=C″，Yld2004屈服準則等同于Yld91屈服準則；當18個各向異性系數ci=1～18全部為1，m=2（或4）時，它也可以退化到Mises各向同性屈服準則。

2 Yld2004屈服準則本構模型的UMAT子程序二次開發

2.1 與Yld2004屈服準則有關的變量推導

構建彈-塑性本構模型的首要條件是確定屈服函數，及與屈服函數有關的變量，它們與硬化法則、應力更新算法是構建彈-塑性本構模型的必要條件。為了將Yld2004屈服準則本構模型通過有限元二次開發--UMAT子程序嵌入ABAQUS軟件，屈服條件、屈服函數對應力分量σij的一次導數、二次導數是必不可少的。它們在整個彈-塑性本構模型發揮著關鍵性的作用，是更新迭代步增量（例如應力增量、應變增量、塑性應變增量、切線模量、塑性參數增量等）、輸出變量的重要環節。

對于Yld2004屈服準則本構模型，結合等式（1），可以得到等效應力σˉ的表達式：

結合等式（1）、（7），等效應力σˉ對應力分量 σij

其中：

等式（9）是Φ對σij的“鏈式求導”，可結合等式（1）-（6）求得。等效應力σˉ對應力分量 σij的二階導為：

其中：

以上這些公式將會在下文應用于UMAT子程序本構模型。

2.2 彈-塑性率本構方程

材料屈服狀態的判斷根據是復雜應力狀態下的等效應力與參考方向的單向屈服應力之間的關系，如此可以建立屈服條件：

從初始加載到屈服之前，材料發生彈性變形，符合廣義Hook定律，彈性響應對應于彈性應變率，滿足如下關系：

式中：

λ*，μ是獨立的材料常數，稱為拉梅（Lame）常數，可以用更接近于物理度量的常數表示它們：彈性模量E、泊松比v，如下式：

塑性應變率由流動法則確定，常表示為塑性流動勢能Ψ的形式：

令：

從幾何角度上，塑性流動方向與屈服面的法線方向相同。當塑性加載時，λ˙＞0，應力保持在屈服表面F=0，也可以用一致性條件表示：F˙=0，通過“鏈規則”擴展，推導得到如下：

2.3 UMAT子程序基本流程的構建

由于ABAQUS不包含高級的Yld2004屈服準則本構模型，所以將高級的Yld2004屈服準則本構模型應用到球形壓痕有限元分析，就必須開發相應的屈服準則本構模型UMAT子程序。

根據ABAQUS提供的二次開發接口要求，用Fortran語言編寫程序，得到.for格式的UMAT子程序，將其嵌入ABAQUS，就可以實現調用。UMAT子程序二次開發使用戶能夠自定義ABAQUS材料模型庫中沒有的材料模型，分析實際中更具體、更復雜的新問題。UMAT子程序的功能非常強悍，可以根據材料特性，自定義新的材料本構；可以實現新問題的有限元分析；它與主程序相互調用，可以用于任何加載或卸載階段的分析；可以在材料單元每個積分點上進行調用等。

上述第2.1、2.2節已經完成了UMAT子程序本構模型需要的主要變量、彈塑性率本構方程的推導過程。為了方便編程，直觀地表示迭代循環的步驟，需要構建UMAT子程序本構模型，將第2.1、2.2節與完全隱式向后Euler圖形返回算法結合，它的簡要過程如下：

（1）初始狀態：

初始值設置：k=0，Δλ=0；

（2）在第k次迭代檢查屈服條件和收斂性：

如果 F（k）≤tolerance 且‖R（k）‖≤tolerance，則儲存變量，此時間步結束；否則，繼續下一步驟；

（3）在第k次迭代，計算中間變量：

（5）更新變量：

返回第（2）步。

2.4 UMAT子程序結構模式及編程

作為獨立的程序模塊，UMAT子程序既能被ABAQUS主程序調用，也能被單獨編譯或存儲。UMAT子程序的編程應該遵循結構化程序設計的思想，即：靈活使用三種基本程序結構（順序結構、循環結構、選擇結構），按照自頂向下的設計思路，分解程序功能，使程序合理的模塊化、分支化，可使子程序足夠簡單，柔性化，具有互換性，更加通用。

本文開發的Yld2004屈服準則本構模型UMAT子程序，從整體上看，都具有以下模塊：①ABAQUS規定的題名聲明語句；②ABAQUS的接口參數聲明語句；③自定義的變量聲明語句；④程序主體，包括主程序和5個嵌套的功能分支子程序；⑤程序結束語和返回語句。在開發的UMAT子程序中，一致切線模量矩陣C的求逆運算、與材料硬化行為有關的硬化函數運算、與屈服函數有關的三個變量：等效應力，一階導，二階導，它們是相互獨立的運算過程，可以將它們分解為獨立的功能模塊——分支子程序，在程序內部實現自調用。結合第2.3節內容，我們提出的UMAT子程序結構模式如圖1所示。等效應力σˉ模塊、一階導模塊、二階導模塊，它們只與選擇的屈服準則有關；屈服應力σY模塊只與加載過程中材料的硬化行為有關；一致切線模量矩陣C的求逆模塊只與矩陣C有直接關系，這5個功能模塊之間相互獨立，運算過程沒有直接關系。針對不同的屈服準則，只需改變與屈服準則有關的模塊（等效應力σˉ模塊、一階導模塊、二模塊）就可以嵌入新的屈服準則；針對不同的材料硬化行為，只需改變與硬化法則有關的模塊（屈服應力σY模塊）就可以研究新的材料硬化行為。因而，上述UMAT子程序結構模式具有互換性，是通用的，它可以將具有新的屈服準則、新的硬化函數的本構模型嵌入ABAQUS.

圖1 通用的UMAT子程序結構模式

表1 Yld2004屈服準則本構模型的UMAT子程序材料常數定義

UMAT子程序的編程需要定義的材料參數有：彈性模量E，泊松比v，屈服準則各向異性系數（Yld2004屈服準則有18個：ci=1～18），硬化模型參數（本文所研究的材料硬化符合Hollomon硬化法則，階導表達式為：σY=K（ε0+ˉp）n，包含 K，ε0，n 三個參數）。Yld2004屈服準則本構模型的UMAT子程序材料常數被定義在表1，同時，UMAT子程序的編程需要定義一個狀態變量數組，用于存放與求解過程有關的的變量，這些變量將在求解過程中不斷地被更新，并傳遞到主程序，然后在下一個增量步開始時再傳遞給UMAT子程序。本文開發的UMAT子程序包含13維的狀態變量數組：彈性應變分量存儲在1～6維，塑性應變分量存儲在7～12維，等效塑性應變存儲在第13維。其他局部變量的名稱定義應該簡潔明了、易辨識。UMAT子程序中數組與數組、矩陣與數組、矩陣與矩陣的矢量乘積編程規則應該特別注意。

UMAT子程序在分析計算過程中會與ABAQUS主程序結合，進行數據的傳輸、更新。為了實現這些過程，UMAT子程序與ABAQUS主程序共享一些變量，即ABAQUS用戶子程序接口，接口定義在UMAT子程序的題名中，對數據的傳輸、更新和相互調用起著橋梁紐帶作用。以下就是UMAT子程序與主程序的接口共享變量：

3 實例有限元分析

為了驗證UMAT子程序開發過程的正確性，必須把UMAT子程序嵌入ABAQUS軟件中進行實例模型有限元分析。根據第1節內容，當Yld2004屈服準則的各向異性系數都為1且指數m=2（或4）時，它可以退化為Mises各向同性屈服準則。通常建立簡單的單拉或者單壓模型，將Yld2004屈服準則的UMAT子程序與單拉或者單壓模型耦合，設置系數全為1、指數m為2或4，將UMAT子程序有限元分析得到的結果，與用ABAQUS自帶的Mises屈服準則計算結果進行對比，可以直接簡潔地驗證Yld2004屈服準則的UMAT子程序是否滿足精度要求，是否在誤差允許范圍內。

本研究建立了簡化的單向壓縮有限元模型，如圖2所示，有限元模型一端完全固定，另一端施加合適的加載力（均布載荷為320MPa），分析步類型為靜力、通用，設置UMAT子程序中各向異性系數都為1即表1材料常數6～23都為1，且指數m=2（或4）即表1材料常數24為2（或4），進行有限元分析，完成有限元分析之后，比較各自的云圖，結果如圖3所示；選取主要的變形路徑提取等效應力、等效塑性應變進行對比，結果如圖4所示。

圖2 簡化的單向壓縮有限元模型

由圖3所示，將Yld2004屈服準則退化為Mises屈服準則模擬的結果與ABAQUS自帶的Mises屈服準則模擬的結果，兩者等效應力分布云圖、等效塑性應變分布云圖整體上都保持一致。

由圖4可知，將Yld2004屈服準則退化為Mises屈服準則模擬的結果與ABAQUS自帶的Mises屈服準則模擬的結果，兩者單元結點上等效應力、等效塑性應變值基本相同，精度保持在99.9%左右。

4 結束語

由于ABAQUS軟件沒有高級的Yld2004屈服準則本構模型，使之受限于精密塑性成形和材料加工有限元分析。因此，本研究推導了與Yld2004-18p屈服準則有關的重要變量，結合完全隱式的Euler圖形返回算法，構建了Yld2004-18p屈服準則本構模型，用Fortran語言編程，通過ABAQUS提供的UMAT子程序二次開發接口，將Yld2004-18p屈服準則本構UMAT子程序嵌入ABAQUS軟件。建立單壓有限元模型，將Yld2004-18p屈服準則UMAT子程序退化到Mises屈服準則，用于實例有限元分析；將分析結果與ABAQUS自帶的Mises屈服準則分析結果作比較，兩者保持高度的一致性，驗證Yld2004-18p屈服準則UMAT子程序二次開發的正確性，同時提出了一個通用的、柔性的屈服準則二次開發結構模式。

圖3 有限元模擬變形圖對比

圖4 變形路徑上節點變量對比