郭雯, 馬玉娥*, Sundar Natarajan, 陳鵬程, 彭帆

（1.西北工業(yè)大學 航空學院，西安 710072；2.Indian Institute of Technology Madras, School of Mechanical Engineering, Chennai 600036；3.長安大學 理學院，西安 710064）

纖維增強復合材料(Fiber-reinforced composites，F(xiàn)RC)由于其優(yōu)異的力學性能及結構的可設計性在航空航天、土木建筑、船舶等工程領域被廣泛應用[1]。但同時，力學行為的各向異性和不均勻性也導致復合材料在斷裂模式上與傳統(tǒng)均質材料有所區(qū)別，即使是單向纖維增強復合材料也會由于載荷(拉伸、壓縮或剪切)大小、方向及作用面的不同而產生不同的失效模式[2-3]，因此研究復合材料的斷裂問題具有一定的挑戰(zhàn)性。

斷裂相場法最早可追溯至上世紀90年代末Francfort等[4]提出的變分斷裂原理，后經Bourdin等[5]完善。斷裂相場法是一種彌散裂紋模型，其通過引入一個序參量來表征材料損傷程度，無需額外的斷裂準則，僅通過相場的演化自動就可以追蹤裂紋擴展路徑，因此在模擬復合材料的復雜斷裂行為時有著明顯優(yōu)勢[6-8]。近年來，斷裂相場法還被廣泛應用于多物理場耦合斷裂[9-11]、疲勞斷裂[12-13]、動態(tài)斷裂[14-15]和超彈性材料斷裂[16]等多種斷裂問題中。

為了解釋纖維增強復合材料在斷裂過程中表現(xiàn)出的各向異性，Clayton等[17]在裂紋面密度函數(shù)中引入了表征材料各向異性的二階結構張量，并多次應用于復合材料斷裂行為的研究中[6]。該結構張量通過一個罰系數(shù)來控制裂紋的演化方向，當罰系數(shù)取較大值時，裂紋只能沿著弱失效方向進行擴展。Natarajan等[18]利用該模型分別對單向和變剛度纖維增強單層板進行了損傷模擬，但該模型并不適用于三維情況。Zhang等[19]引入了一種修正的結構張量，應用到相場模型中對纖維增強復合材料的裂紋擴展路徑和力學響應進行預測，并將其推廣到三維斷裂問題，模擬了纖維增強復合材料層合板的裂紋擴展。Zhang等[20]還進一步結合細觀力學分析了纖維增強復合材料的漸進破壞過程。以上模型均建立在利用帶有罰系數(shù)的結構張量基礎上，但罰系數(shù)的物理含義和取值標準尚不明確，且利用單一相場在模擬復雜裂紋擴展時有裂紋路徑與試驗結果存在偏差的問題[21]。因此，Bleyer等[22]提出了一種各向異性相場模型，該模型區(qū)分了不同的斷裂模式，將標準的單變量相場擴展為多個不同的損傷變量，即利用多個相場變量分別表征纖維和基體的損傷。該研究表明相較于單相場模型，采用多相場模型在模擬各向異性材料裂紋擴展時更加準確。Dean等[23]基于Puck失效理論提出了一種多相場模型，但該模型需要引入額外的失效準則。Zhang等[24]和Pan等[21]則在雙相場模型中綜合考慮了基體的拉伸和剪切臨界能量釋放率，同時分別將應變能分解為纖維主導和基體主導的兩部分，對纖維增強復合材料層板的損傷進行了模擬，結果與試驗均吻合良好。

