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

喬順成，吳建軍，展學鵬

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

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

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

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

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

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

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

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

式中：指數(shù)m與材料的晶體結(jié)構(gòu)類型有關(guān)，對于BCC（體心立方）材料，m=6，對于FCC（面心立方）材料；S?′，S?″是與應(yīng)力偏張量s?′，s?″有關(guān)的主值；s?′，s?″由應(yīng)力偏張量s通過線性變換得到，如下式：

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

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

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

2 Yld2004屈服準則本構(gòu)模型的UMAT子程序二次開發(fā)

2.1 與Yld2004屈服準則有關(guān)的變量推導(dǎo)

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

對于Yld2004屈服準則本構(gòu)模型，結(jié)合等式（1），可以得到等效應(yīng)力σˉ的表達式：

結(jié)合等式（1）、（7），等效應(yīng)力σˉ對應(yīng)力分量 σij

其中：

等式（9）是Φ對σij的“鏈式求導(dǎo)”，可結(jié)合等式（1）-（6）求得。等效應(yīng)力σˉ對應(yīng)力分量 σij的二階導(dǎo)為：

其中：

以上這些公式將會在下文應(yīng)用于UMAT子程序本構(gòu)模型。

2.2 彈-塑性率本構(gòu)方程

材料屈服狀態(tài)的判斷根據(jù)是復(fù)雜應(yīng)力狀態(tài)下的等效應(yīng)力與參考方向的單向屈服應(yīng)力之間的關(guān)系，如此可以建立屈服條件：

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

式中：

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

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

令：

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

2.3 UMAT子程序基本流程的構(gòu)建

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

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

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

（1）初始狀態(tài)：

初始值設(shè)置：k=0，Δλ=0；

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

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

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

（5）更新變量：

返回第（2）步。

2.4 UMAT子程序結(jié)構(gòu)模式及編程

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

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

圖1 通用的UMAT子程序結(jié)構(gòu)模式

表1 Yld2004屈服準則本構(gòu)模型的UMAT子程序材料常數(shù)定義

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

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

3 實例有限元分析

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

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

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

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

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

4 結(jié)束語

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

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

圖4 變形路徑上節(jié)點變量對比