郭榮鑫，索玉霞，牛治亮，夏海廷，顏 峰，張玉波，李一博

(1.昆明理工大學 建筑工程學院 云南省土木工程防災重點實驗室(籌)，云南 昆明 650500；2.昆明理工大學 云南省先進材料力學行為與微結(jié)構(gòu)設計高校重點試驗室，云南 昆明 650500；3.云南農(nóng)業(yè)大學 建筑工程學院，云南 昆明 650201；4.黃河勘測規(guī)劃設計有限公司，河南 鄭州 450000)

1 引 言

功能梯度材料(Functionally Gradient Materials，F(xiàn)GMs)是為了滿足在極限溫度環(huán)境(超高溫、大溫度落差)下正常工作而開發(fā)的一種新型復合材料[1]，Cu/WCpFGMs具備高強度、高硬度、高導電性和優(yōu)異的熱穩(wěn)定性等優(yōu)異力學性能，以及較好的耐電弧侵蝕性、抗熔焊性和耐高溫抗氧化性，被廣泛應用于高壓斷路器、繼電器和開關柜等設備中的高壓電觸頭[2-3]，在電阻焊電極及電火花加工電極等領域也有著廣闊的應用前景[2-3]。

隨著計算機技術的發(fā)展，以及無網(wǎng)格伽遼金法(EFG)[4]和擴展有限元法(XFEM)[5]等方法的提出，國內(nèi)外學者對疲勞裂紋擴展行為的研究取得了重要進展，尤其在一些斷裂參數(shù)(如應力強度因子、J積分等)的計算[6-7]和裂紋擴展路徑的觀測上，但對于裂紋擴展速率方面的研究仍然存在很多問題。楊利等[8]運用分子動力學方法，研究了微觀尺度下鈮元素對單晶γ-TiAl裂紋擴展過程的影響，但該方法針對的是合金材料，僅得到了γ-Ti Al裂紋擴展的軌跡圖和能量演變圖，并未得到裂紋擴展速率的結(jié)果。Harrison[9]提出了一種運用ABAQUS軟件計算普通試件疲勞裂紋擴展速率的方法，該方法通過ABAQUS腳本界面進行模擬，提取出裂紋尖端相關的應力、應變值，用來計算各積分點的疲勞壽命，對裂紋擴展進行定義后，結(jié)合Paris公式，得到了α/β鈦合金的疲勞裂紋擴展速率曲線。然而，該方法需要通過確定裂紋尖端的應力場來確定疲勞壽命，且需要自行定義裂紋尖端網(wǎng)格，通過相關公式計算每一個分析步中的裂紋增量，雖然結(jié)果精確但整個裂紋擴展的分析比較復雜。Nittur等[10]提出了一種采用有限元方法和塑性耗能標準相結(jié)合的數(shù)值計算方法，得到了鉻鎳鐵合金和鈦合金試件在循環(huán)載荷下的疲勞裂紋擴展速率，但該研究僅僅是針對合金等均質(zhì)材料，對 FGMs 的研究還有欠缺。Bhattacharya[11]和王志勇[12]等運用擴展有限元方法(XFEM)，研究了FGMs的疲勞裂紋擴展行為，利用相互作用積分法計算應力強度因子，得到了FGMs的疲勞壽命及裂紋擴展路徑。然而，這種方法在輔助應力場和位移場的確定上較復雜，并且未得出材料的疲勞 裂 紋 擴 展 速 率。JEONG-HO KIM 等[13-14]對FRANC2D 程序進行了改進，對FGMs的裂紋擴展行為進行數(shù)值模擬，但只得出了FGMs的應力強度因子值和裂紋路徑的變化情況，該種方法研究的是材料參數(shù)梯度呈指數(shù)變化的FGMs，對于疊層FGMs而言，在模型的建立、材料參數(shù)的輸入上較復雜。

