王效貴,李曉葉

(浙江工業(yè)大學 機械工程學院，浙江 杭州 310014)

GTN損傷模型的算法研究及試驗驗證

王效貴,李曉葉

(浙江工業(yè)大學 機械工程學院，浙江 杭州 310014)

摘要：基于GTN細觀損傷本構(gòu)理論，引入了塑性極限載荷模型，將完全隱式積分算法與隱式有限元計算相結(jié)合，建立了相應的數(shù)值積分算法.通過編寫用戶自定義材料子程序(UMAT)，實現(xiàn)了修正的GTN模型在有限元軟件ABAQUS中的應用.通過缺口圓棒試樣損傷與斷裂的數(shù)值模擬，對修正的GTN模型進行了驗證，并分析了缺口圓棒試樣的裂紋萌生點位置，裂紋擴展情況及應力三軸度的變化.數(shù)值預測與試驗結(jié)果吻合較好，表明修正的GTN模型可以有效的描述材料的損傷和斷裂現(xiàn)象.

關鍵詞：GTN模型；數(shù)值算法；損傷斷裂；有限元法

Algorithm study and experiment validation of GTN damage model

WANG Xiaogui, LI Xiaoye

(College of Mechanical Engineering, Zhejiang University of Technology, Hangzhou 310014, China)

Abstract:The corresponding numerical integral algorithm is established by adding the plastic limit load model into the GTN mesoscopic damage constitutive model and combining the fully implicit integration algorithm and the implicit finite element. It is implemented in the finite element software ABAQUS by using the user material subroutine (UMAT). The modified GTN model is verified by simulating the damage and fracture of notched round bar. The crack initiation site and propagation path, and the change of stress triaxiality were analyzed. The numerical predications are in good agreement with the experimental observations. It also shows that the modified GTN mesoscopic damage model can effectively reveal the damage behavior and predict the fracture of the material.

Key words:GTN model; numerical algorithm; damage and failure; finite element method

損傷和斷裂是工程材料的常見失效模式.材料中的夾雜物或二相粒子與基體的脫粘或自身開裂，將會引起孔洞的形核、長大及聚合，最終形成一定尺度的宏觀裂紋，這是導致材料失效的根本原因.近年來，很多學者采用不同方法對損傷斷裂問題展開研究.劉俊新[1]等采用三維重建技術研究了泥頁巖蓋層在載荷作用下的裂紋擴展與演化規(guī)律；蔡增伸[2]等通過力控制法對金剛石膜的斷裂韌性進行了研究.蘇彬[3]等自行設計實驗裝置研究了循環(huán)載荷作用下的微動疲勞特性和裂紋位置.

目前應用最為廣泛的損傷模型是GTN模型，它是一種以孔洞體積分數(shù)為損傷變量反映材料內(nèi)部微小缺陷的損傷模型，基于這一模型國內(nèi)外學者取得了諸多研究成果.Vadillo[4]等通過數(shù)值算法研究了GT模型中的損傷參數(shù)q1和q2，發(fā)現(xiàn)這兩個參數(shù)并不是常數(shù)，而是與應力三軸度和初始孔隙率相關的變量；邢佶慧[5]等通過數(shù)值模擬的方法，分析了初始孔洞體積分數(shù)，孔洞形核體積分數(shù)及材料失效時的孔洞破壞體積分數(shù)對斷裂失效的影響.由于原Gurson模型沒有考慮不連續(xù)孔洞的影響，不能預測孔洞之間的頸縮現(xiàn)象.因此，筆者基于GTN細觀損傷理論，引入塑性極限載荷模型，采用Ramberg-Osgood硬化準則，從數(shù)值積分算法著手，模擬高導無氧銅(OFHC)材料的損傷演化和斷裂失效過程，并與試驗作對比，驗證其可行性.

1GTN模型的算法研究及數(shù)值實現(xiàn)

GTN模型的屈服函數(shù)[6]表示為

1-(q1f*)2

(1)

