基于擴(kuò)展有限元法的重組竹Ⅰ型斷裂數(shù)值模擬分析

劉明凱，周愛萍，劉燕燕，盛寶璐

(南京林業(yè)大學(xué)生物質(zhì)材料國家地方聯(lián)合工程研究中心，南京 210037)

重組竹是竹纖維束通過順紋組胚、膠合熱壓制成的材料，其力學(xué)性能優(yōu)于常用結(jié)構(gòu)木材。但由于竹束疏解工藝簡單，難以均勻浸膠，重組竹不可避免會帶裂縫工作。當(dāng)復(fù)合材料存在基體裂紋、孔洞和缺口時，層間斷裂是最常見的失效形式。這種層間分層會導(dǎo)致纖維斷裂，降低復(fù)合材料的使用壽命。目前，重組竹處于基礎(chǔ)性研究階段，由于完整標(biāo)準(zhǔn)體系和設(shè)計計算理論的欠缺，該材料還未大規(guī)模投入實際工程應(yīng)用中。重組竹構(gòu)件在工程應(yīng)用時，體內(nèi)裂紋一般以復(fù)合型裂紋的形式出現(xiàn)，應(yīng)力狀態(tài)復(fù)雜。單一的Ⅰ/Ⅱ型裂紋是研究復(fù)合型裂紋的前提，其中，Ⅰ型層間斷裂一直是研究人員最為關(guān)注的問題，因為與剪切模式相比，拉伸分層所需的能量較低，裂紋更易在該模式下萌生。為防止斷裂失效的發(fā)生，有必要對復(fù)合材料Ⅰ型裂紋的啟裂和擴(kuò)展規(guī)律進(jìn)行研究，這對工程設(shè)計和應(yīng)用指導(dǎo)具有重要意義。

除解析法和試驗法是研究斷裂問題的重要方法，數(shù)值模擬法也是研究裂紋擴(kuò)展問題的有效手段。隨著計算機(jī)技術(shù)的迅猛發(fā)展，運算能力極大提高，斷裂力學(xué)在有限元軟件方面的應(yīng)用更加廣泛。傳統(tǒng)有限元法存在局限性，預(yù)制裂紋面要與劃分的網(wǎng)格重合，并且均要預(yù)先確定裂紋擴(kuò)展的方向。擴(kuò)展有限元法(extended finite element method，XFEM)克服了上述缺陷，給裂紋擴(kuò)展在實際應(yīng)用中帶來了便利。虛擬裂紋閉合技術(shù)(VCCT)和內(nèi)聚區(qū)模型(CZM)是模擬裂紋擴(kuò)展的2種常用方法。VCCT基于線彈性理論對裂縫進(jìn)行分析，忽略了塑性區(qū)能量耗散內(nèi)聚力。CZM的軟化階段有效彌補(bǔ)了VCCT中缺少的塑性效應(yīng)，該方法已被證明是模擬雙材料界面[1-2]和層壓復(fù)合材料[3-4]斷裂的有力工具。Mo?s等[5]最初提出了XFEM與CZM相結(jié)合的方法，得到了與網(wǎng)格無關(guān)的裂紋表示和基于CZM的裂紋擴(kuò)展。Heidari-Rarani等[6]將XFEM與VCCT、CZM相結(jié)合對復(fù)合材料進(jìn)行模擬對比，與XFEM-VCCT相比，XFEM-CZM提供了更準(zhǔn)確的結(jié)果，可以很好地預(yù)測裂紋的萌生和擴(kuò)展。

目前，有關(guān)重組竹的斷裂研究主要集中在斷裂破壞機(jī)理分析和斷裂韌度測定，而重組竹的斷裂仿真模擬以及斷裂韌性與試件尺寸相關(guān)的研究還很少。筆者以ABAQUS有限元軟件為平臺，基于擴(kuò)展有限元法對不同尺寸的重組竹雙懸臂試件進(jìn)行裂紋擴(kuò)展數(shù)值模擬，將模擬結(jié)果和試驗結(jié)果進(jìn)行對比來驗證擴(kuò)展有限元法的準(zhǔn)確性，該內(nèi)聚模型可用于建立構(gòu)件計算模型和確定構(gòu)件承載力；再通過荷載-位移曲線得到裂紋擴(kuò)展阻力曲線(R曲線)，并依此進(jìn)行重組竹斷裂韌度的尺寸效應(yīng)分析。

1 擴(kuò)展有限元法和內(nèi)聚區(qū)模型

