基于均匀网格的复合材料层合板分层损伤建模

钱若力，李梁轶，李顶河，蔡舒妤

（中国民航大学航空工程学院，天津 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]建立复合材料层合板损伤分析模型，扩展逐层理论厚度方向上的位移场包括线性拉格朗日插值函数、一维弱不连续函数和强不连续函数。强不连续函数和弱不连续函数分别模拟由层间分层引起的位移不连续性和由层间界面引起的应变不连续性。基体裂纹则由扩展有限元方法在面内位移离散中进行模拟。扩展逐层理论不仅可以描述复合材料层合板的多层分层和基体裂纹损伤，还可以求得裂纹尖端和分层前缘的位移场和应力场。

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 结语

复合材料结构的分层损伤一般可由无损检测技术探测，形状复杂的裂纹尖端会给建模和分析带来巨大困难。基于扩展逐层理论、多维最大熵原理和分层前缘追踪算法，本文提出了一种基于均匀网格的快速分层建模方法，减小建模的难度，还以含方形、圆形、双圆形分层区域的方形层合板为例，验证了该方法的有效性，从而得到以下结论：

（1）尽管本文方法的分层前缘不够光滑，但其重建的分层区域与真实情况十分接近，该方法在建模难度上远小于精确前缘/三维弹性模型和FUGD/VCCTXLWM 模型，且该方法能模拟分层区域的动态变化过程；

（2）有两种方法可以用来降低本文方法计算能量释放率时的缺点，一种是细化网格来减少振荡，但计算的成本会加大，另一种是利用水平集函数重建更光滑、更精确的分层前缘；

（3）在实际情况中，分层形状不是简单光滑的，而是随分层扩展而改变的，本文方法能模拟任意分层形状，且能有效跟踪分层扩展过程中前缘形状的任意改变，因此，基于非结构化网格和无损检测图像，本文方法能快速有效地得到分层扩展形状；

（4）对于多层分层损伤区域，快速建模和分析在工程运用上意义重大，扩展逐层理论能用于分析多层分层损伤问题，但必须进一步改进多维最大熵原理，使其能识别多层分层损伤区域。