基于均勻網格的復合材料層合板分層損傷建模

錢若力，李梁軼，李頂河，蔡舒妤

（中國民航大學航空工程學院，天津 300300）

預測復合材料結構分層損傷的起始和擴展是目前復合材料結構損傷研究中最棘手的問題之一。近年來，大量學者針對分層損傷研究提出各種分析方法和計算模型，并得到了廣泛應用。目前，擴展有限元方法（XFEM，extended finite element method）與虛擬裂紋閉合技術（VCCT，virtual crack closure technique）或內聚力模型（CZM，cohesive zone model）相結合的思路被廣泛應用于預測和分析分層損傷的擴展問題[1]。

Ashari 等[2-3]利用正交各向異性擴充函數改進了擴展有限元方法，用于模擬分析兩種正交各向異性介質之間的界面裂紋，擴充函數基于無外力作用下界面裂紋的漸進解析位移場得到。Nagashima 等[4-5]利用三維實體元和二維擴展有限元方法研究了層合板結構的分層擴展問題。

Benvenuti 等[6]基于連續-非連續擴展有限元方法提出了纖維增強混凝土板的分層分析模型，把纖維增強混凝土板和粘接層看成線彈性材料，而混凝土的力學性能由連續損傷模型模擬。Borst 等[7]基于擴展有限元方法和內聚力模型研究了圓形初始分層損傷以及正弦型面外損傷問題。Curiel 等[8]利用擴展有限元方法研究了金屬纖維層合板的分層問題。Grogan 等[9]建立了一種用于模擬熱疲勞分層問題的分析方法，預測了分層損傷擴展方向及復合材料低溫燃料箱的滲透率。Jas＇kowiec[10]對含分層損傷的玻璃纖維夾層結構進行三維分析，通過擴展有限元方法引入由于分層引起的位移不連續函數，用于描述有限元中的分層節點。

在基于內聚力模型的有限元方法中，兩層之間的界面采用內聚力單元建模，引入損傷系數后可以用逐漸減小單元剛度模擬分層損傷的漸進擴展。目前，基于內聚力模型的擴展有限元方法也被廣泛用于分析層合板的分層問題。Yazdani 等[11]建立了一個無摩擦線性接觸模型，并利用一階剪切變形理論（FSDT，firstorder shear deformation theory）和擴展有限元方法來模擬層合板的線性和幾何非線性分層問題。Sun 等[12]提出一種將內聚力模型與擴展有限元方法相結合的新方法建立三維分層模型，通過擴展有限元方法的擴充節點與單元劃分對基體裂紋進行有效建模，并通過附著分離法對層間分層進行建模，用于分析失效過程中的層間分層和基體裂紋之間的相互作用。內聚力單元分析分層損傷的漸進擴展存在計算效率低，計算結果難以收斂的問題。Motamedi 等[13]運用內聚力模型對單向復合材料的I 型分層進行數值模擬，只在分層前緣引入內聚力模型，有效提升了計算效率。Wang 等[14]結合離散損傷模型、擴展有限元以及連續損傷力學方法，建立了一種分層損傷的漸進分析模型，用牛頓迭代求解，計算結果可以快速收斂。Zhao 等[15]利用內聚力模型和擴展有限元方法在Abaqus 軟件中建立了分層模型。Jiao 等[16]基于有限元網格重分技術[17-18]提出了一種自適應分層分析法，以解決大規模層合板分層分析的計算效率問題。該方法證明了有限元網格重分技術、擴展有限元方法與逐層離散化在分析層合板分層問題上是等效的[19]。

VCCT 常用來計算分層前緣能量釋放率[20]。近30年中，基于VCCT 的有限元方法被廣泛用于分析復合材料的分層問題[21-24]。Lua 等[25]結合虛擬裂紋閉合技術和有限元方法建立了一種自動分層分析方法，該方法能預測含多分層復合材料粘結接頭處的剩余強度。Bacarreza 等[26]分析了含復雜形狀的混合型分層疲勞擴展問題。Alalade 等[27]應用擴展有限元方法檢測分層和裂紋，提出一種動態耦合問題的逆分析方法，該方法能夠識別裂紋的位置、尺寸和方向。Zhang 等[28]采用奈爾-米德（NM）和擬牛頓（QN）優化方法對彈性動力學中的裂紋識別問題進行了數值求解，研究了各種動力試驗載荷對裂紋識別的影響，與靜態載荷相比，動態載荷對裂紋檢測更有效。VCCT 雖然可以精確分析在靜態和動態載荷下的分層擴展，但是需要在分層前緣進行密集的網格劃分，內聚力單元可以很好地解決這一問題。

