先进复合材料具有比刚度高、比强度高和优秀的抗疲劳特性等出色性能特点，被广泛地应用于各个行业之中。尤其是在航空航天领域中，最新军用或民用飞机的大多部件都是由复合材料制成的^[1]。这些结构在生产、服役以及维修期间可能会遭受到各类冲击事件，例如天气突变落下的冰雹、跑道上溅起的碎石和掉落的维修工具等^[2]。虽然这类低速冲击可能无法对表面造成可视损伤，但是在复合材料的内部已经造成了一定程度的基体裂纹或者层间分层等不可逆的损伤模式，严重降低了其结构强度和性能，造成了一定程度的安全隐患^[3]。这种目视勉强可检损伤(BVID)对复合材料的设计和发展提出了很大的挑战^[4]。鉴于试验测试的高昂成本以及无法精准地实时监控BVID，数值模拟被认为是一种很有效的方法，它能在较短的时间内模拟冲击行为并预测出复合材料内部复杂的破坏机制。因此，国内外很多学者使用有限元数值模拟的方法来研究复合材料层合板在低速冲击载荷作用下的损伤问题^[5-12]。

建立合适的本构模型可以模拟复合材料层合板的冲击行为。在众多发展良好的数值模拟技术中，考虑了损伤起始以及刚度退化过程的渐进损伤模型(PDM)逐渐成为了最适用的一种模拟方案^[10-12]。拓宏亮等^[10]使用最大应变失效准则与三维Puck准则建立了关于复合材料层合板低速冲击的三维数值模型，得到了四级能量冲击下的损伤和响应规律。Li等^[11]评估了失效准则和损伤演化方法在低速冲击下复合材料层合板有限元分析中的适用性。Zhou等^[12]基于三维Hashin失效准则，提出了含有厚度方向应力的损伤演化模型，并且分析了模型中参数对预测结果的影响。通常而言，使用带有评估纤维与基体在拉伸和压缩载荷作用下失效准则是研究人员的首要选择。当损伤开始时，也就是达到失效准则时，复合材料的受损区域由于微裂纹的扩展而引起局部刚度的降低，承载能力逐步减弱，此时需要一个合适的损伤演化模型来模拟受损区域的损伤累积过程。一部分学者在刚度矩阵中加入了预定义的常数来规定刚度退化的过程^[7-9]。显然，这种方法只能用于解决特殊问题，并不适用于一般的情况。因此，为了能够预测大多数复合材料的渐进损伤行为，在损伤演化模型中应用结合材料参数的交互式方法是一项重要的研究内容。

在复合材料层合板的典型损伤模式中，层合板间分层是最重要的模式。共有3种常用的模拟复合材料分层的模型：基于破坏准则的模型、基于断裂力学的模型和内聚区模型(CZM)。第一种模型将分层作为失效准则中的一类模式，例如Hou提出的三维失效准则。虽然这种方法具有良好的数值计算效率，但是由于不同失效模式之间会相互影响，分层的失效判据会跟随不同情况而改变。第二种模型是基于虚拟裂纹闭合技术(VCCT)^[13]发展的。使用VCCT需要预先设置初始裂纹，并且在分层的前端使用自适应网格重划分技术，在处理三维有限元时会遇到巨大的计算量，降低了数值计算效率。值得一提的是，第三种模型CZM^[14]可以克服前2个模型的缺点。CZM使用基于强度和断裂能量的准则来预测分层损伤的萌生与扩展，通常采用双线性牵引力-分离法则来描述单一方向的内聚元的力学行为。