式中：σeq為宏觀Mises等效應力；σm為宏觀靜水應力；σY為流動應力；q1和q2為損傷參數(shù)；f*為孔洞體積分數(shù)函數(shù)，定義為

(2)

式中：f為孔洞體積分數(shù)；fc為開始聚合時的臨界孔洞體積分數(shù)；ff為材料斷裂時的臨界孔洞體積分數(shù).

由于原Gurson模型不能預測孔洞之間的頸縮現(xiàn)象，因此在GTN模型中引入塑性極限載荷模型來描述這一現(xiàn)象.該模型給出的軸對稱球形孔洞聚合公式[7-8]為

(3)

式中：Rr為球形孔洞的當前半徑；L為單胞的當前長度；σ1為最大主應力.

OFHC在塑性屈服后，表現(xiàn)出輕微的硬化現(xiàn)象，為了反映這個現(xiàn)象，采用Ramberg-Osgood硬化準則，其曲線關系為

(4)

式中：σ0為初始屈服應力；ε0為初始屈服應變；k為硬化系數(shù)；ρ為硬化指數(shù).

根據(jù)Aravas[9]提出的向后Euler完全隱式積分算法，更新應力應變關系.n表示增量步，為方便書寫，省略下標n+1.假設給定的應變增量全為彈性應變增量，計算彈性試應力

σtrial=σn+C：Δε

(5)

式中：C為四階彈性模量；Δε為總的應變增量.

根據(jù)關聯(lián)流動準則，更新n+1增量步時的宏觀塑性應變增量為

(6)

(7)

根據(jù)Hook定律，結(jié)合式(5，6)，得到n+1增量步時的宏觀柯西應力為

σ=σtrial-KΔDmI-2GΔDeqN

(8)

式中：G為剪切模量；K為體積模量.將式(8)分別映射到I和N方向，得到當前時刻的宏觀靜水應力和宏觀等效應力為

(9)

因此，宏觀柯西應力又可表示為

(10)

Δf=(1-f)ΔDm+

(11)

(12)

將式(4，9，11，12)代入式(1)，得到n+1增量步時的屈服函數(shù)為

F(ΔDeq,ΔDm)=0

(13)

(14)

根據(jù)式(7，13，14)推導一致切線剛度矩陣T，得到

(15)

式中：Mijkl的表達式見文獻[9]，Jijkl=(δikδjl+δilδjk)/2.

根據(jù)以上推導的離散公式，下面給出修正的GTN模型在用戶子程序(UMAT)中的實現(xiàn)流程圖，如圖1所示.

圖1　流程圖Fig.1　Flow chart

2單軸拉伸試驗

缺口圓棒試樣的單軸拉伸試驗是研究韌性金屬材料損傷斷裂的最好方法，因此，采用如圖2所示的高導無氧銅光滑和缺口圓棒試樣，進行單軸拉伸試驗.光滑試樣標距段的直徑為10mm.缺口試樣有三個重要幾何參數(shù)，分別為缺口寬度d、缺口半徑r和凈斷面直徑D，實測數(shù)據(jù)如表1所示.采用RG4100微機控制電子萬能試驗機，在室溫狀態(tài)下完成單軸拉伸試驗.試驗采用位移控制，引伸計的標距為20mm.光滑試樣的位移加載速率為0.3mm/min.基于不同缺口尺寸，試樣φ6r1，φ6r2，φ8r1，φ8r2的位移加載速率分別為0.05，0.08，0.10，0.15mm/min.試驗完畢后，測量缺口部位的d和D，結(jié)果列于表1.由表1中的斷面收縮率可以看出：凈斷面直徑越大，斷面收縮率越大；同一種凈斷面直徑，缺口半徑越大，斷裂收縮率越大.

圖2　單軸拉伸試樣Fig.2　Uniaxial tensile specimens