由于復合材料層合板的損傷特性和分層前緣的形狀極其復雜，在實際工程領域中，通常采用無損檢測技術來探測分層損傷，通過C 掃描或X 射線掃描得到的分層損傷邊緣是非常模糊的，因此傳統的建模方法不能快速精確地重建分層區域。在基于有限元的建模方法中，單元邊界必須與分層前緣保持一致，在建模與分析過程中很難實現。在損傷區域不規則的情況下，通常需借助有限元軟件完成不規則的網格劃分，需要繁瑣的建模過程，快速建模方法對模型進行均勻網格劃分，只需程序就可以實現，明顯提升了建模效率。在損傷邊緣位置需要一定程度上加密網格，以提高計算精度，勢必會造成計算效率的降低。本文利用C 掃描或X 射線得到初始分層損傷的檢測照片，基于擴展逐層方法[29-32]、多維最大熵原理[33]以及分層前緣追蹤算法[21-22]提出一種基于均勻有限元網格且能夠快速實現分層的建模方法，可以有效提高建模效率和精度。

1 復合材料分層損傷的快速建模方法

基于無損檢測掃描圖像和均勻有限元網格，本文建立的快速建模方法，如圖1 所示。該方法建模過程分為兩步：首先，通過掃描圖像和均勻有限網格，對分層損傷區域進行快速重建；然后，通過虛擬裂紋閉合技術重建分層前緣，并建立分析模型。最初的分層形狀可由C 掃描或X 射線掃描得到，多維最大熵原理用于確定分層區域，找出分層區域的像素點，處于分層區域的有限元節點被標記為分層節點。通過這些已經標記為分層的節點，結合分層前緣追蹤算法[21-22]，可以精確跟蹤分層前緣。

圖1 復合材料層合板分層損傷區域的快速建模方法Fig.1 Fast modeling scheme for delamination damage region in the composite laminates

由快速建模方法得到的最終分層前緣示意圖如圖2 所示，根據Xie 等[21-22]的追蹤算法得到分層前緣是鋸齒狀的，并在一些節點上發生偏移。這可能會導致能量釋放率隨分層前緣的擴展而發生振蕩，精細化的有限元網格可減小這種振蕩，提高對分層擴展路徑的預測精度。

圖2 任意分層前緣分析模型Fig.2 Analysis model of arbitrary delamination front

2 多維最大熵原理

利用多維最大熵原理確定分層區域[33-34]，該方法已經成功解決圖像重建和目標識別等問題。多維最大熵原理表示，在信息量給定且約束條件也確定的條件下，其分布的最佳狀態是熵總和達到最大值。

如圖2 所示，假設掃描圖中有M×N 個像素點，像素點xp的灰度函數定義為

式中：n 表示維度；xp（x，y）為像素點的坐標，1≤x≤M，1≤y≤N；Fm∈[0，1，…，L -1]，m=0，1，…，L -1，L=256 代表灰度。

根據灰度閾值Ti，i=1，2，…，n，可將掃描后得到的像素點分為兩大類：

第一類為目標類

第二類為背景類

其中，目標類處于分層區域內，而背景類則位于分層區域外。

將灰度函數的頻率設定為him（fi=Fm），則聯合概率質量函數可定義如下

利用二維最大熵方法確定復合材料層合結構的分層損傷，找到每個像素點周圍3 × 3 范圍內像素點平均灰度值，則熵函數可定義為

式中，上標1 和2 分別代表目標類和背景類。

利用文獻[35]中的C 掃描圖像來驗證快速建模方法的有效性，如圖3 所示。圖3（c）和圖3（e）以像素點的灰度值分別表示灰度圖像和最終識別的結果。對于3 種不同的有限元網格（81×51，121×81，201×131），分層節點表示的分層區域如圖3（b）、3（d）、3（f）所示?？梢钥闯?，隨著節點數量的增加，分層區域更接近C掃描的真實情況。

