纤维增强树脂基复合材料（fibre reinforced plastic composites，FRP）具有比强度和比模量高、可设计性强、耐疲劳性好等特点，在诸多工程中得到了越来越广泛的应用^[1]。但由于FRP本身对低速冲击比较敏感，使得FRP层合板在受到外物撞击后很容易出现内部损伤。

目前对FRP层合板的数值模拟研究，重点在于对其本构关系的研究。在宏观尺度上，主要有2种本构模型：（1）脆性失效模型（brittle failure model，BFK）^[2-6]，即在未满足失效准则时材料的响应是线弹性的，一旦满足失效准则，材料即刻完全失效，表现为其应力、强度和模量立刻变为零；（2）连续介质损伤力学模型（continuum damage model，CDM）^[7-11]，即一旦满足损伤准则，材料并非立刻丧失承载能力，材料发生损伤，但随着损伤的发展，材料的刚度慢慢减小直到完全失去承载能力。Xin等^[12]提出了一个基于连续损伤介质力学的三维模型，该模型借鉴了Hou等^[2]提出的失效准则并通过引入断裂能的概念和利用应力-位移关系来描述材料的本构方程，从而在一定程度上减少了由于应变集中而造成的网格依赖性。模型考虑了多种破坏形式，如面内纤维的拉伸断裂、压缩破坏、基体开裂、层间剪切，厚度方向上的压溃失效以及分层。在失效准则中，还考虑了沿着厚度方向的剪应力和正应力对不同失效模式的影响，并认为层间分层破坏不仅发生在层间正应力为拉应力的情况下，而且当层间正应力为压应力但剪切应力足够大的情况下也有可能发生。一旦满足失效准则，便认为材料发生损伤。除此之外，该模型还通过引入动态增强因子在统一框架内考虑了FRP的应变率效应。

分层是FRP层合板在冲击荷载作用下主要的损伤模式之一，预测FRP层合板的分层也是目前研究的热点和难点。Olsson^[13]提出了一个分析模型来预测复合材料在受到低速冲击时的分层损伤的产生和发展。Espinosa等^[14]提出了一个针对GFRP(glass fiber reinforced plastic)层合板三维各向异性的黏弹性模型，并结合黏性接触理论，分析了GFRP在动态载荷下的分层情况。Aymerich等^[15]用黏性面单元来研究复合材料的分层并取得了很好的结果。Johnson等^[16]利用二维损伤模型和界面层单元研究了编织GFRP的动态响应。Singh等^[17]提出了一种三维的弹塑性本构模型并结合基于界面间的黏性区域的相互作用来模拟CFRP(carbon fiber reinforced plastic)层合板在低速冲击下的分层行为，模拟结果与实验吻合得较好。

本文中建立有限元模型来预测FRP层合板在低速冲击下的响应和破坏。模型不再将整个层合板当成连续均匀的各向异性材料，而是认为层合板是由单层板通过黏结层粘结而成。单层板的本构模型采用新近提出的渐进损伤本构关系^[12]，粘结层的本构模型采用了改进的黏性区域理论。最后，将有限元模型预测的结果与实验观察结果进行比较和讨论。

1.   本构模型

1.1   单层板连续介质损伤模型

Xin等^[12]提出了一个FRP层合板冲击载荷作用下的连续介质损伤模型，将单层板当成连续均匀的各向异性材料，并不具体区分基体和纤维。该模型包含3个部分即：失效准则、损伤演化、应变率效应。

当以上任何一种失效模式满足时，便认为材料受到不可逆的损伤，直到材料完全失效。为了表征材料的损伤和失效，通过损伤变量来描述材料的刚度折减因子。需要注意的是，单层板的连续损伤模型用基于应力-位移关系的断裂能方法来计算损伤变量。通过引入动态增强因子概念，考虑了材料模量和强度的应变率效应。具体的单层板的连续介质损伤模型的详细介绍可以参考文献[12]。

