基于逐層理論的復(fù)合材料層合結(jié)構(gòu)分層損傷擴展研究

李頂河，曹江濤，郭巧榮

中國民航大學(xué)，天津 300300

廣泛應(yīng)用于航空、航天等領(lǐng)域的復(fù)合材料層合結(jié)構(gòu)存在多種損傷形式，最主要、最常見的損傷形式是分層損傷，它的存在嚴(yán)重影響了復(fù)合材料在服役過程中的強度和剛度。利用試驗方法對復(fù)合材料結(jié)構(gòu)進(jìn)行研究和分析，將消耗大量的時間、人力和物力成本?；跀?shù)值模擬技術(shù)對分層損傷進(jìn)行分析將在很大程度上克服這些問題，成為研究分層損傷問題的主要手段之一[1-4]。

針對復(fù)合材料厚板結(jié)構(gòu)，J. N. Reddy[5]提出了逐層理論，該理論能給出精確的層間應(yīng)力應(yīng)變，也能描述厚板面內(nèi)位移在厚度方向上的Zig-Zag 特征。近年來，國內(nèi)外許多學(xué)者對Reddy 的逐層理論進(jìn)行了優(yōu)化和改進(jìn)，并將此理論廣泛應(yīng)用于復(fù)合材料層合結(jié)構(gòu)、功能梯度層合結(jié)構(gòu)、壓電層合結(jié)構(gòu)等問題的分析中[6]。K.M.Liew 等[7]總結(jié)回顧了逐層理論的最新發(fā)展。M.S.Beg等[8]基于逐層理論提出了一種用于分析弧形彎曲梁問題的梁理論，預(yù)測了功能梯度弧形彎曲梁在靜態(tài)、自由和受迫振動情況下的精確響應(yīng)。王瑞鵬[9]將逐層理論與三維實體元結(jié)合，對夾芯結(jié)構(gòu)進(jìn)行了靜力和自由振動分析，表明逐層理論與三維實體元法進(jìn)行耦合能用于分析復(fù)雜點陣夾芯結(jié)構(gòu)的靜動力問題。

虛擬裂紋閉合技術(shù)被廣泛應(yīng)用于計算應(yīng)變能量釋放率和應(yīng)力強度因子。余毅等[10]基于虛擬裂紋閉合法計算應(yīng)力強度因子，模擬了三維裂紋受拉伸載荷作用下的裂紋擴展情況。魯龍坤等[11]解釋了裂紋尖端張開角的定義，并驗證裂紋尖端張開角是一個有效的斷裂參數(shù)。D.H.Li[12]將擴展逐層理論與虛擬裂紋閉合技術(shù)相結(jié)合，用以預(yù)測復(fù)合材料層合板殼任意形狀分層和橫向裂紋的同時擴展。J.Jokinen 等[13]等利用虛擬裂紋閉合技術(shù)和內(nèi)聚力模型結(jié)合對裂紋萌生和擴展進(jìn)行了分析，萌生階段使用內(nèi)聚力方法模擬，擴展過程使用虛擬裂紋閉合技術(shù)模擬。

內(nèi)聚力模型能有效地與有限元方法相結(jié)合，能克服裂紋尖端應(yīng)力奇異性問題，已被廣泛應(yīng)用于研究裂紋的萌生與擴展問題[14]。L.Xia 等[15]利用連續(xù)損傷力學(xué)和內(nèi)聚定律之間的關(guān)系提出了一種離散損傷模型，考慮了離散界面單元和連續(xù)單元的損傷。Y.Wang 等[16]將擴展有限元法與離散損傷模型結(jié)合，模擬了復(fù)合材料結(jié)構(gòu)的漸進(jìn)分層現(xiàn)象，采用端部開口彎曲試件、雙懸臂梁和混合模式彎曲試樣驗證了漸進(jìn)分層分析的有效性。W.Jia 等[17]建立了一種微觀內(nèi)聚力模型，用于預(yù)測單向纖維復(fù)合材料的力學(xué)性能，模擬了界面強度對復(fù)合材料力學(xué)性能和損傷的影響規(guī)律。