圖3 C 掃描圖像和分層識別結果Fig.3 C-scan image and identification results of delamination

3 基于虛擬裂紋閉合技術的擴展逐層方法

利用擴展逐層理論[29-32]建立復合材料層合板損傷分析模型，擴展逐層理論厚度方向上的位移場包括線性拉格朗日插值函數、一維弱不連續函數和強不連續函數。強不連續函數和弱不連續函數分別模擬由層間分層引起的位移不連續性和由層間界面引起的應變不連續性?；w裂紋則由擴展有限元方法在面內位移離散中進行模擬。擴展逐層理論不僅可以描述復合材料層合板的多層分層和基體裂紋損傷，還可以求得裂紋尖端和分層前緣的位移場和應力場。

3.1 復合材料層合板的擴展逐層理論

對于含多層分層損傷的復合材料層合板，擴展逐層理論的位移場如圖4 所示。其中，hk表示第k 層的厚度，zk代表第k 層和第k-1 層之間厚度方向上的坐標。左側數字1st，2nd，…，（NC+2）th 表示復合材料層合板物理層的中性層，即擴展逐層理論在復合材料層合板厚度方向上建立的位移場節點；右側數字1st，2nd，…，（NC-1）th 表示復合材料層合板的實際分層界面。對于含多層分層損傷的復合材料層合板上任一點（x，y，z）的位移可表示為

圖4 含NC 層的復合材料層合板Fig.4 Composite laminate containing NC layer

式中：α=1，2，3 表示x，y，z 方向上的分量；t 表示時間變量；uαik，uαlk，uαrk分別表示插值點的標準位移自由度、分層損傷引起的位移不連續附加自由度以及層間界面引起的應變不連續附加自由度；下標i，l 和r 分別表示標準自由度、分層的附加自由度和層間界面的附加自由度；NC表示層合板數學層層數；ND表示由于分層而擴充的插值點的數量。從圖4 中可以看出，層間界面處的標準自由度和附加自由度數量分別為NC+2和NC-1。令Φik=?k（z），?k（）是層合板厚度方向上的拉格朗日插值函數；令Φrk=Θk（z），Θk（z）=?k（z）χk（z）是用于模擬分層界面的弱不連續函數，χk（z）是一維符號距離函數；令Φlk=Ξk（z），Ξk（z）=?k（z）Hk（z）是用于模擬分層損傷的強不連續函數，Hk（z）是一維Heaviside 函數。

運用Einstein 求和約定，式（7）可簡化為

式中：重復指標ζ 和k 在取值范圍內求和，k 的取值取決于ζ 的取值，如當ζ=i 時，k 可以取1～NC+2，但當ζ=r 時k 只能取1～NC-1。

擴展逐層理論的基本思想是將復雜的三維斷裂問題分解為一維和二維斷裂問題。對于含多層分層損傷和基體裂紋的復合材料層合結構，將拉格朗日插值函數、裂紋面擴充函數和裂尖擴充形函數線性相加，則每個單元上的節點位移自由度uαik以及附加自由度uαlk和uαrk可表示為

若無基體裂紋，則式（9）可表示為

式中：Uαζkq（t）為分層節點間位移，q=u，s，hb；u=1，2，…，NE，NE為平面內的有限元節點數；s=1，2，…，NEP，NEP為橫向裂紋的面內擴充節點數；hb=1，2，…，NEQ，NEQ為橫向裂紋尖端的面內擴充節點數分別為標準節點自由度、橫向裂紋面自由度和橫向裂紋尖端自由度；ψu（x，y）=?u（x，y）為拉格朗日插值函數；Λs（x，y）=?s（x，y）Hs（x，y）為橫向裂紋的形函數，Hs（x，y）為二維Heaviside 函數；Πhb（x，y）=?hb（x，y）Fhb（x，y）為橫向裂紋尖端的形函數，Fhb（x，y）為分支函數。

Hamilton 原理中的虛應變能、外力虛功和虛動能可以根據式（8）、式（10）的假設以及層合板的應變-位移關系來建立。含多層分層損傷的層合板的控制方程為

式中：u=m，s，hb；v=n，g，fb；m，n=1，2，…，NE；s，g=1，2，…，NEP；hb，fb=1，2，…，NEQ；Kαβζηkeuv為不含橫向裂紋的層合板的剛度矩陣，表示為

