先进复合材料因其优异的性能，在航空航天等领域得到了广泛应用。但在低速冲击载荷作用下，复合材料结构内部很容易出现目视不可检的基体裂纹、分层等损伤，会使层合板的强度削弱35%~40%^[1]，严重影响结构的安全使用。因此，复合材料结构的冲击损伤问题一直是国内外学者研究的一个重要内容。

考虑到冲击试验耗时长、成本高，并且无法呈现冲击损伤的萌生与扩展的全过程，国内外许多学者采用有限元技术研究低速冲击损伤。特别是近年来，基于连续介质损伤力学^[2-3]的渐进损伤分析方法被广泛用于复合材料低速冲击损伤的模拟。Donadon等^[4]结合与各种失效模式相关的能量耗散描述损伤演化过程，提出了一种基于连续损伤力学方法的三维损伤模型，并用于模拟复合材料层合板的动态响应。刘向民等^[5]通过引入再平衡次数m来协调冲击接触时间与冲击波的传播时间，从而考虑冲击损伤分析过程中的“链锁破坏”。李念等^[3]建立了一个将多个铺层损伤用单个实体单元模拟的三维模型，该方法降低了有限元模型的规模，提高了计算效率。以上大部分学者的模型未考虑层合板中子层的面内剪切非线性，但相关分析^[6-7]表明剪切非线性对基体损伤的起始判定及损伤演化和分层损伤有很大的影响。

基于以上考虑，本文基于连续介质损伤力学提出了一种纤维增强复合材料(Fiber reinforced polymer/plastic, FRP)结构冲击损伤预测的三维动力学分析数值模型，包含非线性剪应力应变关系和归一化的混合模式基体损伤演化，来仿真层合板低速冲击过程。模型区分了纤维拉伸/压缩失效、纤维间拉伸/压缩失效4种层内损伤以及层间分层损伤，失效判定时考虑了就位效应对强度的影响。模型中加入单元特征长度以考虑网格依赖性。对层内损伤起始及演化的模拟在ABAQUS软件平台上通过自编材料子程序VUMAT实现，对于分层损伤采用ABAQUS内嵌的cohesive单元进行模拟，完成了相关层合板低速冲击损伤响应参数的预测。

1 冲击损伤模型

连续介质损伤力学方法包括损伤表征、损伤判定、损伤演化三部分。对应复合材料不同的损伤模式，损伤表征通过损伤状态变量的引入从而建立损伤材料与完好材料间的本构关系，由失效准则判定复合材料的损伤萌生。对冲击后层合板截面的微观形貌研究表明，在层合板厚度方向上存在许多明显的基体裂纹及断裂面，基体裂纹的扩展是导致其所处子层与相邻子层间产生分层的根本原因，因此如何在损伤模型中模拟基体裂纹的产生及不同断裂面角度对材料损伤的影响是合理模拟冲击伤的关键。基于单层板破坏机理的Puck失效准则^[8]考虑了断裂面角度对基体损伤的影响并能确定损伤后的断裂面角度值，该准则对大量复合材料典型试验^[9-10]的预测结果良好。损伤演化过程实际上是应变能释放的过程，当材料点应变能释放密度等于其断裂能密度时，意味着该材料点完全失效。在应变能释放过程中，材料会软化，宏观表现为弹性模量的退化和承载能力的下降。常见分析中，材料软化形式有线性^{[2, 11-12]}和指数^[13-14]形式。早期材料的损伤演化是在材料坐标系中，通过引入损伤张量对材料的刚度进行退化。当利用Puck准则作为失效准则时，断裂面的出现使以往的方法不再适用。因此本文提出了一种在断裂面坐标系进行损伤演化的方法。

为了准确地模拟层内损伤，在冲击损伤模型中将层合板的每个子层等效为正交各向异性体，采用3D实体单元模拟。通过ABAQUS动力学求解器求解应变增量，利用VUMAT子程序对材料进行损伤分析与损伤演化，并将计算得到的应力张量传回求解器继续运算。模型的总体框架见图 1，模型的关键要素在1.1~1.4节中阐述。

1.1 本构模型

在连续介质损伤力学中，将微观上不连续的缺陷或者损伤看作是宏观的连续损伤变量，用其描述材料的损伤程度。类似于文献[4]中的说明和讨论，损伤连续均匀地分布在有限单元体中。有效应力张量即完好应力张量σ，与真实应力张量σ通过损伤张量D进行关联，即采用损伤张量D直接退化本构关系式中的应力分量