但是相比于一般的斷裂數(shù)值方法，相場模型在進行斷裂模擬時需要在裂紋路徑處設置較細的網格。由于復合材料的斷裂模式復雜多樣，難以預知裂紋的擴展路徑，因此計算前處理過程會十分復雜，甚至需要對網格全局細化來保證計算精度，這就導致了極高的計算成本[25]。而自適應網格細化的方式可以較好地解決這一問題[26]，Gupta等[27]將自適應有限元方法與斷裂相場模型相結合，應用于脆性斷裂、內聚斷裂和動態(tài)斷裂問題中，與傳統(tǒng)非自適應網格細化方式相比，計算時間縮短了5～50倍，極大地提高了計算效率。Noii等[28]提出了一種高效的全局-局部自適應相場方法，并應用于各向異性脆性斷裂相場模型中，并輔以預測校正，通過計算橫觀各向同性材料的位移-載荷曲線驗證了模型的可行性。Chen等[29]在各向異性脆性斷裂相場中提出了兩種基于等幾何方法的自適應網格細化方案，與全局細化相比，兩種自適應方案均合理地減少了計算工作量。Hirshikesh等[30]采用了基于四叉樹分解的網格自適應策略，結合相場模型對準靜態(tài)加載下的脆性斷裂進行模擬，該方法的優(yōu)勢在于易于實現(xiàn)且能保留分層的網格結構。

目前基于四叉樹結構的自適應算法在復合材料斷裂問題的研究中應用較少。本文結合該方法的簡單高效性和雙相場模型模擬復材斷裂的優(yōu)越性，建立了基于四叉樹分解數(shù)據(jù)結構的自適應網格細化算法，在裂紋擴展路徑上進行局部網格細化。同時采用各向異性雙相場模型分別模擬了不同纖維角度下單向纖維增強復合材料單層板和變剛度纖維增強復合材料單層板的拉伸斷裂。

1 雙相場斷裂模型

1.1 雙相場基本斷裂模型

尖銳裂紋模型和雙相場模型如圖1所示，考慮一個內部存在纖維斷裂面Γf和基體斷裂面Γm的復合材料控制體單元??Rδ,δ∈[1,2,3]，其中Rδ表示維度為δ的空間，本文中僅考慮二維情況，因此δ=2?？刂企w外邊界為???Rδ-1，外邊界的法向量為n。其中，Dirichlet邊界??u和Neumann邊界??t滿足以下條件：??u∪??t=??，??u∩??t=?(?為空集)。相場法通過引入輔助場變量di(x)∈[0,1]，i=f,m及其空間梯度?di來構造裂紋演化狀態(tài)，di=0和di=1分別表示材料未發(fā)生斷裂和完全斷裂狀態(tài)。

圖1 (a) 尖銳裂紋模型；(b) 雙相場模型Fig.1 (a) Sharp crack model; (b) Double phase-field model

根據(jù)Francfort等[4]提出的斷裂變分原理，正則化后的系統(tǒng)總勢能由彈性應變能、斷裂表面能和外力功3部分構成，可表示為

式中：ψe(ε)為彈性應變能密度函數(shù)；Gc,i為材料的臨界能量釋放率；u為位移向量；b為體力；t為面力；γi(di,?di)為裂紋面密度函數(shù)，在標準的各向同性材料的AT2相場模型中[31]，γi(di,?di)為

其中，li為裂紋特征長度。引入各向異性結構張量，將裂紋面密度函數(shù)推廣至纖維增強復合材料，

式中，ωf、ωm分別為表征纖維失效方向和基體失效方向的二階結構張量。本文沿用Pan等[21]提出的結構張量形式

式中，θ為纖維方向與全局坐標系下x軸的夾角。

彈性應變能密度由纖維主導斷裂和基體主導斷裂兩部分組成，

為了區(qū)分應變能的不同分量對裂紋擴展的影響，將主應變張量與應變能進行拉壓分解，應變分解的向量形式為

式中：〈?〉為麥考利符號；γ12為切應變。

裂紋驅動力只由正應變所對應的應變能部分構成，因此彈性應變能密度函數(shù)寫為如下形式：

其中，g(di)=(1-di)2+k為能量退化的經驗公式。k是為避免di=1時數(shù)值奇異所取的附加參數(shù)。一般而言，k?1，在本文中均取k=1×10-6。