本文結(jié)合逐層理論、虛擬裂紋閉合技術(shù)和離散損傷模型，結(jié)合三種方法的優(yōu)勢，建立了含分層層合結(jié)構(gòu)的三維分析模型，開發(fā)數(shù)值模擬程序，對復(fù)合材料結(jié)構(gòu)層間分層的萌生和漸進(jìn)擴展行為進(jìn)行數(shù)值模擬，研究了復(fù)合材料雙懸臂梁模型的純I型分層漸進(jìn)擴展行為，討論了離散損傷模型對網(wǎng)格尺寸的敏感性，并對復(fù)合材料雙懸臂梁分層損傷擴展機理進(jìn)行探討。

1 分層損傷的逐層理論模型

1.1 逐層理論

復(fù)合材料層合板逐層理論的位移假設(shè)模式如圖1所示，x3方向上的位移插值點位于層合板的上下表面和層間界面上，層合板結(jié)構(gòu)上任意點P(x,y,z)上的位移模式可假設(shè)為

式中:a=1,2,3 分別為x,y,z方向上的位移分量；φk為沿板厚方向的一元拉格朗日插值函數(shù)。如圖1 所示，層合板均勻的分為N層，插值點數(shù)為N+ 1，uαk為厚度方向上的節(jié)點位移值，k為厚度方向上插值點的編號。

圖1 復(fù)合材料層合板逐層理論的位移假設(shè)模式Fig.1 Displacement fields of layerwise theory of composite laminate

將式（1）代入幾何方程中[5]，通過虛功原理可得到復(fù)合材料層合板結(jié)構(gòu)逐層理論的運動方程為

式中:H為復(fù)合材料層合板的厚度。

復(fù)合材料層合板第μ層本構(gòu)方程為

由于式（2）中的位移依舊是x和y的連續(xù)函數(shù)，還需在面內(nèi)對位移進(jìn)行數(shù)值離散才能對逐層理論的運動方程進(jìn)行求解。通過有限元方法求解逐層法的運動方程，厚度方向上的位移節(jié)點值u1k，u2k和u3k在第k個平面內(nèi)可離散為

式中:m= 1,2,…,nelm為二元單元的節(jié)點數(shù)，Ψm(x,y)為二元單元的形函數(shù)，uαkm為第m個平面內(nèi)節(jié)點的位移值。

將式（13）代入虛功原理，經(jīng)變分并分部積分可得到層合板結(jié)構(gòu)逐層理論運動方程的有限元列式為

其中:

1.2 分層損傷模型

在逐層理論中，因為控制方程中包含了層合板上、下表面的位移自由度，使逐層法能方便地通過層合板上/下表面的位移協(xié)調(diào)條件和內(nèi)力平衡條件與其他方法耦合，建立分析復(fù)雜層合結(jié)構(gòu)的方法。本文通過此方法，將含分層損傷復(fù)合材料層合板分解成上下兩個子板，如圖2所示，未分層連接界面處，利用連續(xù)條件將兩子板耦合，得到耦合后含分層損傷模型的整體控制方程。

圖2 復(fù)合材料層合板兩子板的耦合Fig.2 Coupling of two sub‐plates of composite laminates

對于上下子板，式（14）中的位移矢量U可以分為兩部分:未分層區(qū)域位移自由度UP1和分層區(qū)域位移自由度UP2。因此，復(fù)合材料層合板上子板總控制方程可寫為

對層合板下子板總控制方程相同處理，可得

同理，將式（16）和式（17）第二行與式（18）進(jìn)行聯(lián)合，可得到含分層復(fù)合材料層合板的總體控制方程為

1.3 基于逐層理論的虛擬裂紋閉合技術(shù)

基于線彈性斷裂力學(xué)概念發(fā)展而來的VCCT被廣泛用于計算二元和三元連續(xù)體的應(yīng)變能量釋放率[18]，本文參考D.H.Li[12]的方法，基于逐層理論，計算純I型分層前緣處的應(yīng)變能量釋放率，通過斷裂準(zhǔn)則，預(yù)測復(fù)合材料層合板的分層擴展行為。

