基于逐层理论的复合材料层合结构分层损伤扩展研究

李顶河，曹江涛，郭巧荣

中国民航大学，天津 300300

广泛应用于航空、航天等领域的复合材料层合结构存在多种损伤形式，最主要、最常见的损伤形式是分层损伤，它的存在严重影响了复合材料在服役过程中的强度和刚度。利用试验方法对复合材料结构进行研究和分析，将消耗大量的时间、人力和物力成本。基于数值模拟技术对分层损伤进行分析将在很大程度上克服这些问题，成为研究分层损伤问题的主要手段之一[1-4]。

针对复合材料厚板结构，J. N. Reddy[5]提出了逐层理论，该理论能给出精确的层间应力应变，也能描述厚板面内位移在厚度方向上的Zig-Zag 特征。近年来，国内外许多学者对Reddy 的逐层理论进行了优化和改进，并将此理论广泛应用于复合材料层合结构、功能梯度层合结构、压电层合结构等问题的分析中[6]。K.M.Liew 等[7]总结回顾了逐层理论的最新发展。M.S.Beg等[8]基于逐层理论提出了一种用于分析弧形弯曲梁问题的梁理论，预测了功能梯度弧形弯曲梁在静态、自由和受迫振动情况下的精确响应。王瑞鹏[9]将逐层理论与三维实体元结合，对夹芯结构进行了静力和自由振动分析，表明逐层理论与三维实体元法进行耦合能用于分析复杂点阵夹芯结构的静动力问题。

虚拟裂纹闭合技术被广泛应用于计算应变能量释放率和应力强度因子。余毅等[10]基于虚拟裂纹闭合法计算应力强度因子，模拟了三维裂纹受拉伸载荷作用下的裂纹扩展情况。鲁龙坤等[11]解释了裂纹尖端张开角的定义，并验证裂纹尖端张开角是一个有效的断裂参数。D.H.Li[12]将扩展逐层理论与虚拟裂纹闭合技术相结合，用以预测复合材料层合板壳任意形状分层和横向裂纹的同时扩展。J.Jokinen 等[13]等利用虚拟裂纹闭合技术和内聚力模型结合对裂纹萌生和扩展进行了分析，萌生阶段使用内聚力方法模拟，扩展过程使用虚拟裂纹闭合技术模拟。

内聚力模型能有效地与有限元方法相结合，能克服裂纹尖端应力奇异性问题，已被广泛应用于研究裂纹的萌生与扩展问题[14]。L.Xia 等[15]利用连续损伤力学和内聚定律之间的关系提出了一种离散损伤模型，考虑了离散界面单元和连续单元的损伤。Y.Wang 等[16]将扩展有限元法与离散损伤模型结合，模拟了复合材料结构的渐进分层现象，采用端部开口弯曲试件、双悬臂梁和混合模式弯曲试样验证了渐进分层分析的有效性。W.Jia 等[17]建立了一种微观内聚力模型，用于预测单向纤维复合材料的力学性能，模拟了界面强度对复合材料力学性能和损伤的影响规律。

本文结合逐层理论、虚拟裂纹闭合技术和离散损伤模型，结合三种方法的优势，建立了含分层层合结构的三维分析模型，开发数值模拟程序，对复合材料结构层间分层的萌生和渐进扩展行为进行数值模拟，研究了复合材料双悬臂梁模型的纯I型分层渐进扩展行为，讨论了离散损伤模型对网格尺寸的敏感性，并对复合材料双悬臂梁分层损伤扩展机理进行探讨。

1 分层损伤的逐层理论模型

1.1 逐层理论

复合材料层合板逐层理论的位移假设模式如图1所示，x3方向上的位移插值点位于层合板的上下表面和层间界面上，层合板结构上任意点P(x,y,z)上的位移模式可假设为

式中:a=1,2,3 分别为x,y,z方向上的位移分量；φk为沿板厚方向的一元拉格朗日插值函数。如图1 所示，层合板均匀的分为N层，插值点数为N+ 1，uαk为厚度方向上的节点位移值，k为厚度方向上插值点的编号。

图1 复合材料层合板逐层理论的位移假设模式Fig.1 Displacement fields of layerwise theory of composite laminate

将式（1）代入几何方程中[5]，通过虚功原理可得到复合材料层合板结构逐层理论的运动方程为

