李云飛，曾祥國

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

0 引 言

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

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

(1)

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

(2)

對于fcc金屬材料：

(3)

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

(4)

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

(5)

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

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

1 ABAQUS簡介與子程序算法

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

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

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

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

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

(6)

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

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

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

(7)

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

(8)

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

(9)

(10)

(11)

(12)

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

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

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

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

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

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

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

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

2 結(jié)果與分析

圖3 VUMAT子程序計(jì)算流程圖Fig.3 Flowchart of VUMAT subroutine

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

圖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

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

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

以上為對單個單元本構(gòu)行為的驗(yàn)證，對于結(jié)構(gòu)的準(zhǔn)確度以及計(jì)算效率的驗(yàn)證，本研究選擇平面應(yīng)變懸臂梁的實(shí)例，調(diào)用VUMAT子程序進(jìn)行數(shù)值模擬，然后與采用ABAQUS自帶的Johnson-Cook模型模擬結(jié)果進(jìn)行對比。懸臂梁長度為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自帶模型的計(jì)算結(jié)果吻合較好。然而，子程序的計(jì)算時間為21 s，自帶Johnson-Cook模型的計(jì)算時間為230 s，可見VUMAT子程序的計(jì)算效率明顯優(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子程序進(jìn)行了數(shù)值程序?qū)崿F(xiàn)，解決了商業(yè)軟件中自帶模型無法描述材料在高速切削等條件下力學(xué)性能的問題。與實(shí)驗(yàn)數(shù)據(jù)對比發(fā)現(xiàn)，子程序能夠較好地描述TC4鈦合金或其它金屬的應(yīng)變率強(qiáng)化與溫度軟化效應(yīng)。同時通過與ABAQUS軟件自帶的Johnson-Cook本構(gòu)模型數(shù)值模擬結(jié)果對比發(fā)現(xiàn)兩者結(jié)果吻合良好，子程序在計(jì)算效率方面與前者相比有較大提高。該子程序研究可為鈦合金等金屬高速碰撞沖擊、切削或爆炸變形等情況的數(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] 楊揚(yáng), 曾毅, 汪冰峰. 基于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ìn)展[J]. 機(jī)械工程材料, 2012, 36(3): 5-10.

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

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

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

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