本文在ABAQUS / Explicit中建立了一种三维有限元模型用来模拟并分析碳纤维增强复合材料层合板的低速冲击损伤情况；结合Hashin失效准则和Hou失效准则的损伤判据，基于等效位移法的线性退化方案作为损伤演化模型，在层间插入了零厚度的内聚元，结合双线性牵引力-分离法则与Benzeggagh and Kenane(B-K)准则预测了分层损伤情况，并将此渐进损伤模性编入了用户自定义子程序(VUMAT)。应用了一般接触算法模拟冲击头与层合板表面的接触问题，同时按照试验施加约束条件并赋予冲头预定义初速度，与Hongkarnjanakul等^[15]做的冲击试验结果进行了比较，验证了数值模型，详细地讨论了复合材料层合板在低速冲击作用下局部损伤情况。

1 低速冲击的渐进损伤模型

一般而言，低速冲击对复合材料层合板造成的损伤可以分为2大类：①单向板内发生的层内损伤，包括纤维的拉伸/压缩断裂和基体的拉伸/压缩破坏；②2个相邻单向板之间发生的层间损伤，也就是分层破坏。低速冲击损伤的过程可以通过渐进损伤模型来描述，包含了连续损伤模型、损伤起始判据和损伤演化模型。

1.1 连续损伤模型

因为复合材料结构会积累损伤直至完全破坏，所以不能应用简单的破坏准则来预测结构的完整性。连续损伤模型就是适用于评估这类损伤累积的破坏模式，其本构方程如(1)式所示

(1)

式中：σ_ij, ε_kl和C_ijkl(i, j, k, l=1, 2, 3)分别代表工程应力、工程应变和刚度系数; 下角标1, 2, 3分别指面内纵向, 面内横向和面外方向。对于复合材料层压板来说, 每一层单向板都可被视为横观各向同性材料, 可以用5个独立常数的本构方程来描述未受损的单向板。通过引入不同的损伤变量来规定复合材料在加载过程中的退化, 使用修正的本构方程来描述其应力应变关系。含有损伤变量的柔度退化矩阵S_d如(2)式所示

(2)

式中, d_f和d_m分别代表纤维和基体的损伤变量。由柔度退化矩阵S_d可以推出刚度退化矩阵C_d如式(3)所示

(3)

式中：

(4)

d_ft和d_fc表示纤维拉伸失效的损伤变量与纤维压缩失效的损伤变量; d_mt和d_mc表示基体拉伸失效的损伤变量与基体压缩失效的损伤变量。

刚度退化矩阵C_d中含有损伤变量使得数值计算很难收敛。因此, 采用Duvaut and Lions^[16]正则化模型使得刚度退化过程平滑而容易收敛, 如(5)式所示

(5)

式中: d_I^ν表示失效模式I的正则化损伤变量; 符号上一点表示时间的一阶导数; λ_I表示失效模式I的黏性系数。

1.2 损伤起始判据

采用Hashin失效准则作为纤维损伤的起始判据, 如下所示:

纤维拉伸失效(σ₁₁≥0)

(6)

纤维压缩失效(σ₁₁ < 0)

(7)

采用Hou失效准则作为基体损伤的起始判据, 如下所示:

基体拉伸失效(σ₂₂≥0)

(8)

基体压缩失效(σ₂₂ < 0)

(9)

式中: σ_ij(i, j=1, 2, 3)代表各方向的应力; X_T和X_C代表纵向拉伸强度与压缩强度; Y_T和Y_C代表横向拉伸强度与压缩强度; S₁₂, S₁₃和S₂₃代表纵向和横向的剪切强度。

1.3 损伤演化模型

当满足损伤起始判据后, 持续地加载会使得刚度系数降低, 结构承载能力逐步减弱, 此时需要定义一个合适的损伤演化模型来模拟受损区域的损伤累积过程。在有限元模拟中, 材料的损伤会表现出局部化特性, 细化的网格会使得能量消耗减少, 这意味着网格尺寸大小将影响数值结果。因此, 采用Bažant and Oh^[17]研究的裂纹带模型来减少数值结果对单元尺寸的依赖性, 也就是将单元特征长度加到损伤演化模型的公式中。每个单元消耗能量可用(10)式表示

(10)