在逐層法中，應(yīng)變能量釋放率的計算過程如圖3（b）所示。以四節(jié)點單元為例，復(fù)合材料兩子板分別由4個單元構(gòu)成，上子板下表面與下子板上表面構(gòu)成分層前緣的接觸界面。u1,l u2,l u3,l和u1,l*u2,l*u3,l*分別為層合板上、下子板x，y，z三個方向的節(jié)點張開位移，裂紋在整個分層前緣的距離Δa1內(nèi)可以虛擬閉合，虛擬閉合區(qū)域ΔA由ΔA1和ΔA2兩部分組成。s為沿分層前緣的局部坐標(biāo)，q為垂直于分層前緣的法線局部坐標(biāo)。純I 型分層的應(yīng)變能量釋放率分量GI為

圖3 逐層法計算應(yīng)變能量釋放率Fig.3 Calculation of SERR based on layerwise method

對式（20）求解，可得GI簡化式為

式中:Δa2/Δa1是用于矯正單元寬度Δai沿分層前緣變化時的相對位移系數(shù)，u3lk為全局坐標(biāo)x3方向的分層位移張開量，F(xiàn)3為使接觸節(jié)點張開的節(jié)點力。

應(yīng)變能量釋放率的分層擴展準(zhǔn)則為

對于純I型裂紋，式（22）中的GII和GIII項省略不計。計算得到GI后，通過式（22）進(jìn)行判定，若成立，則分層擴展[19]。Ed為分層擴展參數(shù)，GIC為模式I的臨界應(yīng)變能量釋放率，指數(shù)a,b,c在本文中取為1[20]。

2 基于逐層理論的離散損傷模型

2.1 離散損傷模型

基于線彈性斷裂力學(xué)的VCCT 雖然能用于預(yù)測分層，但無法在裂紋萌生階段對復(fù)合材料進(jìn)行模擬和預(yù)測[15]。B.N.Cox 等[21]指出線彈性斷裂力學(xué)概念的關(guān)鍵限制因素在于:復(fù)合材料中的損傷涉及極其復(fù)雜的非線性過程，與VCCT相比，內(nèi)聚力理論可以有效預(yù)測裂紋的萌生和擴展。因此，本文將內(nèi)聚力理論下的離散損傷模型引入逐層理論中，用以精確模擬和預(yù)測含損傷層合結(jié)構(gòu)的分層漸進(jìn)擴展行為。

結(jié)合逐層理論和內(nèi)聚力理論的三維離散損傷模型如圖4所示，模型由兩個復(fù)合材料層合板子板構(gòu)成，界面節(jié)點通過內(nèi)聚力單元（也稱彈簧單元）建立黏結(jié)關(guān)系。

圖4 三維離散損傷模型Fig.4 Three‐dimensional discrete damage model

2.2 控制方程

圖5 為含損傷Γc求解區(qū)域Ω的示意圖，Ω1+Ω2=Ω，Ω?RN(N= 1,2,3)，包括位移邊界Γu、載荷邊界Γt和損傷區(qū)域Γc，uˉ和tˉ分別為Γu和Γt上的已知位移和載荷，上子板下表面Γtf1和下子板上表面Γtf2組成表面力區(qū)域Γtf，內(nèi)部不連續(xù)損傷區(qū)域由正在發(fā)生損傷擴展區(qū)域的內(nèi)聚力區(qū)域Γcoh和表面力區(qū)域Γtf組成。作用于Γcoh上的ttop為上子板下表面內(nèi)聚力，tbot為下子板上表面內(nèi)聚力。

圖5 含離散損傷模型的求解區(qū)域Fig.5 Solving area with DDZM

裂紋尖端區(qū)域平衡方程為

式中:b為作用在求解區(qū)域的體積力；σ為應(yīng)力；?為柯西梯度算子。

Γt和Γcoh的自然邊界條件為

式中:n為邊界上的外法線單位矢量；ntop和nbot分別表示上子板下表面和下子板上表面Γcoh上的內(nèi)法線單位矢量。

內(nèi)聚力t=ttop=-tbot，通常與損傷區(qū)域上子板下表面與下子板上表面的張開位移有關(guān)，張開位移表示為

求解區(qū)域能量平衡的弱形式為

離散損傷模型的內(nèi)聚力單元總勢能可表示為

式中:NSP為沿內(nèi)聚力區(qū)域的內(nèi)聚力單元總數(shù)；tC為內(nèi)聚力，xC為空間坐標(biāo)。

