王永剛，徐東明

(寧波大學(xué) 力學(xué)與材料科學(xué)研究中心，浙江 寧波315211)

0 引言

延性金屬材料在高應(yīng)變率加載下發(fā)生拉伸斷裂主要包括3 個(gè)物理過程:微孔洞的成核、長(zhǎng)大和貫通［1］。迄今，人們對(duì)于微孔洞的成核和長(zhǎng)大規(guī)律進(jìn)行了大量的實(shí)驗(yàn)和理論分析工作［2-4］，但微孔洞相互貫通行為的研究工作相對(duì)很少。已有學(xué)者研究指出微孔洞的貫通在動(dòng)態(tài)拉伸斷裂過程中起著很關(guān)鍵的作用［5］，因此，研究微孔洞貫通機(jī)理對(duì)揭示材料動(dòng)態(tài)拉伸斷裂的整個(gè)物理過程和本質(zhì)有重要意義。

Curran 等［6］基于實(shí)驗(yàn)觀測(cè)結(jié)果認(rèn)為有兩類微孔洞貫通機(jī)理，即微孔洞之間直接貫穿和微孔洞之間的變形局部化。事實(shí)上，微孔洞相互作用是通過它們周圍的塑性應(yīng)變場(chǎng)來實(shí)現(xiàn)的。當(dāng)兩個(gè)微孔洞之間距離達(dá)到一個(gè)臨界值時(shí)，它們周圍的塑性應(yīng)變場(chǎng)發(fā)生交接，從而改變微孔洞的長(zhǎng)大速率。一些學(xué)者曾對(duì)微孔洞之間開始貫通臨界條件進(jìn)行了理論研究。Needleman 和Tvergarrd［7］認(rèn)為材料內(nèi)部微孔洞體積分?jǐn)?shù)f 達(dá)到一個(gè)臨界值fc時(shí)微孔洞貫通開始。Brock和Sun［8］研究指出當(dāng)材料中初始微孔洞體積分?jǐn)?shù)較小時(shí)，臨界微孔洞體積分?jǐn)?shù)fc對(duì)應(yīng)力三軸度和加載應(yīng)變率不敏感。而Pardoen 和Hutchinson［9］認(rèn)為fc依賴于微孔洞初始體積分?jǐn)?shù)、加載應(yīng)力三軸度及微孔洞形狀。Brown 和Embury［10］提出了一個(gè)微孔洞開始貫通的臨界韌帶距離(ILDc)準(zhǔn)則，指出當(dāng)微孔洞間韌帶長(zhǎng)度等于微孔洞的直徑時(shí)，微孔洞從獨(dú)立長(zhǎng)大開始向貫通轉(zhuǎn)變。Horstemeyer 等［11］采用微力學(xué)有限元計(jì)算方法研究準(zhǔn)靜態(tài)加載條件下304L不繡鋼和6061T6 鋁中兩個(gè)微孔洞的貫通行為，結(jié)果表明:ILDc的值依賴于邊界條件和材料的特性，約為2～8 倍微孔洞直徑。最近，Sepp?l? 等［12-13］采用分子動(dòng)力學(xué)方法對(duì)金屬銅動(dòng)態(tài)拉伸斷裂的微孔洞長(zhǎng)大和聚集行為進(jìn)行了數(shù)值模擬，計(jì)算結(jié)果顯示:當(dāng)ILDc等于0.5 倍微孔洞直徑時(shí)，微孔洞周邊塑性形變區(qū)相遇，微孔洞長(zhǎng)大的速率突然增加，最終導(dǎo)致了微孔洞之間發(fā)生了直接貫穿。

無論是采用fc還是ILDc作為微孔洞貫通開始的準(zhǔn)則，人們對(duì)微孔洞貫通行為的認(rèn)識(shí)還是比較初步的，存在一些爭(zhēng)議。本文擬采用動(dòng)力學(xué)有限元計(jì)算模型，對(duì)高應(yīng)變率加載條件下純鋁材料中雙孔洞貫通行為進(jìn)行數(shù)值分析，重點(diǎn)討論初始孔洞間韌帶距離和加載應(yīng)變速率對(duì)雙孔洞貫通行為的影響規(guī)律，為合理建立孔洞貫通的理論模型奠定一定的基礎(chǔ)，從而進(jìn)一步揭示延性金屬拉伸型斷裂的物理機(jī)制。

