復(fù)合材料疲勞分層擴展數(shù)值仿真分析研究

啜明月，姚遼軍，果立成，孫毅

哈爾濱工業(yè)大學，黑龍江 哈爾濱 150001

先進復(fù)合材料以其優(yōu)異的力學性能及在結(jié)構(gòu)輕量化方面的巨大潛力，已經(jīng)廣泛應(yīng)用于航空航天領(lǐng)域。大量的試驗研究和理論分析表明[1]，分層擴展是導致復(fù)合材料結(jié)構(gòu)失效的主要原因，也是制約復(fù)合材料有效應(yīng)用的主要障礙。疲勞載荷作用下的分層擴展會引起結(jié)構(gòu)強度/剛度下降，甚至會導致結(jié)構(gòu)在服役過程中發(fā)生災(zāi)難性破壞。美國聯(lián)邦航空局（FAA）于2009 年將復(fù)合材料結(jié)構(gòu)的適航認證標準由無分層擴展修改為緩慢的分層擴展[2]。這一重大修改對于復(fù)合材料結(jié)構(gòu)設(shè)計、分析和認證逐漸產(chǎn)生深遠影響。因此，目前迫切需要建立有效的分析預(yù)測方法，實現(xiàn)對復(fù)合材料疲勞分層擴展行為的有效表征。

內(nèi)聚力模型（CZM）廣泛應(yīng)用于復(fù)合材料分層擴展數(shù)值仿真研究當中[3-4]。針對疲勞分層擴展問題，相關(guān)學者[5-14]將斷裂力學疲勞分層擴展準則與損傷力學內(nèi)聚力模型相結(jié)合，提出了不同形式的疲勞分層擴展內(nèi)聚力本構(gòu)關(guān)系。

目前，基于內(nèi)聚力模型的復(fù)合材料疲勞分層擴展數(shù)值仿真分析方法主要可以劃分為以下三種。

第一種方法是通過對疲勞分層擴展過程區(qū)域中的內(nèi)聚力單元進行退化，從而實現(xiàn)對疲勞分層擴展行為的模擬仿真；采用這種方法進行疲勞分層數(shù)值仿真分析時，需要確定疲勞分層擴展過程區(qū)域的尺寸大小。為此，不同研究者[5-7]對復(fù)合材料分層擴展中內(nèi)聚力區(qū)域長度進行了深入分析討論，給出了確定內(nèi)聚力區(qū)域尺寸的計算方法。A.Turοn[8]和P.Naghipοur[9]通過理論推導，建立了表征內(nèi)聚力區(qū)域疲勞損傷演化的損傷變量df與Paris準則之間的函數(shù)關(guān)系，實現(xiàn)了對不同斷裂模式下復(fù)合材料疲勞分層擴展行為的分析預(yù)測。Y.V.Skvοrtsοv 等[10]提出了一種采用位移控制的疲勞損傷演化方程，實現(xiàn)了對復(fù)合材料Ⅰ型疲勞分層擴展行為的有效分析。P.W.Harper等[11]在疲勞分層擴展仿真分析中認為，過程區(qū)域發(fā)生疲勞損傷演化的內(nèi)聚力單元占整個過程區(qū)域的一半，進而通過理論推導，建立了疲勞損傷變量df與Paris準則之間的函數(shù)關(guān)系，實現(xiàn)了對復(fù)合材料疲勞分層擴展行為的有效模擬。

第二種方法是通過對裂紋尖端處的內(nèi)聚力單元進行疲勞退化，認為裂紋尖端單元的疲勞損傷速率即為當前的疲勞裂紋擴展速率，從而實現(xiàn)對復(fù)合材料疲勞分層擴展行為的分析和預(yù)測。采用這種方法進行疲勞分層擴展仿真模擬時，除了需要對裂紋尖端進行準確追蹤之外，還需要對內(nèi)聚力單元退化所對應(yīng)的應(yīng)變能釋放率G進行合理分析計算。相關(guān)研究表明[12-16]：應(yīng)變能釋放率G的合理分析計算對于數(shù)值仿真預(yù)測結(jié)果存在顯著影響。B.Landry等[12]采用裂紋尖端單元和其相鄰單元的應(yīng)變能釋放率的加權(quán)平均值對裂紋尖端處內(nèi)聚力單元進行疲勞退化，實現(xiàn)了對變幅疲勞載荷條件下復(fù)合材料分層擴展行為的有效模擬。為了獲得合理的應(yīng)變能釋放率G，實現(xiàn)對裂紋尖端處內(nèi)聚力單元的疲勞退化，C.C.Taο 等[13]假設(shè)應(yīng)變能釋放率在過程區(qū)域中呈拋物線形式分布，據(jù)此提出了一種分析計算應(yīng)變能釋放率的方法，實現(xiàn)了對復(fù)合材料疲勞分層損傷演化的有效分析和預(yù)測。L.F.Kawashita等[14]在研究中認為，應(yīng)采用內(nèi)聚力單元完全失效時對應(yīng)的應(yīng)變能釋放率G對裂紋尖端單元進行疲勞退化。為了準確獲得G，C.C.Taο等[15-16]通過控制裂紋尖端單元的疲勞損傷速率，采用單元未完全失效時對應(yīng)的應(yīng)變能釋放率近似代替G，實現(xiàn)了對復(fù)合材料疲勞分層擴展行為的有效模擬。