因此，如何快速建立疊層FGMs的有限元模型，快速計算出疊層FGMs在裂紋擴展中的應力強度因子，從而得出疲勞裂紋擴展速率曲線，是亟待解決的問題之一?；贔RANC2D(Fracture Analysis Code in 2 Dimensions)有限元斷裂分析軟件，建立了一種FRANC2D 與ABAQUS相結(jié)合的FGMs裂紋擴展模擬方法，對Cu/WCp雙層及多層FGMs標準拉伸試樣的疲勞裂紋擴展行為進行數(shù)值模擬，該方法可以快速地建立疊層FGMs的有限元模型，迅速地計算出裂紋擴展中的應力強度因子，從而繪制出疲勞裂紋擴展速率曲線，在比較復雜的疊層FGMs裂紋擴展模擬研究中有很大優(yōu)勢。

2 FGMs裂紋擴展速率的數(shù)值模擬方法

FRANC2D 斷裂有限元分析軟件是由美國Cornell大學斷裂力學小組開發(fā)的，是用于裂紋分析的專業(yè)軟件，其最大優(yōu)點在于在模擬自動擴展裂紋的同時，可以自動重置網(wǎng)格，并計算裂紋尖端的應力強度因子(SIF)，節(jié)約了時間又提高了計算精度[15-16]，在模擬裂紋的起裂和擴展方面有很大優(yōu)勢。

2.1 應力強度因子

在FRANC2D 中，對裂紋擴展過程中應力強度因子的計算主要基于三種方法：位移相關法(DCT)、虛裂紋閉合法(MD CRK-CLOS)和J 積分法(JINTEGRAL)。Ⅰ型裂紋的應力強度因子可以使用位移相關法來進行計算[17]，這種方法基于節(jié)點位移來計算用應力強度因子，且在有限元方法計算中，基于位移公式的計算方法比基于應力的計算方法更為精確。

采用DCT 方法對層間應力強度因子進行計算，當設置裂紋和裂紋擴展增量Δa后，F(xiàn)RANC2D 將對裂紋尖端和裂紋尖端附近的單元重新進行單元結(jié)構(gòu)劃分，由裂尖附近裂紋面的節(jié)點位移，可以實現(xiàn)每一步裂紋擴展后應力強度因子的計算。通過DCT 方法得到的應力強度因子為[18]：

式中，G=E/(1+ν)，為剪切模量；ν為泊松比；對于平面應變狀態(tài)，к=3-4ν；對于平面應力狀態(tài)，к=(3-4ν)/(1+ν)；a、b、c、d分別為裂紋尖端有限元單元上的節(jié)點，如圖1所示；V a、V b、V c、V d分別為裂紋尖端相對于節(jié)點a、b、c、d處的y方向位移；U a、U b、U c、U d分別為裂紋尖端相對于節(jié)點a、b、c、d處的x方向位移；Δa為裂紋尖端的單元長度。由DCT 方法計算出的應力強度因子將被代入所選定的斷裂判據(jù)中，斷裂判據(jù)用來判斷局部裂紋的穩(wěn)定性及預測裂紋擴展角度，新裂紋尖端的位置則完全取決于所設定的裂紋擴展增量Δa。

圖1 裂紋尖端附近單元劃分Fig.1 Element division near crack tip

2.2 裂紋擴展速率

如果在應力循環(huán)ΔN次后，裂紋擴展量為Δa，則應力每循環(huán)一周，裂紋擴展量為Δa/ΔN(mm/周)，這稱為裂紋擴展速率。對于疲勞裂紋擴展速率的研究，主要在于尋求裂紋擴展速率與有關各種力學參量之間的數(shù)學表達式。斷裂力學認為裂紋的擴展速率同裂紋尖端應力強度因子幅值ΔK有關，P.C.Paris提出的描述這種關系的裂紋擴展公式為：

式中，a為裂紋長度，N為循環(huán)次數(shù)，C和m是與試驗條件(環(huán)境、頻率、溫度和應力比R 等)有關的材料常數(shù)。ΔK稱為應力強度因子幅，是每個周期內(nèi)的最大與最小應力強度因子之差：