試樣編號實驗前缺口半徑r/mm缺口寬度d/mm凈斷面直徑D/mm試驗后缺口寬度d'/mm凈斷面直徑D'/mm斷面收縮率/%?6r1-11.02.106.053.263.8060.55?6r1-21.02.085.953.091)4.1052.52?6r1-31.02.126.073.061)4.0655.26?6r2-12.03.845.955.793.6063.39?6r2-22.03.865.935.733.6362.53?6r2-32.03.835.925.521)3.8058.80?8r1-11.02.157.924.594.4568.43?8r1-21.02.107.944.564.4768.31?8r1-31.02.137.954.524.5067.76?8r2-12.03.237.935.744.3070.60?6r2-22.03.257.945.774.2771.08?6r2-32.03.207.925.794.3470.00

注：1) 為在單軸拉伸試驗中材料沒有拉斷.

由光滑試樣單軸拉伸試驗得到的真實應力—應變曲線，如圖3所示.由彈性階段的斜率，得到OFHC的彈性模量為E=108 GPa.基于國標GB/T 228.1—2010，確定初始屈服應變?yōu)棣?=0.026 85，且滿足σ0=Eε0，屈服極限為σy=290 MPa，強度極限為σu=314 MPa.OFHC在塑性屈服后，表現(xiàn)出輕微的硬化現(xiàn)象.采用公式(4)給出的Ramberg-Os-good硬化準則，曲線擬合后得到硬化系數(shù)為k=0.75和硬化指數(shù)為ρ=27.

圖3　光滑圓棒試樣的真實應力—應變曲線Fig.3　True stress-strain curves of smooth round bar specimens

由缺口試樣單軸拉伸試驗得到的載荷—位移曲線如圖4所示.對每種缺口試樣進行3次重復試驗，獲得的載荷—位移曲線具有很好的一致性.試驗曲線分成三個階段：第一，試驗開始到載荷最大值，這一期間載荷快速增大，而材料的變形程度非常小；第二，載荷最大值到裂紋萌生點(拐點)，當達到最大載荷時孔洞開始聚合，材料承載能力逐漸減弱，同時缺口半徑處出現(xiàn)頸縮現(xiàn)象，縱向變形快速增大；第三，裂紋萌生點到材料斷裂失效，隨著裂紋的萌生和擴展，損傷演化程度加快，載荷急劇降低，材料承載能力下降，最終導致斷裂破壞.圖4可以看出：缺口試樣的凈斷面直徑越小，材料的斷裂失效出現(xiàn)的越早；同一種凈斷面直徑，缺口半徑越小，斷裂破壞產(chǎn)生的越快；凈斷面直徑越大，載荷—位移曲線中峰值越高.

圖4　缺口試樣的載荷—位移曲線Fig.4　Load-Displacement curves of notched specimens

3算例及討論

3.1有限元模型

建四分之一軸對稱有限元模型，如圖5所示.邊界條件的施加分為3部分：在左端面，施加軸向位移載荷，試樣φ6r1，φ6r2，φ8r1，φ8r2的位移大小分別為0.6，0.75，1.5，1.5 mm；在右端面，施加y向?qū)ΨQ位移約束；在x=0對稱面上，施加x向?qū)ΨQ位移約束.為了能夠很好地模擬缺口根部區(qū)域較大的塑性應力應變梯度，將該區(qū)域的網(wǎng)格細化，最小網(wǎng)格尺寸為0.1 mm.為了提高計算速度，遠離缺口區(qū)域的網(wǎng)格尺寸為0.2 mm.采用嵌入GTN模型的有限元軟件ABAQUS，選用八節(jié)點軸對稱縮減積分單元(CAX8R)，進行數(shù)值模擬斷裂預測.

圖5　缺口圓棒試樣的有限元模型Fig.5　FEA model of notched specimens

3.2GTN模型的參數(shù)確定