擴(kuò)展有限元法和常規(guī)有限元法的區(qū)別在于，擴(kuò)展有限元法在單位分解法的基礎(chǔ)上增加了階躍函數(shù)和漸進(jìn)場函數(shù)來對位移進(jìn)行具體描述，實現(xiàn)了裂紋和網(wǎng)格獨立存在。因此，在模擬裂紋擴(kuò)展時無須重新劃分網(wǎng)格，能方便地分析不連續(xù)問題，特別是不連續(xù)邊界演化問題。求解域上的近似位移場[uh(x)]為：

(1)

在黏性裂紋中，裂紋的擴(kuò)展受尖端附近裂紋面上牽引位移關(guān)系的控制。該模型是由Dugdale[7]和Barenblatt[8]引入金屬領(lǐng)域的，后來演化為內(nèi)聚區(qū)模型，能夠消除裂紋尖端的應(yīng)力奇異性，可將裂紋尖端處的應(yīng)力幅度限制在有物理意義的水平。

圖1 內(nèi)聚區(qū)模型Fig. 1 Cohesive zone model

2 試驗概述和仿真模擬

2.1 試驗概述

斷裂韌性測定試驗對試件的形狀和尺寸有嚴(yán)格的要求，雙懸臂梁(DCB)是測試復(fù)合材料斷裂性能最常用的結(jié)構(gòu)形式。梁高1/2處的預(yù)制裂紋將整梁劃分為上下兩部分，通過拉伸試驗研究材料的界面強(qiáng)度。Ⅰ型層間韌性能量釋放率試驗基于ASTM D5528-13“Standard test method for Mode I interlaminar fracture toughness of unidirectional fiber-reinforced polymer matrix composites” 測定，試件尺寸如圖2所示。

a0為初始裂紋長度；a為裂紋長度；l為試件長度；b為試件厚度；h為試件1/2高度；P為荷載。圖2 DCB測試試件尺寸Fig. 2 DCB test specimen dimension

分析不同厚度、不同初始裂紋長(同向)的幾何相似試樣，確定適當(dāng)?shù)暮穸群统跏剂鸭y長，從而評估有效的斷裂能，試件具體尺寸見表1。其中，厚度的取值是滿足平面應(yīng)變的條件之一；裂紋長度的取值是滿足小范圍屈服(塑性區(qū)尺寸遠(yuǎn)小于裂紋尺寸)的條件，只有小范圍屈服，才能按線彈性計算出滿足工程精度要求的雙懸臂梁能量釋放率(GIC)。本試驗分為A、B兩組，其中，A組試件初始裂紋長度相同，厚度不同；B組試件厚度相同，初始裂紋長度不同。試驗共7組，其中每組至少包含3個試件。準(zhǔn)靜態(tài)DCB試驗在萬能試驗機(jī)上進(jìn)行，采用0.5 mm/min 恒定速率加載，低速率可以避免裂紋的不穩(wěn)定傳播。在裂紋擴(kuò)展過程中，載荷-位移(P-δ)曲線是唯一需要記錄的數(shù)據(jù)，是獲得R曲線的基礎(chǔ)。最后，將R曲線的GIC與模擬得到的GIC進(jìn)行比較，評估該方法的可行性。

表1 試件尺寸Table 1 Specimen sizes mm

2.2 仿真模擬

模型是仿真的基礎(chǔ)，模型的還原度直接決定了分析結(jié)果的可靠度?；谟邢拊浖嗀BAQUS，建立帶預(yù)制裂紋的雙懸臂模型。在對模型進(jìn)行荷載設(shè)定時，為消除耦合帶來的局部扭轉(zhuǎn)效應(yīng)，建立剛體并將其貫穿于雙懸臂預(yù)制的銷孔內(nèi)，設(shè)置2個剛體的相反位移。

網(wǎng)格劃分時，在保證計算精度的前提下為減少相對所需的計算時間，采用網(wǎng)格局部加密。劃分部件并在裂紋擴(kuò)展區(qū)附近細(xì)化局部網(wǎng)格，得到精確的模擬結(jié)果(高度方向每1 mm劃分1次)；而在遠(yuǎn)離裂紋的區(qū)域，由于對裂紋尖端附近的應(yīng)力狀態(tài)影響較小，使用相對較大的網(wǎng)格(高度方向每10 mm劃分1次)，并設(shè)置順紋方向為坐標(biāo)軸的1方向，其余2個方向為2和3方向，如圖3所示。進(jìn)一步的網(wǎng)格加密研究表明，更細(xì)的網(wǎng)格劃分對結(jié)果精度影響較小。