2.3 本構(gòu)關(guān)系

在離散損傷模型中，復(fù)合材料分層漸進(jìn)擴展的實質(zhì)是損傷沿材料界面不斷累積的過程?；趹?yīng)變能量釋放率和離散損傷模型損傷演化規(guī)律，損傷系數(shù)Dn為[16]

式中:G1=δ1/I1，其中I1為應(yīng)變能量釋放率與張開量的轉(zhuǎn)化系數(shù)，δ1為內(nèi)聚力單元張開量。Dn為內(nèi)聚力單元損傷系數(shù)，取值范圍為0～1，Dn值從0到1表示內(nèi)聚力單元受到載荷后，從未損傷狀態(tài)到完全損傷，內(nèi)聚力單元逐漸失效的過程。

式（28）中的n代表I 型分層，Bn=1/,為內(nèi)聚力單元臨界張開量。

內(nèi)聚力單元承載力Fn由張開量δn和內(nèi)聚力剛度Kn計算得到，表達(dá)式為

將式（28）代入式（29），離散損傷模型本構(gòu)關(guān)系最終形式為

其中:

式中:B為模型寬度；Ny為寬度方向的網(wǎng)格數(shù)；ls為離散損傷模型的有限元網(wǎng)格長度。

離散損傷模型的損傷演化規(guī)律如圖6所示，圖6分別為本構(gòu)和損傷狀態(tài)圖，數(shù)字表示內(nèi)聚力單元加載狀態(tài)，其中1～2 表示內(nèi)聚力單元張開量隨位移載荷呈線性增加趨勢，不斷接近臨界張開量，此時內(nèi)聚力單元無損傷；2～4表示張開量超過臨界張開量δcr后，內(nèi)聚力單元處于非線性退化狀態(tài)，內(nèi)聚力單元剛度不斷降低，損傷系數(shù)呈非線性增長趨勢；4點后，內(nèi)聚力單元失效并完全損傷，如圖6中三維模型圖所示，內(nèi)聚力單元斷裂失效，節(jié)點失去粘結(jié)關(guān)系，分層擴展。由于位移載荷分步施加在模型上，載荷需要通過加載、卸載、再加載的循環(huán)過程，使內(nèi)聚力單元不斷張開、退化、失效，分層得以漸進(jìn)擴展。因此，位移載荷的加載-卸載-再加載過程描述如圖6所示。

圖6 離散損傷模型的損傷演化關(guān)系Fig.6 Damage evolution relationship for DDZM

2.4 有限元方程

對于逐層理論，考慮離散損傷模型的虛功原理為

式中:δKD為分層區(qū)域內(nèi)聚力的虛勢能，具體表達(dá)式如下

因此，

式中:內(nèi)聚力合力矩具體表達(dá)式為

式（33）可進(jìn)一步推導(dǎo)為

含離散損傷模型的逐層理論最終控制方程為

本文方法的整體框架如圖7 所示。首先，基于逐層理論建立分層損傷模型，得到模型所有節(jié)點位移值；其次，引入虛擬裂紋閉合技術(shù)，計算分層前緣應(yīng)變能量釋放率，依據(jù)式（22），研究分層損傷模型分層擴展行為；最后，引入離散損傷模型，以應(yīng)變能量釋放率作為內(nèi)聚力單元剛度退化的依據(jù)，模擬和預(yù)測二元和三元雙懸臂梁純I型分層的漸進(jìn)萌生和擴展行為。

圖7 本文方法的整體框架Fig.7 Overall framework of the method in this paper

2.5 程序?qū)崿F(xiàn)

為了便于理解建立的分層損傷模型，給出了相應(yīng)的程序流程圖，如圖8所示。

圖8 程序?qū)崿F(xiàn)流程圖Fig.8 Flowchart of computational procedure

首先，分別讀取復(fù)合材料上/下子板的輸入文件，隨后處理邊界條件，確定兩子板的稀疏矩陣維度，并將內(nèi)聚力單元嵌入至上下兩子板層間界面處。根據(jù)式（17）確定兩子板的接觸和非接觸自由度，并進(jìn)行行列變換，得到兩子板耦合后的總剛度矩陣。位移載荷做分步施加處理，獲得每一步兩子板節(jié)點位移矩陣，并根據(jù)式（21），計算分層前緣應(yīng)變能量釋放率值。計算離散損傷模型損傷系數(shù)，判定內(nèi)聚力單元是否完全損傷，若不滿足條件，則內(nèi)聚力剛度退化，直至達(dá)到完全損傷條件，內(nèi)聚力單元完全失效，分層得到漸進(jìn)擴展。輸出結(jié)果，程序停止運算。