對于現(xiàn)有的有限元軟件來說，一般很難直接輸出疲勞裂紋擴展速率曲線，同樣，F(xiàn)RANC2D 也不能直接給出疲勞裂紋擴展速率曲線。因此，本研究根據(jù)Cornell大學斷裂力學小組所提供的表格，通過FRANC2D 計算出每一步裂紋擴展增量所對應的的應力強度因子值，從而確定出應力強度因子幅(ΔK)，結(jié)合Paris公式以及由試驗確定的C和m值，可以繪制出疲勞裂紋擴展速率曲線(da/d N-ΔK)；再根據(jù)裂紋尖端距界面的距離，繪制出隨裂尖距界面距離變化的裂紋擴展速率曲線(da/d N-l)。

本研究結(jié)合ABAQUS和FRANC2D 兩個有限元軟件，對疊層FGMs的裂紋擴展行為進行數(shù)值模擬，最終得到疲勞裂紋擴展速率曲線，過程流程圖如圖2所示。

3 FGMs疲勞裂紋擴展行為數(shù)值模擬過程

為研究Cu/WCp疊層FGMs的疲勞裂紋擴展行為，驗證本研究中所采用數(shù)值模擬方法的準確性及可靠性，文中基于FRANC2D 斷裂分析軟件，分別對疲勞試驗所采用的Cu/WCp疊層FGMs拉伸試件進行疲勞裂紋擴展數(shù)值模擬。

3.1 有限元模型的建立

疲勞試驗研究所采用的試件為雙層及多層Cu/WCpFGMs啞鈴狀單邊缺口拉伸試件，經(jīng)粉末冶金熱壓燒結(jié)制備而成，試件平行段尺為30×12×5mm，厚度為2mm，如圖3所示。初始裂紋分別在不同WCp含量方向。試件分層布置為：雙層啞鈴狀單邊缺口拉伸試件中，梯度層的厚度從左到右分別為3和9mm，WCp體積分數(shù)變化為3%～15%；多層啞鈴狀單邊缺口拉伸試件中，梯度層的厚度從左到右依次為2、1和9mm，WCp體積分數(shù)變化為3%～9%～15%。

因此，本研究使用ABAQUS 有限元軟件，完成FGMs幾何模型的建立、材料參數(shù)的輸入(彈性模量E和泊松比ν)以及有限元網(wǎng)格的劃分，再經(jīng)MATLAB程序?qū)⑺Ｐ娃D(zhuǎn)化為適用于FRANC2D 斷裂模擬的輸入程序。這種建模方法簡單方便，能夠得到已經(jīng)賦予了材料屬性并進行了網(wǎng)格劃分的有限元模型，在FRAN2D 中只需對模型的邊界條件(荷載和約束)進行設置，即可進行斷裂模擬。

在ABAQUS有限元軟件中建立試件的結(jié)構(gòu)模型和有限元模型，如圖4所示。采用六節(jié)點三角形單元進行自由網(wǎng)格劃分，在模型上下兩邊和弧形部分所布的種子長度為0.2mm，考慮到應力應變在裂紋尖端具有奇異性，將裂紋附近的網(wǎng)格進行加密，所布種子長度為0.1mm；將相應的材料參數(shù)E和ν輸入有限元軟件中，具體數(shù)值見表1。隨后，在FRANC2D中對所建模型施加荷載和約束，通過在模型下邊界定義約束邊界的起始點、起始點的相鄰點和約束邊界的截止點，能夠?qū)κ┘铀胶拓Q直兩個方向的約束；同時，在模型上邊界施加一個均勻拉力，參考試驗，控制最大力不變，輸入的拉力大小分別為40 MPa(當裂紋從WCp體積分數(shù)為3%方向開始擴展)和22 MPa(當裂紋從WCp體積分數(shù)為15%方向開始擴展)。