1 計(jì)算模型建立

高應(yīng)變率拉伸載荷下延性金屬中出現(xiàn)的損傷主要是大量的隨機(jī)分布的球形微孔洞。因此，在一定假設(shè)條件下，可以把含損傷材料離散成大量含孔洞的圓柱體胞的集合體。這里，為了研究孔洞之間的貫通行為，在數(shù)值模擬中設(shè)定每一個(gè)圓柱體胞含有兩個(gè)大小相同的球形孔洞，二分之一模型見圖1(a)所示。在對(duì)稱載荷條件下，我們采用軸對(duì)稱的二維單元來進(jìn)行計(jì)算分析，有限元模型見圖1(b)所示。微孔洞的直徑d0為50 μm，圓柱體胞的直徑D 和高度H 相等，均取為1 mm.圓柱體胞的直徑和高度是球形微孔洞直徑的20 倍。

圖1 二分之一體胞模型及軸對(duì)稱有限元分析模型Fig.1 Half of the axisymmetric unit cell and FEM model

計(jì)算中，邊界條件采用速度邊界，如圖1所示，vz、vr分別表示在軸向和徑向施加的速度函數(shù)。速度從0 開始經(jīng)歷0.1 μs 線性增加到恒值，隨后保持不變。選取的模擬材料為純鋁，在高應(yīng)變率加載下，采用考慮塑性硬化的流體彈塑性本構(gòu)關(guān)系和Grüneisen 狀態(tài)方程進(jìn)行數(shù)值分析。

式中:σij為應(yīng)力張量;p 為壓力;sij為偏應(yīng)力張量。

式中:σy為屈服應(yīng)力考慮塑性硬化時(shí)，

式中:σ0為初始屈服應(yīng)力;為等效塑性應(yīng)變;Eh塑性硬化模量。

在材料處于壓縮狀態(tài)下，Grüneisen 狀態(tài)方程的壓力計(jì)算:

當(dāng)材料處于拉伸狀態(tài)下，

計(jì)算參數(shù)見表1.

表1 純鋁的本構(gòu)和狀態(tài)方程參數(shù)Tab.1 Constitutive and equation of state parameters for pure aluminium

2 計(jì)算結(jié)果討論

在軸向速度場(chǎng)和徑向速度場(chǎng)相同條件下，圖2顯示軸對(duì)稱體胞模型邊界單元上宏觀軸向應(yīng)力和徑向應(yīng)力幾乎完全相等，應(yīng)力開始時(shí)基本呈線性增加，隨后由于微孔洞長(zhǎng)大的影響，出現(xiàn)非線性增長(zhǎng)，很快達(dá)到極值點(diǎn)，隨后快速下降。由此可見，體胞模型中微孔洞承受的是一個(gè)準(zhǔn)靜水拉伸載荷條件，這與一維應(yīng)變層裂實(shí)驗(yàn)中樣品中微孔洞所承受了應(yīng)力狀態(tài)基本一致。因此，后續(xù)的有關(guān)微孔洞之間貫通行為的研究結(jié)果可以直接應(yīng)用于層裂實(shí)驗(yàn)的理論分析和數(shù)值模擬中。另外，在準(zhǔn)靜水拉伸載荷下，圖3顯示不同時(shí)刻兩個(gè)微孔洞周圍的塑性應(yīng)變場(chǎng)的演化云圖。圖3(a)顯示兩個(gè)微孔洞各自形成幾乎均勻的球?qū)ΨQ塑性應(yīng)變場(chǎng)，微孔洞呈現(xiàn)相互獨(dú)立球形快速長(zhǎng)大，而隨著微孔洞不斷長(zhǎng)大，兩個(gè)微孔洞之間韌帶區(qū)域塑性應(yīng)變場(chǎng)相互交疊，如圖3(b)所示，當(dāng)此區(qū)域的應(yīng)變場(chǎng)達(dá)到一個(gè)臨界值時(shí)，兩個(gè)微孔洞即發(fā)生貫通行為，直至直接貫穿，如圖3(c)和圖3(d)所示。