單元剛度矩陣子項的具體表達式參考文獻[32]。在式（12）中，（m，n），（s，g）和（hb，fb）分別對應形狀函數ψ（），Λ（），Π（）。對于復合材料層合板，可令1/R=0，其中R 是層合板的內徑。

3.2 分層前緣跟蹤算法

利用Xie 等[21-22]提出的算法來追蹤形狀上任意改變的移動分層前緣，如圖5 所示。

圖5 分層前緣和封閉裂紋區域的定義Fig.5 Definition of the delamination front and closed crack area

定義從中心節點到周圍節點的8 個基本向量Ri（i=1，2，...，8），分層由變量mi（i=1，2，…，8）定義，mi=0 表示完好節點，mi=1 則表示分層節點。Rb和Re兩向量分別表示分層區域的初始端和末端。完好區可由Rb和Re逆時針轉動得到，因此分層前緣也可以由這兩個向量來定義。在已知損傷狀態變量mi的情況下，分層前緣（Rb和R）e可由基向量的線性組合得到，即

中心節點處的單位向量以及切向量定義為

式中：nb=∑mi-1（1 -m）i（-Riyi+Rixj）是初始端向量Rb的法方向向量；ne=∑（1-m）imi+1（Riyi-Rixj）是末端向量Re的法方向向量；Rix，Riy分別表示Ri在x，y 方向上的分量；i，j 和k 是單位向量。

3 個向量（n0，q，k）構成局部坐標系的一組基，能量釋放率也是在這個坐標系下計算的。在自然坐標系（ξ，η）中，如圖5（b）所示，閉合區域可以分為A1和A2兩部分，可由高斯積分計算得到。每4 個點2；j=1，2，3，4）分別圍成兩個子區域，這些點在自然坐標系下的具體表達式見參考文獻[21-22]。為了得到兩個未知點P21和P42的具體位置，首先需要得到兩個臨時點的位置。點是平行于向量n^ 過點P21的線（1）與平行于向量Re過點P31的線（2）的交點。

一旦計算出能量釋放率分量，可利用復合型斷裂判據來預測分層是否擴展，即

式中：Ed是分層擴展參數；GⅠ、GⅡ和GⅢ分別對應Ⅰ、Ⅱ和Ⅲ型分層前緣能量釋放率；GⅠC、GⅡC和GⅢC分別對應Ⅰ、Ⅱ和Ⅲ型分層臨界能量釋放率；指數a，b，c 取1[21]。

本文研究的分層擴展是由mi表示的已知分層區域開始的。首先，在每一步分析過程中，根據多維熵原理確定所有節點上的mi值。然后，用Xie 等[21-22]的算法追蹤分層前緣，用虛擬裂紋閉合技術計算分層前緣處節點上（mi=0）的能量釋放率。最后，若分層前緣上節點的能量釋放率滿足混合型斷裂判據（式（16）），則令該點的mi=1，并繼續下一步計算。

4 數值算例

4.1 矩形分層區域

考慮一塊含正方形分層區域（10 m × 10 m）的方板（20 m×20 m），厚度為ht=0.1 m，如圖6（a）所示。采用10 201 單元10 404 節點的均勻有限元網格對求解區域進行劃分，如圖6（b）所示。在厚度方向上，層合板被均勻地分為8 個子層，分層位于第4 和第5 層之間，分層區域中心與板的中心重合。材料參數取為：彈性模量E=3.585×1011Pa，泊松比μ=0.3。板的上表面受P=120 kN/m2的均布載荷作用，邊界條件為四邊固支。

圖6 含矩形分層區域的方板Fig.6 The square plate with a square delamination region

矩形分層區域的識別結果如圖7 所示，快速建模/擴展逐層法（FUGD/VCCT-XLWM）重建的分層區域比實際的分層區域大。但從圖7 可以看出，重建分層區域與實際分層區域之間的差別由前緣處劃分的有限單元的大小控制。因此，分層前緣附近采用更密的有限元網格可提高識別精度。

圖7 矩形分層區域的識別結果Fig.7 The identification result of the square delamination region