圖2 模擬計算疲勞裂紋擴展速率曲線的流程圖Fig.2 Flowchart of simulating fatigue crack growth rate curves

圖3 雙層3%～15%(a)、三層3%～9%～15%(b)、雙層15%～3%(c)和三層15%～9%～3%(d)的Cu/WCp FGMs試驗疲勞試樣Fig.3 Material specimen in experiment of 3%～15%bi-layered(a)，3%～9%～15% multilayer(b)，15%～3%bi-layered(c)，and 15%～9%～3% multilayer(d)

圖4 雙層(a)、三層(b)FGMs幾何模型和Cu/WCp FGMs幾何模型及有限元模型(c)Fig.4 Geometric model of bi-layered(a)，geometric model of multilayer(b)and finite element model(c)in Cu/WCp FGMs

表1 材料參數(shù)Table 1 Material parameters

3.2 分析過程

FRANC2D 對裂紋擴展行為的分析主要分為前處理、模擬計算和后處理三個步驟，有平面應力(Plane Stress)、平 面 應 變 (Plane Strain)、軸 對 稱(Axisymmetric)和彎曲(Bending)四種分析類型，本研究所模擬的對象為薄板狀拉伸試件，因此選擇的分析類型為平面應力狀態(tài)。

FRANC2D 在普通有限元的等參數(shù)單元基礎上，引入奇異單元，采用標準的具有二次形函數(shù)的六節(jié)點三角形單元(Bilinear3 side)和八節(jié)點四邊形單元(Bilinear4 side)。在進行網(wǎng)格劃分時，采用裂紋尖端的K主導區(qū)適當加密網(wǎng)格，對裂紋尖端徑向和環(huán)向的網(wǎng)格進行細化。這種特殊梯度單元，對各向異性材料(如FGMs)裂紋尖端的網(wǎng)格劃分具有很大的優(yōu)勢，即使采用較稀疏的網(wǎng)格劃分也能得到滿足模擬要求的結(jié)果。由用戶定義的裂紋增量確定新的裂紋尖端位置后，沿著增量上裂紋擴展路徑的單元將被消除，并自動進行局部網(wǎng)格的重構(gòu)。通過這種網(wǎng)格局部重新劃分方法以及特殊網(wǎng)格的移動，可以進行試件疲勞裂紋擴展的模擬，如圖5所示。

圖5 裂紋尖端局部網(wǎng)格重新劃分 (a)用戶在初始幾何構(gòu)型中自定義裂紋幾何形狀；(b)刪除新裂紋尖端位置附近的網(wǎng)格并進行網(wǎng)格更新與重構(gòu)；(c)在裂紋尖端建六節(jié)點三角形單元；(d)增加裂紋尖端附近徑向和環(huán)向單元數(shù)量并進行局部網(wǎng)格細化Fig.5 Mesh local area around the crack tip (a)user defines the crack geometry in the initial geometry；(b)mesh near the tip of the new crack is removed and the mesh is updated and reconstructed；(c)bilinear3 side are constructed at the crack tip；(d)number of radial and ring elements near the crack tip is increased and local mesh is refined

根據(jù)試驗中試件的情況，在FRANC2D 中先引入缺口，長度為1mm；通過選取裂紋的中心點坐標和裂紋尖端位置，定義一個線彈性裂紋，長度為0.5mm；設置裂紋擴展量Δa為65μm，分別對不同的Cu/WCp疊層FGMs模型進行裂紋擴展，圖6為雙層FGMs拉伸試件的裂紋擴展路徑。

圖6 裂紋擴展 (a)裂紋開始擴展；(b)裂紋尖端到達WCp 含量3%(體積分數(shù))和15%(體積分數(shù))界面；(c)裂紋擴展結(jié)束；(d)試件變形情況Fig.6 Crack propagation (a)start propagating；(b)crack tip reaches the interfaces between two contents of WCp 3%(Vol)and 15%(Vol)；(c)stop propagating；(d)deformation of specimen

