高應變率下延性金屬中微孔洞貫通行為的數值分析

王永剛，徐東明

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

0 引言

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

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

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

1 計算模型建立

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

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

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

式中:σij為應力張量;p 為壓力;sij為偏應力張量。

式中:σy為屈服應力考慮塑性硬化時，

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

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

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

計算參數見表1.

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

2 計算結果討論

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

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

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

圖3 不同時刻微孔洞周圍等效塑性應變云圖Fig.3 Effective plastic strain contours around void

2.1 ILD0的影響

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

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

2.2 加載應變速率的影響

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

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

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

圖6 ILD/d 時程曲線和A 和B 點軸向位移時程曲線Fig.6 Evolution of ILD/d profiles and A and B points’axial displacement profiles

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

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

3 關于ILDc和fc的討論

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

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

圖9 臨界孔洞間韌帶距離和臨界孔洞體積分數的幾何說明Fig.9 Geometry relationship between ILDc and fc

4 結論

采用LS-DYNA 動力學有限元程序，對軸對稱體胞模型中雙孔洞的貫通行為進行了數值分析，討論了初始孔洞間韌帶距離和加載應變速率的影響，得到了如下幾點認識:

1)提出了一個微孔洞開始貫通的判斷標準，即通過比較微孔洞上頂點和下頂點的軸向位移曲線的明顯偏離來判斷微孔洞開始貫通，數值結果顯示該方法具有一定的合理性和有效性。

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

3)微孔洞間貫通行為是通過塑性應變場相互作用實現的，數值計算結果顯示:在孔洞間韌帶區(qū)域的塑性應變場需要達到一個臨界閾值才能驅動微孔洞開始貫通，該臨界塑性應變閾值基本保持在0.5左右。

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

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.