圖3給出了微孔洞貫通行為的一些定性描述，為了定量討論微孔洞之間的貫通行為，這里采用模擬計(jì)算中實(shí)時(shí)跟蹤微孔洞長(zhǎng)大過程中形狀改變來判斷，即通過比較圖1中標(biāo)識(shí)的微孔洞上頂點(diǎn)A 和下頂點(diǎn)B 的軸向位移。如果微孔洞處于獨(dú)立長(zhǎng)大階段，A 點(diǎn)和B 點(diǎn)軸向位移是一致的，如果A 和B 點(diǎn)軸向位移發(fā)生明顯的改變，表明微孔洞之間韌帶區(qū)域塑性應(yīng)變場(chǎng)交疊改變了B 點(diǎn)的軸向位移增長(zhǎng)速率，此時(shí)刻即可判斷微孔洞貫通開始，所對(duì)應(yīng)的圖1中標(biāo)識(shí)的B 點(diǎn)和C 點(diǎn)之間距離即為微孔洞開始貫通的臨界韌帶距離(ILDc).下面首先討論初始孔洞間韌帶距離ILD0對(duì)微孔洞貫通動(dòng)力學(xué)行為的影響。

圖2 宏觀軸向應(yīng)力和徑向應(yīng)力的比較Fig.2 Comparison of global axial stress and radial stress

圖3 不同時(shí)刻微孔洞周圍等效塑性應(yīng)變?cè)茍DFig.3 Effective plastic strain contours around void

2.1 ILD0的影響

在加載應(yīng)變率104s-1條件下，考察不同初始微孔洞間韌帶距離對(duì)孔洞貫通動(dòng)力學(xué)行為的影響，選擇ILD0/d0為2、3、4 和5，即ILD0分別為100 μm、150 μm、200 μm、250 μm(上面圖2和圖3的數(shù)據(jù)對(duì)應(yīng)的ILD0/d0=3).圖4(a)～圖4(d)左邊縱軸給出了不同ILD0/d0條件下實(shí)時(shí)的微孔洞間韌帶距離ILD 與實(shí)時(shí)微孔洞直徑d 之比隨時(shí)間演化過程，而其右邊縱軸同時(shí)給出了A 點(diǎn)和B 點(diǎn)的軸向位移時(shí)程。這里需要說明一下，為什么沒有直接來討論臨界ILD，而是討論它與實(shí)時(shí)微孔洞直徑d 之比，這主要是因?yàn)槲⒖锥簇炌ㄊ峭ㄟ^微孔洞周圍的塑性應(yīng)變場(chǎng)來實(shí)現(xiàn)的，而微孔洞的直徑直接影響著其周圍的塑性應(yīng)變場(chǎng)。對(duì)比A 點(diǎn)和B 點(diǎn)的軸向位移時(shí)程，發(fā)現(xiàn):微孔洞貫通開始前A 點(diǎn)和B 點(diǎn)軸向位移是相同的，兩個(gè)微孔洞之間沒有發(fā)生相互作用，微孔洞處于獨(dú)立長(zhǎng)大階段，隨著微孔洞不斷長(zhǎng)大，它們周圍的塑性應(yīng)變場(chǎng)發(fā)生相互影響，導(dǎo)致隨后A 點(diǎn)和B 點(diǎn)軸向位移明顯發(fā)生偏離，利用此點(diǎn)恰好可以判斷微孔洞開始貫通時(shí)刻，因此，這里的計(jì)算結(jié)果也檢驗(yàn)了上面提出的通過微孔洞形狀改變來判斷微孔洞貫通開始的有效性。而從ILD/d 的時(shí)程曲線上來看，在0.3 μs 之前，ILD/d 減小速率緩慢，此階段微孔洞基本還處于彈性變形階段，隨后進(jìn)入微孔洞獨(dú)立快速長(zhǎng)大階段，這時(shí)ILD/d 迅速減小，緊接著微孔洞間發(fā)生貫通行為，導(dǎo)致ILD/d 減小速率再次變緩，最后慢慢趨于0，即微孔洞發(fā)生了直接貫穿，如上面圖3(d)所示。對(duì)比圖4(a)～圖4(d)，由A 點(diǎn)和B 點(diǎn)軸向位移明顯發(fā)生偏離點(diǎn)容易確定微孔洞間貫通開始的臨界韌帶距離，在圖中用黑點(diǎn)標(biāo)識(shí)，結(jié)果顯示:微孔洞從獨(dú)立長(zhǎng)大向微孔洞貫通轉(zhuǎn)變開始的臨界韌帶距離約為0.5 倍(ILDc=0.5 ±0.1d)的微孔洞直徑，并且它對(duì)初始的ILD0/d0不敏感。