考慮到纖維增強復合材料的各向異性，在二維情況下，將纖維和基體主導的不同斷裂模式所對應的應變能密度函數(shù)定義如下：

式中，

其中坐標轉換矩陣T為[32]

1.2 能量耗散不等式

假定在恒溫狀態(tài)且沒有其他外力作用的條件下，控制體的能量耗散總是大于等于0，可得Clausius-Duhem不等式如下：

式中，Yf、Ym、Hf、Hm分別為損傷df、dm及其梯度的共軛廣義力[33]。而自由能密度函數(shù)ψ(ε,df,dm,?df,?dm)的時間導數(shù)為

可將式(15)簡化為

為滿足任意狀態(tài)下的能量耗散不等式，可得

式中，Cg為全局坐標系下材料的彈性張量。

1.3 控制方程

對位移和相場變量作一階變分：

外力功的變分形式為

根據(jù)虛功原理，在任意的δu和δdi下，應滿足

經推導可得本文各向異性雙相場斷裂模型中各變量的控制方程及邊界條件為

2 網格自適應

2.1 網格自適應過程

為保證相場模型的求解精度，一般需要在裂紋擴展路徑上覆蓋極為細化的網格。相比于全局細化，采用自適應算法進行網格局部細化可以避免繁瑣的網格細化前處理過程，同時在達到相同計算精度的前提下降低網格數(shù)量，節(jié)省計算成本，因此在提高計算效率方面具有很大優(yōu)勢。

首先對試件進行粗略的網格劃分，對可能產生裂紋的位置進行局部細化以避免在裂紋萌生的初始階段出現(xiàn)較大誤差。為了識別裂紋面，需要引入相應的準則來確定有限元網格細化的區(qū)域，而本文同時將相場閾值dc及最小網格尺寸s作為網格細化的誤差指示。在每個時間步的初始網格上求解單元的相場變量、位移和歷史最大變量，當單元的相場值達到di>dc且s>li,i=f,m時，該單元將會被標記。結合四叉樹分解的數(shù)據(jù)結構對標記的單元進行細化并編號，這種數(shù)據(jù)結構可以減小存儲空間。具體的細化方式如圖2所示，將每個單元遞歸的劃分為4個子單元。每次細化后，相場等變量都需要進行重新計算以判斷是否繼續(xù)進行細化，直到滿足s≤li/2,i=f,m后，輸出最新的網格及變量值用于下一步計算。在生成新網格的過程中，會不斷產生新的節(jié)點，而在原始未細化和細化單元的過渡區(qū)域，可能出現(xiàn)節(jié)點不適配的單元，即帶有懸點的單元。未被細化的單元因懸點的存在不再是傳統(tǒng)的四節(jié)點單元，簡單的雙線性形函數(shù)也不再適用。本文采用多邊形單元解決上述問題，即將帶有懸點的單元視為n邊多邊形，以Gauss積分點為極點，連接極點與各個節(jié)點和懸點，將單元劃分為i個三角形，節(jié)點坐標均采用極坐標表示。利用Floater[34]提出的二維凸多邊形均值坐標理論定義此類單元內部的形函數(shù)為

圖2 四叉樹網格及帶懸點單元的形函數(shù)構造Fig.2 Quadtree mesh and the construction of shape function of elements with hanging nodes

其中：wi(x)為各個節(jié)點的權重系數(shù)；αi為內角；xi為第i個節(jié)點的極坐標。

2.2 相場方程的數(shù)值離散

在構造好自適應多邊形網格形函數(shù)后，再對位移場、相場變量及其一階變分進行數(shù)值離散：

式中：和Ni分別為位移場和相場的形狀函數(shù)矩陣；和Bi為相應形函數(shù)的一階偏導，其形式如下：

對勢能方程的弱形式進行離散化，可得位移場及兩個相場的殘差如下：