$ \mathit{\boldsymbol{\bar \sigma }} = \mathit{\boldsymbol{D\sigma }} $

(1)

式中完好应力张量σ通过将层合板的每个子层等效为正交各向异性体，根据经典层合板理论求得。依据Puck准则，纤维间断裂面(图 2)与纤维平行，由此可知纤维损伤变量不会随断裂角的变化而变化。也就是说, 可以将损伤张量D分为独立的纤维损伤变量和基体损伤变量两部分来分别求解。

由于剪切非线性的存在, 本文将基体损伤具体分为两个阶段：

(1) 由Puck准则判断未出现纤维间损伤，但此时在宏观上应力-应变表现为剪切非线性，本文将此非线性当量为线性与“唯象损伤”之和，其损伤张量为D=diag[1-d₁₁ 1 1 1-d₁₂ 1-d₂₃ 1-d₃₁]，其中基体损伤变量d_ij，i≠j为唯象损伤，用来描述割线刚度的变化。

(2) 用Puck准则判断出现纤维间损伤时，由归一化的混合模式损伤变量d_m表征基体损伤，此时应力分量的退化在基体断裂面坐标系完成

$ \mathit{\boldsymbol{\bar \sigma }} = {\mathit{\boldsymbol{T}}_{ - {\theta _{{{\rm{f}}_{\rm{p}}}}}}}\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{T}}_{{\theta _{{{\rm{f}}_{\rm{p}}}}}}}\mathit{\boldsymbol{\sigma }} $

(2)

式中：θ_fp为断裂面角度, T_θfp为应力转轴坐标转换矩阵。断裂面坐标系损伤张量定义如下, D=diag[1-d₁₁ 1-d_n 1-d_t 1-d_ln 1-d_nt 1-d_lt]，其中d_n=d_ln=d_nt=d_m，d_t=d_lt=0。

1.2 纤维损伤判定及演化

本文采用最大应变准则(图 3)作为纤维失效的判定准则，考虑到载荷反向，刚度回复，纤维损伤d₁₁表示为

$ {d_{11}} = \max \left( {d_{11}^t\frac{{\left\langle {\sigma _{11}^t} \right\rangle }}{{\left| {\sigma _{11}^t} \right|}},d_{11}^c} \right) $

(3)

式中

$ \begin{array}{l} d_{11}^t\left( {{\varepsilon _{11}}} \right) = \frac{{\varepsilon _{11}^{ft}\left( {{\varepsilon _{11}} - \varepsilon _{11}^{0t}} \right)}}{{{\varepsilon _{11}}\left( {\varepsilon _{11}^{ft} - \varepsilon _{11}^{0t}} \right)}}\\ d_{11}^c\left( {{\varepsilon _{11}}} \right) = \frac{{\varepsilon _{11}^{fc}\left( {{\varepsilon _{11}} - \varepsilon _{11}^{0c}} \right)}}{{{\varepsilon _{11}}\left( {\varepsilon _{11}^{fc} - \varepsilon _{11}^{0c}} \right)}} \end{array} $

(4)

式中：d₁₁^t和d₁₁^c分别对应纤维拉伸和压缩损伤变量，〈x〉=max(0, x)为McCauley算子。

考虑到在损伤区域，压缩时碎块的相互作用，假设存在一很小的剩余强度等于基体压缩强度^[4]。

1.3 基体损伤起始及演化 1.3.1 剪切非线性

基体损伤起始与断裂面上的应力状态有关，包括线性的正应力-应变关系与非线性的剪应力-应变关系，见图 4。本文采用三次多项式拟合试验数据来描述剪应力的非线性，表达式为

$ {\tau _{ij}} = {c_1}\gamma _{ij}^3 - {\rm{sign}}\left( {{\gamma _{ij}}} \right){c_2}\gamma _{ij}^2 + {c_3}{\gamma _{ij}} $

(5)

式中c_i(i = 1, 2, 3)为待定系数。

本文将此非线性当量为线性与“唯象损伤”之和，在这一阶段损伤变量可定义为

$ {d_{ij}} = \alpha {\gamma _{ij}} $

(6)

式中：α为一材料常数；总剪应变γ_ij为弹性剪应变γ_ij^e和不可逆损伤应变γ_ij^d之和。在卸载时，τ_ij=G₀(1-d_ij)(γ_ij-γ_ij^d)，其中G₀取为三次曲线在原点的切线斜率。

1.3.2 基体损伤起始判据

