复合材料层板冰雹高速冲击损伤预测及失效分析

张超，方鑫，刘建春

（江苏大学机械工程学院，镇江 212013）

碳纤维增强复合材料层板具有比强度高、比模量大、可设计性强等性能优势，在现代飞机结构中应用非常广泛。飞机在飞行、停靠期间常会遇到冰雹等恶劣天气，冰雹冲击属于软体冲击事件，对复合材料层板往往会产生肉眼难以观察的隐形损伤，如层间分层、基体开裂等，这对飞机关键部位的结构承载非常危险。因此，需要深入开展碳纤维增强复合材料层板冰雹高速冲击力学响应及损伤机理研究。

复合材料层板冰雹冲击是一个高度非线性的瞬态过程，冲击过程中涉及接触分析、冰雹破裂、层板损伤等问题，研究方法主要是实验研究和数值模拟。实验研究成本高、周期长且受限于具体的实验条件；而有限元数值模拟可以较好地解决冰雹冲击这一瞬态动力学问题，详细获取冲击过程中冰雹破裂、流动特性和层板损伤演化过程，具有边界适应性强、求解精度及求解效率高等优点，已成为研究复合材料高速冲击问题的首选方法［1-2］。

在实验研究方面，学者们通过开展冰球高速冲击复合材料层板实验［3-6］，发现冰球冲击动能和冲击力峰值之间呈线性关系，层间损伤能量阈值和冰雹直径与层板厚度的比值呈线性关系，并利用微观扫描技术获得复合材料层板的损伤形貌及破坏阈值。数值模拟方面，冰雹数值模型主要采用传统的拉格朗日法和光滑粒子流体动力学（smoothed particle hydrodynamics，SPH）法。早期，学者Kim和Kedward［3］使用拉格朗日网格并采用弹塑性模型对冰雹进行建模，模型中未考虑冰雹的应变率效应及压力与体积的非线性关系，在冰雹破碎后的物理特征模拟方面不够理想。SPH法可以有效解决高速冲击下网格大变形的问题，冰雹单元失效后转化为粒子，可以有效再现冰破碎、流动的整个力学行为。

Coles等［7］对金属弹丸和冰球冲击复合材料层板进行了数值对比研究。Pernas-Sanchez等［8］采用SPH法建立了冰雹高速冲击复合材料层板的有限元模型，采用界面接触（cohesive contact）来模拟层间分层行为，预测了不同冰雹直径和层板厚度下分层面积和冲击速度之间的定量关系。文献［5，9-10］在冰雹冲击复合材料层板有限元模拟中，采用内聚力界面单元（cohesive element）结合Traction-Separation本构关系模拟层间分层现象，发现层板的分层开裂损伤形貌与实验结果较为吻合，说明界面接触和界面单元均可以有效模拟冰雹冲击下复合材料层板的层间分层现象，但界面接触只是一种接触属性，不是材料属性，没有创建额外单元代表界面这种组分实体，难以直接表征材料的界面力学行为。

张晓晴等［11］建立了冰雹冲击复合材料加筋结构的有限元模型，预测了冲击位置、冲击能量及冲击入射角对复合材料加筋壁板冲击动力学响应的影响。Dolati等［12］从能量吸收和破坏模式角度，对复合材料／波纹铝芯夹芯板开展了冰雹冲击数值模拟，讨论了冰雹弹形状、面板堆叠顺序和加固类型对夹芯板冲击响应的影响。

综上，国内外学者对冰雹冲击复合材料结构开展了一定的实验和模拟工作。但现有数值模型在冰雹破碎后的物理行为模拟方面还不理想，对冰雹高速冲击下纤维增强复合材料应变率效应的认识不足，相关的数值模拟和分析工作还需进一步深入。本文建立冰雹高速冲击复合材料层板有限元模型，模拟冰雹高速冲击复合材料层板的瞬态过程，分析材料的损伤失效机理并探讨冲击速度、角度对层板损伤性能的影响。

1 损伤本构模型