3 數(shù)值算例

3.1 直線型分層前緣應(yīng)變能量釋放率驗證

為驗證本文方法的正確性，利用Li[12]中的算例進(jìn)行對比，含分層雙懸臂梁模型如圖9所示，每一單層材料特性為

在圖9中，a為預(yù)置分層區(qū)域長度，區(qū)域L-a通過位移協(xié)調(diào)條件和內(nèi)力平衡條件約束，分層前緣區(qū)域網(wǎng)格局部加密處理，雙懸臂梁張口端施加P= 0.5064N/m的線分布均布載荷，總厚度H為3mm。

圖9 具有分層的復(fù)合材料DCB模型Fig.9 Composite DCB model with delamination

Kruger[22]基于非線性三維殼單元和VCCT 研究了同樣的模型，對GI值進(jìn)行了如下?lián)Q算

式中:h為梁總厚度的一半，即單個子模型總厚度；B為梁的寬度。

本文方法所得結(jié)果與文獻(xiàn)結(jié)果如圖10 所示，計算出的應(yīng)變能量釋放率與參考文獻(xiàn)給出的結(jié)果吻合較好。由于做了歸一化處理，以GIC值為分界點，大于等于GIC值時，分層將漸進(jìn)擴展，小于GIC值時，分層將不擴展。由圖10 歸一化GI值沿寬度方向的分布可知，在分層擴展前，GI值呈兩邊低，中間高的分布趨勢，即中間分層快于兩邊，這是因為DCB模型受到均布載荷后，寬度方向兩側(cè)邊附近（W=0為模型最左側(cè)，W=1為最右側(cè)）受到長度方向的壓縮、拉伸應(yīng)力，DCB 模型寬度方向的壓縮與拉伸由泊松效應(yīng)引起，在W=0和W=1 附近產(chǎn)生正、負(fù)應(yīng)變，使DCB 模型的寬度方向上呈現(xiàn)圖10 所示的中間高、兩邊低分布。此分布狀態(tài)顯示了分層漸進(jìn)擴展優(yōu)先發(fā)生在分層前緣中間區(qū)域，隨后以中間節(jié)點為基準(zhǔn)向兩側(cè)邊節(jié)點逐個擴展，使分層前緣產(chǎn)生C 形彎曲形態(tài)。

圖10 分層前緣應(yīng)變能量釋放率結(jié)果對比Fig.10 Comparison between strain energy release rate results at the front edge of the delamination

DCB模型厚度方向應(yīng)力云圖如圖11所示，從圖11中可以看出，厚度方向上的應(yīng)力集中現(xiàn)象主要出現(xiàn)在分層前緣附近，并且隨著分層的漸進(jìn)擴展，分層前緣和應(yīng)力集中現(xiàn)象不斷向固定端移動。

圖11 分層前緣橫向應(yīng)力的分布Fig.11 Distribution of transverse stress at the leading edge of the delamination

3.2 內(nèi)聚力單元受壓縮載荷的驗證

復(fù)合材料結(jié)構(gòu)實際服役過程中，不僅會受到拉伸載荷，也會受到壓縮載荷的作用[23]。因此，在不考慮損傷擴展引起材料性能退化的情況下，驗證復(fù)合材料結(jié)構(gòu)受壓時內(nèi)聚力單元傳遞載荷的行為。內(nèi)聚力單元受壓，內(nèi)聚力張開量δn為負(fù)值，此時懲罰剛度Kp代替拉伸剛度Kn，防止含損傷復(fù)合材料結(jié)構(gòu)受壓后，分層區(qū)域發(fā)生網(wǎng)格嵌入現(xiàn)象。為避免計算中出現(xiàn)數(shù)值病態(tài)和結(jié)果無法收斂或收斂困難等問題，懲罰剛度Kp通常取模型整體剛度矩陣最大值的104～107倍[24]，本文取104倍。