基体失效，采用Puck准则^[8]判定。将子层绕着纤维方向旋转作180°划分，在θ∈[-90°, 90°]的角度范围内对应力危险系数进行搜索，寻找应力危险系数f_{E, IFF}最大时所对应的角度θ，此角度所对应的载荷作用面即为潜在纤维间失效断裂面。当f_{E, IFF}≥1时，认为基体裂纹出现。

考虑到层合板中子层的就位效应，按照Pinho等^[15]提出的公式确定各子层的横向就位拉伸强度Y_{T, is}与面内就位剪切强度S_{12, is}。

1.3.3 基体混合模式损伤演化

复合材料的损伤过程实质为应变能释放的过程(图 5)，材料点在损伤萌生后(点O)开始释放能量，在损伤过程中(点C)的应变能释放密度大小为红色三角形的面积，当总应变能释放密度(为唯象损伤的黑色部分面积与红色部分之和)等于弥散分布在特征长度l_m内的临界断裂能密度时，表示该材料点完全失效。

如上所述，完全失效时有

$ \frac{1}{2}\sigma _{\rm{r}}^0\left( {\varepsilon _{\rm{r}}^{\rm{f}} - \varepsilon _{\rm{r}}^0} \right) + {g_0} = \frac{{G_{\rm{r}}^{\rm{c}}}}{{{l_{\rm{m}}}}} $

(7)

由此可得, 完全失效时的应变

$ \varepsilon _{\rm{r}}^{\rm{f}} = \frac{2}{{\sigma _{\rm{r}}^0}}\left( {\frac{{G_{\rm{r}}^{\rm{c}}}}{{{l_{\rm{m}}}}} - {g_0}} \right) + \varepsilon _{\rm{r}}^0 $

(8)

式中σ_r⁰为损伤起始时的等效应变

$ \sigma _{\rm{r}}^0 = \sqrt {{{\left\langle {\sigma _n^0} \right\rangle }^2} + {{\left\langle {\tau _{nt}^0} \right\rangle }^2} + {{\left\langle {\tau _{nl}^0} \right\rangle }^2}} $

(9)

ε_r⁰为损伤起始时的等效应变

$ \varepsilon _{\rm{r}}^0 = \sqrt {{{\left\langle {\varepsilon _n^0} \right\rangle }^2} + {{\left\langle {\gamma _{nt}^0} \right\rangle }^2} + {{\left\langle {\gamma _{nl}^0} \right\rangle }^2}} $

(10)

式中：σ_n⁰，σ_nl⁰和σ_nt⁰为损伤起始时断裂面上的正应力与剪应力；ε_n⁰, γ_nt⁰和γ_nl⁰为相应的应变。

g₀为损伤起始时的等效体积应变能，由断裂面上的应力分量按式(11)计算

$ {g_0} = g_n^0{\left( {\frac{{\left\langle {\sigma _n^0} \right\rangle }}{{\sigma _{\rm{r}}^0}}} \right)^2} + g_{\mathit{ln}}^0{\left( {\frac{{\tau _{\mathit{ln}}^0}}{{\sigma _{\rm{r}}^0}}} \right)^2} + g_{nt}^0{\left( {\frac{{\tau _{nt}^0}}{{\sigma _{\rm{r}}^0}}} \right)^2} $

(11)

与各个应力分量相关的体积应变能为

$ g_i^0 = \int_0^{\varepsilon _i^0} {{\sigma _i}{\rm{d}}{\varepsilon _i}} \;\;\;\;\;i = n,\mathit{ln},nt $

(12)

G_r^c为混合模式的临界能量释放率，由式(13)给出

$ G_{\rm{r}}^{\rm{c}} = G_{{\rm{22}}}^{\rm{t}}{\left( {\frac{{\left\langle {\sigma _n^0} \right\rangle }}{{\sigma _{\rm{r}}^0}}} \right)^2} + G_{{\rm{12}}}^{\rm{c}}{\left( {\frac{{\tau _{\mathit{ln}}^0}}{{\sigma _{\rm{r}}^0}}} \right)^2} + G_{{\rm{23}}}^{\rm{c}}{\left( {\frac{{\tau _{nt}^0}}{{\sigma _{\rm{r}}^0}}} \right)^2} $

(13)

式中G_ij^c(t)(ij = 22, 12, 23)为各应力分量对应的临界能量释放率。

进而可以得到混合模式的损伤参量为