1.2   粘结层损伤模型

黏性区域方法（cohesive zone method, CZM）^[18-19]把基于应力的初始损伤判断准则和基于断裂力学的损伤演化规律结合起来。然而，黏性区域方法和其后相关文献[14, 17]都没有考虑应变率效应。在下面的章节中，较为详细地介绍一种考虑应变率效应的修正的黏性区域方法（modified cohesive zone method)。

黏性区域方法借鉴断裂力学的概念，认为分层的破坏有3种基本形式，即：张开型分层（Mode Ⅰ，open mode）、剪切型分层（Mode Ⅱ，shear mode）、撕裂型分层（Mode Ⅲ，tearing mode）。由于在实验中环氧树脂撕裂断裂韧度很难测量，因此在本文中将剪切型分层和撕裂型分层当成相同形式，将都等效为剪切型分层，推导出相关的理论解。

1.2.1   牵引力位移准则

如图1所示，牵引力-位移(traction-separation)准则包括3个部分，即：初始线性响应阶段、分层开始发生阶段、分层扩展阶段。

图  1  牵引力-位移准则示意图

Figure  1.  Schematic diagram of traction-separation law

下载: 全尺寸图片 幻灯片

在线性阶段，在不考虑耦合的情况下，可以把应力-牵引准则写成如下的形式：

式中，${t_{{n}}}$、${t_{{s}}}$、${t_{{t}}}$分别代表牵引t在沿着法向和2个切向的分量即${t_{{n}}}$、${t_{{s}}}$、${t_{{t}}}$分别为法向和2个切向的应力，${\delta _{{n}}}$、${\delta _{{s}}}$、${\delta _{{t}}}$是指在相应方向上分离的距离，${K_{{n}}}$、${K_{{s}}}$、${K_{{t}}}$则分别是在相应方向上的刚度值。

随着载荷的进一步加载，一旦应力场满足以下二次应力损伤准则，则认为初始的分层开始发生。二次应力损伤准则可以写成：

式中：$\tau _{{n}}^{\left( 0 \right)}$、$\tau _{{s}}^{\left( 0 \right)}$、$\tau _{{t}}^{\left( 0 \right)}$分别为法向和2个切向的极限应力。将Mode Ⅱ和Mode Ⅲ当成都是剪切分层，则式（2）可以重新写为：

式中：${t_{{\rm{shear}}}} = \sqrt {{t_{{s}}}^2 + {t_{{t}}}^2}$，$\tau _{{\rm{shear}}}^{\left( 0 \right)} = \tau _{{s}}^{\left( 0 \right)} = \tau _{{t}}^{\left( 0 \right)}$。其中＜ ＞ 是Macaulay括号，定义为：

等效位移${\delta _{\rm{m}}}$可以写成如下表达式:

式中：${\delta _{{\rm{shear}}}} = \sqrt {{\delta _{{s}}}^2 + {\delta _{{t}}}^2}$。当开裂位移$\delta _{{n}}^{\left( 0 \right)}$大于0时，可以定义参数$\beta$

将式（6）代入方程（5）中，得到：

进一步假设3种模式的刚度相同，即$K = {K_{{n}}} = {K_{{s}}} = {K_{{t}}}$，将${t_{{i}}} = {K_{{i}}}{\delta _{{i}}}$代入式（3）中，通过相关推导，并考虑$\delta _{{n}}^{\left( 0 \right)} {\text{≤}} 0$时的情况，可以得出在发生初始分层时对应的等效位移：

式中：$\delta _{{n}}^{\left( 0 \right)}$、$\delta _{{s}}^{\left( 0 \right)}$、$\delta _{{t}}^{\left( 0 \right)}$是对应牵引力极限值$\tau _{{n}}^{\left( 0 \right)}$、$\tau _{{s}}^{\left( 0 \right)}$、$\tau _{{t}}^{\left( 0 \right)}$时的位移，其可以写成