臨界韌帶距離作為微孔洞貫通開始的準(zhǔn)則，比較直觀，但相對(duì)缺乏物理本質(zhì)的認(rèn)識(shí)和提煉。分子動(dòng)力學(xué)研究成果［12-13］已經(jīng)在微觀層次上揭示微孔洞間韌帶區(qū)域位錯(cuò)密度迅速增長(zhǎng)導(dǎo)致微孔洞發(fā)生貫通。位錯(cuò)運(yùn)動(dòng)宏觀上對(duì)應(yīng)于塑性應(yīng)變，因此為了進(jìn)一步理解微孔洞貫通動(dòng)力學(xué)行為，圖5給出了不同ILD0/d0條件下微孔洞間韌帶距離中點(diǎn)處典型單元的等效塑性應(yīng)變時(shí)程曲線。利用微孔洞開始貫通的時(shí)間對(duì)應(yīng)條件，可以在圖5中找出微孔洞貫通時(shí)刻所對(duì)應(yīng)的塑性應(yīng)變，圖中用黑點(diǎn)標(biāo)識(shí)，結(jié)果顯示:臨界塑性應(yīng)變閾值約為0.5，隨著ILD0/d0增大，臨界塑性應(yīng)變閾值緩慢減小。

2.2 加載應(yīng)變速率的影響

通過改變邊界速度函數(shù)，可以調(diào)整加載應(yīng)變速率。這里以ILD0/d0=3 模型為例，討論加載應(yīng)變速率對(duì)微孔洞貫通行為的影響。圖6給出了不同加載應(yīng)變速率下ILD/d 時(shí)程曲線和A 和B 點(diǎn)軸向位移時(shí)程曲線，結(jié)果顯示:隨著加載應(yīng)變速率增大，微孔洞長(zhǎng)大速率提高，微孔洞間貫通時(shí)間提前，但對(duì)微孔洞間開始貫通臨界韌帶距離沒有影響，還保持在0.5 倍微孔洞直徑。圖7給出了不同加載應(yīng)變速率條件下微孔洞間韌帶距離中點(diǎn)處典型單元的等效塑性應(yīng)變時(shí)程曲線，結(jié)果也同樣顯示微孔洞開始貫通的臨界塑性應(yīng)變閾值為0.5，圖中黑點(diǎn)標(biāo)識(shí)。圖8給出了不同應(yīng)變速率條件下體胞模型邊界單元上的宏觀平均應(yīng)力隨時(shí)間的變化曲線，從中可見:隨著加載應(yīng)變速率增大，在微孔洞開始貫通時(shí)宏觀平均應(yīng)力(黑點(diǎn)標(biāo)識(shí)處)在增大，不對(duì)應(yīng)于應(yīng)力峰值，也就是說微孔洞間貫通開始后，宏觀平均應(yīng)力場(chǎng)還在增大，材料仍然可以繼續(xù)承載。以前有很多學(xué)者為了回避對(duì)微孔洞貫通過程的理論描述，常常假定微孔洞貫通開始后材料迅速失去承載能力，這個(gè)觀點(diǎn)是值得商榷的。以后在建立合理的微孔洞貫通理論模型時(shí)必須考慮微孔洞貫通也是需要一個(gè)過程的，它不是瞬時(shí)事件。