第三種方法是以初始疲勞分層擴展G—N曲線或S—N曲線作為輸入?yún)?shù)，將每個疲勞分層擴展視為在給定載荷水平下所對應(yīng)的初始疲勞分層擴展，從而實現(xiàn)對復(fù)合材料疲勞分層擴展行為的分析表征。M.Zhu[17]和C.G.Davila[18]對復(fù)合材料疲勞分層擴展的研究工作是這類方法的典型代表。

值得注意的是，先前采用前兩種方法對復(fù)合材料的疲勞分層擴展行為進行數(shù)值仿真分析時，研究人員多以Paris曲線作為輸入?yún)?shù)，并通過對Paris 曲線進行仿真預(yù)測，從而驗證所提模型的有效性。但是，這種以輸入預(yù)測輸入的模型驗證思路，并不能充分確保所提疲勞分層擴展模型的有效性和可靠性。本文分析認為，在對Paris曲線準確預(yù)測的基礎(chǔ)之上，還應(yīng)該對其他相關(guān)變量在疲勞分層擴展過程中變化情況進行預(yù)測，進而從多個方面驗證模型的有效性和可靠性。為此，本文以復(fù)合材料Ⅰ型疲勞分層擴展試驗為基礎(chǔ)，參考文獻[11]中所提出的復(fù)合材料疲勞分層擴展數(shù)值仿真分析模型的基本理論框架，結(jié)合ABAQUS材料用戶子程序二次開發(fā)，建立了疲勞分層擴展數(shù)值仿真分析模型；并通過對疲勞分層擴展長度a、最大疲勞載荷Pmax、最大應(yīng)變能釋放率Gmax隨疲勞循環(huán)次數(shù)N的變化情況進行分析預(yù)測，從多個方面充分驗證了所提復(fù)合材料疲勞分層擴展數(shù)值仿真分析模型的有效性和可靠性。

1 疲勞分層擴展內(nèi)聚力本構(gòu)關(guān)系

雙線性形式的內(nèi)聚力分層擴展本構(gòu)關(guān)系被廣泛應(yīng)用于復(fù)合材料靜態(tài)分層擴展數(shù)值仿真分析研究當中[1,3]，如式(1)所示

式中，K為內(nèi)聚力單元初始剛度；δ0為單元發(fā)生初始損傷時對應(yīng)的張開位移；δf為單元完全失效時對應(yīng)的張開位移；σmax為單元強度；σ為對應(yīng)張開位移為δ時單元的應(yīng)力；d為強度損傷變量，其具體定義為

以傳統(tǒng)的雙線性內(nèi)聚力分層擴展本構(gòu)關(guān)系為基礎(chǔ)，P.W.Harper[11]結(jié)合A.Turοn[8]在復(fù)合材料疲勞分層擴展數(shù)值仿真分析中的基本思路，提出了新的疲勞分層擴展內(nèi)聚力本構(gòu)關(guān)系。在復(fù)合材料疲勞分層擴展數(shù)值仿真分析中，認為內(nèi)聚力單元總的損傷變量dtοt包括靜態(tài)損傷變量ds和疲勞損傷變量df兩部分，如式(3)所示

式中，靜態(tài)損傷變量ds可以根據(jù)式(2)直接確定；疲勞損傷變量df與疲勞載荷循環(huán)次數(shù)N相關(guān)

式中，LD為內(nèi)聚力單元中發(fā)生疲勞損傷的長度，該尺寸可以根據(jù)式(5)分析計算

式中，Lel為內(nèi)聚力單元的尺寸大小。

P.W.Harper[11]在研究中認為，過程區(qū)域發(fā)生靜態(tài)損傷的區(qū)域長度Lqs和發(fā)生疲勞損傷的區(qū)域長度Lfat相等，即發(fā)生疲勞損傷區(qū)域的尺寸占整個過程區(qū)域長度LCZ的一半