為了驗證方法的正確性，建立了三維彈性模型，并采用8 節點實體單元對該問題進行分析。在三維彈性模型中，分層界面上的一對重復節點用于模擬位移的不連續性，且與單元的邊界重合，因此，三維彈性模型中的分層前緣是精確的，該模型被稱為精確前緣/三維彈性模型。此外，直接利用VCCT-XLWM 建立了分析模型，采用10 201 單元10 404 節點進行計算，平面內網格與三維彈性模型一樣，因此分層前緣也是精確的，該模型稱為精確前緣/擴展逐層法（EF/VCCT-XLWM）模型。

由FUGD/VCCT-XLWM、精確前緣/三維彈性和EF/VCCT-XLWM 模型計算得到的最大位移（u1和u3）和應力（σ11和σ12）如表1 所示，能量釋放率（ERR，energy relense rate）GI如圖8 所示，其中，θ 是分層區域的幅角（極點位于分層區域的中心），θ=0°代表軸x1的方向。

表1 由3 種模型計算出的矩形分層板的最大位移和應力Tab.1 The maximum displacements and stresses of the square delamination plate calculated by three kinds of models

圖8 由FUGD/VCCT-XLWM 和EF/VCCT-XLWM得到的能量釋放率G1Fig.8 ERR G1 obtained by FUGD/VCCT-XLWM and XLWM

從表1 和圖8 中可以看出，FUGD/VCCT-XLWM計算結果足夠精確。對于各向同性板上的方形分層，GⅠ的最大值出現在邊緣中心點，而在頂點處GⅠ幾乎為0。因此，4 條邊的中心點是最先出現分層擴展的位置。

下面考慮厚度方向上分為8 個子層（[0]8）的層合板。該板的幾何形狀、邊界條件與各向同性板相同。上表面受P=140 kN/m2載荷，每單層的材料參數為：E11=1.09×1011Pa，E22=E33=8.82×109Pa，G12=G13=4.32 × 109Pa，G23=3.2 × 109Pa，μ32=0.52，μ12=μ13=0.342，GⅠC=0.306 N/mm，GⅡC=0.632 N/mm，GⅢC=0.817 N/mm。

對比算例采用EF/VCCT-XLWM 方法分析該問題。圖9 是這兩種方法得到的最終擴展區域，可以看出，由FUGD/VCCT-XLWM 預測的擴展路徑與EF/VCCTXLWM 模型的吻合得非常好。在均勻受載的情況下，層合板的方形分層區域最終擴展成橢圓，其長軸沿材料主方向。由于裂紋前緣為鋸齒狀的，因此G 隨分層前緣的擴展而振蕩。從圖9 中可以看出，如果主要關注分層擴展的形狀，這種振蕩并不會影響分層擴展路徑的預測。在大多數實際情況下，分層不會是簡單光滑的形狀，而是會隨擴展而變化。FUGD/VCCT-XLWM 可以近似模擬出分層形狀并且使分層任意擴展。因此，本文提出的模擬方案可從非結構化網格和X 射線圖像得到合理的近似的分層擴展形狀。

圖9 含方形分層層合板的擴展路徑Fig.9 Extended path of the laminate with square delamination

4.2 圓形分層區域

以一個含圓形分層區域（半徑0.5 m）的方板為例，分層的位置、材料參數、幾何尺寸以及邊界條件都與4.1 節的方板相同。該板在厚度方向上分為8 個單層，板的上表面受P=140 kN/m2的均布載荷作用。

EF/VCCT-XLWM 模型的有限單元網格和FUGD/VCCT-XLWM 的網格如圖10 所示，其中，紅色代表分層區域。

圖10 含圓形分層區域的方形板Fig.10 The square plate with a circular delamination region

從圖10 中可以看出，盡管FUGD/VCCT-XLWM方法重建的分層區域的前緣不是很光滑，但與EF/VCCT-XLWM 模型相近。在EF/VCCT-XLWM 模型中，邊緣處的單元與分層前緣一致，因此在建模方面比較困難。而FUGD/VCCT-XLWM 的分析模型基于均勻有限單元網格，在建模方面的難度要低于精確前緣/三維彈性模型和EF/VCCT-XLWM 模型，而且可以很大程度提高計算效率。

