鈦合金熱粘塑性本構關系的子程序開發及有限元驗證

李云飛，曾祥國

(1.中國工程物理研究院總體工程研究所，四川 綿陽 621999)(2.四川大學，四川 成都 610065)

0 引 言

鈦合金具有優良的力學性能以及耐腐蝕性能，在航空航天、海洋工程、石油化工等領域得到廣泛應用。近年來，隨著各工業領域對金屬材料在極端服役條件下的性能要求愈加嚴苛，因此需要對鈦合金等材料在高應變率、高溫等載荷下的特性進行探究。

對于鈦合金在高速沖擊、碰撞等極端載荷下的力學行為，研究者提出了諸多熱粘塑性本構模型，這些模型主要可分為2類。一類是通過大量離散實驗數據擬合得到的唯象經驗型本構模型，其中Johnson-Cook唯象本構方程的應用最為廣泛，其流變應力-塑性應變關系式為[1-2]：

(1)

此外，針對Johnson-Cook模型在高應變率下無法準確描述材料流變應力下降趨勢的缺陷[3]，一些學者提出了改良的物理本構模型。Zerilli與Armstrong[4]基于晶體材料塑性變形過程中位錯運動的熱激活機制，提出了可描述晶體結構分別為體心立方(bcc)與面心立方(fcc)的金屬材料力學行為的本構模型，對于bcc金屬材料：

(2)

對于fcc金屬材料：

(3)

Gao 等人[5]基于金屬塑性變形細觀位錯理論建立了bcc金屬材料的塑性本構方程：

(4)

綜上所述，實際上任何基于物理的或復雜的熱粘塑性唯象經驗本構模型，在同時考慮應變硬化、應變率強化以及溫度軟化效應時，其流變應力表達式都可寫成以下形式[6]：

(5)

近年來隨著計算機模擬技術的快速發展，利用有限元軟件進行分析計算已成為現代科學研究中不可或缺的部分。相比ANSYS、MSC與AIDINA等通用有限元軟件，ABAQUS對于計算不同材料、復雜載荷以及變化接觸條件的非線性組合問題，尤其是非線性力學的分析求解功能處于世界領先水平。目前，雖然ABAQUS等通用有限元軟件的材料庫中自帶了多種本構模型，但都無法準確描述材料在高速切削下的行為。而ABAQUS軟件具備強大的自我擴展開發能力，為各專業領域的用戶提供了若干子程序接口，允許將用戶自定義的材料本構模型導入到軟件的主程序中。

本研究則引入了一種顯式積分算法，介紹了基于VUMAT子程序接口的將通用形式如式(5)的金屬熱粘塑性本構模型進行數值實現的具體方法，并基于Johnson-Cook模型編寫VUMAT子程序對TC4鈦合金的熱粘塑性力學行為進行有限元模擬，對子程序的準確性與計算效率進行驗證分析。

1 ABAQUS簡介與子程序算法

1.1 ABAQUS動態分析技術

ABAQUS于1978年由有限元分析軟件公司ABAQUS推出，被譽為國際上功能最強大的有限元分析軟件之一，不僅可以進行靜態分析，還可以準確分析碰撞、沖擊、爆炸與斷裂等瞬態問題，特別是在非線性分析領域可以解決復雜的工程力學問題。在結構、傳熱學、流體以及流固耦合、熱固耦合等方面具有龐大求解規模的能力，在機械、土木、船舶、汽車、航空航天等各工程領域均發揮了巨大作用。

ABAQUS由ABAQUS/Standard、ABAQUS/Explicit和ABAQUS/CFD 3個主要分析模塊組成，其產品模塊如圖1所示。在這3個模塊中，Explicit可進行顯式動態分析，它使用的顯式求解算法特別適用于求解復雜非線性動力學問題。對于受沖擊載荷并隨后在結構內部發生復雜相互作用的結構瞬態響應問題可以很好的模擬。因此本研究選用ABAQUS/Explicit模塊中的VUMAT接口進行子程序的開發與應用。

圖1 ABAQUS產品模塊示意圖Fig.1 Schematic diagram of ABAQUS product modules

1.2 應力更新算法及VUMAT子程序