1.1 冰雹本构模型

在高速冲击条件下，冰的弹塑性模型与应变率具有明显的相关性，低应变率下冰材料表现为韧性，而高应变率下表现为脆性。其中，冰的拉伸性能受应变率影响不大，而压缩屈服强度受应变率影响较大［13］。基于Tippmann［14］建立的考虑应变率效应的弹塑性模型，通过定义张力失效准则来描述冰雹破裂后的流动特性。达到拉伸失效应力时，偏应力分量设为零且材料只承受压缩应力，失效后单元不删除。冰雹本构的偏应力失效类型设置为脆性，压应力失效类型设置为韧性。

引入水的Gruneisen状态方程来控制冰破碎后压力与体积之间的非线性关系，该方程描述了压力、密度和能量之间的关系。

压缩状态时（μ＞0），冰雹模型的状态方程可以表示为［12］

1.2 单向复合材料损伤本构模型

碳纤维增强复合材料由纤维和基体2种组分材料构成。纤维通常被认为是线弹性、横观各向同性材料，可以用线弹性单元来描述；基体可以看作为黏弹性、均匀各向同性材料，可以用黏弹性单元来描述。它们的应力-应变关系分别为

式中：Ef1为纤维轴向弹性模量；Vf为单向复合材料的纤维体积含量。

对于冰雹高速冲击下的复合材料层板，层内具有纤维拉伸、纤维压缩、基体拉伸、基体压缩4种典型失效模式。三维Hashin失效准则［16］能够很好地预测这4种失效模式，在此作为复合材料单层的失效判据，具体表达式如下。

纤维拉伸失效（σ11≥0）：

式中：S0为当前应变率下的强度；S0r为参考应变率˙ε0下的强度；ξ为应变率修正系数。

上述4种典型失效模式贯穿于冰雹高速冲击复合材料层板的整个过程。当复合材料层板出现局部破坏后，层板不会直接失去全部承载能力，而是材料局部性能出现下降。因此，当复合材料层板开始发生局部损伤后，根据损伤模式的不同，将纤维和基体的弹性常数分别进行相应的折减，折减方案如表1所示。

表1 材料刚度折减方案Table 1 Material stiffness reduction scheme

1.3 层间损伤模型

为了模拟冰雹冲击下复合材料层板的层间分层现象，引入界面单元结合双线性内聚力模型来描述界面的力学行为。

界面初始为线弹性阶段，本构关系为

2 有限元模型及验证

2.1 冰雹冲击复合材料层板有限元模型

为了分析冰雹高速冲击下复合材料层板的损伤特性和失效机理，选取文献［6］的实验结果作为对比验证数据，复合材料层板夹持方式为四周固支，通过拧紧螺栓压紧底板和压板对实验件起到稳固作用（见图1（a））。窗口尺寸为200 mm×200 mm，压板和底板尺寸为300 mm×300 mm，层板尺寸为250 mm×250 mm。复合材料层板的铺层方式为正交铺层［0°／90°］4S和准各向同性铺层［0°／45°／90°／-45°］2S，共16层，每层厚度为0.125 mm。冰雹直径为25.4 mm，采用具有沙漏控制的八节点C3D8R 单元离散，单元尺寸为1 mm。达到失效应力时自动转化为SPH 粒子。复合材料层板每层均采用C3D8R单元离散，并且建立局部坐标系来定义每层纤维的局部方向。

图1 冰雹高速冲击复合材料层板实验装置及有限元模型Fig.1 Experimental device and finite element model of composite laminates under high-velocity hailstone impact

为模拟冰雹冲击下层板分层现象，在复合材料各层之间插入零厚度内聚力界面单元COH3D8。在冰雹高速冲击过程中，层板的损伤主要集中在冲击中心区域。在层板网格划分时，在中心冲击区域细化网格，单元尺寸为1.74 mm×1.74 mm。从中心向外网格密度由密到疏，采用了bias功能从内到外布置种子，比例为2，且该网格尺度下计算结果收敛，这样既保证了计算精度又提高了运算效率。