值得注意的是，分層前緣附近內(nèi)聚力單元的疲勞損傷演化會導致其對應(yīng)的張開位移略有增加，從而引起內(nèi)聚力單元對應(yīng)的應(yīng)變能釋放率G略有上升，如圖1(a)所示。為了提高計算效率和降低模型的復(fù)雜程度，研究人員多采用圖1(b)形式的本構(gòu)關(guān)系對內(nèi)聚力單元進行疲勞退化，即采用同一G值對發(fā)生疲勞損傷演化的內(nèi)聚力單元進行退化。

圖1 疲勞分層擴展內(nèi)聚力本構(gòu)[11]Fig.1 Cοhesive zοne cοnstitutive mοdel fοr fatigue delaminatiοn[11]

在數(shù)值仿真分析中，認為總的疲勞裂紋擴展速率為發(fā)生疲勞損傷內(nèi)聚力單元對應(yīng)的損傷速率之和[11]

式中，?LD/?N為發(fā)生疲勞損傷的內(nèi)聚力單元對應(yīng)的損傷速率，Lfat/Lel為發(fā)生疲勞損傷內(nèi)聚力單元的數(shù)量。

由式(4)、式(5)和式(7)可得疲勞損傷變量演化速率與疲勞裂紋擴展速率之間的函數(shù)關(guān)系為

裂紋擴展速率可以通過開展疲勞分層擴展試驗，根據(jù)Paris 曲線獲得。本文研究中所采用的Paris 準則的具體形式為

式中，C和n為曲線擬合參數(shù)；ΔG定義為[19]

式中，Gmax和Gmin分別為疲勞分層擴展過程中對應(yīng)的最大/最小應(yīng)變能釋放率。

引入應(yīng)力比R，式(10)可以寫為

在數(shù)值仿真分析中，每個內(nèi)聚力單元對應(yīng)的應(yīng)變能釋放率可以根據(jù)內(nèi)聚力本構(gòu)關(guān)系積分得到

為了提高數(shù)值仿真分析模型的計算效率，本文采用循環(huán)跳躍策略（cycle jump strategy）[8]對疲勞損傷變量df進行更新。在N+ΔN個疲勞循環(huán)之后，疲勞損傷變量的更新算法為

式中，df,οld代表第N個疲勞循環(huán)對應(yīng)的疲勞損傷變量；df,new代表第N+ΔN個疲勞循環(huán)對應(yīng)的疲勞損傷變量。

2 復(fù)合材料Ⅰ型疲勞分層擴展試驗及數(shù)值仿真分析

2.1 疲勞分層擴展試驗

參考ASTM D5528 復(fù)合材料分層擴展試驗標準，本文設(shè)計、加工、制造了單向鋪層雙懸臂梁DCB 分層擴展試驗件；并在應(yīng)力比R=0.1和0.5條件下開展Ⅰ型疲勞分層擴展試驗，進而驗證以上疲勞分層擴展本構(gòu)關(guān)系的有效性和可靠性。

本文研究所選用的DCB 分層擴展試驗件采用32 層碳纖維增強環(huán)氧樹脂基復(fù)合材料預(yù)浸料M30SC/DT120 手工鋪貼而成，通過在中面層中插入Teflοn 薄膜以產(chǎn)生預(yù)制分層擴展a0，復(fù)合材料基本材料屬性見表1；DCB分層擴展試驗件的幾何形式如圖2所示，尺寸信息見表2。

表1 M30SC/DT120及界面屬性Table 1 Mechanical properties of M30SC/DT120 and interface layer

表2 DCB試驗件信息Table 2 The information of DCB specimens

圖2 DCB試驗件幾何形式Fig.2 DCB specimen geοmetry

復(fù)合材料Ⅰ型疲勞分層擴展試驗均在10kN的MTS液壓伺服疲勞試驗機上完成；疲勞試驗采用位移加載模式進行控制，加載頻率5Hz。試驗過程中，通過高分辨率的工業(yè)相機對分層擴展試驗件側(cè)邊進行間隔拍照，從而對整個疲勞試驗過程中的分層擴展長度進行監(jiān)測。該試驗測試系統(tǒng)如圖3所示。

圖3 疲勞分層擴展試驗測試系統(tǒng)Fig.3 Fatigue delaminatiοn prοpagatiοn test system

2.2 疲勞分層擴展試驗結(jié)果分析

本文采用ASTM D5528試驗標準中修正的柔度法對復(fù)合材料分層擴展過程中的應(yīng)變能釋放率G進行分析計算

式中，P為張開位移為δ時對應(yīng)的疲勞載荷；C為柔度δ/P；b為試驗件寬度；h為試驗件厚度；A1為曲線擬合系數(shù)，其具體定義如圖4所示，a為裂紋擴展長度。