式中:H为复合材料层合板的厚度。

复合材料层合板第μ层本构方程为

由于式（2）中的位移依旧是x和y的连续函数，还需在面内对位移进行数值离散才能对逐层理论的运动方程进行求解。通过有限元方法求解逐层法的运动方程，厚度方向上的位移节点值u1k，u2k和u3k在第k个平面内可离散为

式中:m= 1,2,…,nelm为二元单元的节点数，Ψm(x,y)为二元单元的形函数，uαkm为第m个平面内节点的位移值。

将式（13）代入虚功原理，经变分并分部积分可得到层合板结构逐层理论运动方程的有限元列式为

其中:

1.2 分层损伤模型

在逐层理论中，因为控制方程中包含了层合板上、下表面的位移自由度，使逐层法能方便地通过层合板上/下表面的位移协调条件和内力平衡条件与其他方法耦合，建立分析复杂层合结构的方法。本文通过此方法，将含分层损伤复合材料层合板分解成上下两个子板，如图2所示，未分层连接界面处，利用连续条件将两子板耦合，得到耦合后含分层损伤模型的整体控制方程。

图2 复合材料层合板两子板的耦合Fig.2 Coupling of two sub‐plates of composite laminates

对于上下子板，式（14）中的位移矢量U可以分为两部分:未分层区域位移自由度UP1和分层区域位移自由度UP2。因此，复合材料层合板上子板总控制方程可写为

对层合板下子板总控制方程相同处理，可得

同理，将式（16）和式（17）第二行与式（18）进行联合，可得到含分层复合材料层合板的总体控制方程为

1.3 基于逐层理论的虚拟裂纹闭合技术

基于线弹性断裂力学概念发展而来的VCCT被广泛用于计算二元和三元连续体的应变能量释放率[18]，本文参考D.H.Li[12]的方法，基于逐层理论，计算纯I型分层前缘处的应变能量释放率，通过断裂准则，预测复合材料层合板的分层扩展行为。

在逐层法中，应变能量释放率的计算过程如图3（b）所示。以四节点单元为例，复合材料两子板分别由4个单元构成，上子板下表面与下子板上表面构成分层前缘的接触界面。u1,l u2,l u3,l和u1,l*u2,l*u3,l*分别为层合板上、下子板x，y，z三个方向的节点张开位移，裂纹在整个分层前缘的距离Δa1内可以虚拟闭合，虚拟闭合区域ΔA由ΔA1和ΔA2两部分组成。s为沿分层前缘的局部坐标，q为垂直于分层前缘的法线局部坐标。纯I 型分层的应变能量释放率分量GI为

图3 逐层法计算应变能量释放率Fig.3 Calculation of SERR based on layerwise method

对式（20）求解，可得GI简化式为

式中:Δa2/Δa1是用于矫正单元宽度Δai沿分层前缘变化时的相对位移系数，u3lk为全局坐标x3方向的分层位移张开量，F3为使接触节点张开的节点力。

应变能量释放率的分层扩展准则为

对于纯I型裂纹，式（22）中的GII和GIII项省略不计。计算得到GI后，通过式（22）进行判定，若成立，则分层扩展[19]。Ed为分层扩展参数，GIC为模式I的临界应变能量释放率，指数a,b,c在本文中取为1[20]。

2 基于逐层理论的离散损伤模型

2.1 离散损伤模型

基于线弹性断裂力学的VCCT 虽然能用于预测分层，但无法在裂纹萌生阶段对复合材料进行模拟和预测[15]。B.N.Cox 等[21]指出线弹性断裂力学概念的关键限制因素在于:复合材料中的损伤涉及极其复杂的非线性过程，与VCCT相比，内聚力理论可以有效预测裂纹的萌生和扩展。因此，本文将内聚力理论下的离散损伤模型引入逐层理论中，用以精确模拟和预测含损伤层合结构的分层渐进扩展行为。

结合逐层理论和内聚力理论的三维离散损伤模型如图4所示，模型由两个复合材料层合板子板构成，界面节点通过内聚力单元（也称弹簧单元）建立黏结关系。

图4 三维离散损伤模型Fig.4 Three‐dimensional discrete damage model

2.2 控制方程

