平面狀態下EWK延性斷裂準則與K準則對比研究

楊鋒平，孫 秦，羅金恒，張 華，張 奕

（1中國石油集團石油管工程技術研究院，石油管力學和環境行為重點實驗室，西安 710065；2西北工業大學 航空學院，西安 710072）

平面狀態下EWK延性斷裂準則與K準則對比研究

楊鋒平1，孫 秦2，羅金恒1，張 華1，張 奕1

（1中國石油集團石油管工程技術研究院，石油管力學和環境行為重點實驗室，西安 710065；2西北工業大學 航空學院，西安 710072）

為研究EWK延性斷裂準則對金屬起裂及擴展的預測能力，將EWK斷裂算法嵌入ABAQUS主程序。針對金屬I型裂紋，在平面應變條件下，通過試驗與有限元計算相互對比的方法，研究了該準則與應力強度因子K準則對緊湊拉伸CT試樣起裂的預測結果，發現兩者有近似等效性。在平面應力條件下，以緊湊拉伸CS試樣為模型，通過有限元計算，提取由EWK控制斷裂的結構KR-曲線，將該曲線與試驗資料對比，發現兩者吻合良好。

EWK延性斷裂準則；應力強度因子；KR-曲線；有限元方法

1 引 言

金屬材料斷裂往往是一種彈塑性斷裂，相關研究工作非常多見。按理論歸類，這些工作主要集中在以下三類：第一、以古典強度理論為依據。這部分內容非常豐富，中文可參考俞茂宏教授的論文[1]或

專著[2]；第二、以斷裂力學的應力強度因子和J積分等為理論依據；第三、以損傷力學為理論依據。三者特別是前兩者為許多工程問題的解決提供了強有力的支持，而第三者則正在成為更好的理論。三者各自的缺陷在于，古典強度理論對材料進行了理想化的假設，對斷裂的研究大部分停留在經驗上。斷裂力學有比較完善的理論體系，但其研究對象是含裂體，且對裂紋作了理想化的假設。對于實際結構從無初始裂紋開始斷裂或含復雜形狀裂紋的問題則無能為力。損傷力學由50年代前蘇聯學者Kachanov提出雛形，其發展特別是連續損傷力學（CDM）的發展是在70年代由Lemaitre教授[3]領導完成。宏細觀結合的損傷力學被認為是研究斷裂最有前途的方法，不過就目前而言，損傷因子的演化方程一直是理論與試驗亟待突破的瓶頸。

基于理論一或三的一類重要工程應用是各種金屬延性斷裂準則，主要解決無初始裂紋結構的開裂問題，在金屬成型等領域起著重要作用。如何評定某種準則的精確程度成為其控制斷裂的關鍵，目前較為多見的方法是通過幾類不同的試驗[4-10]，包括不同尺寸的鐓粗、拉伸試驗，然后建立試驗的有限元模型，將基于該延性斷裂準則的計算結果與試驗對比，以此確定這種延性斷裂準則的優劣，本文作者也研究過此類方法。該方法的不足之處在于其判定的依據僅來自有限的幾類試驗，無法從理論上獲得更多的支持。

由于斷裂力學有比較完善的理論基礎，若某種延性斷裂準則的判定結果與斷裂力學的判定結果一致，則說明該準則與斷裂力學K或J準則等效，從而獲得了一定的理論支持。在等效關系確立后，可進一步認為對于K準則或J準則無法解決的問題（例如無初始裂紋結構開裂問題），可以試圖用該延性斷裂準則解決，當然其最終結果需要試驗評定。這樣做的實質是把基于古典強度理論或損傷力學理論的金屬延性斷裂準則與斷裂力學理論結合在一起考慮，充分發揮了各自的優勢，使得到的結論更加合理。然而很少有學者在同樣的結構中將兩者一起對比，本文的研究工作即在于此。通過兩種緊湊拉伸試樣（兩類模型均滿足小范圍屈服假設，金屬小范圍屈服僅表示塑性區范圍較小，但塑性區內變形仍非常劇烈，材料斷裂屬于韌性斷裂），將Kamoulakos博士提出的EWK延性斷裂準則[11]與K準則對比。結果證實了EWK準則與K準則有較好的互換性。