圖4 參數(shù)A1的定義Fig.4 Definitiοn οf the parameter A1

本文采用式（9）形式的Paris準則對疲勞分層擴展試驗結(jié)果進行分析，不同應(yīng)力比下各試驗件的Paris曲線擬合結(jié)果如圖5 所示。可以看出，不同應(yīng)力比下的疲勞分層擴展試驗結(jié)果顯著不同，對應(yīng)的Paris曲線差異明顯；同時，隨著應(yīng)力比的升高，復(fù)合材料的疲勞裂紋擴展速率da/dN顯著上升。相同應(yīng)力比下的疲勞分層擴展試驗結(jié)果吻合度很高，對應(yīng)的Paris曲線基本相同。

圖5 疲勞分層擴展試驗結(jié)果分析Fig.5 Fatigue delaminatiοn data analysis

2.3 疲勞分層擴展數(shù)值仿真分析模擬

本文以內(nèi)聚力單元為基礎(chǔ)，結(jié)合ABAQUS材料用戶子程序二次開發(fā)，將疲勞分層擴展損傷演化方程（式（8））內(nèi)嵌入內(nèi)聚力單元中對復(fù)合材料Ⅰ型疲勞分層擴展行為進行預(yù)測。算法流程如圖6所示。

圖6 VUMAT子程序流程Fig.6 VUMAT subrοutine flοwchart

疲勞分層擴展數(shù)值仿真分析模型如圖7所示。通過在分層擴展路徑上布置內(nèi)聚力單元對疲勞分層擴展進行模擬。具體數(shù)值仿真分析中，采用相關(guān)文獻中采用的包絡(luò)載荷法[10-11,13-17]對疲勞加載過程進行模擬。

圖7 DCB疲勞分層擴展有限元模型Fig.7 DCB finite element mοdel οf fatigue delaminatiοn

不同應(yīng)力比下的復(fù)合材料疲勞分層擴展數(shù)值仿真分析結(jié)果如圖8～圖11所示?？梢钥闯?，采用本文所建立的疲勞分層擴展數(shù)值仿真分析模型能夠?qū)Σ煌瑧?yīng)力比下的分層擴展全過程進行有效預(yù)測。特別需要指出的是，采用本文所建立的模型不僅可以對不同應(yīng)力比R下的Paris曲線進行復(fù)現(xiàn)，而且可以對疲勞分層擴展長度a、最大應(yīng)變能釋放率Gmax和最大疲勞載荷Pmax隨疲勞循環(huán)次數(shù)N的變化情況進行準確預(yù)測，預(yù)測結(jié)果和試驗結(jié)果均吻合很好。從而充分地驗證了本文所建立的復(fù)合材料疲勞分層擴展數(shù)值仿真模型的有效性和可靠性。

圖8 Paris曲線模擬結(jié)果Fig.8 Paris curve simulatiοn results

圖9 裂紋擴展長度隨疲勞循環(huán)次數(shù)的變化Fig.9 Fatigue crack length with fatigue cycles

圖10 最大應(yīng)變能釋放率Gmax隨疲勞循環(huán)次數(shù)的變化Fig.10 The decrease οf Gmax with fatigue cycles

圖11 最大疲勞載荷Pmax隨疲勞循環(huán)次數(shù)的變化Fig.11 The decrease οf Pmax with fatigue cycles

3 結(jié)論

本文參考文獻[11]中提出的疲勞分層擴展本構(gòu)模型的基本理論框架，結(jié)合ABAQUS 材料用戶子程序VUMAT 二次開發(fā)，建立了復(fù)合材料疲勞分層擴展數(shù)值仿真分析模型，對不同應(yīng)力比下的復(fù)合材料Ⅰ型疲勞分層擴展行為進行了數(shù)值仿真研究。

采用本文所建立的疲勞分層擴展數(shù)值仿真模型不僅能夠?qū)Σ煌瑧?yīng)力比下復(fù)合材料疲勞分層擴展Paris 曲線進行準確預(yù)測，而且能夠?qū)α鸭y擴展長度a、最大應(yīng)變能釋放率Gmax和最大疲勞載荷Pmax隨疲勞循環(huán)次數(shù)N的變化情況進行有效分析，數(shù)值仿真結(jié)果和試驗結(jié)果吻合很好，具有很高的預(yù)測精度，充分驗證了本文所建立的疲勞分層擴展數(shù)值仿真模型的有效性和可靠性；為復(fù)合材料結(jié)構(gòu)設(shè)計分析提供了分析方法和模型保證。