图5 为含损伤Γc求解区域Ω的示意图，Ω1+Ω2=Ω，Ω⊂RN(N= 1,2,3)，包括位移边界Γu、载荷边界Γt和损伤区域Γc，uˉ和tˉ分别为Γu和Γt上的已知位移和载荷，上子板下表面Γtf1和下子板上表面Γtf2组成表面力区域Γtf，内部不连续损伤区域由正在发生损伤扩展区域的内聚力区域Γcoh和表面力区域Γtf组成。作用于Γcoh上的ttop为上子板下表面内聚力，tbot为下子板上表面内聚力。

图5 含离散损伤模型的求解区域Fig.5 Solving area with DDZM

裂纹尖端区域平衡方程为

式中:b为作用在求解区域的体积力；σ为应力；∇为柯西梯度算子。

Γt和Γcoh的自然边界条件为

式中:n为边界上的外法线单位矢量；ntop和nbot分别表示上子板下表面和下子板上表面Γcoh上的内法线单位矢量。

内聚力t=ttop=-tbot，通常与损伤区域上子板下表面与下子板上表面的张开位移有关，张开位移表示为

求解区域能量平衡的弱形式为

离散损伤模型的内聚力单元总势能可表示为

式中:NSP为沿内聚力区域的内聚力单元总数；tC为内聚力，xC为空间坐标。

2.3 本构关系

在离散损伤模型中，复合材料分层渐进扩展的实质是损伤沿材料界面不断累积的过程。基于应变能量释放率和离散损伤模型损伤演化规律，损伤系数Dn为[16]

式中:G1=δ1/I1，其中I1为应变能量释放率与张开量的转化系数，δ1为内聚力单元张开量。Dn为内聚力单元损伤系数，取值范围为0～1，Dn值从0到1表示内聚力单元受到载荷后，从未损伤状态到完全损伤，内聚力单元逐渐失效的过程。

式（28）中的n代表I 型分层，Bn=1/,为内聚力单元临界张开量。

内聚力单元承载力Fn由张开量δn和内聚力刚度Kn计算得到，表达式为

将式（28）代入式（29），离散损伤模型本构关系最终形式为

其中:

式中:B为模型宽度；Ny为宽度方向的网格数；ls为离散损伤模型的有限元网格长度。

离散损伤模型的损伤演化规律如图6所示，图6分别为本构和损伤状态图，数字表示内聚力单元加载状态，其中1～2 表示内聚力单元张开量随位移载荷呈线性增加趋势，不断接近临界张开量，此时内聚力单元无损伤；2～4表示张开量超过临界张开量δcr后，内聚力单元处于非线性退化状态，内聚力单元刚度不断降低，损伤系数呈非线性增长趋势；4点后，内聚力单元失效并完全损伤，如图6中三维模型图所示，内聚力单元断裂失效，节点失去粘结关系，分层扩展。由于位移载荷分步施加在模型上，载荷需要通过加载、卸载、再加载的循环过程，使内聚力单元不断张开、退化、失效，分层得以渐进扩展。因此，位移载荷的加载-卸载-再加载过程描述如图6所示。

图6 离散损伤模型的损伤演化关系Fig.6 Damage evolution relationship for DDZM

2.4 有限元方程

对于逐层理论，考虑离散损伤模型的虚功原理为

式中:δKD为分层区域内聚力的虚势能，具体表达式如下

因此，

式中:内聚力合力矩具体表达式为

式（33）可进一步推导为

含离散损伤模型的逐层理论最终控制方程为

本文方法的整体框架如图7 所示。首先，基于逐层理论建立分层损伤模型，得到模型所有节点位移值；其次，引入虚拟裂纹闭合技术，计算分层前缘应变能量释放率，依据式（22），研究分层损伤模型分层扩展行为；最后，引入离散损伤模型，以应变能量释放率作为内聚力单元刚度退化的依据，模拟和预测二元和三元双悬臂梁纯I型分层的渐进萌生和扩展行为。

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

2.5 程序实现

为了便于理解建立的分层损伤模型，给出了相应的程序流程图，如图8所示。

图8 程序实现流程图Fig.8 Flowchart of computational procedure