圖4 ILD/d 演化時(shí)程曲線及A 點(diǎn)和B 點(diǎn)軸向位移時(shí)程曲線Fig.4 Evolution of ILD/d profiles and A and B points’axial displacement profiles

圖5 微孔洞間韌帶中心處典型單元的等效塑性應(yīng)變時(shí)程Fig.5 Effective plastic strain profiles of representative element in the center of intervoid ligament distance

圖6 ILD/d 時(shí)程曲線和A 和B 點(diǎn)軸向位移時(shí)程曲線Fig.6 Evolution of ILD/d profiles and A and B points’axial displacement profiles

圖7 微孔洞間韌帶中心處典型單元的等效塑性應(yīng)變時(shí)程Fig.7 Effective plastic strain profiles of representative element in the center of intervoid ligament distance

圖8 微孔洞間韌帶中心處典型單元的等效塑性應(yīng)變時(shí)程Fig.8 Effective plastic strain profiles of representative element in the center of intervoid ligament distance

3 關(guān)于ILDc和fc的討論

目前，關(guān)于微孔洞開始貫通的準(zhǔn)則主要有兩種觀點(diǎn)，有些學(xué)者采用微孔洞間臨界韌帶距離ILDc，還有些學(xué)者常采用臨界微孔洞體積分?jǐn)?shù)fc.如果基于微孔洞在材料中統(tǒng)計(jì)均勻分布的假定，設(shè)用邊長(zhǎng)為L(zhǎng) 的立方體來代表受損傷材料的一個(gè)特征體積，假設(shè)有8 個(gè)直徑為d 的微孔洞均勻分布于立方體的8 個(gè)角上，如圖9所示。每個(gè)微孔洞只有八分之一屬于此立方體，因此微孔洞開始貫通時(shí)的臨界微孔洞體積分?jǐn)?shù):

由此可見，ILDc和fc是可以相互轉(zhuǎn)換的，在物理本質(zhì)也是一致的。這里近似假定圖9中微孔洞兩兩之間的相互貫通行為是相互獨(dú)立的，當(dāng)然這個(gè)假定具有一定的合理性，因?yàn)槲⒖锥粗g相互貫通是通過塑性應(yīng)變場(chǎng)交疊來實(shí)現(xiàn)的，圖9中其他6 個(gè)微孔洞由于距離相對(duì)比較遠(yuǎn)，對(duì)兩個(gè)微孔洞之間的塑性應(yīng)變場(chǎng)的影響相對(duì)較弱。因此，在微孔洞兩兩之間的相互貫通行為是相互獨(dú)立的假定下，可以通過含兩個(gè)微孔洞圓柱體胞組合排列形成空間的8 個(gè)微孔洞，同時(shí)體胞模型中微孔洞間臨界韌帶距離ILDc在物理本質(zhì)上可以近似替代含8 個(gè)微孔洞的立方體損傷單元中ILDc.對(duì)于純鋁，上述數(shù)值模擬得到ILDc為0.5 倍微孔洞直徑，通過(7)式，計(jì)算得到的fc約為0.15.文獻(xiàn)［14］曾對(duì)純鋁開展過系列層裂實(shí)驗(yàn)研究，基于對(duì)回收樣品中含有的微孔洞宏觀統(tǒng)計(jì)分析，獲得了微孔洞開始貫通的微孔洞體積分?jǐn)?shù)fc約為0.11，與數(shù)值分析結(jié)果非常接近，這表明通過數(shù)值分析獲得微孔洞開始貫通臨界韌帶距離或臨界孔洞體積分?jǐn)?shù)具有一定的合理性。

圖9 臨界孔洞間韌帶距離和臨界孔洞體積分?jǐn)?shù)的幾何說明Fig.9 Geometry relationship between ILDc and fc

4 結(jié)論

采用LS-DYNA 動(dòng)力學(xué)有限元程序，對(duì)軸對(duì)稱體胞模型中雙孔洞的貫通行為進(jìn)行了數(shù)值分析，討論了初始孔洞間韌帶距離和加載應(yīng)變速率的影響，得到了如下幾點(diǎn)認(rèn)識(shí):