圖3 有限元模型網(wǎng)格劃分Fig. 3 Mesh generation of finite element model

在材料模塊中定義重組竹的材料特性(E為彈性模量，E1=11 212 MPa、E2=E3=2 561 MPa；ν為泊松比，ν12=ν13=0.304、ν23=0.054；G為剪切模量，G12=G13=1 418 MPa、G23=749 MPa)。重組竹為纖維定向增強(qiáng)復(fù)合材料，順紋與橫紋的力學(xué)性能差異較大，順紋方向(1方向)的強(qiáng)度遠(yuǎn)高于另外2個方向(2和3方向)。實體部件使用C3D8R單元，采用靜態(tài)分析并打開幾何大變形。在XFEM裂紋模塊中，定義裂紋面所處位置并允許裂紋沿順紋方向擴(kuò)展；“場輸出”勾選PHILSM、PSILSM可觀察裂紋的擴(kuò)展?fàn)顟B(tài)，勾選STATUSXFEM可在結(jié)果中查看擴(kuò)展有限單元的失效狀態(tài)。

損傷起始準(zhǔn)則采用最大名義應(yīng)力準(zhǔn)則，當(dāng)分層過程中最大接觸應(yīng)力比值達(dá)到1時，材料開始發(fā)生損傷，該準(zhǔn)則可表示為：

(2)

損傷演化準(zhǔn)則為基于位移的、線性軟化的損傷演化規(guī)律，一旦滿足損傷起始準(zhǔn)則，材料將根據(jù)損傷演化規(guī)律發(fā)生損傷，該規(guī)律描述了內(nèi)聚區(qū)模型剛度退化的速率。令內(nèi)聚剛度損傷因子為D，在起始階段D=0，表示材料沒有發(fā)生損傷，仍處于彈性狀態(tài)。損傷開始后，D逐漸增加，直到D=1時材料完全失效。線性軟化的內(nèi)聚剛度損傷因子可表示為：

(3)

3 模擬結(jié)果及分析

3.1 結(jié)果對比

模擬結(jié)果的應(yīng)力分布如圖4所示。隨著加載點相對位移的增大，裂紋沿著纖維方向向后擴(kuò)展，與試驗觀察到的現(xiàn)象一致。通過對模擬數(shù)據(jù)的后處理獲得P-δ曲線，進(jìn)而獲得R曲線，將模擬曲線和試驗曲線進(jìn)行對比。

圖4 雙懸臂梁仿真模擬應(yīng)力圖Fig. 4 Simulated stress diagram of double cantilever beam

3.1.1 荷載-位移曲線對比與分析

數(shù)值模擬得到的P-δ曲線與試驗曲線吻合良好，如圖5和6所示。峰值荷載的誤差也在容許范圍內(nèi)，具體誤差值見表2。

對于彈性階段，由A組曲線對比可知，試件厚度較大時，模擬曲線中彈性階段的直線斜率和試驗斜率吻合度較好。其中，b=100 mm的試件(B-1)模擬曲線中的直線段斜率和峰值荷載均與試驗曲線吻合較好。相反，較薄的試件承載能力模擬值略低于試驗值，模擬曲線彈性階段斜率也小于試驗曲線的斜率。

最薄試件荷載差異(絕對值)為12.39%，其他厚度試件模擬差異(絕對值)為1.37%～8.40%。該結(jié)果與Manshadi等[13]對雙懸臂梁模擬的結(jié)果一致，認(rèn)為該差異是由于薄試件缺乏大規(guī)模纖維橋連而導(dǎo)致的。模擬曲線的下降段均與試驗曲線較吻合，體現(xiàn)出ABAQUS中的損傷演化準(zhǔn)則有效地模擬了材料軟化階段。

B組曲線進(jìn)一步驗證了厚試件模擬結(jié)果的準(zhǔn)確性，3組曲線彈性階段均吻合，不同的初始裂紋長度對初始剛度的影響較小。初始裂紋越長，雙懸臂梁剛度越小，能承受的荷載越小。

綜上所述，ABAQUS在對重組竹仿真模擬時，薄試件DCB模擬結(jié)果的準(zhǔn)確度不及厚試件，雙懸臂梁厚度越大，模擬結(jié)果越吻合。在實際應(yīng)用中，為節(jié)省財力和物力，厚度大于某一值時，任意尺寸雙懸臂梁可以通過該仿真模擬方法得到其可靠剛度。

圖5 A組試件荷載-位移曲線對比Fig. 5 Comparison of load-displacement curves of specimens in Group A

圖6 B組試件荷載-位移曲線對比Fig. 6 Comparison of load-displacement curves of specimens in Group B