幾何尺寸如圖12（a）所示，未分層區(qū)域I 和II 內(nèi)的節(jié)點滿足位移協(xié)調(diào)和內(nèi)力平衡條件。網(wǎng)格尺寸如圖12（b）所示，長a= 1.0m；寬b= 1.0m；高h(yuǎn)= 0.1m，面內(nèi)x、y方向均采用36個網(wǎng)格，分層區(qū)域x方向網(wǎng)格數(shù)和未分層區(qū)域x方向網(wǎng)格數(shù)均為12，單個子板厚度方向均分為4 層，鋪層順序為[0/90/0/90]s，在分層區(qū)域施加垂直于z軸向下的均布載荷P= 1 × 108N，內(nèi)聚力單元內(nèi)嵌于兩板接觸表面分層區(qū)域節(jié)點上，懲罰剛度為Kp= 1 × 1014N/m，邊界條件為:x= 0，x=1兩邊固支。材料參數(shù)為

模型考慮了三種情況，情況1:無分層損傷；情況2:兩子板接觸面含分層區(qū)域，區(qū)域內(nèi)嵌入內(nèi)聚力單元，即圖12（a）中灰色分層區(qū)域；情況3:兩子板接觸面之間含分層區(qū)域，區(qū)域不嵌入內(nèi)聚力單元。

圖12 內(nèi)聚力單元受壓縮載荷模型Fig.12 Cohesive element subjected to compressive loading model

上板中心點z方向第10步最大位移計算結(jié)果見表1，整體變形如圖13所示。

表1 上板中心點z方向最大位移（單位：mm）Table 1 Maximum displacement in the z-direction at the center of upper subplate

對比情況1 和情況2，z方向最大位移和變形云圖幾乎一致，說明在不考慮損傷導(dǎo)致層合板材料性能退化的前提下，內(nèi)聚力單元受壓而不受損，并傳遞壓縮載荷至下子板，引起下子板發(fā)生變形，并且厚度方向兩子板接觸界面附近的位移、應(yīng)力、應(yīng)變值連續(xù)。對比情況2和情況3，情況3沒有嵌入內(nèi)聚力單元，兩子板分層區(qū)域不受位移協(xié)調(diào)和內(nèi)力平衡條件的約束。因此，上子板受壓后，與下子板發(fā)生網(wǎng)格嵌入現(xiàn)象，導(dǎo)致厚度方向上的位移、應(yīng)力、應(yīng)變值不連續(xù)，如圖13（c）所示。

圖13 三種情況的位移云圖Fig.13 Displacement clouds for three cases

因此，本節(jié)驗證了數(shù)值模擬程序的正確性，顯示了內(nèi)聚力單元能有效阻止兩板接觸界面與分層區(qū)域發(fā)生相嵌入的現(xiàn)象。

3.3 雙懸臂梁分層漸進(jìn)擴展

本節(jié)對雙懸臂梁模型進(jìn)行研究，分析分層的萌生和分層漸進(jìn)擴展行為,模型如圖14所示，尺寸為:長a=100.0mm，寬b=1.0mm，高2h=3.0mm，位移載荷P=12.0mm，臨界應(yīng)變能量釋放率為GIC=1.0N/mm，材料力學(xué)參數(shù)為:E=135.0×109Pa，G=54.0×109Pa，ν=0.25。

圖14 離散損傷雙懸臂梁模型Fig.14 Discrete damage double cantilever beam model

模型長度方向，考慮6 種網(wǎng)格劃分:40、50、100、200、400、800，寬度方向網(wǎng)格被均劃分為1，厚度方向兩子板均劃分為4個網(wǎng)格，探究網(wǎng)格密度對純I型雙懸臂梁分層擴展結(jié)果的敏感性?？紤]5.7MPa、20MPa、57MPa三種不同界面強度[16]，與數(shù)值解[16]和解析解[25]進(jìn)行對比。載荷位移曲線如圖15 所示，圖中分別給出了界面強度分別為5.7MPa、20MPa、57MPa的載荷位移曲線結(jié)果。