冰雹和层板之间采用通用接触，接触属性为无摩擦的硬接触。层板与压板、底板之间建立绑定约束，压板和底板均为离散刚体，并通过固定压板与底板的参考点来模拟实验中的夹持条件。所建立的冰雹冲击复合材料层板有限元模型如图1（b）所示，该模型共由165 331个节点、6 176个R3D4单元、87 200个C3D8R单元和67 500个COH3D8单元组成。

水的Gruneisen状态方程参数、冰雹基本力学性能参数及基于应变率的冰雹屈服强度参数如表2［12］、表3［18］和表4［14］所示，界面单元材料参数和单向复合材料的材料参数如表5［19-20］和表6［6，15，21］所示。

表2 Gruneisen状态方程参数［12］Table 2 Parameters of Gruneisen state equation［12］

表3 冰雹基本力学性能参数［18］Table 3 Basic mechanical performance parameters of hailstone［18］

表4 基于应变率的冰雹屈服强度［14］Table 4 Strain rate based yield strengths of hailstone［14］

表5 界面单元材料参数Table 5 Material parameters of interface element［19-20］

表6 单向复合材料的材料参数［6，15，21］Table 6 Material parameters of unidirectional composite［6，15，21］

2.2 有限元模型验证

为了验证所建有限元模型的准确性和可靠性，将冰雹高速冲击下2种不同铺层方式层板的数值模拟结果与文献［6］中编号为P010101和P020101的实验结果进行对比。

图2为铺层方式［0°／90°］4S的层板在冰雹以206 m／s速度冲击下宏观损伤的实验和模拟对比情况。从图2（a）、（c）可以发现，层板正面出现少量纤维断裂和基体开裂现象；从图2（b）、（d）可以发现，层板背面出现鼓包现象。实验与模拟的宏观损伤形貌比较吻合，初步验证了所建有限元模型的有效性。

图2 ［0°／90°］4S层板实验［6］与模拟宏观损伤对比Fig.2 Comparison of macroscopic damage of［0°／90°］4S laminates between experiment［6］and simulation

图3为铺层方式［0°／90°］4S的层板在冰雹以206 m／s速度冲击后层间分层的C扫描和数值模拟结果对比。选中心200 mm×200 mm区域进行查看，如图3（a）所示，中间为发生分层区域；如图3（b）所示，其为每一层的分层区域投影轮廓。可以发现，分层主要从撞击中心处沿0°纤维方向展开，故层板在206 m／s冰雹冲击下损伤主要是层间分层这种肉眼难以观察到的内部损伤。

图3 ［0°／90°］4S层板分层情况实验［6］与模拟对比Fig.3 Comparison of delamination of［0°／90°］4S laminates between experiment［6］and simulation

图4为铺层方式［0°／45°／90°／-45°］2S的准各向异性层板在冰雹206 m／s速度冲击后层间分层的C扫描和数值模拟结果对比。可以看出，其分层现象较正交铺层更为严重，主要从冲击中心沿着45°、-45°方向扩展。整体的分层形状为梅花形，与实验较为吻合，进一步验证了该有限元模型的有效性。

图4 ［0°／45°／90°／-45°］2S层板分层情况实验［6］与模拟对比Fig.4 Comparison of delamination of［0°／45°／90°／-45°］2S laminates between experiment［6］and simulation

3 数值结果分析与讨论

3.1 冰雹冲击过程分析

冰雹高速冲击复合材料层板是一个瞬时的非线性动态过程，图5为初速度206 m／s冰雹冲击［0°／90°］4S复合材料层板的动态过程。结合层板中心处位移历程图可以看出，冰雹高速冲击层板过程中层板变形主要分为3个阶段：①冲击初期，层板的快速弯曲变形阶段（0～0.2 ms）；②达到最大位移后，层板的快速回弹阶段（0.2～0.4 ms）；③冰雹大面积沿层板扩散的回弹阶段（0.4～0.6 ms）。