表2 中列出了用3 種模型計算出的最大位移及應力。對于精確前緣/三維彈性模型，有限元網格如圖10（a）所示。FUGD/VCCT-XLWM 用到了兩種均勻網格：網格1，10 201 單元10 404 節點；網格2，22 500 單元22 801 節點。EF/VCCT-XLWM 模型中有10 900 個四節點單元和11 041 個節點。

表2 3 種模型計算出的含圓形分層板的最大位移和應力Tab.2 The maximum displacements and stresses of the plate with circular delamination calculated by three kinds of models

從表2 可以看出，本文方法的網格1 能夠得到足夠精確的靜力結果。盡管使用網格1 的FUGD/VCCTXLWM 和EF/VCCT-XLWM 模型的計算成本一樣，但前者在最大位移、應力上更接近精確前緣/三維彈性模型的計算結果。

接下來考慮含圓形分層區域的復合材料層合板。其分層的位置、材料參數、幾何尺寸以及邊界條件與4.1 節中的層合板相同。板的上表面受P=180 kN/m2的均勻載荷作用，鋪層順序為[0]8。

圖11 給出了分層前緣的最終擴展區域，FUGD/VCCT-XLWM 預測得到的擴展路徑與EF/VCCT-XLWM模型得到的很接近。與方形分層區域一樣，層合板上的圓形分層區域最終會擴展成橢圓形，其長軸方向為材料主方向。

圖11 含圓形分層層合板的擴展路徑Fig.11 Extended path of the laminate with circular delamination

4.3 雙圓形分層區域

以含雙圓形分層區域的方形層合板為例，兩個圓形分層區域的半徑為0.5 m，且兩圓心之間的距離為5 m。該層合板分層的位置、材料參數、幾何尺寸以及邊界條件都與4.1 節中的層合板一致。板的上表面受P=140 kN/m2的均勻載荷作用。圖12 給出了兩種方法的有限元網格劃分。

圖12 含雙圓形分層區域的方形板有限元網格Fig.12 The finite element meshing of the square plate with doubly circular delamination region

從圖12 中可以看出，重建的分層區域與真實情況非常接近，尤其是在尖角處。FUGD/VCCT-XLWM 也使用4.2 節中模型的2 種網格，EF/VCCT-XLWM 模型和精確前緣/三維彈性模型中都有10 201 個四節點單元和10 404 個節點。

表3 給出了這3 種分析模型的最大位移和應力，從表3 中可以看出，采用網格1 的FUGD/VCCT-XLWM最適合雙圓形分層區域。圖13 給出了3 種方法得到的能量釋放率GI，最大值出現在兩個圓形分層的交點處（θ=90°和270°）。因此，雙圓形分層最先在這2 個交點處擴展。

表3 含雙圓形分層板的最大位移和應力Tab.3 The maximum displacements and stresses of the plate with doubly circular delamination

圖13 含雙圓形分層區域方形板的能量釋放率GIFig.13 ERR GI of the square plate with doubly circular delamination region

5 結語

復合材料結構的分層損傷一般可由無損檢測技術探測，形狀復雜的裂紋尖端會給建模和分析帶來巨大困難?；跀U展逐層理論、多維最大熵原理和分層前緣追蹤算法，本文提出了一種基于均勻網格的快速分層建模方法，減小建模的難度，還以含方形、圓形、雙圓形分層區域的方形層合板為例，驗證了該方法的有效性，從而得到以下結論：

（1）盡管本文方法的分層前緣不夠光滑，但其重建的分層區域與真實情況十分接近，該方法在建模難度上遠小于精確前緣/三維彈性模型和FUGD/VCCTXLWM 模型，且該方法能模擬分層區域的動態變化過程；

（2）有兩種方法可以用來降低本文方法計算能量釋放率時的缺點，一種是細化網格來減少振蕩，但計算的成本會加大，另一種是利用水平集函數重建更光滑、更精確的分層前緣；

（3）在實際情況中，分層形狀不是簡單光滑的，而是隨分層擴展而改變的，本文方法能模擬任意分層形狀，且能有效跟蹤分層擴展過程中前緣形狀的任意改變，因此，基于非結構化網格和無損檢測圖像，本文方法能快速有效地得到分層擴展形狀；

（4）對于多層分層損傷區域，快速建模和分析在工程運用上意義重大，擴展逐層理論能用于分析多層分層損傷問題，但必須進一步改進多維最大熵原理，使其能識別多層分層損傷區域。