由結(jié)果可以看出，隨著界面強度的增大，載荷位移曲線結(jié)果不斷向解析解靠近；界面強度的不同，影響了界面分層起裂成程度的不同，即相同網(wǎng)格密度下，較大界面強度需要較大的位移載荷，使模型開始發(fā)生分層行為。

圖15 中載荷位移曲線可分為三段:AB和CD段中，載荷隨位移的循環(huán)加載呈線性增長趨勢，循環(huán)載荷加載到B處時，模型的GI值等于GIC，預(yù)置分層起裂；加載B處后，GI值大于GIC，分層失穩(wěn)擴展，載荷下降；BC段隨著載荷的增加，內(nèi)聚力單元開始退化，內(nèi)聚力單元剛度值不斷下降，分層漸進(jìn)擴展。由此可見，臨界應(yīng)變能量釋放率決定了界面層間分層是否可以發(fā)生擴展行為，并且，載荷位移曲線反映了界面層間內(nèi)聚力單元剛度的變化，與預(yù)置分層的起裂和擴展相關(guān)，層間的分層擴展導(dǎo)致了內(nèi)聚力單元剛度的下降。

圖15 不同網(wǎng)格密度和界面強度的位移載荷曲線Fig.15 Load‐displacement curves obtained with different mesh densities and interface strengths

隨著DCB長度方向網(wǎng)格的不斷加密，載荷位移曲線不斷逼近解析解，網(wǎng)格密度為40 × 1、50 × 1計算結(jié)果與解析解曲線偏差較遠(yuǎn)，較容易發(fā)生數(shù)值震蕩現(xiàn)象，導(dǎo)致結(jié)果收斂困難。這是因為小網(wǎng)格數(shù)模型，內(nèi)聚力單元數(shù)量比大網(wǎng)格數(shù)模型少，單個內(nèi)聚力單元的剛度增大，此時需要更大的位移載荷使內(nèi)聚力單元張開量增大，所以得到更高的位移載荷值。在加載—卸載—再加載過程中，由于內(nèi)聚力剛度的變化，出現(xiàn)了數(shù)值震蕩現(xiàn)象。網(wǎng)格800 × 1 得到的曲線與Wang 的結(jié)果吻合較好；結(jié)果表明，當(dāng)網(wǎng)格密度為200～800時，離散損傷模型對網(wǎng)格變化不敏感。同時也驗證了逐層理論、虛擬裂紋閉合技術(shù)、內(nèi)聚力理論三種方法結(jié)合建立分層損傷模型的正確性。

考慮厚度方向網(wǎng)格數(shù)對計算結(jié)果的影響，界面強度20MPa，臨界應(yīng)變能量釋放率為GIC= 1.0N/mm，網(wǎng)格數(shù)為400 × 1，對厚度方向網(wǎng)格數(shù)為4、8、12、16進(jìn)行了分析，結(jié)果如圖16所示。

由圖16 可推斷離散損傷模型在厚度方向上的網(wǎng)格敏感性沒有面內(nèi)長度方向網(wǎng)格敏感性強，厚度方向網(wǎng)格數(shù)量幾乎對計算出來的結(jié)果無明顯影響。

圖16 厚度方向分層數(shù)收斂性驗證Fig.16 Verification of convergence of stratification number in thickness direction

圖17 給出了DCB 模型的應(yīng)力云圖，Δ 為位移載荷值。隨著位移載荷的增加，分層處于漸進(jìn)擴展?fàn)顟B(tài)，并且分層前緣處出現(xiàn)很大的應(yīng)力集中現(xiàn)象。

圖17 離散損傷模型下的DCB模型σzz方向分層擴展過程Fig.17 DCB model under DDZM model σzz directional hierarchical delamination process

3.4 三元雙懸臂梁的分層漸進(jìn)擴展

本節(jié)分析了三元雙懸臂梁的分層漸進(jìn)擴展行為，模型尺寸如圖18 所示。長L= 150.0mm，寬度B分為10mm、20mm 和30mm 三種情況。對不同寬度，選取相同y方向網(wǎng)格數(shù)進(jìn)行分析計算，高2h= 3.0mm，位移載荷P= 3.0mm，臨界應(yīng)變能量釋放率為GIC= 1.0N/mm，兩子模型厚度方向均分為4 層，初始分層長度L0= 35.0mm，各層材料參數(shù)為