图5 206 m／s冰雹冲击速度下［0°／90°］4S复合材料层板动态过程Fig.5 Dynamic process of［0°／90°］4S composite laminates under hailstone impact at velocity of 206 m／s

由图5可以看出，该有限元模型很好地描述了冰雹的破碎流动特征，冲击过程中，应力波主要沿纤维铺层方向和层板厚度方向扩散而形成明显的圆形鼓包。冲击初期，层板中心处快速弯曲变形，t＝0.2 ms时，层板背面中心处产生最大拉应力而出现局部纤维拉伸断裂现象，并伴随中心区域基体开裂，层间分层面积快速增长；然后，冰雹完全失效破碎，以粒子形式沿层板四周继续冲击，冲击力急剧减小，层板快速回弹变形，层间分层面积逐渐增加到最大，层内出现了不连续的基体拉伸损伤；冲击后期，冰雹继续以粒子形式沿层板大范围扩散，该阶段层板缓慢回弹变形，层间分层面积基本不再变化，层内不连续的基体拉伸损伤继续小幅扩展。冲击结束后，复合材料层板呈轻微凸起状态，出现不可恢复的凹坑。

3.2 复合材料层板损伤分析

3.2.1 层间损伤

冰雹等类柔性体冲击复合材料层板时，与其他损伤模式相比，层间损伤最为严重，这是考虑冲击后层板剩余力学性能的一个重要因素。本节选取正交铺层和准各向同性层板为研究对象，在冰雹冲击速度为206 m／s时进行层间分层损伤分析。

图6为2种不同铺层方式层板分层面积随冲击时间的变化。数值模拟中，当界面的损伤变量达到0.97时，认为界面单元完全损坏即发生层间分层。可以发现，层间分层大面积扩散主要集中于冲击的起始阶段（0～0.27 ms），在0～0.2 ms之间，分层扩展速率较快，且2种铺层方式层板分层总面积随时间变化趋势较为一致。冲击结束后，正交铺层层板的分层总面积约为36 506 mm2，准各向异性层板分层总面积约为38 008 mm2，准各向异性层板分层总面积稍大于正交铺层层板。

图6 两种典型层板分层面积随时间变化历程Fig.6 Evolution of delamination area of two typical laminates with time

3.2.2 层内损伤

本节选取正交铺层层板为研究对象，在冰雹冲击速度为206 m／s时进行层内损伤分析。图7和图8为冲击后复合材料层板中8个单层基体拉伸开裂和纤维拉伸断裂的损伤分布情况，其中，白色部分为纤维和基体均完全失效，相关单元已从模型中删除。可以发现，在冰雹高速冲击过程中，基体开裂失效面积较大，基体开裂区域主要沿着层板中心以矩形向外扩展；而纤维拉伸断裂失效分布于层板中心处，面积相对较小。纤维主要沿0°和90°方向发生断裂，部分层纤维沿45°和-45°方向发生断裂。随着冲击的进行，纤维损伤面积不再扩展，但由于冰雹冲击过程中的流动特征，在冲击的后段每层的边缘处继续产生了一些基体拉伸开裂现象。

图7 206 m／s冰雹冲击速度下层板典型层基体拉伸损伤分布Fig.7 Distribution of matrix tensile damage in typical layers of laminates under hailstone impact at velocity of 206 m／s

图8 206 m／s冰雹冲击速度下层板典型层纤维断裂损伤分布Fig.8 Distribution of fiber breaking damage in typical layers of laminates under hailstone impact at velocity of 206 m／s

3.3 冲击速度、冲击角度对复合材料损伤的影响

为了详细研究不同冰雹冲击载荷工况下碳纤维复合材料层板的力学响应和损伤特性，本节选取25.4 mm直径冰雹，以不同冲击速度、不同冲击角度冲击［0°／90°］4S正交铺层层板，开展了多组数值模拟工作。

