复合材料冲击损伤数值仿真模型评估

2019-12-31张嘉睿吴富强姚卫星

航空工程进展 2019年6期

张嘉睿，吴富强，姚卫星,

(1.南京航空航天大学机械结构力学及控制国家重点实验室，南京 210016) (2.南京航空航天大学飞行器先进设计技术国防重点学科实验室，南京 210016)

0 引言

纤维增强树脂基复合材料(FRP)已被广泛应用于航空航天等先进结构。FRP抗冲击能力较弱而易产生冲击损伤，对结构的承载能力和使用安全造成了严重威胁。因此，FRP低速冲击损伤问题受到了工程界和学术界的广泛关注。

目前,研究FRP低速冲击问题主要有实物试验和数值仿真两种方法。贾建东[1]从试验角度对冲击损伤进行了总结，综述了FRP层合板的损伤形式、损伤检测等方面的研究进展;张铱芬等[2]从数值仿真角度对破坏准则、冲击损伤产生和扩展进行了综述；张颖军等[3]从试验和数值方法角度对FRP层合板冲击损伤的表征进行了综述；崔海坡等[4]对FRP数值模拟模型进行了分类和介绍。

由于数值模型的复杂性，前人对于数值模型的分类方法有很多种，但都不尽完备。本文认为一个FRP低速冲击数值模型应包括两部分：第一部分描述了冲击物和被冲击FRP结构之间的相互作用及其响应，第二部分描述了FRP结构内部损伤变化。本文将按此分类，通过两个算例对相关模型进行评估。

1 冲击过程模型分类

FRP层合板与冲击物接触时，涉及冲击物、接触过程、被冲击物等问题的描述，因此冲击模型可以按不同角度分类。按外载荷的处理方式，可分为准静态模型和动态模型；按载荷-位移关系，可分为线性模型和非线性模型；按分析模型中是否利用了试验数据，可分为经验、半经验和全数值模型；按模型的维数可分为二维模型和三维模型。

本文认为数值模型是物理模型的映射，数值模型应该在物理机制上再现冲击损伤试验。从冲击过程看，冲击损伤数值模型应包含冲击载荷、冲击物模型、接触过程描述、控制方程求解等要素，每个要素都可以有不同的处理方法，如表1所示。

表1 冲击过程模型分类Table 1 Classification of impact process models

从理论上讲，上述四个要素的任意组合都可构成一个FRP结构冲击损伤数值仿真模型，即可能有54类模型。但常见的数值模型主要有以下四种，如表2所示。

表2 冲击过程典型模型Table 2 Typical models of impact process

从表2可以看出：冲击物目前鲜有采用柔性体模型的，主要原因是计算成本太高；准静态求解也主要是为了节省计算时间；这样处理对仿真精度有一定的影响。

近几年随着计算机运算能力的提高，动态有限元算法越来越多地被人们所采用，冲击物柔性体模型也开始有所尝试。

2 冲击载荷下材料损伤演化模型

冲击载荷下FRP结构状态变化的数值模型主要包括初始损伤场、材料本构、就位特征、破坏准则、损伤演化、网格效应六个要素，如表3所示。

表3 材料损伤演化数值模型要素和取值Table 3 Elements and values of numerical model for damage evolution of materials

上述六个要素的每个都可以有多种处理方法，不同组合可以构成数百类模型，对它们一一进行评估的工作量十分巨大。下面对这六个要素略作介绍，并指出目前较主流的组合方法。

(1) 初始损伤

FRP组分材料本身就具有一定的分散性和不确定性。同时，在FRP结构的加工成型制造过程中，会出现气泡、夹杂、基体富聚、孔隙、基体和纤维界面粘接不匀称、纤维局部弯曲等工艺缺陷。因此FRP结构一定存在初始的物理损伤，目前的数值模型中均未直接考虑初始损伤问题，一般将其均匀化到宏观力学性能中，也有少量采用随机变量等不确定性方法来描述这种初始损伤。

(2) 材料本构

大量研究表明单层FRP的纵向与横向应力-应变呈现近线性关系，而剪切应力-应变呈现明显的非线性特征，因此考虑FRP的非线性本构关系是有必要的。目前处理FRP剪切非线性行为的模型主要有单系数模型和三系数模型。H.T.Hahn等[26]单系数模型为

(1)

式中：σij和εij分别为剪应力和剪应变，ij=13,23;β为材料剪切非线性系数；Gij为起始剪切模量。

文献[27]对式(1)进行了稳定性分析，并引入折减系数dnon对剪切应力应变进行线性化处理

(2)

在纤维间损伤出现前，用折减系数dnon对剪切割线模量在起始剪切模量基础上进行折减。损伤出现后，剪切模量退化在采用dnon退化后的剪切模量上进行。

D.Feng等[28]和段永照等[29]采用三系数剪切非线性模型：

(3)

图1 三系数剪切非线性模型Fig.1 Three-coefficient shear nonlinear model