2 EWK斷裂模型

EWK模型是Kamoulakos博士在McClintock、Wilkins（主要貢獻）等人的工作基礎上加以總結，在假設裂紋形成和擴展是材料一個連續特性的前提下，認為結構斷裂是材料應變損傷的累積結果。裂紋的起始、擴展和結構的斷裂主要取決于危險區域現時和歷史的受力情況，獨立于危險區域的形狀、邊界條件，除非危險區域的形狀和邊界條件影響著它的受力情況。他給出公式如下：

當材料中某個點的損傷變量Dp值達到臨界值時，以該點為中心的某個局部區域（半徑為Rc的圓）斷裂失效。式中εp表示等效塑性應變，p代表該點所受的靜水壓力，Si表示應力偏量；α、β為材料常數，分別代表拉、剪對材料微孔洞生長的影響因子，plim表示不考慮孔洞效應時的理想極限靜水壓力，為材料常數。該模型與損傷力學概念的區別之處在于，認為損傷變量Dp的演化不影響材料的本構關系。因為從表達式看，Dp是由代表外力作用因素的量（w1代表拉、w2代表剪）構成，而不是由諸如材料微孔洞面積等內部因素構成，因此它不影響本構曲線。微孔洞的生長等材料內部因素則由α、β確定，認為是材料常數。從斷裂力學角度考慮，Dp具有應變能密度性質，一定程度上和斷裂力學中的能量釋放率理論吻合。從古典強度理論體系出發，也可將其僅僅理解為與第二強度理論等相似的一種理論，區別在于EWK表達式是一種積分。該理論涉及的4個材料常數Dp、α、β、plim可由4種斷裂實驗確定。本文以鋁合金2A12-CZ（LY12-CZ）為材料，分別進行如圖1所示5種試驗以確定材料常數（α試驗為了得到材料本構關系）。 由試驗確定的該材料常數分別為：Dp=0.234 6，α=2.311 0，β=-0.385 7，plim=1 717MPa。

在有限元計算時，EWK算法的實現是基于單元積分點計算的。將自行編制的EWK斷裂子程序嵌入ABAQUS主程序中。對于每個單元的每個積分點，主程序每完成一次增量步計算，子程序提取一次EWK判斷準則需要的參數，并進行判斷，若達到斷裂門檻值，積分點楊氏模量置零，否則不作任何變化。遍歷所有單元的積分點后，程序進行下一個增量步計算。這樣保證了主、子程序之間的實時通信，直到外載全部添加完畢為止。算法流程圖如圖2所示，該算法將使斷裂發生在單元上，故單元越小，其計算所得的結果越精確，但它對單元的形狀沒有任何要求。另外，只要計算收斂，該算法可一直進行到裂紋擴展完畢。

3 緊湊拉伸（CT）試樣模擬

3.1 方法

該模型目的在于考察K準則與EWK準則的關系。K準則對于彈塑性斷裂的有效性僅僅限于小范圍屈服。因此，為了保證裂紋頂端處于小范圍屈服條件，標準CT試樣（其試樣設計以及結果評定都需要核實小范圍屈服及平面應變條件）是可靠的。用標準CT試樣可以得到材料的斷裂韌性KIC。反過來說，如果用測得的KIC值（一般需要用有效裂紋長度修正后的K準則）來預測CT試驗的P-V曲線，其預測結果應當與試驗曲線一致。若采用EWK準則預測，如果其預測所得的P-V曲線也與試驗曲線一致，則證明了在小范圍屈服條件下，EWK與K準則有等效性；如果不一致，則說明EWK準則的預測能力在材料小范圍屈服條件下沒有很好的精度。因此，本節通過試驗與基于EWK準則的有限元模擬對比的方法來確認上述思想。