首先，分别读取复合材料上/下子板的输入文件，随后处理边界条件，确定两子板的稀疏矩阵维度，并将内聚力单元嵌入至上下两子板层间界面处。根据式（17）确定两子板的接触和非接触自由度，并进行行列变换，得到两子板耦合后的总刚度矩阵。位移载荷做分步施加处理，获得每一步两子板节点位移矩阵，并根据式（21），计算分层前缘应变能量释放率值。计算离散损伤模型损伤系数，判定内聚力单元是否完全损伤，若不满足条件，则内聚力刚度退化，直至达到完全损伤条件，内聚力单元完全失效，分层得到渐进扩展。输出结果，程序停止运算。

3 数值算例

3.1 直线型分层前缘应变能量释放率验证

为验证本文方法的正确性，利用Li[12]中的算例进行对比，含分层双悬臂梁模型如图9所示，每一单层材料特性为

在图9中，a为预置分层区域长度，区域L-a通过位移协调条件和内力平衡条件约束，分层前缘区域网格局部加密处理，双悬臂梁张口端施加P= 0.5064N/m的线分布均布载荷，总厚度H为3mm。

图9 具有分层的复合材料DCB模型Fig.9 Composite DCB model with delamination

Kruger[22]基于非线性三维壳单元和VCCT 研究了同样的模型，对GI值进行了如下换算

式中:h为梁总厚度的一半，即单个子模型总厚度；B为梁的宽度。

本文方法所得结果与文献结果如图10 所示，计算出的应变能量释放率与参考文献给出的结果吻合较好。由于做了归一化处理，以GIC值为分界点，大于等于GIC值时，分层将渐进扩展，小于GIC值时，分层将不扩展。由图10 归一化GI值沿宽度方向的分布可知，在分层扩展前，GI值呈两边低，中间高的分布趋势，即中间分层快于两边，这是因为DCB模型受到均布载荷后，宽度方向两侧边附近（W=0为模型最左侧，W=1为最右侧）受到长度方向的压缩、拉伸应力，DCB 模型宽度方向的压缩与拉伸由泊松效应引起，在W=0和W=1 附近产生正、负应变，使DCB 模型的宽度方向上呈现图10 所示的中间高、两边低分布。此分布状态显示了分层渐进扩展优先发生在分层前缘中间区域，随后以中间节点为基准向两侧边节点逐个扩展，使分层前缘产生C 形弯曲形态。

图10 分层前缘应变能量释放率结果对比Fig.10 Comparison between strain energy release rate results at the front edge of the delamination

DCB模型厚度方向应力云图如图11所示，从图11中可以看出，厚度方向上的应力集中现象主要出现在分层前缘附近，并且随着分层的渐进扩展，分层前缘和应力集中现象不断向固定端移动。

图11 分层前缘横向应力的分布Fig.11 Distribution of transverse stress at the leading edge of the delamination

3.2 内聚力单元受压缩载荷的验证

复合材料结构实际服役过程中，不仅会受到拉伸载荷，也会受到压缩载荷的作用[23]。因此，在不考虑损伤扩展引起材料性能退化的情况下，验证复合材料结构受压时内聚力单元传递载荷的行为。内聚力单元受压，内聚力张开量δn为负值，此时惩罚刚度Kp代替拉伸刚度Kn，防止含损伤复合材料结构受压后，分层区域发生网格嵌入现象。为避免计算中出现数值病态和结果无法收敛或收敛困难等问题，惩罚刚度Kp通常取模型整体刚度矩阵最大值的104～107倍[24]，本文取104倍。

几何尺寸如图12（a）所示，未分层区域I 和II 内的节点满足位移协调和内力平衡条件。网格尺寸如图12（b）所示，长a= 1.0m；宽b= 1.0m；高h= 0.1m，面内x、y方向均采用36个网格，分层区域x方向网格数和未分层区域x方向网格数均为12，单个子板厚度方向均分为4 层，铺层顺序为[0/90/0/90]s，在分层区域施加垂直于z轴向下的均布载荷P= 1 × 108N，内聚力单元内嵌于两板接触表面分层区域节点上，惩罚刚度为Kp= 1 × 1014N/m，边界条件为:x= 0，x=1两边固支。材料参数为

模型考虑了三种情况，情况1:无分层损伤；情况2:两子板接触面含分层区域，区域内嵌入内聚力单元，即图12（a）中灰色分层区域；情况3:两子板接触面之间含分层区域，区域不嵌入内聚力单元。