从式(8)可以清楚地看出，当$\beta = 0$时是单纯的拉伸，即Mode Ⅰ模式；而当$\beta \to \infty$时便是单纯的剪切模式。

当等效位移${\delta _{\rm{m}}}$达到损伤发生的临界位移$\delta _{\rm{m}}^{\left( 0 \right)}$时，便认为损伤开始发生了。此时，随着载荷的增加，便开始进入损伤演化阶段。B-K准则^[20]是应变率释放率的函数，认为B-K准则一旦满足，便发生了完全的损伤，分层便发生。B-K准则表达式如下：

式中：${G_{\rm{I}}}\text{、}{G_{{\rm{II}}}}\text{、}{G_{{\rm{III}}}}$分别是Mode Ⅰ、Mode Ⅱ和Mode Ⅲ的应变能释放率，${G_{{\rm{shear}}}} = {G_{\rm I}} + {G_{{\rm{II}}}}$是剪切模式的应变能释放率，${G_{\rm{T}}} = {G_{\rm I}} + {G_{{\rm{shear}}}}$是总的应变能释放率，${G_{{\rm{IC}}}}\text{、}{G_{{\rm{IIC}}}}$分别是Mode Ⅰ、Mode Ⅱ的断裂韧度，${G_{\rm{C}}}$是总体的断裂韧度，$\eta$是从混合模式弯曲实验中拟合得到的一个经验参数，本文中取1。

如图1所示，可以把断裂韧度${G_{\rm{C}}}$写成：

将等式（11）代入式（10），并且考虑$\delta _n^{\left( 0 \right)} {\text{≤}} 0$时的情况，得出最后完全分层的失效位移$\delta _{\rm{m}}^{\left( {\rm{f}} \right)}$，表达式为：

从式(11)中可以清楚地看出，当$\beta = 0$时是单纯的拉伸，即Mode Ⅰ模式；而当$\beta \to \infty$时便是单纯的剪切模式。

损伤的形式表述为单元刚度的退化。通过等效位移${\delta _{\rm{m}}}$，开始发生损伤的等效位移$\delta _{\rm{m}}^{\left( 0 \right)}$，破坏完全发生时的失效位移$\delta _{\rm{m}}^{\left( {\rm{f}} \right)}$，来定义一个单元的损伤情况，单元的刚度退化系数D表达式如下所示：

为了描述该损伤的不可逆性，引入一个状态变量$\delta _{{\rm{m}},\max }^{}$

式中：${\delta _{\rm{m}}}(T)$是当前分析步长的等效位移，${\delta _{\rm{m}}}(T - \Delta T)$是前一个分析步的等效位移。这样，就可以把牵引力-位移准则写成统一的形式：

1.2.2   应变率效应

与单层板连续介质损伤模型相类似，粘结层的应变率效应可表示为^[12]：

式中：$\psi$为动态增强因子，A、B、C是描述FRP层合板应变率效应的经验常数。动态的强度和模量等于静态的强度和模量乘以动态增强因子，其表达式如下：

式中：$\tau _{{i}}^{\left( 0 \right)}\text{、}K_{{i}}^{}$是动态的强度和模量，$\tau _{{i}}^{\left( {\rm{s}} \right)}\text{、}K_{{i}}^{\left( {\rm{s}} \right)}$是静态的强度和模量。

2.   有限元模型

FRP层合板的有限元模型包括两个部分，即单层板和粘结层，如图2所示。在实际建模中，先利用Truegrid软件生成单层板模型，然后再在ABAQUS/CAE中通过偏置方法在指定位置生成8节点零厚度的粘结层单元，粘结单元和实体单元是共节点的，相当于在单层板之间插入粘结单元，如图3所示。单层板用C3D8R单元来创建，而粘结层用COH3D8单元（cohesive单元）来创建。弹体和靶板间的接触是通用接触，其摩擦因数设为0.1。

图  2  FRP层合板组成示意图