表2 雙懸臂梁的峰值荷載誤差Table 2 Peak load errors of double cantilever beams

3.1.2 模型準(zhǔn)確性驗證

以裂紋擴(kuò)展量(Δa)為函數(shù)的裂紋擴(kuò)展阻力曲線稱為R曲線，可由荷載-位移曲線得到。R曲線是一種幾何相關(guān)性質(zhì)[14-15]，即初始裂紋長度、厚度、高度、曲率等幾何參數(shù)均可能影響R曲線。本研究對初始裂紋長度和厚度2種參數(shù)進(jìn)行分析。

對于含有穿透裂紋的DCB試件，本研究采用柔度法計算應(yīng)變能釋放率。在線彈性斷裂力學(xué)中，柔度法是測定材料斷裂韌度的方法之一。

由經(jīng)典梁理論，加載點的位移δ，即雙懸臂梁的撓度可表示為：

(4)

式中：ω為單懸臂的撓度；E為彈性模量，此處使用順紋方向的彈性模量E1。則柔度C可表示為：

(5)

考慮DCB為線彈性的情況下，Irwin[16]給出了能量釋放率：

(6)

式中，柔度的立方根(C1/3)作為a的函數(shù)。

在實際應(yīng)用中，由于雙懸臂梁端并非處于完全固定邊界條件下，且裂紋前緣可能因旋轉(zhuǎn)偏離中間平面，該方法計算得到的G偏大，因此，須對G進(jìn)行修正。糾正裂紋尖端旋轉(zhuǎn)的一種方法是引入校正參數(shù)Δ，即雙懸臂梁的有效裂紋長度為a+Δ。對于半厚度的DCB試樣，Williams[17]提出了Δ的計算公式：

(7)

以a+Δ代替式(6)中的a得到修正后的G。A組和B組試件R曲線對比分別如圖7和8所示。由圖7和8可知，修正后的R曲線結(jié)果較準(zhǔn)確，a0符合實際尺寸，且模擬曲線的水平段和試驗曲線相近。

圖7 A組試件R曲線對比Fig. 7 Comparison of R curves of specimens in Group A

圖8 B組試件R曲線對比Fig. 8 Comparison of R curves of specimens in Group B

R曲線的試驗結(jié)果與模擬結(jié)果在趨勢上存在差異，試驗結(jié)果在轉(zhuǎn)折后呈上升趨勢，而模擬結(jié)果在轉(zhuǎn)折后呈下降趨勢。該偏差可從兩方面進(jìn)行分析：一方面，仿真模擬時采用理想化模型的裂紋面為理想平面，而實際復(fù)合材料的非均質(zhì)特性導(dǎo)致裂紋面呈現(xiàn)曲面，纖維橋連現(xiàn)象的增多使得斷裂所需的能量相應(yīng)增多，因此試驗曲線呈上升趨勢，仿真模型和實際試驗存在差異；另一方面，牽引參數(shù)通常被認(rèn)為是材料參數(shù)，文獻(xiàn)[12]中的內(nèi)聚參數(shù)是基于一定尺寸構(gòu)件通過試驗反推得到的，跟尺寸有一定的關(guān)聯(lián)。該文獻(xiàn)中采用的構(gòu)件厚度為80和100 mm，初始裂紋長度為100和134 mm，本研究使用該內(nèi)聚參數(shù)得到的模擬結(jié)果中，厚度為100 mm且初始裂紋長度為102和136 mm的DCB模擬得到的R曲線吻合度較高，R曲線趨于水平。對于其他尺寸的試件，模擬結(jié)果趨勢存在一定偏差。

3.2 尺寸效應(yīng)

3.2.1 厚 度

匯總不同厚度試件的臨界能量釋放率，對比試驗值和模擬值，如圖9所示。當(dāng)b≥40 mm時，臨界能量釋放率均保持在2.0 J/m2左右，即該厚度范圍均可看成平面應(yīng)變問題；但當(dāng)厚度較小且b=20 mm時，試驗測得的GIC高于模擬得到的GIC，較薄的試樣具有較大的斷裂韌度。材料斷裂的臨界能量釋放率隨厚度的增加而逐漸減小，最終趨于一個恒定的較低值。這一結(jié)果驗證了斷裂韌度隨試樣厚度的變化關(guān)系，如圖10所示。

圖9 不同厚度試件的臨界能量釋放率Fig. 9 Critical energy release rate of specimens with different thicknesses

圖10 斷裂韌度與厚度關(guān)系圖Fig. 10 Relationship between fracture toughness and thickness