图12 内聚力单元受压缩载荷模型Fig.12 Cohesive element subjected to compressive loading model

上板中心点z方向第10步最大位移计算结果见表1，整体变形如图13所示。

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

对比情况1 和情况2，z方向最大位移和变形云图几乎一致，说明在不考虑损伤导致层合板材料性能退化的前提下，内聚力单元受压而不受损，并传递压缩载荷至下子板，引起下子板发生变形，并且厚度方向两子板接触界面附近的位移、应力、应变值连续。对比情况2和情况3，情况3没有嵌入内聚力单元，两子板分层区域不受位移协调和内力平衡条件的约束。因此，上子板受压后，与下子板发生网格嵌入现象，导致厚度方向上的位移、应力、应变值不连续，如图13（c）所示。

图13 三种情况的位移云图Fig.13 Displacement clouds for three cases

因此，本节验证了数值模拟程序的正确性，显示了内聚力单元能有效阻止两板接触界面与分层区域发生相嵌入的现象。

3.3 双悬臂梁分层渐进扩展

本节对双悬臂梁模型进行研究，分析分层的萌生和分层渐进扩展行为,模型如图14所示，尺寸为:长a=100.0mm，宽b=1.0mm，高2h=3.0mm，位移载荷P=12.0mm，临界应变能量释放率为GIC=1.0N/mm，材料力学参数为:E=135.0×109Pa，G=54.0×109Pa，ν=0.25。

图14 离散损伤双悬臂梁模型Fig.14 Discrete damage double cantilever beam model

模型长度方向，考虑6 种网格划分:40、50、100、200、400、800，宽度方向网格被均划分为1，厚度方向两子板均划分为4个网格，探究网格密度对纯I型双悬臂梁分层扩展结果的敏感性。考虑5.7MPa、20MPa、57MPa三种不同界面强度[16]，与数值解[16]和解析解[25]进行对比。载荷位移曲线如图15 所示，图中分别给出了界面强度分别为5.7MPa、20MPa、57MPa的载荷位移曲线结果。

由结果可以看出，随着界面强度的增大，载荷位移曲线结果不断向解析解靠近；界面强度的不同，影响了界面分层起裂成程度的不同，即相同网格密度下，较大界面强度需要较大的位移载荷，使模型开始发生分层行为。

图15 中载荷位移曲线可分为三段:AB和CD段中，载荷随位移的循环加载呈线性增长趋势，循环载荷加载到B处时，模型的GI值等于GIC，预置分层起裂；加载B处后，GI值大于GIC，分层失稳扩展，载荷下降；BC段随着载荷的增加，内聚力单元开始退化，内聚力单元刚度值不断下降，分层渐进扩展。由此可见，临界应变能量释放率决定了界面层间分层是否可以发生扩展行为，并且，载荷位移曲线反映了界面层间内聚力单元刚度的变化，与预置分层的起裂和扩展相关，层间的分层扩展导致了内聚力单元刚度的下降。

图15 不同网格密度和界面强度的位移载荷曲线Fig.15 Load‐displacement curves obtained with different mesh densities and interface strengths

随着DCB长度方向网格的不断加密，载荷位移曲线不断逼近解析解，网格密度为40 × 1、50 × 1计算结果与解析解曲线偏差较远，较容易发生数值震荡现象，导致结果收敛困难。这是因为小网格数模型，内聚力单元数量比大网格数模型少，单个内聚力单元的刚度增大，此时需要更大的位移载荷使内聚力单元张开量增大，所以得到更高的位移载荷值。在加载—卸载—再加载过程中，由于内聚力刚度的变化，出现了数值震荡现象。网格800 × 1 得到的曲线与Wang 的结果吻合较好；结果表明，当网格密度为200～800时，离散损伤模型对网格变化不敏感。同时也验证了逐层理论、虚拟裂纹闭合技术、内聚力理论三种方法结合建立分层损伤模型的正确性。

考虑厚度方向网格数对计算结果的影响，界面强度20MPa，临界应变能量释放率为GIC= 1.0N/mm，网格数为400 × 1，对厚度方向网格数为4、8、12、16进行了分析，结果如图16所示。