試驗依據HB5142-96，試樣厚度B取25mm，寬厚比W/B=2。有限元計算時材料模型使用圖1(a)所得的本構關系，其真實應力應變曲線特征數據如下：彈性模量71 200MPa，泊松比0.33，屈服應力380MPa，極限應力590MPa，最大應變0.15。裂尖模型建成半徑為0.08mm的半圓，裂尖局部細化網格，使最小網格尺度為0.1mm。采用8節點平面應變單元。分別在模型的兩個銷軸孔處添加1mm的位移載荷，并在模型對稱軸處約束好1方向的位移，如圖3所示（線切割及疲勞裂紋建成了開口裂紋）。試樣斷口及測得的P-V曲線取部分點如圖4“Experiment”所示。有限元計算在斷裂處理上采用1節中的辦法，單元變剛度引起的收斂性問題可參看文獻[12]。計算一直持續到載荷加載完畢或計算不收斂，其結果如圖4的“Simulation”所示。

3.2 結果及討論

對于平面應變的應力條件（厚板）而言，一旦起裂，則迅速達到失穩擴展，而不像薄板那樣起裂后還可繼續承載。因此對于圖4中的曲線，以最大載荷為界限，只有前半段曲線對KIC計算有意義。對于前半段曲線的線性段，有限元模擬比試驗略顯剛硬，約2%。原因與計算采用平面應變單元有關。因為真正試驗中，試樣只有沿厚度方向正中間一層才處于純粹平面應變狀態，其上下兩個表面均處于純粹平面應力狀態，沿厚度方向中心向兩邊是一種平面應變狀態向平面應力狀態的過渡。當然這種過渡劇烈，使得絕大部分處于近似的平面應變狀態。而有限元模擬將這“絕大部分”變成了“全部”，于是顯得有些剛硬。有限元計算所得的最大載荷則顯示了該延性斷裂準則的精度。由曲線可以看出，三件試樣所得最大載荷均值為21.41kN，有限元計算所得最大載荷為20.90kN，相對誤差2.38%。這種考察精度的缺點在于忽略了實際裂紋長度的影響。有限元建模時取裂紋長度a′=0.5W=25mm，實際試件在疲勞裂紋預制時不可避免會有長有短，本試驗三個試樣的裂紋長度分別為25.29mm、24.00mm、24.57mm。

為了避免裂紋長度的影響，還可從另一個種更加客觀的角度去評價其精度，即：真實試驗可以得到這種材料的平面應變斷裂韌性KIC，而基于EWK斷裂子程序的有限元計算可看成一種虛擬試驗，通過這種虛擬試驗也可計算得到這種材料的平面應變斷裂韌性KIC′。于是可通過KIC與KIC′的接近程度來考察EWK延性斷裂準則的精度。根據HB5142-96數據處理方法，可得三件試樣所得KIC分別為34.15MPa·m1/2、31.82MPa·m1/2、34.85MPa·m1/2，且三組數據均有效，均值 KIC為 33.61 MPa·m1/2。 對于虛擬試驗，條件載荷 PQ′=20.05kN，最大載荷 Pmax′=21.41kN，載荷比 Pmax′/PQ′=1.04<1.1，裂紋長度 a′=25mm，計算得 KIC′的條件值 KQ′=34.65MPa·m1/2，計算 2.5(KQ′/σs)2=20.78

4 緊湊拉伸（CS）試樣模擬

4.1 方法