圖18 三維DCB模型Fig.18 3D DCB model

根據(jù)式（31），將模型擴展至三維后，內(nèi)聚力初始剛度通過2B/(Ny+1)項修正，使y方向網(wǎng)格數(shù)對結(jié)果產(chǎn)生的影響降至最低。

圖19 給出了三種寬度計算的載荷位移曲線，與解析解[25]進(jìn)行了對比。結(jié)果表明，雙懸臂梁模型寬度的不同，引起支反力最大值發(fā)生變化，大寬度的模型需要較大的支反力使分層漸進(jìn)擴展；小寬度載荷位移曲線更接近解析解，因為y方向網(wǎng)格數(shù)通過公式修正后，對離散損傷模型結(jié)果的影響不大，但本文通過應(yīng)變能量釋放率作為內(nèi)聚力剛度退化的依據(jù)，雙懸臂梁模型y方向的網(wǎng)格少時，計算出的應(yīng)變能量釋放率與文獻(xiàn)解[12]存在誤差，引起離散損傷模型中損傷系數(shù)Dn發(fā)生微小變化，進(jìn)一步影響載荷位移曲線偏離解析解。因此，較小寬度模型計算的應(yīng)變能量釋放率結(jié)果更準(zhǔn)確，載荷位移曲線更接近解析解。

圖19 不同寬度雙懸臂梁模型載荷—位移曲線Fig.19 Load‐displacement curves of laminates with different widths

圖20（a）分別為分層起裂和分層漸進(jìn)擴展時，層間界面分層前緣輪廓線的移動過程圖，圖20（b）分別為分層起裂和分層漸進(jìn)擴展時的應(yīng)力云圖。隨著位移載荷加載—卸載—再加載的過程，模型的內(nèi)聚力單元剛度逐漸降低并退化，分層前緣內(nèi)聚力單元斷裂失效，分層前緣向固定端移動，層間分層漸進(jìn)擴展。

由圖20（a）中分層前緣輪廓線可知，層間分層起裂后，分層前緣呈中間高，兩邊低的擴展形態(tài)，與本文3.1 節(jié)虛擬閉合裂紋技術(shù)模擬分層漸進(jìn)擴展的結(jié)果相符合。從三元角度對分層漸進(jìn)擴展過程進(jìn)行模擬和預(yù)測，驗證本文方法的正確性。

圖20 復(fù)合材料層合板分層擴展過程圖Fig.20 Delamination extension process of composite laminates

4 結(jié)論

本文結(jié)合逐層理論、虛擬裂紋閉合技術(shù)和離散損傷模型，建立了用于分析復(fù)合材料分層損傷的漸近擴展分析模型。該模型將含分層損傷復(fù)合材料層合板分解為上下兩子板，通過位移協(xié)調(diào)和內(nèi)力平衡條件，建立兩子板界面之間連接關(guān)系。通過數(shù)值算例，研究了復(fù)合材料雙懸臂梁和三維雙懸臂梁分層損傷的漸近擴展機理，得到以下幾點結(jié)論。

（1）虛擬裂紋閉合技術(shù)可以有效與離散損傷模型結(jié)合，通過應(yīng)變能量釋放率作為內(nèi)聚力單元剛度退化的依據(jù)，能準(zhǔn)確地描述復(fù)合材料分層的萌生與漸進(jìn)擴展過程。

（2）雙懸臂梁和三元雙懸臂梁厚度方向上的網(wǎng)格密度不影響計算精度；面內(nèi)網(wǎng)格數(shù)控制在合理區(qū)間內(nèi)，計算結(jié)果對網(wǎng)格數(shù)不敏感，若超出合理的網(wǎng)格區(qū)間，計算結(jié)果會出現(xiàn)結(jié)果收斂困難和震蕩等現(xiàn)象。

（3）界面強度對載荷位移曲線結(jié)果有較大影響，為了得到更精確的結(jié)果可選用較大的界面強度。

（4）模型寬度方向上的網(wǎng)格密度通過修正后幾乎對結(jié)果無明顯影響，寬度方向網(wǎng)格數(shù)主要影響應(yīng)變能釋放率結(jié)果，進(jìn)而影響載荷位移曲線結(jié)果。