李云飛，曾祥國

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

0 引 言

鈦合金具有優(yōu)良的力學(xué)性能以及耐腐蝕性能，在航空航天、海洋工程、石油化工等領(lǐng)域得到廣泛應(yīng)用。近年來，隨著各工業(yè)領(lǐng)域?qū)饘俨牧显跇O端服役條件下的性能要求愈加嚴苛，因此需要對鈦合金等材料在高應(yīng)變率、高溫等載荷下的特性進行探究。

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

(1)

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

(2)

對于fcc金屬材料：

(3)

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

(4)

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

(5)

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

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

1 ABAQUS簡介與子程序算法

1.1 ABAQUS動態(tài)分析技術(shù)

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

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

圖1 ABAQUS產(chǎn)品模塊示意圖Fig.1 Schematic diagram of ABAQUS product modules

1.2 應(yīng)力更新算法及VUMAT子程序

(6)

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

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

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

(7)

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

(8)

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

(9)

(10)

(11)

(12)

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

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

(3)調(diào)用子程序，計算初始流變應(yīng)力σf；

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

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

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

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

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

2 結(jié)果與分析

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

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

表1 Johnson-Cook動態(tài)本構(gòu)參數(shù)優(yōu)化值

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

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

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

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

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

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

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

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

3 結(jié) 論

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

[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鈦合金動態(tài)本構(gòu)關(guān)系[J].中國有色金屬學(xué)報, 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] 彭鴻博, 張宏建. 金屬材料本構(gòu)模型的研究進展[J]. 機械工程材料, 2012, 36(3): 5-10.

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

[8] 李宏偉.宏細觀本構(gòu)關(guān)系數(shù)值化及其在有限元模擬中的應(yīng)用[D].西安:西北工業(yè)大學(xué), 2007.

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

[10] 劉洋, 王秀梅, 王幸敏，等.基于VUMAT的鈦合金切削仿真[J]. 工業(yè)控制計算機, 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.