采用分步式算法，利用Newton-Raphson方法交錯迭代求解位移u和相場變量df、dm：

相應的剛度矩陣通過對殘差求偏導可得：

3 數(shù)值算例

為驗證各向異性雙相場模型模擬纖維增強復合材料斷裂問題的可行性，研究了含單邊裂紋及中心含孔的纖維增強復合材料單層板在拉伸位移載荷下的裂紋擴展，并將計算結果與試驗進行了對比。利用Matlab開發(fā)了自適應有限元相場求解程序，分別對構件進行全局網格細化和自適應網格細化，通過對比兩種不同細化方式下的計算結果、網格數(shù)和計算時間來驗證自適應算法的高效性。另外，研究了自適應雙相場模型在模擬周期性變剛度纖維增強復合材料單層板斷裂問題的可行性，討論了不同的纖維路徑設計參數(shù)對結構等效強度的影響。

3.1 單邊裂紋FRC單層板拉伸試驗

含單邊裂紋的FRC單層板模型的幾何尺寸與邊界條件如圖3所示，材料參數(shù)見表1。取s=0.02mm的網格尺寸劃分初始網格，設置網格尺寸閾值sc=0.005mm，位移增量步取?u=1×10-4mm，纖維相場和基體相場的特征長度取lf=lm=0.01mm。

表1 單邊裂紋纖維增強復合材料試樣材料參數(shù)[35]Table 1 Material parameters of single edge notched fiber reinforced composite specimen[35]

圖3 單邊裂紋纖維增強復合材料(FRC)試樣：幾何尺寸及邊界條件Fig.3 Single edge notched fiber-reinforced composites (FRC) specimen:Geometric size and boundary conditions

為研究相場閾值dc對計算結果的影響，以纖維角度θ=30o為例，分別取相場閾值dc=0.1、0.3、0.5、0.7、0.9進行計算，計算得到的位移-載荷曲線和最終網格數(shù)如圖4和表2所示。

表2 不同相場閾值對應的網格數(shù)Table 2 Number of elements with different phase-field thresholds

圖4 不同相場閾值下單邊裂紋FRC試樣的位移-載荷曲線Fig.4 Displacement-load curves of single edge notched FRC specimen with different phase-field thresholds

由圖4可知，相場閾值取0.1～0.5時，閾值大小對計算結果的影響較??；閾值取0.7～0.9時，計算結果偏差較大且極不穩(wěn)定。根據(jù)表2，相場閾值取0.1～0.5時最終網格數(shù)隨閾值的增大而減少。為節(jié)約計算成本，本文中自適應網格細化的相場閾值均取dc=0.5。

對纖維角度θ為30°、45°、60°、90°的FRC單層板施加拉伸位移載荷，計算和試驗所得的裂紋路徑及對應的自適應網格分別如圖5所示。在本文算例中，由于纖維損傷較小，因此只給出了基體損傷變量的云圖。根據(jù)試驗所得的真實裂紋路徑可知，對于在圖3所示加載條件下的單向纖維增強復合材料單層板，基體裂紋會沿平行于纖維角度的方向進行擴展。對比圖5(a)、5(b)可見，雙相場模型計算得到的裂紋路徑與試驗結果相同。圖5(c)中自適應網格細化的區(qū)域也與裂紋擴展路徑保持一致。

圖5 單邊裂紋FRC試樣沿不同纖維方向下的裂紋擴展：(a)計算結果；(b)試驗結果[35]；(c)自適應網格Fig.5 Crack propagation paths of single edge notched FRC specimen with different fiber orientation: (a) Simulation results;(b) Experimental results[35]; (c) Adaptive elements

同時，為驗證自適應算法的高效性，對試件以0.005 mm的網格尺寸進行全局細化，并施加相同的拉伸位移載荷。計算結果如圖6和表3所示。

表3 計算時間和網格數(shù)量對比Table 3 Comparison of computation time and number of elements