式中: G_I是失效模式I的断裂能, 也是材料属性之一; σ_eq^f和ε_eq^f是失效时等效应力和等效应变; l_c是单元特征长度, 在本文中使用单元体积的立方根作为单元特征长度。

失效模式I的损伤变量d_I可用(11)式表示

(11)

式中

(12)

(13)

(14)

式中: δ_{eq, I}表示失效模式I的等效位移; σ_{eq, I}表示失效模式I的等效应力; 上角标0和f分别表示损伤起始时刻与最终失效时刻; F_I表示失效模式I的损伤起始值。由(11)式可以看出, 损伤变量的值取决于等效位移δ_{eq, I}。

4种失效模式的等效位移以及等效应力的计算方法如下所示:

纤维拉伸失效模式

(15)

(16)

纤维压缩失效模式

(17)

(18)

基体拉伸失效模式

(19)

(20)

基体压缩失效模式

(21)

(22)

式中, 〈〉为Macaulay算子。通过以上公式可以看出, 单元的等效位移和等效应力均与单元特征长度以及局部应力应变有关。

1.4 内聚区模型

分层损伤是复合材料层合板的一种典型损伤模式。本文将零厚度的内聚单元引入有限元模型, 用来模拟复合材料层合板中相邻层之间的分层损伤。

从图 1中可以看出, 内聚单元遵循双线性牵引力-分离法则, 其中下角标n, s, t分别代表牵引力的法线方向和面内剪切力的第一、第二方向; 上角标0表示单元损伤的初始时刻, f表示最终失效时刻。在牵引力达到峰值之前, 内聚元素的初始响应被假设为线性关系。如果单元形变完全沿法线方向或者第一、第二剪切方向, 达到A点之后, 单元将开始产生损伤。

在实际应用中, 复合材料结构不仅仅承受单一方向的载荷, 而是受到多个方向的载荷, 也就是混合模式加载。图 2是混合模式加载下的内聚区模型, 其中阴影三角形表示内聚单元在混合模式加载下的响应; 2个无阴影的三角形分别表示内聚单元沿纯法线方向的响应和沿纯剪切方向的响应; A点和B点代表混合模式加载下内聚单元损伤的开始时刻与完全失效时刻。因此, 采用二次方破坏准则作为分层损伤起始判据, 如下所示

(23)

式中: t_n, t_s和t_t分别为法向牵引力和第一、第二切向牵引力; N和S分别表示法向强度和切向强度。当满足起始判据时, 复合材料层合板开始分层, 引入B-K准则计算混合模式加载下的能量耗散值G^C, 如下所示

(24)

式中: G_n^C和G_s^C分别为沿法线方向的断裂能和沿剪切方向的断裂能; G_S为剪切方向的能量耗散值; G_T为总共的能量耗散值; η为B-K准则中相关系数, 一般取值为1.45。

2 有限元模型

在ABAQUS中, 建立了T700CG/M21碳纤维环氧树脂基复合材料层合板在低速冲击下的有限元模型。该模型为界观尺度模型, 将CFRP层合板的每个单向纤维层板作为最小研究单位, 并不考虑纤维取向或纤维编织方式。使用VUMAT子程序编写相应的失效准则与损伤演化方法, 用来模拟并分析该层合板在低速冲击下的力学响应以及损伤情况。

2.1 几何尺寸与边界条件

根据Hongkarnjanakul^[15]的试验, 建立了T700CG/M21碳纤维环氧树脂基复合材料层合板的有限元模型, 铺层顺序为[0₂, 45₂, 90₂, -45₂]_S, 几何尺寸为150 mm×100 mm×4 mm。使用局部坐标系定义材料铺层方向。冲击器为头部直径16 mm的半球型圆柱体, 质量2 kg, 沿Z轴方向初始速度5 m/s, 得到冲击能量25 J。冲击器外表面刚好接触到层合板上表面, 距离为0 mm。将复合材料层合板放置在中心镂空125 mm×75 mm的支座上, 并且约束了层合板4个角的Z轴方向位移, 来模拟试验的边界条件。图 3展示了T700CG/M21碳纤维环氧树脂基复合材料层合板的低速冲击有限元模型。