$ {d_{\rm{m}}} = \frac{{\varepsilon _{\rm{r}}^f - \varepsilon _{{\rm{r,d}}}^0}}{{\varepsilon _{\rm{r}}^{\rm{f}} - \varepsilon _{\rm{r}}^0}}\left( {\frac{{\varepsilon _{\rm{r}}^0 - {\varepsilon _{\rm{r}}}}}{{{\varepsilon _{\rm{r}}} - \varepsilon _{{\rm{r,d}}}^0}}} \right) $

(14)

基体损伤变量d_m为ε_r的函数，用来定义断裂面上的应力σ_n，σ_ln，σ_nt的退化。ε_r为当前增量步的等效应变，表达式为

$ {\varepsilon _{\rm{r}}} = \sqrt {{{\left\langle {{\varepsilon _n}} \right\rangle }^2} + {{\left( {{\gamma _{nt}}} \right)}^2} + {{\left( {{\gamma _{nl}}} \right)}^2}} $

(15)

ε_{r, d}⁰为不可逆等效损伤应变

$ \varepsilon _{{\rm{r,d}}}^0 = \sqrt {{{\left\langle {\varepsilon _{n{\rm{,d}}}^0} \right\rangle }^2} + {{\left( {\gamma _{nt,{\rm{d}}}^0} \right)}^2} + {{\left( {\gamma _{nl,{\rm{d}}}^0} \right)}^2}} $

(16)

式中：ε_{n, d}⁰, γ_{nt, d}⁰, γ_{nl, d}⁰为材料坐标系中, 损伤起始时弹性剪应变γ_ij^e和损伤应变γ_ij^d转换到断裂面上的分量。

断裂面上的应力按照以下方式进行退化