Figure  2.  Schematic diagram of FRP laminates

下载: 全尺寸图片 幻灯片

图  3  将粘结单元插入实体单元示意图

Figure  3.  Illustration of cohesive elements inserted in solid elements

下载: 全尺寸图片 幻灯片

对于实体单元失效问题，认为当单元等效应变达到极限等效应变时单元便被删除，具体细节可以参考文献[12]。而对于粘结单元，通过在用户子程序中定义一个状态变量令其值等于粘结单元的损伤变量D的值，来控制单元的状态。当D=0时，表示单元未发生损伤，而当D=1时，认为单元完全破坏，单元将被删除。粘结单元被删除表示此处发生了分层破坏。

对于各种不同铺层方向的单层板，在总体坐标系下，对每层单层板来建立一个局部坐标系，来确定单层板的主方向，用来表征单层板的不同铺层方向。粘结单元和实体单元是共节点的，两者面内单元尺寸保持一致，在文献[7]中详细讨论了网格收敛性，结果表明单元面内尺寸为1 mm×1 mm时数值计算结果具有良好的精度，因此本文中撞击区域的单元面内尺寸都选为1 mm×1 mm。

从几何上来看，2个单层板通过零厚度的粘结层连结在一起，但从计算的角度来说，由于要将应变增量变成等效位移的增量，必须对粘结单元赋予一个厚度。在这里，计算厚度取单位长度，粘结单元的计算厚度的选取关系到粘结单元的刚度参数的单位的选取。在这里，粘结单元的计算厚度取1，则粘结单元的刚度的单位与模量的刚度一致。

将粘结层当成一种与环氧树脂相接近的材料，其相关材料性质如断裂韧度、模量等都可以通过设计相关实验测出，并认为粘结层的应变率效应与单层板的应变率效应保持一致。表1给出了粘结层的相关材料参数。其中${A_{\rm{i}}}\text{、}{B_{\rm{i}}}\text{、}{C_{\rm{i}}}$是强度的应变率参数${A_{\rm{m}}}\text{、}{B_{\rm{m}}}\text{、}{C_{\rm{m}}}$是模量的应变率参数。

表  1  模型中用到的粘结层参数

Table  1.  cohesive interface properties used in present model

材料参数	数值		材料参数	数值		材料参数	数值
${\tau _{n}^{\left( {\rm{s}} \right)}}$	11 MPa		${\tau _{s}^{\left( {\rm{s}} \right)}}$	17 MPa		${\tau _{t}^{\left( {\rm{s}} \right)}}$	17 MPa
${G_{{\rm{IC}}}}$	0.3 N/mm[17]		${G_{{\rm{IIC}}}}$	0.8 N/mm[17]		${G_{{\rm{IIIC}}}}$	0.8 N/mm[17]
${K_{n}^{\left( {\rm{s}} \right)}}$	850 MPa[21]		${K_{s}^{\left( {\rm{s}} \right)}}$	850 MPa[21]		${K_{t}^{\left( {\rm{s}} \right)}}$	850 MPa[21]
${A_{\rm{i}}}$	2.5[22]		${B_{\rm{i}}}$	0.9[22]		${C_{\rm{i}}}$	3.7[22]
${A_{\rm{m}}}$	1.85[22]		${B_{\rm{m}}}$	0.5[22]		${C_{\rm{m}}}$	1.3[22]

下载: 导出CSV 

| 显示表格