2.2 单元选择与网格划分

在文献[21]的试验中, 使用2层同向纤维的层板是为了增加同一纤维方向的层板厚度, 在无损检测时更容易观测到不同铺层角度层板间的分层损伤, 因而并未记录同向纤维层板间的分层情况。为了提高计算效率, 本文建立的有限元模型中将相邻并且纤维方向相同的2层单向板合并为1个层板。复合材料层合板的有限元模型可以分为2个区域: 采用精密网格的冲击区域和采用较为粗糙网格的其他区域。采用单元尺寸为1.2 mm×0.9 mm×0.5 mm的带有减缩积分的八节点实体单元C3D8R对72 mm×36 mm的冲击区域划分网格。选取松弛刚度的沙漏控制方法以避免零能量模式。此外, 每层间插入零厚度的八节点内聚单元COH3D8以模拟分层损伤, 与相邻的实体单元共享节点来确保位移的连续性。采用刚体来模拟冲击器与支座。总的来说, 本文的复合材料层合板低速冲击有限元模型由70 440个C3D8R单元和61 635个COH3D8单元组成。

2.3 接触设置与材料属性

考虑到试验过程中, 冲击器在25 J的低速冲击能量下不能完全穿透层合板, 因此应用ABAQUS的通用接触算法来模拟复杂的接触过程。设置硬接触与摩擦因数(0.3), 分别描述法向与切向接触属性。T700GC/M21复合材料参数采用文献[21]中的数值模拟参数, 其中包括单向板与层间界面的材料参数, 列于表 1。

3 数值结果与分析 3.1 模型验证

为了验证模型的正确性, 将数值仿真结果与Hongkarnjanakul^[21]的试验结果进行比较。对比结果表明了数值模型预测值与试验测量值吻合良好。图 4是在25 J冲击能量下, 试验测量与数值预测得到的T700CG/M21复合材料层合板接触力-时间曲线。从图 4中可以看出, 在冲击阶段的初期, 接触力-时间曲线呈现出震荡上升的趋势, 这可能是由于冲击层合板后振动造成的。随着时间的推移, 接触力逐渐增大, 损伤也逐步积累。当接触力接近峰值时, 会产生高频率震荡, 并且接触力增加的速度也逐渐缓慢。峰值过后, 冲击器开始反弹, 接触力逐渐减小, 直至0值, 也就是冲击器完全离开复合材料层合板表面。仿真预测得到接触力的波动过程反映了渐进损伤模型中不同损伤模式的萌生与扩展, 它与试验记录值吻合较好, 除了峰值略低于试验值, 需要稍长一些时间回到0值。

图 5是T700CG/M21复合材料层合板接触力-时间曲线。显而易见得是, 数值模型预测的最大位移略大于试验测量的位移。当冲击速度变为0值时, 复合材料层合板达到最大中心位移处。此时, 层合板中心位移远大于层压板的厚度, 这是由于层压板受冲击导致弯曲变形所造成的。最后, 当冲击器回弹并离开层合板表面时, 会在层合板表面留下凹坑, 使得中心位移并没有回到0值。

表 2给出了模型预测的冲击响应与试验数据的误差分析。从中可以看出，模型对最大接触力与冲击持续时间的预测较为准确，误差均在7.8%左右；预测最大中心位移时，仿真结果略大于试验数据，误差为12.82%；而在预测最终中心位移时，模型表现得更加良好，误差为6.12%。考虑到复合材料试验数据具有一定的分散性，可以认为本文建立的模型对碳纤维增强复合材料低速冲击响应的预测结果合理有效。