由图16 可推断离散损伤模型在厚度方向上的网格敏感性没有面内长度方向网格敏感性强，厚度方向网格数量几乎对计算出来的结果无明显影响。

图16 厚度方向分层数收敛性验证Fig.16 Verification of convergence of stratification number in thickness direction

图17 给出了DCB 模型的应力云图，Δ 为位移载荷值。随着位移载荷的增加，分层处于渐进扩展状态，并且分层前缘处出现很大的应力集中现象。

图17 离散损伤模型下的DCB模型σzz方向分层扩展过程Fig.17 DCB model under DDZM model σzz directional hierarchical delamination process

3.4 三元双悬臂梁的分层渐进扩展

本节分析了三元双悬臂梁的分层渐进扩展行为，模型尺寸如图18 所示。长L= 150.0mm，宽度B分为10mm、20mm 和30mm 三种情况。对不同宽度，选取相同y方向网格数进行分析计算，高2h= 3.0mm，位移载荷P= 3.0mm，临界应变能量释放率为GIC= 1.0N/mm，两子模型厚度方向均分为4 层，初始分层长度L0= 35.0mm，各层材料参数为

图18 三维DCB模型Fig.18 3D DCB model

根据式（31），将模型扩展至三维后，内聚力初始刚度通过2B/(Ny+1)项修正，使y方向网格数对结果产生的影响降至最低。

图19 给出了三种宽度计算的载荷位移曲线，与解析解[25]进行了对比。结果表明，双悬臂梁模型宽度的不同，引起支反力最大值发生变化，大宽度的模型需要较大的支反力使分层渐进扩展；小宽度载荷位移曲线更接近解析解，因为y方向网格数通过公式修正后，对离散损伤模型结果的影响不大，但本文通过应变能量释放率作为内聚力刚度退化的依据，双悬臂梁模型y方向的网格少时，计算出的应变能量释放率与文献解[12]存在误差，引起离散损伤模型中损伤系数Dn发生微小变化，进一步影响载荷位移曲线偏离解析解。因此，较小宽度模型计算的应变能量释放率结果更准确，载荷位移曲线更接近解析解。

图19 不同宽度双悬臂梁模型载荷—位移曲线Fig.19 Load‐displacement curves of laminates with different widths

图20（a）分别为分层起裂和分层渐进扩展时，层间界面分层前缘轮廓线的移动过程图，图20（b）分别为分层起裂和分层渐进扩展时的应力云图。随着位移载荷加载—卸载—再加载的过程，模型的内聚力单元刚度逐渐降低并退化，分层前缘内聚力单元断裂失效，分层前缘向固定端移动，层间分层渐进扩展。

由图20（a）中分层前缘轮廓线可知，层间分层起裂后，分层前缘呈中间高，两边低的扩展形态，与本文3.1 节虚拟闭合裂纹技术模拟分层渐进扩展的结果相符合。从三元角度对分层渐进扩展过程进行模拟和预测，验证本文方法的正确性。

图20 复合材料层合板分层扩展过程图Fig.20 Delamination extension process of composite laminates

4 结论

本文结合逐层理论、虚拟裂纹闭合技术和离散损伤模型，建立了用于分析复合材料分层损伤的渐近扩展分析模型。该模型将含分层损伤复合材料层合板分解为上下两子板，通过位移协调和内力平衡条件，建立两子板界面之间连接关系。通过数值算例，研究了复合材料双悬臂梁和三维双悬臂梁分层损伤的渐近扩展机理，得到以下几点结论。

（1）虚拟裂纹闭合技术可以有效与离散损伤模型结合，通过应变能量释放率作为内聚力单元刚度退化的依据，能准确地描述复合材料分层的萌生与渐进扩展过程。

（2）双悬臂梁和三元双悬臂梁厚度方向上的网格密度不影响计算精度；面内网格数控制在合理区间内，计算结果对网格数不敏感，若超出合理的网格区间，计算结果会出现结果收敛困难和震荡等现象。

（3）界面强度对载荷位移曲线结果有较大影响，为了得到更精确的结果可选用较大的界面强度。

（4）模型宽度方向上的网格密度通过修正后几乎对结果无明显影响，宽度方向网格数主要影响应变能释放率结果，进而影响载荷位移曲线结果。