如图4所示，建立了Aymerich等^[15]的低速撞击实验的有限元模型。根据实验中的描述，建立了尺寸为65 mm×87.5 mm×2 mm的CFRP（carbon fiber reinforced plastic）层合板模型，并将落锤系统简化为质量为2.3 kg，直径为12.5 mm的半球头刚弹。模型中，在开口区域45 mm×67.5 mm处单元面内尺寸约为1 mm×1 mm，在开口区域外单元尺寸约为2 mm×2 mm，中间用一个蝴蝶网格进行过渡。在厚度方向上，层合板的实际铺层方式为${[{0_3}{90_3}]_{\rm{s}}}$，如图4所示，层合板从上到下依次是纤维取向为${0^{\circ} }$单层板3层、${90^{\circ} }$单层板3层、${90^{\circ} }$单层板3层、${0^{\circ} }$单层板3层，下标s代表对称。其中0°方向沿着长边的方向。在实际建模中，将相同铺层方向的三层板当作单层板处理，每层单层板厚度为0.5 mm，模型中的铺层方式为$[0{\rm{/}}90{\rm{/}}90{\rm{/}}0]$。在不同铺层方向的单层板之间插入零厚度的粘结单元，因此共有2层粘结单元，如图4所示。实验中，层合板是放在开口区域为45 mm×67.5 mm的垫板上，并用夹具将4个角固定住。

图  4  低速冲击下CFRP层合板的有限元模型

Figure  4.  Finite element model for CFRP laminates under low velocity impact

下载: 全尺寸图片 幻灯片

3.   结果和讨论

利用前面介绍的模型，来对落锤低速撞击CFRP层合板的实验进行数值模拟，并于相关实验进行了比较。落锤分别以2.1 J（半球头弹撞击速度为1.35 m/s）和4.9 J（半球头弹撞击速度为2.06 m/s）的能量撞击CFRP板。单层板的相关参数如表2所示。其中E是弹性模量，参数下标1、2、3分别表示面内纵向，面内横向和面外方向；G是剪切模量，参数下标1、2、3分别表示面内纵向，面内横向和面外方向；${X_{\rm{t}}}$、${Y_{\rm{t}}}$、${Z_{\rm{t}}}$分别表示1、2、3方向的拉伸强度；${X_{\rm{c}}}$、${Y_{\rm{c}}}$、${Z_{\rm{c}}}$分别表示1、2、3方向的压缩强度；$S$是剪切强度，参数下标1、2、3分别表示面内纵向，面内横向和面外方向；$\nu$是泊松比，参数下标1、2、3分别表示面内纵向，面内横向和面外方向；${{\delta }_{1}^{\left( \rm{f} \right)}}\text{、}{{\delta }_{2}^{\left( \rm{f} \right)}}\text{、}{{\delta }_{3}^{\left( \rm{f} \right)}}\text{、}{{\delta }_{4}^{{{\left( \rm{f} \right)}}}}\text{、}{{\delta }_{5}^{\left( \rm{f} \right)}}\text{、}{{\delta }_{5}^{\left( \rm{f} \right)}}\text{、}{{\delta }_{7}^{\left( \rm{f} \right)}}$分别表示不同失效准则对应的断裂位移。

表  2  Graphite/epoxy单层板的参数值

Table  2.  Parameters for graphite/epoxy laminates

材料参数	数值		材料参数	数值		材料参数	数值
${{E}_{11}}$	143.4 GPa[17]		${Y_{\rm{t}}}$	54 MPa[17]		${{\nu }_{23}}$.	0.52[17]
${{E}_{22}}$	9.27 GPa[17]		${{X}_{\rm{c}}}$	1 650 MPa[17]		${{\nu }_{13}},{{\nu }_{12}}$	0.31[17]
${{E}_{33}}$	9.27 GPa[17]		${{Y}_{\rm{c}}}$	240 MPa[17]		${{\delta }_{1}^{\left( \rm{f} \right)}},\rm{ }{{\delta }_{2}^{\left( \rm{f} \right)}}$	5.7×10−2 mm[12]
${G_{12}}$	3.8 GPa[17]		${Z_{\rm{t}}}$	54 MPa[17]		${\delta _7^{\left( {\rm{f}} \right)}}$	2.7×10−2 mm[12]
${G_{23}}$	3.2 GPa[17]		${Z_{\rm{c}}}$	240 MPa[17]		${\delta _3^{\left( {\rm{f}} \right)}},{\rm{ }}{\delta _4^{\left( {\rm{f}} \right)}}$	7.65×10−2 mm[12]
${G_{31}}$	3.8 GPa[17]		${S_{12}}$	100 MPa[17]		${\delta _5^{\left( {\rm{f}} \right)}}$	5×10−3 mm[12]
${X_{\rm{t}}}$	2 945 MPa[17]		${S_{31}},{S_{23}}$	100 MPa[17]		${\delta _6^{\left( {\rm{f}} \right)}}$	5×10−2 mm[12]