圖6 不同纖維方向的單邊裂紋FRC試樣位移-載荷曲線Fig.6 Displacement-load curves of single edge notched FRC specimen with different fiber directions

從圖6可以看出，不同纖維角度下，兩種網格細化方式所得到的位移-載荷曲線基本保持一致。表3給出了兩種網格細化方式下的計算時間和最終的網格數(shù)量，經過計算可知自適應網格細化比全局網格細化的計算時間平均減少了56.17%，網格數(shù)量平均減少了89.47%。

以上結論證明了自適應雙相場模型用于模擬纖維增強復合材料斷裂行為的可行性和高效性。在準確模擬裂紋擴展路徑的同時，該模型可以減少計算所需的網格數(shù)量，縮短計算時間，提高計算效率。

3.2 中心含孔FRC板拉伸試驗

研究不同纖維角度下中心含孔的FRC單層板在拉伸載荷下的斷裂問題。模型的幾何尺寸及邊界條件如圖7(a)所示，材料參數(shù)在表4中給出。對試件分別進行自適應網格細化和全局細化，其中自適應細化的初始網格尺寸取0.4 mm，網格尺寸閾值取sc=0.2mm，相場閾值取dc=0.5。全局細化的網格尺寸取0.2 mm。位移增量步為?u=5×10-3mm，基體和纖維相場的特征長度分別為lm=lf=0.4mm。根據(jù)圖7(b)給出的部分計算結果可知，對于圖7(a)所示的中心含孔FRC單層板，試件受拉時基體裂紋會從孔邊萌生，并同樣沿著纖維角度方向進行擴展。

表4 中心含孔FRC試樣材料參數(shù)[36]Table 4 Material parameters of open hole FRC specimen[36]

圖7 中心含孔FRC試樣：(a)幾何尺寸及邊界條件；(b)部分計算結果Fig.7 Open hole FRC specimen: (a) Geometric size and boundary conditions; (b) Some calculation results

將兩種網格細化策略下計算所得的結構等效強度與試驗[36]進行對比，曲線如圖8所示。由圖可知，在兩種網格細化方式下，計算結果均與試驗結果吻合良好。這表明自適應細化策略可以達到與全局細化相同的計算效果。

圖8 不同纖維角度對應下中心含孔FRC試樣的結構等效強度Fig.8 Effective strengths of the central circle FRC specimen with different fiber orientations

圖9對比了計算結果基本一致時采用兩種網格細化策略所需要的最終網格數(shù)量。自適應網格數(shù)量平均約為全局細化網格數(shù)量的70%，有效地降低了計算成本。

圖9 中心含孔FRC試樣的自適應和非自適應的網格數(shù)量對比Fig.9 Comparison of the number of elements between adaptive and unadaptive schemes for central circle FRC specimen

3.3 變剛度FRC板拉伸試驗

變剛度纖維增強復合材料的特征是其纖維角度可以根據(jù)實際應用需求進行設計，與傳統(tǒng)單向纖維增強復合材料相比有力學性能可控的優(yōu)勢。一般變剛度FRC單層板的纖維角度變化函數(shù)選取周期性函數(shù)，具體函數(shù)形式見下式：

式中：以圖10中做出標注的纖維為例，φ為纖維路徑繞幾何坐標系的旋轉角度；纖維在試件中心處與x軸的夾角定義為初始角度T0；T1為距離中心a/2處的纖維終止角度；a表示一個周期的特征尺寸；i為周期數(shù)。纖維路徑函數(shù)可表示為[φ〈T0|T1〉]。由于θ是x的函數(shù)，因此同一單元內的不同節(jié)點所對應的材料參數(shù)不同。本文在高斯積分點處對同一單元內各節(jié)點的纖維角度進行插值，從而保證一個單元只對應一個纖維角度。