在有限元中用来描述分层的主要方法有两种，一种是基于断裂力学的直接分析方法，该方法认为能量释放率分量组合大于或等于临界能量释放率时发生分层，虽然该方法能够预测分层的起始情况，但是对于分层的扩展需要通过网格重构才能实现。另一种方法是粘聚区模型(CZM)，该方法不仅能够预测新裂纹的萌生，同时还可以避免在裂纹扩展过程中的网格重构，是目前使用较多的一种分层模拟方法。

粘聚区模型的本构关系不考虑轴向耦合效应，层间本构关系为

(4)

式中:Knn、Kss和Ktt为界面元在三个方向上的刚度；δn、δs和δt为界面元在三个方向上的相对位移；tn、ts和tt分别为界面元法向应力及两个方向上的剪应力。

(3) 就位特征

对FRP就位效应的研究主要集中在就位层厚度效应和邻近层约束效应两个方面，处理就位效应目前有三种方法：能量准则法，经验公式法和统计分析方法。F.K.Chang等[30]提出经验公式确定各子层的横向就位拉伸强度YT,is和面内就位剪切强度S12,is，

(5)

式中：YT和S12分别为单层板的横向拉伸强度与面内剪切强度；M为子层的铺层数；Δθa和Δθb分别为子层与上下层的夹角；A、B、C和D为材料的常数。

S.T.Pinho等[31]提出了薄板就位强度计算公式：

(6)

(4) 破坏准则

在层内破坏准则方面已有诸多研究人员进行了深入研究，其中Chang-Chang准则[32]，Hashin准则[33]和Puck准则[34]受到广泛的认可并被采用在各种有限元模型中，来进行FRP单层板损伤起始的判断。

Chang-Chang准则为

(7)

Hashin准则为

(8)

A.Puck等通过大量破坏试验发现：对于FRP存在一个平行于纤维方向的潜在断裂面，如图2所示，断裂面是特定应力状态下最容易发生失效的作用面，当危险应力系数fE大于1时，纤维间失效发生。

(9)

图2 子层应力状态与断裂面Fig.2 Sublayer stress state and fracture surface

在层间起始损伤判定准则方面，常见的准则为ABAQUS自带的基于应力的二次准则：

(10)

式中：N和S分别为层间界面的拉伸与剪切强度；tn、ts和tt分别为粘结元法向应力与两个方向上的剪应力。

(5) 损伤演化

FRP单层损伤演化可分为瞬间卸载模型、逐渐卸载模型和承载力不变模型，如图3所示。目前大多采用基于系数折减的瞬间卸载模型与基于能量释放率的逐渐卸载模型。

图3 FRP损伤演化模型Fig.3 FRP damage evolution model

层间损伤演化方式目前大多采用混合模式的B-K能量释放率来预测分层发生后的扩展，当材料的释放率满足损伤扩展准则时材料完全失效。

(11)

式中：GⅠC、GⅡC、GⅢC为Ⅰ、Ⅱ、Ⅲ型临界能量释放率；η为B-K准则系数。

(6) 网格效应

在仿真模型中，损伤的最小尺度是有限元的单元。通过引入单元的特征长度和修改单元失效应变，可以降低有限元分析过程中能量释放率对网格的依赖性，从而削弱数值计算结果对网格尺度的依赖[35]。若在单轴加载下的材料满足线性软化行为，则修改后的单元材料失效应变不再维持不变，因为

(12)

式中：GC为临界能量释放率；gC为单元体在单位特征长度内的临界应变能密度;Ri为损伤起始时的最大应力；l*为单元特征长度。

(13)

式中：ly、lz为单元在整体坐标系中沿坐标轴y、z方向的尺寸。

特征单元尺寸如图4所示，常见的数值模型搭配如表4所示。

图4 特征单元尺寸Fig.4 Characteristic unit size

表4 材料损伤演化典型模型Table 4 Typical models for damage evolution of materials

3 评估

由于目前冲击模型的种类很多，仅对典型的冲击模型进行分析评估。本文将整个冲击模型分为冲击物和被冲击FRP结构之间相互作用及其响应，以及FRP结构内部损伤变化。对两个部分分别选取典型的方法开展分析与评估。

3.1 模型宏观评价

冲击过程模型中Collombet模型被广泛认可，刘洪权模型比较新，这两个模型的主要区别为冲击载荷的施加方式不同且控制方程的解法不同。考虑到冲击载荷还有一种集中力施加方式，本文增加一种组合方式：将刘洪权模型的位移施加方式改为最大接触力施加方式。在数值分析中，采用ABAQUS中的冲击接触来描述冲头与FRP层合板的动态冲击接触。

在FRP材料损伤演化模型中，Camanho模型和Faggiani模型[14]被认可程度较高。段永照等[29]在Faggiani模型基础上采用三系数非线性本构进行了完善。本文选用Camanho和段永照这两个模型，并同时考虑网格效应和就位效应，具体模型形式和参数详见文献[36]和文献[29]。

将冲击过程模型与材料损伤演化模型组合，构成六种模型，如表5所示。

表5 评估采用模型Table 5 Assessment adoption model

3.2 评估算例一：正交铺设层合板

3.2.1 原始数据