下载: 导出CSV 

| 显示表格

图5给出2种不同低速落锤撞击能量下的数值模拟得到的载荷-时间曲线与实验观察结果的比较。从图5可以看出，与实验结果相比，数值模拟结果的峰值载荷和作用时间都吻合得很好。

图  5  数值模拟得到的载荷-时间历程与实验观察^[15]的比较

Figure  5.  Comparison of the numerically predicted load-time histories with the experimental observations^[15]

下载: 全尺寸图片 幻灯片

图6给出2种能量下的数值模拟得到的载荷-位移曲线与实验观察结果的比较。从图6可以看出，与实验结果相比，数值模拟预测的中心最大位移和最终残余位移都吻合得很好。

图  6  数值模拟得到的载荷-位移曲线与实验观察^[15]的比较

Figure  6.  Comparison of the numerically predicted load-displacement cures with the experimental observations^[15]

下载: 全尺寸图片 幻灯片

图7给出了graphite/epoxy层合板在不同能量（2.1 J、4.9 J）打击下底部[90/0]两层单层板之间分层形貌的数值预测结果与实验观察的比较。在图7(a)和(b)中，左边是数值模拟得到的结果，通过[90/0]中间的粘结单元的破坏来表征分层形貌；右边是实验中通过X光技术扫描得到底部[90/0]间形貌，阴影部分表明发生了分层破坏。界面层的粘结单元的损伤通过损伤因子D来控制。当损伤因子D=0时，说明粘结单元未有任何损伤，云图中颜色为蓝色；而当损伤因子D=1时，说明粘结单元已经完全失效，说明此处发生了分层破坏；云图中为中间过渡的颜色表明粘结单元受到损伤却还未完全失效。不难看出，本模型得到的分层的形貌和面积与实验观察吻合得较好，在[90/0]层间分层的形貌是双火焰形状，分层的长边都是沿着底层单层板纤维的铺层方向。

图  7  数值模拟得到的分层形貌与实验结果^[15]的比较

Figure  7.  Comparison of the numerically predicted delamination areas with the experimental observations^[15]

下载: 全尺寸图片 幻灯片

通过与实验的观察比较，验证了本文模型的可靠性，本文的模型能较好地预测CFRP层合板在落锤低速撞击下的动态响应，并成功地预测了CFRP层合板在横向冲击载荷下的分层行为。本模型考虑了应变率效应，使得模型能够处理动态载荷问题，相比不考虑应变率效应，对于不同能量冲击层合板，本模型计算得到的分层损伤结果与实验结果更加吻合。另外，由于考虑了应变率效应，该模型不仅可用于低速撞击问题而且还可用于中、高速撞击问题。

4.   结束语

本文建立了一个有限元模型来预测纤维增强树脂基复合材料（FRP）层合板在低速冲击下的响应和破坏。模型不再将整个层合板当成连续均匀的各项异性材料，而是认为层合板是由单层板通过粘结层粘结而成。提出考虑应变率效应的基于粘性区域方法的粘结层的损伤模型，并将提出的本构方程和失效准则通过子程序VUMAT嵌入到有限元软件ABAQUS/EXPLICIT中；将FRP层合板当成由单层板通过粘结层粘结而成，并建立了有限元模型；在低速撞击时，有限元模型能够很好地预测FRP层合板的载荷-时间历程、载荷-位移曲线和分层效应。