$ \left\{ \begin{array}{l} {\sigma _n} = {{\bar \sigma }_n} - {d_{\rm{m}}}\left\langle {{{\bar \sigma }_n}} \right\rangle \\ {\sigma _{\mathit{ln}}} = \left( {1 - {d_{\rm{m}}}} \right){{\bar \sigma }_{\mathit{ln}}}\\ {\sigma _{\mathit{nt}}} = \left( {1 - {d_{\rm{m}}}} \right){{\bar \sigma }_{\mathit{nt}}} \end{array} \right. $

(17)

退化完成后，将应力分量转换到材料坐标系得到损伤单元体的应力张量。

1.4 单元特征长度的计算

在ABAQUS商业有限元软件中，实体单元的特征长度默认等于该单元体积的三次方根。显然该方法未考虑单元长宽比，对于长宽高尺寸差别较大的单元，消除网格依赖性的效果有限。本文采用式(18)求得单元特征长度l_m

$ {l_{\rm{m}}} = \left\{ \begin{array}{l} {l_y}\cos {\theta _{{\rm{fp}}}}\;\;\;\;\;\;\;\left| {{\theta _{{\rm{fp}}}}} \right| < \theta \\ {l_z}\sin \left| {{\theta _{{\rm{fp}}}}} \right|\;\;\;\;\;\;\;\left| {{\theta _{{\rm{fp}}}}} \right| \ge \theta \end{array} \right. $

(18)

式中：θ_fp为基体断裂面角度；l_y，l_z为单元体在全局坐标系下沿坐标轴y，z方向的尺寸(图 6)。

2 算例分析

下面通过引用2个文献算例和作者做的1个试验来检验本文模型的有效性。

2.1 [45₄/-45₈/45₄]层合板低速冲击损伤 2.1.1 分析对象

试验源自文献[16]，模型根据试验时试件的真实夹持条件，将两端固支，有限元建模时层合板采用8节点减缩积分实体单元(C3D8R)模拟，厚度方向上将每4个子层离散为一个单元(图 7)。为减小计算时间，网格在冲击中心区域附近局部细化，局部细化网格大小为0.5 mm×0.5 mm，并采用沙漏控制技术。钢制冲头简化为解析刚体，质量为0.16 kg，冲击速度为5.89 m/s。考虑钢制冲头及夹具表面与试件表面的摩擦，滑动摩擦因数取为0.3。在不同角度铺层之间插入零厚度界面元来模拟层间分层，损伤演化采用B-K准则。模型所需材料参数见表 1。

2.1.2 结果与讨论

图 8给出了模型计算结果与试验结果的对比，其中图 8(a, b)为层合板分层损伤投影形貌的对比。可以看到，界面分层形状为花生状，损伤形状的长轴与该界面下方铺层方向一致。图 8(c)为基体背面由层合板弯曲引起的拉伸裂纹, 图 8(d)中间铺层由于剪切应力引起的裂纹，都与试验值吻合得很好。

表 2给出了分层长度、宽度及面积的试验值与本文计算仿真结果对比，同时给出并未考虑剪切非线性的仿真结果(文献[16])。相比之下，可以看出本文的仿真结果与试验值吻合得更好。

2.2 [0₃/45/-45]_S层合板低速冲击损伤 2.2.1 分析对象

试验源自文献[17]，模型根据试验时试件的真实夹持条件，将矩形开口简支夹具固支，同时限制冲头除冲击方向外的其余自由度，有限元建模时层合板采用8节点减缩积分实体单元(C3D8R)模拟，厚度方向上将每一子层离散为一个单元。为减小计算时间，网格在冲击中心区域附近局部细化，局部细化网格大小为0.5 mm×0.5 mm，并采用沙漏控制技术。考虑钢制冲头及夹具表面与试件表面的摩擦，滑动摩擦因数取为0.3。有限元模型见图 9，数值仿真模型所需材料参数见表 3。

2.2.2 结果与讨论

不同冲击能量下有限元预测与试验实测的冲击力-时间历程曲线见图 10。由图 10可以看出，对应于2，4，8J的冲击能量，试验测得的最大冲击力分别为2.3，3.1，4.7kN，预测值分别为2230，3194，4577N，有限元预测的误差分别为3.0%，3.0%，2.6%。两者吻合较好，尽管有限元仿真结果随着冲击能量的增加，冲击历程微长于试验值。同时冲击能量为8J时，可以看到达到最大载荷后，由于纤维断裂出现了突降现象。

图 11为不同冲击能量下层合板内部分层损伤投影形貌的对比。每一界面分层形状为花生状，损伤形状的长轴与该界面下方铺层方向一致。损伤特征与试验观测结果基本一致。

表 4给出了不同冲击能量下分层投影面积试验值与计算结果，表中同时给出了文献[17]的计算结果。可以看出在充分考虑剪切非线性之后，本文的计算结果与试验值吻合得更好。

2.3 碳纤维层合板低速冲击试验与损伤仿真分析 2.3.1 冲击试验

本文设计并完成了一个低速冲击损伤试验。试验件材料为U3160/5284RTM，铺层顺序为[45/-45/0₂/90/45/0₂/-45/45]₃，共30层，材料基本参数见文献[18]。试验件几何尺寸见图 12。

试验方法采用ASTM D7136，试验机为ZCJ9162全自动落锤冲击试验机。试验时将试件置于125 mm×75 mm的矩形开口简支支持夹具上，上部用盖板压紧。试验中通过改变落锤高度来调节冲击能量，落锤质量为5.5 kg，冲头为直径16 mm的半球形钢球，试件的分层面积由IUCS-Ⅱ型便携式C扫描系统测得。

2.3.2 仿真分析

模型根据试验时试件的真实夹持条件，将矩形开口简支夹具和上盖板固支，同时限制冲头除冲击方向外的其余自由度。为减少网格数量，网格划分将相同取向铺层离散为同一层网格，局部网格细化为1.5 mm×1.5 mm。在不同角度铺层之间插入零厚度界面元来模拟层间分层，损伤演化采用B-K准则，有限元模型见图 13，数值仿真模型所需材料参数见表 5。

2.3.3 试验结果与仿真计算结果对比

图 14为冲击能量18 J时，层间分层形状的预测值与试验结果的比较，图中截取的正方形区域大小为75 mm×70 mm。从中可见，有限元预测的分层投影形状与C扫描结果接近，分层投影形状大致呈圆形。表 6给出了损伤大小的详细尺寸比较，从表中可以看出仿真结果与试验值吻合得很好。

图 15依次给出了冲击能量为18J时，各界面层间分层损伤情况。每层分层形状大致为花生状，损伤形状的长轴大部分与该界面下方铺层方向一致。并且所有的分层都发生在单层厚度改变时，在单层厚度相同的±45°界面和90°/45°界面没有分层情况出现。

3 结束语

本文基于连续损伤力学提出了一种适用于复合材料结构低速冲击损伤的数值预测模型，模型考虑了层内损伤和层间损伤。模型包含非线性剪应力应变关系、就位效应、归一化的混合模式基体损伤演化和一种计算特征长度的方法。采用该模型研究了两种不同材料和铺层的复合材料层合板的低速冲击响应问题，得到与试验较为吻合的冲击载荷-时间历程、损伤形状和分层面积，验证了模型的有效性。