本算例以H.Y.Choi等[6]的试验为原始数据。材料为T300/976，层合板铺层为[454/-454]S，材料属性如表6所示[6]。试验件的有效冲击尺寸为100 mm×76 mm×2.3 mm。冲击能量为2.78 J，冲击后分层损伤如图5所示[6]。

表6 T300/976材料属性Table 6 Material properties of T300/976

图5 分层损伤结果Fig.5 Delamination damage results

3.2.2 模型与结果

整个模型分为冲头与FRP层合板两部分。冲头部分将钢质冲头简化为解析刚体，半径为6.35 mm，质量为0.16 kg。FRP层合板部分，按照试验条件对层合板两端进行固支约束，单层板采用8节点减缩积分实体单元(C3D8R)进行模拟，厚度方向上将每四个子层离散为一个单元，在不同角度铺层之间插入零厚度界面元来模拟层间分层，损伤演化采用B-K准则，有限元图如图6所示，层间界面元参数如表7所示[73]。模型分层结果如图7所示，分层面积具体尺寸如表8所示。

图6 算例一有限元模型Fig.6 Finite element model of case 1 表7 层间界面元参数(算例一)Table 7 Interlayer interface element parameters(case 1)

参数数值参数数值GⅠC/(J·m-2)230Knn/(GPa·mm-1)156(GⅡC=GⅢC)/(J·m-2)450(Kss=Ktt)/(GPa·mm-1)70N/MPa45η1.45(S=T)/MPa115

(a) 模型Ⅰ

(b) 模型Ⅱ

(d) 模型Ⅳ

(e) 模型Ⅴ

(f) 模型Ⅵ图7 各模型仿真结果(算例一)Fig.7 Simulation results of each model(case 1)

表8 分层损伤特征尺寸Table 8 Layered damage characteristic size

3.2.3 讨论

从图5和图7可以看出：各仿真模型得到的分层形状都接近花生状，损伤形状的长轴与该界面下方铺层方向一致，这与试验观测到的现象符合；Puck准则材料损伤演化模型得到的分层形状较Hashin准则的材料损伤演化模型更接近于试验，这是因为Puck准则材料损伤演化模型考虑了剪切非线性的影响，但是从计算效率方面，由于Puck准则需要寻找隐藏断裂面，在计算时间上明显高于Hashin准则材料损伤演化模型，Hashin材料损伤演化模型更适用于工程实践；两种材料损伤演化模型均可以较准确的预测冲击后的分层面积。

在冲击过程模型方面，采用位移输入的冲击过程模型得到的分层面积与试验相比存在损伤面积偏小的情况，这是由于位移输入忽略了冲头的惯性效应以及冲头回弹的影响，导致损伤面积偏小。采用最大接触力作为输入时，未考虑在冲击过程中冲击力随时间的变化，得到的分层面积与试验值相比偏大。

3.3 评估算例二：角铺设层合板

3.3.1 原始数据

以D.Feng等[28]的试验为原始数据。材料为HS300/ET223，层合板铺层为[03/454/-454]S，材料属性如表9所示[28]。试验件的有效冲击尺寸为65 mm×87.5 mm×3.2 mm。冲击能量为4 J，冲击后分层损伤如图8所示[28]。

表9 HS300/ET223材料属性Table 9 Material properties of HS300/ET223

图8 试验结果(算例一)Fig.8 The test result(case 1)

3.3.2 模型与结果

整个模型分为冲头与FRP层合板两部分。冲头部分将钢质冲头简化为解析刚体，直径为12.5 mm，质量为2.34 kg。FRP层合板部分，按照试验条件对矩形开口简支夹具固支，单层板采用8节点减缩积分实体单元(C3D8R)进行模拟，厚度方向上将每一个子层离散为一个单元，在不同角度铺层之间插入零厚度界面元来模拟层间分层，损伤演化采用B-K准则，有限元图如图9所示，层间界面元参数如表10所示[28]。模型分层结果如图10所示，分层面积具体尺寸如表11所示。

图9 算例二有限元模型Fig.9 Finite element model of case 2 表10 层间界面元参数(算例二)Table 10 Interlayer interface element parameters(case 2)

参数数值参数数值GⅠC/(J·m-2)520Knn/(GPa·mm-1)120(GⅡC=GⅢC)/(J·m-2)920(Kss=Ktt)/(GPa·mm-1)43N/MPa30η1.45(S=T)/MPa40

(a) 模型Ⅰ

(d) 模型Ⅳ

(e) 模型Ⅴ

(f) 模型Ⅵ图10 各模型仿真结果(算例二)Fig.10 The test result(case 2)

表11 特征尺寸Table 11 Characteristic size

3.3.3 讨论

从图8和图10可以看出：各仿真模型得到的分层形状都接近花生状，损伤形状的长轴与该界面下方铺层方向一致，这与试验观测到的现象符合；但是Hashin准则得到的分层面积较Puck准则明显偏大，对于角铺设层合板，考虑剪切非线性的Puck准则更为合适；计算效率方面，由于Puck准则需要寻找隐藏断裂面，在计算时间上明显高于Hashin准则材料损伤演化模型。