圖10 變剛度FRC試樣的幾何尺寸及邊界條件Fig.10 Geometric size and boundary conditions of variable stiffness FRC specimen

首先研究自適應雙相場模型在模擬變剛度纖維增強復合材料單層板斷裂問題時的可行性。試件的幾何尺寸及邊界條件在圖10中給出，材料參數(shù)與表1相同。劃分的初始網格尺寸為0.2 mm，自適應網格尺寸閾值為dc=0.5，位移增量步取?u=1×10-4mm，纖維相場和基體相場的特征長度分別取lf=lm=0.01mm。

定義特征尺寸a=0.4mm ，纖維路徑為[φ〈-30o|30o〉]。分別計算旋轉角度φ為0o、30o、45o、60o的變剛度FRC單層板在拉伸位移載荷下的斷裂問題。計算得到的裂紋擴展路徑和部分自適應細化網格如圖11所示。

圖11 變剛度FRC試件部分模擬結果：((a)～(c)), ((e)～(g))裂紋擴展路徑；((d), (h))自適應網格Fig.11 Some simulation results of variable stiffness FRC specimen：((a)-(c)), ((e)-(g)) Crack propagation paths; ((d),(h)) Adaptive elements

由圖可以看出，對于拉伸載荷下的變剛度FRC單層板，盡管纖維角度連續(xù)變化，纖維路徑較曲折，基體裂紋也會沿著平行于纖維角度的方向進行擴展。這與前兩個算例所得到的結論相同。同時，自適應方法仍可以捕捉不規(guī)則的裂紋擴展路徑來進行網格細化。對于與變剛度復合材料類似的非均質材料，裂紋擴展路徑難以預測，因此網格細化的前處理過程將會十分復雜，而自適應方法只需對試樣劃分較粗糙的初始網格，有效地簡化了前處理過程，提高工作效率。

其次討論部分纖維路徑參數(shù)對于結構等效強度的影響。設置纖維鋪設路徑為[0〈-30o|30o〉]，計算在不同特征尺寸a下的結構等效強度。計算結果如圖12所示。

圖12 變剛度FRC試樣等效強度對比：不同旋轉角度和不同特征尺寸Fig.12 Comparison of the effective strength of variable stiffness FRC specimen: Different rotation angles and different reference length

圖12采用了雙坐標軸的形式，三角形符號的曲線表示不同旋轉角度下的等效強度(對應于左、下坐標軸)，而圓形符號的曲線則表示不同特征尺寸下的等效強度(對應于右、上坐標軸)。可以看出，對于圖10所示的變剛度纖維增強復合材料單層板，試件受拉時纖維路徑的設計參數(shù)均對結構等效強度存在一定的影響。其中旋轉角度φ產生的影響較大，而特征尺寸a的影響并不顯著。

4 結 論

采用基于四叉樹的網格自適應策略，利用Matlab開發(fā)了網格自適應程序，結合各向異性雙相場模型對纖維增強復合材料單層板的拉伸斷裂行為進行了研究。通過對比數(shù)值計算結果和試驗結果，分析不同網格細化方式下的計算時間和網格數(shù)量，驗證了自適應雙相場模型的可行性和高效性。得到以下結論：

(1) 自適應雙相場模型可以準確地預測纖維增強復合材料單層板的裂紋擴展路徑。在拉伸載荷下，對于纖維體積分數(shù)較小的單向和變剛度纖維增強復合材料單層板，相場裂紋均沿著平行于纖維角度的方向擴展，且纖維鋪設角度對結構等效強度存在一定的影響；

(2) 自適應網格細化策略能夠準確捕捉裂紋路徑進行局部細化。在不影響計算精度的情況下，自適應算法可以減少70%～90%的網格數(shù)量，有效提高計算效率，節(jié)約計算成本；

(3) 用于確定自適應網格細化區(qū)域的相場閾值對計算精度有一定影響。閾值過大會導致求解不穩(wěn)定，閾值越小計算結果越精確。