由于GTN損傷模型中的體胞是假設的，因此f0的取值需根據(jù)模擬結(jié)果和試驗數(shù)據(jù)對比而確定.以試樣φ6r1-1為研究對象，當給定的f0能夠使數(shù)值模擬結(jié)果和試驗數(shù)據(jù)基本相吻合時，則認為此時的f0即為要確定的初始孔隙率，最終確定f0=0.001.文獻[4]給出了損傷參數(shù)q1和q2隨應力三軸度的變化圖，依據(jù)給定的關系圖，可以初步確定損傷參數(shù)q1和q2，由于不同材料的實際損傷情況并不一致，因此根據(jù)數(shù)值預測結(jié)果和試驗數(shù)據(jù)進一步修正，最終得到q1=1.2，q2=0.8.根據(jù)本試驗的實際情況，結(jié)合文獻[5]中介紹的方法，最終確定fN=0.004 5，ff=0.4.

3.3結(jié)果與討論

根據(jù)確定的損傷參數(shù)，將修正的GTN模型嵌入ABAQUS軟件進行斷裂預測，分析缺口圓棒試樣裂紋萌生點位置和擴展情況，以及應力三軸度的變化.

3.3.1裂紋萌生點位置及擴展路線的預測

圖4展示了缺口試樣裂紋萌生點和裂紋擴展路線，圖4中的黑點代表裂紋萌生點位置.可以清楚的看到：裂紋萌生最早產(chǎn)生于缺口試樣的中心位置附近，隨后向缺口試樣的尖端方向逐漸發(fā)展，隨著塑性變形的增大，裂紋快速擴展并貫通缺口根部，最終造成材料的斷裂失效.這是因為，由于塑性變形的增大，孔洞體積分數(shù)也隨之增加，當達到塑性極限載荷時，微孔洞開始聚合，孔洞體積分數(shù)逐漸由f變成f*，并急劇增大，當f*增加到一定程度時，缺口試樣的中間部位裂紋開始萌生.隨著裂紋的快速長大和擴展，載荷急劇下降，材料的承載能力降低，當f*增加到材料斷裂的臨界值時損傷失效產(chǎn)生.

從圖4中的試驗和數(shù)值模擬曲線看到：裂紋萌生點位置的數(shù)值預測與試驗結(jié)果較為接近，裂紋擴展路線的預測結(jié)果與試驗曲線趨勢基本一致.表明修正的GTN模型能夠模擬缺口圓棒試樣的斷裂失效特征，驗證了該模型的有效性，該方法的合理性和可行性.

3.3.2應力三軸度的變化

應力三軸度定義為平均應力與等效應力的比值，它是影響韌性材料斷裂的一個重要因素，通常根據(jù)試樣最小橫截面來評估.圖6給出最大載荷條件下應力三軸度的變化過程，圖7給出裂紋萌生點附近應力三軸度的變化過程.

以上四幅圖顯示，凈斷面直徑相同，缺口半徑越小，最小橫截面上的應力三軸度越大；應力三軸度通常在缺口試樣的中間部位最高，在自由邊附近最低.由3.3.1節(jié)知道：裂紋萌生最早發(fā)生于缺口試樣中心部位附近，隨著塑性變形的增大，加劇了缺口試樣中心位置的裂紋萌生擴展，從而導致缺口根部的孔洞體積分數(shù)快速增大.由GTN模型的屈服函數(shù)可以知道，當孔洞體積分數(shù)增大到一定程度時，平均應力遠大于等效應力，因此出現(xiàn)了上圖所述規(guī)律.比較圖6,7看到：同一種缺口結(jié)構(gòu)，圖7中的應力三軸度要比圖6中的偏大.這是因為，當達到塑性極限載荷時微孔洞剛開始聚合，裂紋萌生并沒有完全開始，孔洞體積分數(shù)相對較小，此時等效應力所起的作用比平均應力大，應力三軸度較小.而在裂紋萌生點附近，由于塑性變形的增大，加劇了裂紋萌生和快速擴展，導致此時的孔洞體積分數(shù)比最大載荷工況下更大，平均應力的影響遠遠大于等效應力，所以應力三軸度較大.同時還可以發(fā)現(xiàn)，同一種缺口半徑，凈斷面直徑越小，應力三軸度越大.這是因為凈斷面直徑越小，此處的應力集中對損傷斷裂的影響更加敏感而造成.

圖6　最大載荷下應力三軸度沿最小截面的變化Fig.6　Variations in stress triaxiality along minimum section corresponding to maximum load