(6)

由式(6)可知該算法沒有迭代積分，計算中只需要恒定的彈性張量De，可以顯著減小計算量，提高計算效率。

圖2 應力補償更新算法示意圖Fig.2 Schematic diagram of stress compensation updating algorithm

顯式積分算法基本控制方程及主要步驟如下，彈性變形過程，將應力寫成應變率形式，則廣義胡克定律為：

(7)

對于各向同性硬化的塑性流動采用Mises屈服準則，屈服方程為：

(8)

在某一增量步開始時，假定所有應變增量△ε均為彈性應變，則試應力σtrial寫為：

(9)

(10)

(11)

(12)

為了滿足實際工程應用，ABAQUS為用戶提供了若干用戶子程序(User Subroutines)接口，與命令行形式的程序格式相比，用戶在子程序的限制少得多，功能強大，更加靈活方便。用戶可以利用用戶子程序UMAT與VUMAT接口自行定義材料的本構模型和有限元算法，其中，顯式用戶程序VUMAT適用于ABAQUS/Explicit中。VUMAT主要由以下幾部分組成：子程序初始變量定義、調用ABAQUS外部材料參數、應力應變等參數更新主體程序、結束語句。

考慮到ABAQUS軟件的材料庫中自帶Johnson-Cook模型，為了便于對比驗證子程序的可行性與準確性，本研究以采用Johnson-Cook模型為例，編制其相應的VUMAT用戶子程序。根據上述的應力補償更新算法，定義了3個需要不斷更新的狀態變量，即等效塑性應變、等效塑性應變率與流變應力。實現Johnson-Cook本構方程的子程序主要計算步驟如下：

(3)調用子程序，計算初始流變應力σf；

(4)將試探應力代入屈服判斷準則，判斷是否發生屈服；

(6)若屈服，計算本增量步的塑性應變增量△εp，利用應力補償更新算法更新本增量步結束時的應力；

(7)更新內能、消耗的無彈性能、各狀態變量的值；

(8)結束，返回主程序。

VUMAT子程序的計算流程如圖3所示，子程序在每一個材料積分點被ABAQUS主程序調用[9-10]。

2 結果與分析

圖3 VUMAT子程序計算流程圖Fig.3 Flowchart of VUMAT subroutine

根據上述顯式積分算法以及簡化的應力補償更新算法，按照ABAQUS用戶子程序接口規范，基于FORTRAN語言編寫Johnson-Cook本構模型的VUMAT子程序。數值模擬中的材料選用TC4鈦合金，其本構參數如表1所示。考慮到材料行為與結構形態無關，在ABAQUS中采用單個八節點六面體等參單元(C3D8R)對TC4鈦合金在單軸加載條件下的應力應變響應進行了數值模擬。測試單元的邊界條件及載荷如圖4所示。

表1 Johnson-Cook動態本構參數優化值

圖4 測試單元的邊界條件與載荷Fig.4 Boundary and loading condition of testing unit

在ABAQUS/Explicit中調用VUMAT用戶子程序，分別得到TC4鈦合金在不同應變率(1、10、100、1 000 s-1)及不同溫度條件(293、323、353、403、503 K)下的應力-應變響應曲線如圖5、圖6所示。由圖可知，TC4鈦合金的塑性流變應力隨應變率的增大而逐漸增大，隨初始溫度的升高而逐漸降低，可見VUMAT子程序能夠較好地描述鈦合金或其他金屬的應變率強化與溫度軟化效應，驗證了子程序顯式積分算法的可行性。

圖5 不同應變率下TC4鈦合金數值模擬應力-應變曲線Fig.5 Simulation results of TC4 titanium alloy under different strain rates

圖6 不同初始溫度下TC4鈦合金數值模擬應力-應變曲線Fig.6 Simulation results of TC4 titanium alloy under different initial temperatures

為了驗證VUMAT子程序對TC4鈦合金本構行為預測的準確性，在ABAQUS/Explicit中調用子程序，得到不同應變率、不同初始溫度條件下的應力-應變曲線，其與文獻[11]的實驗數據對比結果如圖7所示。由圖可知，子程序數值模擬結果與實驗數據吻合良好，可見該子程序可以較好地描述TC4鈦合金在不同載荷條件下的本構行為，該顯式積分算法可以推廣至其他用戶定義的金屬熱粘塑性本構模型的數值應用中。