Hu等[18]認(rèn)為與尺寸效應(yīng)相關(guān)的準(zhǔn)脆性斷裂行為實際上是由于裂紋尖端處的FPZ與最近的結(jié)構(gòu)邊界相互作用的結(jié)果，斷裂過程區(qū)(如橋連區(qū))到最近邊界的距離控制著斷裂行為。FPZ在小試件中相對較大，因此，斷裂以強(qiáng)度準(zhǔn)則為主。在大試樣中，F(xiàn)PZ相對于試樣尺寸及其到最近試樣邊界的距離較小(FPZ的絕對尺寸沒有意義)，因此，斷裂以斷裂韌性準(zhǔn)則為主。

實際上，當(dāng)構(gòu)件厚度與塑性區(qū)尺寸的比值較小時，厚度方向的約束較弱，材料的塑性變形在該方向上不受限制，厚度方向可以自由地發(fā)生屈服，易產(chǎn)生45°斜斷口；當(dāng)構(gòu)件厚度與塑性區(qū)尺寸的比值較大時，厚度方向不能自由產(chǎn)生屈服，周圍彈性材料的限制使厚度方向應(yīng)變?yōu)?，此時裂紋呈平斷口。對于較厚的構(gòu)件，兩自由表面在厚度方向的應(yīng)力為0而處于平面應(yīng)力狀態(tài)，但構(gòu)件內(nèi)部絕大部分處于平面應(yīng)變狀態(tài)，如圖11所示。厚度越大，平面應(yīng)力區(qū)域所占比例越小，因而可將構(gòu)件整體視為完全處于平面應(yīng)變區(qū)域。有限元分析結(jié)果中，雙懸臂梁的剖面應(yīng)力圖驗證了上述關(guān)于厚度的尺寸效應(yīng)，如圖12所示。該結(jié)果和陳濤等[19]研究鋁合金材料R曲線的三維效應(yīng)類似，薄試件有較大的斷裂韌性，隨著厚度的增加，阻力曲線趨于一條水平直線。

圖11 裂紋尖端三維塑性區(qū)Fig. 11 3D plastic zone at the crack tip

圖12 雙懸臂梁模擬剖面應(yīng)力圖Fig. 12 Stress diagram of simulated section of double cantilever beam

3.2.2 初始裂紋長度

ASTM D5528-13建議裂紋分層前沿應(yīng)距加載線50 mm，而ESIS TC4 Protocol “Determination of the mode I delamination resistance of unidirectional fibre-reinforced polymer laminates using the double cantilever beam specimen”建議從加載塊(或鋼鉸鏈)的前邊緣到薄膜前沿的距離應(yīng)至少為45 mm。

隨著初始裂紋的增長，裂紋啟裂所需要的能量減小。初始裂紋較短時，橫向剪切變形占主導(dǎo)地位，裂紋擴(kuò)展所需的能量較大；初始裂紋較長時，雙懸臂梁會產(chǎn)生大撓度變形，裂紋前端接近DCB試件的末端，臨界能量釋放率較小。此時為消除端部效應(yīng)，應(yīng)制備更長的試件。不同初始裂紋長度試件的臨界能量釋放率如圖13所示，當(dāng)試件處于平面應(yīng)變條件時，初始裂紋長度對模擬結(jié)果的吻合度影響較小。

圖13 不同初始裂紋長度試件的臨界能量釋放率Fig. 13 Critical energy release rate of specimens with different initial crack lengths

4 結(jié) 論

本研究運用ABAQUS軟件對重組竹的裂紋擴(kuò)展進(jìn)行了三維仿真模擬，得到了和試驗吻合度較好的模擬結(jié)果，主要結(jié)論如下：

1)本研究得到了基于內(nèi)聚區(qū)模型的擴(kuò)展有限元分析模型，在可靠內(nèi)聚參數(shù)設(shè)定下，模擬裂紋擴(kuò)展過程，得到較吻合的荷載-位移曲線和裂紋擴(kuò)展阻力曲線。

2)試件厚度較大時，P-δ曲線吻合度較好；修正后的R曲線較準(zhǔn)確，雙懸臂梁的能量釋放率隨厚度的增加先下降后趨于平穩(wěn)，并且當(dāng)b≥40 mm時趨于穩(wěn)定。因此，厚度達(dá)到40 mm后，模擬結(jié)果和試驗結(jié)果吻合。

3)初始裂紋長度越大，雙懸臂梁承載能力越低，能量釋放率越小，且初始裂紋長度對模擬的吻合度影響較小。