4 結(jié)果分析與對比驗證

采用FRANC2D 與ABAQUS結(jié)合的數(shù)值模擬方法，計算出了Cu/WCp疊層FGMs單邊缺口拉伸試件的應力強度因子值，結(jié)合Paris公式，得出了疲勞裂紋擴展速率曲線。為了對該數(shù)值模擬方法的可行性進行驗證，在MTS810疲勞實驗機上進行疲勞裂紋擴展試驗。用電火花線切割機分別在WCp體積分數(shù)3%(或15%)一側(cè)預制1mm 的缺口，根據(jù)規(guī)范要求[19]，預制長度約為0.5mm 的疲勞裂紋，采用長焦距顯微鏡M75進行裂紋長度監(jiān)測。疲勞裂紋擴展過程中保持最大循環(huán)載荷Pmax不變，采用正弦波加載，頻率為10 Hz，應力比為0.1，同時記錄裂紋擴展量Δa i以及相鄰兩個裂紋尖端所對應的疲勞載荷循環(huán)次數(shù)，利用傳統(tǒng)均質(zhì)材料的計算公式來確定應力強度因子值，繪制出裂紋擴展速率曲線。研究結(jié)果顯示，該方法可以快速地計算出Cu/WCp疊層FGMs的應力強度因子值，并可用于模擬疊層FGMs的疲勞裂紋擴展速率。

4.1 裂紋擴展速率和應力強度因子幅度的關系(da/d N-ΔK)

圖7為數(shù)值模擬方法及試驗方法下，雙層及多層Cu/WCpFGMs裂紋擴展速率da/d N和應力強度因子幅度ΔK的關系對比曲線圖。從圖可見，梯度層界面兩端的da/d N有很大的差異。由7(a)～(b)可知，雙層Cu/WCpFGMs中，當控制裂紋從低WCp體積分數(shù)(3%)層向高WCp體積分數(shù)(15%)層擴展時，da/d N隨著ΔK的增大而增大，在靠近界面前變化趨勢有所減慢，穿過界面后迅速增大直至試件斷裂；當控制裂紋從高WCp體積分數(shù)(15%)層開始擴展時，在靠近界面前da/d N變化趨勢減緩甚至下降，穿過界面后迅速增大直至試件斷裂。由圖7(c)～(d)可知，多層Cu/WCpFGMs中裂紋擴展速率也存在和雙層Cu/WCpFGMs類似的變化規(guī)律，當控制裂紋從3%(體積分數(shù))WCp層開始擴展時，隨著ΔK的增大，在靠近第一個界面(WCp體積分數(shù)為3%與WCp層與體積分數(shù)為9%)時，裂紋擴展速率減小，穿過界面后速率迅速上升；當靠近第二個界面(WCp體積分數(shù)為9%與 WCp層與體積分數(shù)為15%)時，裂紋擴展速率減小，穿過界面后迅速上升直至試件斷裂；當控制裂紋從15%(體積分數(shù))WCp層開始擴展時，隨著ΔK的增大，裂紋擴展速率分別在靠近第一個界面(WCp體積分數(shù)為15%與WCp層與體積分數(shù)為9%)和第二個界面(WCp體積分數(shù)為9%與WCp層與體積分數(shù)為3%)時，都有降低的趨勢。從數(shù)值模擬方法的結(jié)果中同樣可以發(fā)現(xiàn)，無論裂紋從高WCp體積分數(shù)向低WCp體積分數(shù)擴展還是從低WCp體積分數(shù)向高WCp體積分數(shù)擴展，梯度層界面兩端的da/d N都有很大的差異，裂紋擴展速率在靠近界面時均減緩或降低。因此，梯度層間界面對功能梯度材料裂紋擴展的整個階段都存在很大影響，它可以起到減緩裂紋擴展的作用。