圖7 不同應變率下TC4鈦合金數值模擬結果 與實驗結果的對比Fig.7 Comparison between simulation results and experiment results of TC4 titanium alloy under different strain rates

以上為對單個單元本構行為的驗證，對于結構的準確度以及計算效率的驗證，本研究選擇平面應變懸臂梁的實例，調用VUMAT子程序進行數值模擬，然后與采用ABAQUS自帶的Johnson-Cook模型模擬結果進行對比。懸臂梁長度為1 m，高度為0.2 m，左端固支約束，右端載荷為一垂直向下的高速沖擊位移載荷0.1 m，載荷時間為0.01 s。懸臂梁材料選用TC4鈦合金。子程序與軟件自帶Johnson-Cook模型的數值模擬Mises等效應力結果對比如圖8所示。兩者最大等效應力的誤差約為2.4%，可見子程序與ABAQUS自帶模型的計算結果吻合較好。然而，子程序的計算時間為21 s，自帶Johnson-Cook模型的計算時間為230 s，可見VUMAT子程序的計算效率明顯優于自帶模型。

圖8 子程序與自帶模型等效應力數值模擬結果對比Fig.8 Comparison of Mises stress calculated with subroutine(a) and ABAQUS’s own model(b)

3 結 論

本研究基于顯式應力積分算法將用戶自定義的金屬熱粘塑性本構模型通過VUMAT子程序進行了數值程序實現，解決了商業軟件中自帶模型無法描述材料在高速切削等條件下力學性能的問題。與實驗數據對比發現，子程序能夠較好地描述TC4鈦合金或其它金屬的應變率強化與溫度軟化效應。同時通過與ABAQUS軟件自帶的Johnson-Cook本構模型數值模擬結果對比發現兩者結果吻合良好，子程序在計算效率方面與前者相比有較大提高。該子程序研究可為鈦合金等金屬高速碰撞沖擊、切削或爆炸變形等情況的數值模擬提供技術支撐。

[1] Johnson G R, Cook W H. Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures[J]. Engineering Fracture Mechanics, 1985, 21(1): 31-48.

[2] 楊揚, 曾毅, 汪冰峰. 基于Johnson-Cook模型的TC16鈦合金動態本構關系[J].中國有色金屬學報, 2008, 18(3): 505-510.

[3] Arsecularatne J A, Zhang L C. Assessment of constitutive equations used in machining[J]. Key Engineering Materials, 2004, 274-276:277-282.

[4] Zeriili F J, Armstrong R W. Dislocation-mechanics-based constitutive relations for material dynamics calculations[J]. Journal of Applied Physics, 1987, 61(5):1816-1825.

[5] Gao C Y, Zhang L C, Yan H X. A new constitutive model for HCP metals[J]. Materials Science and Engineering A, 2011, 528(13/14): 4445-4452.

[6] 彭鴻博, 張宏建. 金屬材料本構模型的研究進展[J]. 機械工程材料, 2012, 36(3): 5-10.

[7] 李云飛, 曾祥國, 盛鷹,等. 基于實驗的鈦合金優化動態本構模型與有限元模擬[J]. 材料導報B:研究篇, 2016, 30(12): 137-142.

[8] 李宏偉.宏細觀本構關系數值化及其在有限元模擬中的應用[D].西安:西北工業大學, 2007.

[9] 蘭彬, 葉獻輝, 宋順成, 等. 304NG不銹鋼高應變率材料模型在ABAQUS 中的實現技術[J]. 應用數學和力學, 2015, 36(2): 167-177.

[10] 劉洋, 王秀梅, 王幸敏，等.基于VUMAT的鈦合金切削仿真[J]. 工業控制計算機, 2016, 29(8): 19-21.

[11] Khan A S, Suh Y S, Kazmi R. Quasi-static and dynamic loading responses and constitutive modeling of titanium alloys[J]. International Journal of Plasticity, 2004, 20(12): 2233-2248.