本節目的在于考察平面應力條件下K準則與EWK準則的關系。當試樣厚度足夠大，三維約束已經達到平面應變條件時，斷裂韌度KIC是與厚度B和裂紋擴展量Δa無關的材料特性。但當試樣厚度很薄，塑性區尺寸與厚度相比不再是很小，則裂紋尖端應力應變場將處于平面應力狀態。與平面應變斷裂問題相比，平面應力斷裂問題更為復雜。其主要原因在于材料起裂后，結構并未馬上進入失穩，而是可以繼續承受載荷而出現緩慢、穩態的擴展。在這種情況下，材料的抗裂阻力（R）既與結構厚度有關，又與裂紋擴展量Δa有關。對于給定厚度而言，測量K-Δa曲線成為衡量結構斷裂韌性的一種辦法。

于是為了評價在平面應力條件下EWK斷裂算法的準確性，可將給定材料給定厚度的標準試樣用基于EWK斷裂準則的有限元算法計算其斷裂過程，將計算得到的K-Δa曲線與真實試驗所得曲線對比，若兩者符合較好，則說明EWK斷裂算法在平面應力條件下有較好精度。為此，采用HB5261-83金屬板材KR-曲線試驗方法中的標準緊湊拉伸試樣，取厚度B=2mm，寬度W=400mm，初始裂紋a0=160mm，如圖5所示。裂紋尖端建成半徑為0.08mm的半圓模型，此舉既模擬了裂尖的鈍化效果，又可使劃分網格時避免應力奇異。有限元計算時取其一半模型，將對稱線上沿2方向的位移約束，銷釘孔加拉力載荷，同時為約束1方向運動，將對稱線上遠離裂尖的地方添加1方向約束。模型采用4節點平面應力單元，在裂尖以及裂紋可能擴展的路徑上細化單元，使其最小尺寸為0.1mm（基本為肉眼可見最小裂紋的尺寸），以保證獲得較為精確的K-Δa，如圖6所示。材料的本構關系采用2節中的數據。

4.2 結果及分析

有限元計算時，有效裂紋增量的計算根據HB5261-83“有形裂紋長度加塑性區修正量”方法可得如下

其中da表示當前可見裂紋擴展量，即模型中單元丟失的長度，ry表示裂紋前沿塑性區修正量，可用以下公式計算：

其中KI指當前的應力強度因子。KR的獲取也可根據HB5261-83中提供的方法，即

其中P表示當前所加的載荷量，可由約束端沿2方向的支反力求出；F （a/W） 表示形狀系數，可由以下公式計算：

其中a表示有效裂紋長度，其計算表達式如下

由此可知，此方法中Δa的計算依賴于當前的KR值，而KR的計算又依賴于當前的Δa值，這就形成一對矛盾。解決該矛盾可采用迭代求解，即取很小的增量步，利用上一步的Δai-1值，根據公式（8）算得當前步 KR(i)值，然后根據公式（7）、（6）求解當前步的 Δai，再利用當前步 Δai，求解下一步 KR(i+1)。 在第一步 KR(1)的計算中，可認為上一步 Δa0=0。 將得到的數據點（Δai，KR(i)）每十個取一個如圖 7“simulation”所示。并將其擬合成函數，外推至Δa0=0，可得圖7曲線“Fitting of simulation”。真實試驗擬合所得的KRΔa 曲線如圖 7“Data-1”、“DATA-2”所示，其數據源自文獻《中國航空材料手冊（第三卷）》[13]，列出 6 個數據點的信息見表1。

表1 試驗及模擬所得某些裂紋擴展阻力值KR/MPa·m1/2Tab.1 Some resistance values of crack propagation given by tests and simulation,KR/MPa·m1/2

在認為試驗可靠的前提下，將模擬數據與試驗所得兩組數據的均值做相對誤差分析，其值可見表1的最后一行。最大相對誤差出現在Δa=20mm處，即出現在曲線的中部，為6.22%，最小相對誤差為0.96%，出現在裂紋擴展量Δa=10mm處，起裂時，其相對誤差為1.29%，6個數據點平均誤差為3.70%。這些數據說明，基于EWK斷裂準則對平面應力條件下薄板斷裂韌性的預測有很好的精度。