圖7　裂紋萌生下應力三軸度沿最小截面的變化Fig.7　Variations in stress triaxiality along minimum section corresponding to crack initation

圖4中的試驗曲線顯示：試樣從早到晚的斷裂情況依次是φ6r1，φ6r2，φ8r1，φ8r2.觀察圖6，7發(fā)現(xiàn)：應力三軸度由大到小的順序是φ6r1，φ6r2，φ8r1，φ8r2.而應力三軸度越大，孔洞體積分數(shù)的增加就越急劇，裂紋萌生和擴展速度越快，斷裂破壞產(chǎn)生的就越早.因此GTN模型的模擬結(jié)果和試驗結(jié)果吻合良好.

4結(jié)論

采用嵌入修正的GTN模型的有限元軟件ABAQUS進行斷裂模擬，得出以下幾點結(jié)論：1) 通過對缺口圓棒試樣裂紋萌生點位置和擴展路徑的分析發(fā)現(xiàn)，裂紋萌生首先產(chǎn)生于缺口試樣的中部，而不是缺口試樣的尖端.一旦裂紋萌生，載荷急劇下降，材料的承載能力快速降低.2) 應力三軸度的分析表明，缺口試樣的中間位置應力三軸度最大，自由邊附近最低.隨著塑性變形的增加，孔洞體積分數(shù)急劇增大，加劇了裂紋的萌生和長大，從而導致缺口試樣中間部位的塑性應變最大，應力三軸度也最大，因此斷裂失效最早發(fā)生.3) 缺口試樣的模擬結(jié)果與試驗結(jié)果吻合較好，驗證了修正的GTN模型具有預測材料斷裂行為的可行性.

參考文獻：

[1]劉俊新，楊春和，冒海軍，等．基于CT圖像處理的泥頁巖裂紋擴展與演化研究[J]．浙江工業(yè)大學學報，2015，43(1)：66-72．

[2]蔡增伸，蔣政，王從賢，等．測定金剛石膜斷裂強度及斷裂韌性的力控制法研究[J]．浙江工業(yè)大學學報，2000，28(3)：216-219．

[3]蘇彬，邢海軍，趙穎娣，等．循環(huán)載荷應力比對微動疲勞特性的影響[J]．浙江工業(yè)大學學報，2011，39(4)：445-448．

[4]VADILLO G, FERNANDEZ-SAEZ J. An analysis of Gurson model with parameters dependent on triaxiality based on unitary cell[J]. European Journal of Mechanics. A: Solids，2009，28(3)：417-427．

[5]邢佶慧，史一劍，劉文濤，等．建筑鋼材GTN損傷模型參數(shù)識別[J]．建筑結(jié)構(gòu)學報，2014，35(4)：149-154．

[6]TVERGAARD V, NEEDLEMAN A. Analysis of the cup-cone fracture in a round tensile bar[J]．Acta Metallurgica，1984，32(1)：157-169．

[7]THOMASON P F. Ductile fracture of metals[M]．Oxford：Pergamon Press，1990．

[8]PARDOEN T, HUTCHINSON J W. An extended model for void growth and coalescence[J]．Journal of Mechanics and Physics of Solids，2000，48(12)：2467-2512．

[9]ARAVAS N. On the numerical integration of a class of pressure-dependent plasticity models[J]．International Journal for numerical methods in engineering，1987，24(7)：1397-1416．

[10]CHU C C, NEEDLEMAN A. Void nucleation effects in biaxially stretched sheets[J]．Journal of Engineering Materials and Technology，1980，102(3)：249-256．

(責任編輯：陳石平)

文章編號：1006-4303(2015)06-0660-06

中圖分類號：TG146.1

文獻標志碼：A

作者簡介：王效貴(1973—)，男，山東青州人，教授，研究方向為結(jié)構(gòu)完整性，E-mail：hpcwxg@zjut.edu.cn.

基金項目：國家自然科學基金資助項目(51175469)

收稿日期：2015-05-04