應用FRANC2D 與ABAQUS結(jié)合的數(shù)值模擬方法，可以更為快速的計算出FGMs在疲勞載荷下裂紋尖端的應力強度因子，并用于分析FGMs中梯度層的存在對疲勞裂紋擴展速率的影響規(guī)律。然而，數(shù)值模擬結(jié)果中梯度層界面的應力強度因子值較試驗中的低，這主要是因為WCp在銅基體中的分散并不能保證是完全均勻的，試驗中采用傳統(tǒng)均質(zhì)材料的計算方法得到的應力強度因子值與有限元模型建立中對材料的均質(zhì)化假設，均會造成試驗和模擬結(jié)果曲線的不完全吻合。

圖7 裂紋擴展速率曲線(da/d N-ΔK)Fig.7 Crack growth rate curve(da/d N-ΔK)(a)3%～15%Bi-layered FGMS；(b)15%～3%Bi-layered FGMS；(c)3%～9%～15% Multilayer FGMS；(d)15%～9%～3% Multilayer FGMS

4.2 裂紋擴展速率和裂紋尖端距試件界面距離的關系(da/d N-l)

圖8 裂紋擴展速率曲線(da/d N-l)Fig.8 Crack growth rate curve(da/d N-l)

圖8 為數(shù)值模擬方法下，雙層及多層Cu/WCpFGMs裂紋擴展速率da/d N和裂紋尖端距試件界面距離l的關系曲線。從圖可見，雙層梯度材料的裂紋擴展速率明顯快于多層功能梯度材料的裂紋擴展速率；在雙層功能梯度材料中，裂紋從高含量梯度(WCp體積分數(shù)為15%)層開始擴展，裂紋擴展速率快于裂紋從低含量梯度(WCp體積分數(shù)為3%)層開始擴展的裂紋擴展速率；在多層功能梯度材料中也存在相同的變化規(guī)律，裂紋從高含量梯度(WCp體積分數(shù)為15%)層開始擴展時，裂紋擴展速率也要稍快一些。說明梯度含量的變化對功能梯度材料的疲勞裂紋擴展影響很大；同時，梯度層的數(shù)量對功能梯度材料的疲勞裂紋擴展速率也有一定的影響。

5 結(jié) 論

本研究運用二維斷裂分析有限元軟件FRANC2D和ABAQUS相結(jié)合的方法，對Cu/WCp疊層FGMs的疲勞裂紋擴展行為進行數(shù)值模擬，分別對四種不同模型進行了數(shù)值分析，并與試驗計算結(jié)果進行對比，得到以下結(jié)論：

1.所建立的FRANC2D 與ABAQUS 相結(jié)合的FGMs裂紋擴展模擬方法，可以較快速建立有限元模型，迅速計算出裂紋擴展中的應力強度因子，有利于比較復雜的疊層FGMs的裂紋擴展模擬研究。

2.應用該方法對Cu/WCp疊層FGMs單邊缺口拉伸試件進行疲勞裂紋擴展數(shù)值模擬，得到了裂紋從不同WCp含量(體積分數(shù)為3%或15%)層開始擴展時的疲勞裂紋擴展速率曲線。發(fā)現(xiàn)曲線吻合度較高且變化規(guī)律及趨勢基本相同，證明使用FRANC2D 計算FGMs疲勞裂紋擴展速率是可行的。

3.無論裂紋從低含量(WCp體積分數(shù)為3%)層還是從高含量層(WCp體積分數(shù)為15%)開始擴展，裂紋擴展速率曲線在界面處均明顯減緩甚至下降，裂紋尖端穿過界面后，裂紋擴展速率迅速增大直至試件破壞，說明梯度變化及梯度層間界面的存在對功能梯度材料裂紋擴展的整個階段都存在很大影響。

4.梯度含量的變化對功能梯度材料的疲勞裂紋擴展影響很大；同時，梯度層的數(shù)量對功能梯度材料的疲勞裂紋擴展速率也有一定影響。