图 6是复合材料层合板通过C扫描与数值模型预测得到的分层损伤情况。仿真结果与试验结果的每层间分层损伤分布基本相似，除了预测损伤面积略大于C扫描面积。各层间具体情况将在下节详细讨论。综上所述，本文建立的有限元模型与相应的渐进损伤模型能够较好地预测复合材料层合板受低速冲击载荷作用下的力学响应与分层情况。

3.2 损伤分析

低速冲击载荷作用下的复合材料层合板会出现纤维断裂、基体开裂和分层等典型损伤模式。尽管维修人员几乎不能通过肉眼观察出这类损伤，但是它们对复合材料结构的机械性能有重要的影响。因此，损伤分析就显得格外重要。图 7是预测的T700CG/M21复合材料层合板损伤分布，包括了各层的纤维拉伸损伤、基体拉伸/压缩损伤与分层。从图中可以得到，每层板中出现纤维拉伸损伤的单元很少，基本位于板中央。从仿真结果中发现并没有单元发生纤维压缩失效。随着冲击过程的进展，层合板中心位置产生了沿冲击方向的位移，纤维与基体承受拉伸应力，引起拉伸损伤。在冲头与层合板的接触区域，少许纤维断裂使得接触力出现震荡，材料刚度退化程度加大，导致层合板承载能力下降。基体损伤是低速冲击损伤的主要模式之一。从图 7中可以看出每一层都有大量单元出现了基体拉伸损伤，尤其是靠近底部两层。很明显，这2层的损伤面积远大于其他层板。层合板受到冲击产生弯曲变形，使得底层先满足拉伸失效，并逐渐向受冲击面扩展。底层0°板的基体拉伸损伤沿纤维方向扩展，其它层板并无明显扩展规律。至于基体压缩损伤，靠近受冲击的上2层板有明显损伤，其他层板损伤面积较小，甚至没有单元出现基体压缩破坏。冲击过程中，随着冲头与层合板表面的接触，靠近接触面的层板中基体处于压缩状态，会率先出现压缩损伤，并沿着纤维方向扩展。其他层板受弯曲影响，拉伸应力处于主导地位，导致基体拉伸损伤面积大于基体压缩损伤面积。

图 7显示了在25 J冲击能量下T700CG/M21复合材料层合板的分层情况。从图中可以看出，分层在每个层间都有发生。分层损伤从每层间中心向周围扩展，呈不规则形状分布，靠近冲击面的分层面积略小于背面分层面积。除了第一层，其余每层的损伤形状长轴方向与相邻下层板的纤维方向一致。预测分层损伤分布取决于单元尺寸和层间界面材料属性，在数值模型中更改参数会影响预测结果。从图 6可以观察出，数值仿真预测的分层面积要比C扫描面积略大。这说明研究人员使用该数值方法设计复合材料许用值时会得到更为保守的结果，能够更好地保证复合材料的安全性。

4 结论

1) 建立了一种适用于碳纤维增强复合材料层合板的低速冲击的渐进损伤模型。该模型使用Hashin和Hou失效准则作为损伤起始判据，结合等效位移法的线性退化方案进行刚度折减，采用内聚区模型来模拟层间损伤情况。对比试验结果，使用该模型预测的力学响应与损伤情况吻合较好，验证了模型的正确性与有效性。

2) 板中央有少量的纤维拉伸损伤，并没有发现纤维压缩失效；每层单向板，尤其是远离接触面的层板，受到层合板弯曲变形所产生的拉应力影响，出现不同程度的基体拉伸失效；冲头在厚度方向上的接触力使得靠近冲击面的层板出现基体压缩失效。

3) 每个层间都有不同程度的分层损伤，呈不规则形状分布；远离冲击面的分层面积略大于靠近冲击面的分层面积；多数分层损伤形状长轴方向与相邻下层板的纤维方向一致。

4) 给出碳纤维增强复合材料层合板的材料参数，使用该模型能够分析冲击力学响应，预测层内与层间的冲击损伤，预期可为研究其他条件的碳纤维增强复合材料层合板低速冲击问题提供一种数值模拟方案。