5 結 論

（1）平面應變條件下，通過數值計算與試驗對比的方法確立了EWK延性斷裂準則與K準則之間存在近似等效關系。

（2）平面應力條件下，EWK準則預測裂紋穩態擴展所獲得的KR-曲線與現有試驗數據吻合較好。

[1]俞茂宏.強度理論百年總結[J].力學進展,2004,34(4):529-560.

[2]俞茂宏.雙剪理論及其應用[M].北京:科學出版社,1998.

[3]Lemaitre J,Desmorat R.Engineering Damage Mechanics[M].Berlin:Springer,2005.

[4]虞 松,陳 軍,阮雪榆.韌性斷裂準則的試驗與理論研究[J].中國機械工程,2006,17(19):2049-2052.

[5]Venugopal Rao A,Ramakrishnan N,Krishna kumar R.A comparative evaluation of the theoretical failure criteria for workability in cold forging[J].Journal of Materials Processing Technology,2003,142:29-42.

[6]Ozturk F,Lee D.Analysis of forming limits using ductile fracture criteria[J].Journal of Materials Processing Technology.2004,147:397-404.

[7]Clift SE,Hartley P,Sturgess C E N,Rowe G W.Fracture prediction in plastic deformation process[J].Int.J Mech.Sci.1990,32:1-17.

[8]Bao Y B,Wierzbicki T.A comparative study on various ductile crack formation criteria[J].Journal of Engineering Materials and Technology,2004,126:314-324.

[9]Narayana Murty S V S,Nageswara Rao B,Kashyap B P.Improved ductile fracture criterion for cold forming of spheroidised steel[J].Journal of Materials Processing Technology,2004,147:94-101.

[10]Gouveia B P P A,Rodrigues J M C,Martins P A F.Ductile fracture in metalworking:experimental and theoretical research[J].Journal of Materials Processing Technology,2000,101:52-63.

[11]Kamoulakos A.The ESI-Wilkins-Kamoulakos(EWK)rupture model[M].In:Raabe D,ed.editors.Continuum Scale Simulation of Engineering Materials.Berlin:Wiley-VCH,2004:796-804.

[12]楊鋒平,孫 秦.韌性金屬材料漸進斷裂的有限元算法研究[J].金屬學報,2008,44(4):489-494.

[13]顏鳴皋,劉伯操,李金桂等.中國航空材料手冊（第三卷）[M].北京:中國標準出版社,2002:39-83.

Studies on comparison between EWK ductile criterion and K criterion in plane state

YANG Feng-ping1,SUN Qin2,LUO Jin-heng1,ZHANG Hua1,ZHANG Yi1
(1 The Key Laboratory for Mechanical and Environmental Behavior of Tubular Goods,CNPC Tubular Goods Research Institute,Xi’an 710065,China;2 School of Aeronautics,Northwestern Polytechnical University,Xi’an 710072,China)

In order to study the precision of EWK ductile criterion for mode I cracking of metallic material,a subroutine of EWK fracture criterion was programmed and added to ABAQUS’main program.For 2A12-T4 aluminum alloy,under the stress condition of plane strain state,this criterion was compared with K-Criterion in the initial cracking prediction of compact tension specimen.Through comparison of experimental result and FEM result,it is found that EWK criterion and K-Criterion are approximatly equivalent.Under the stress condition of plane stress state,another compact tension specimen was taken as computational model whose cracking was judged by EWK criterion.KR-curve was derived after this computation was complete.It shows that the KR-curve is close to the published experimental data.

EWK ductile criteria;stress intensity factor(SIF);KR-curve;finite element method(FEM)

TB125

A

1007-7294（2011）05-0506-07

2010-05-10

楊鋒平（1982-），男，博士，中國石油集團石油管工程技術研究院工程師；

孫 秦（1956-），男，西北工業大學教授，博士生導師。