1)提出了一個(gè)微孔洞開始貫通的判斷標(biāo)準(zhǔn)，即通過比較微孔洞上頂點(diǎn)和下頂點(diǎn)的軸向位移曲線的明顯偏離來判斷微孔洞開始貫通，數(shù)值結(jié)果顯示該方法具有一定的合理性和有效性。

2)微孔洞從獨(dú)立長(zhǎng)大向微孔洞貫通轉(zhuǎn)變開始的臨界韌帶距離約為0.5 倍微孔洞實(shí)時(shí)直徑，并且不依賴于初始微孔洞間韌帶距離和加載應(yīng)變速率。

3)微孔洞間貫通行為是通過塑性應(yīng)變場(chǎng)相互作用實(shí)現(xiàn)的，數(shù)值計(jì)算結(jié)果顯示:在孔洞間韌帶區(qū)域的塑性應(yīng)變場(chǎng)需要達(dá)到一個(gè)臨界閾值才能驅(qū)動(dòng)微孔洞開始貫通，該臨界塑性應(yīng)變閾值基本保持在0.5左右。

4)微孔洞在材料中統(tǒng)計(jì)均勻分布和微孔洞之間的相互貫通行為是相互獨(dú)立的假定下，臨界孔洞間韌帶距離ILDc判據(jù)和臨界孔洞體積分?jǐn)?shù)判據(jù)在物理本質(zhì)也是一致的，可以相互轉(zhuǎn)換。

References)

［1］ Antoun T，Seaman L，Curran D R，et al.Spall fracture［M］.Berlin:Springer，2003.

［2］ Johnson J N.Dynamic fracture and spallation in ductile solids［J］.Journal of Applied Physics，1981，52:2812 -2825.

［3］ Molinari A，Wright T W.A physical model for nucleation and early growth of voids in ductile materials under dynamic loading［J］.Journal of the Mechanics and Physics of Solids，2005，53:1476-1504.

［4］ Chen D，Tan H，Yu Y，et al.A void coalescence-based spall model［J］.International Journal of Impact Engineering，2006，32(11):1752 -1767.

［5］ Czarnota C，Mercier S，Molinari A.Modeling of nucleation and void growth in dynamic pressure loading，application to spall test in tantalum［J］.International Journal of Fracture，2006，141(1 -2):177 -194.

［6］ Curran D R，Seaman L，Shockey D A.Dynamic failure of solids［J］.Physics Reports，1987，147:253 -388.

［7］ Tvergarrd V，Needleman A.Analysis of the cup-cone fracture in a round tensile bar［J］.Acta Metal.，1984，27:1 -15.

［8］ Brock W，Sun D Z，Honign A.Verification of the transferability of micromechanical parameters by cell model calculations with visco-plastic materials［J］.International Journal Plasticity，1995，11(8):971 -989.

［9］ Pardoen T，Hutchinson J W.An extended model for void growth and coalescence［J］.Journal of the Mechanics and Physics of Solids，2000，48:2467 -2512.

［10］ Brown L M，Embury J D.The initiation and growth of void at second phase particles［C］∥Proceedings 3rd International Conference on Strength of Metals and Alloys.London:Institute of Metals，1973.

［11］ Horstemeyer M F，Matalanis M M，Sieber A M，et al.Micromechanical finite element calculations of temperature and void configuration effects on void growth and coalescence［J］.International Journal Plasticity，2000，16:979 -1015.

［12］ Sepp?l? E T，Belak J，Rudd R E.Onset of void coalescence during dynamic fracture of ductile metals［J］.Physics Review Letters，2004，93(24):245503.

［13］ Sepp?l? E T，Belak J，Rudd R E.Three-dimensional molecular dynamic simulation of void coalescence during dynamic fracture of ductile metals［J］.Physics Review B，2005，71:064112.

［14］ Qi M L，He H L.Simulation of critical behavior on damage evolution［J］.Chinese Physics B，2010，19(3):036201-5.