图9为不同冲击速度、冲击角度下层板所受到的最大冲击力变化曲线。可以发现，最大冲击力在冲击角度为45°、60°时随着冲击速度的增大呈线性增大趋势。在同一冲击速度下，最大冲击力随冲击角度的增大而减小，很大程度上是由于在带有偏角的冲击条件下，层板所受的法向压应力小于垂直冲击条件下，同时受到一部分切应力作用。冲击角度为0°（垂直冲击层板）时，在160 m／s和200 m／s的冲击速度下层板所受到的最大冲击力比较接近，此时层板也承载着最大的法向压应力及弯曲变形，一般认为层板发生纤维断裂现象后力学性能大幅降低，故可以将200 m／s冲击速度时对应的能量作为此层板可承受的最大冲击能量阈值（154.44 J）。同时，根据数值模拟得出的冲击力时间曲线可以发现，最大冲击力出现在冲击的初始阶段，在冲击角度为45°时，最大冲击力出现的时间最早；在同一冲击速度下，0°冲击角度下最大冲击力出现的最早。

图9 冲击速度和冲击角度对最大冲击力的影响Fig.9 Influence of impact velocity and impact angle on maximum impact force

图10为不同冲击速度、冲击角度下层板分层总面积的变化曲线。可以发现，随着冲击角度的增加，分层面积逐渐减小，即冰雹垂直冲击层板时分层总面积最大，法向冲击力是造成层间分层现象的主要因素。垂直冲击下，分层面积随着冲击速度的增大而近似呈现出线性增加趋势。其中冲击角度为0°和30°时，分层面积的差值相对较小，冲击角度从45°增加到60°时，分层总面积减小幅度较大。此外，在较低的冲击速度下（低于120 m／s），冲击角度对分层面积的影响较小，分层面积低于5 000 mm。

图10 冲击速度和冲击角度对分层面积的影响Fig.10 Influence of impact velocity and impact angle on delamination area

为了有效分析冰雹冲击速度对层板分层面积的影响，本文进行了冰雹垂直冲击层板情况下，不同冲击速度下的多组数值模拟工作。结果发现，层板分层面积在较高范围内（大于120 m／s）随冲击速度的增大而近似线性增大；当冲击速度为80 m／s时，层板开始出现层间分层损伤。由此可以推测，此种复合材料层板在直径25.4 mm冰雹的冲击下，分层损伤能量阈值约为24.71 J。

4 结 论

1）本文考虑了冰雹破裂后的流动特性和冰雹高速冲击下复合材料层板的应变率效应，建立了复合材料层板冰雹高速冲击有限元分析模型。模拟了冰雹高速冲击下复合材料层板的损伤过程，分析了材料的损伤特性和失效机理，通过与已有实验结果对比分析，验证了模型的有效性。

2）层间分层是冰雹高速冲击下层板的主要损伤形式。在25.4 mm冰雹200 m／s垂直冲击速度下，层板出现少量的纤维拉伸断裂，大面积基体裂纹及层间分层现象，分层总面积达到54 457 mm2。该冲击速度对应能量可以作为层板此种工况下可承受的最大冲击能量。冰雹80 m／s冲击速度时，层板开始出现分层损伤，故此复合材料层板的分层损伤能量阈值约为24.71 J。

3）层板分层总面积随着冲击角度的增加而减少，在200 m／s速度垂直冲击下，分层面积随冲击速度增大而线性增大，而在一定冲击角度下（45°、60°），最大冲击力与冲击速度呈线性关系。冲击速度较低时，冲击角度对分层总面积的影响较小。相同冲击条件下，准各向异性铺层层板分层面积大于正交铺层层板，分层损伤主要发生在冲击的初始阶段，冰雹破碎后，层板卸载回弹阶段分层面积基本不变。

4）所建立的有限元模型能够有效获取高速冲击过程中冰雹的破裂、流动特性，模拟复合材料层板的应变率效应和层内、层间损伤过程，且所采用的损伤模型和损伤准则具有良好的通用性，研究结果可以为其他复合材料结构冰雹冲击问题数值分析提供参考。