纤维增强树脂基复合材料作为新型高性能结构材料，由于其高比强度、高比模量及各向刚度及强度可设计性等优点，被广泛应用于土木工程、航空航天工程、车辆工程和海洋工程等领域，具有良好的应用前景。而大量试验研究^[1-5]表明，由于大部分基体材料具有塑性力学行为，大多数纤维增强树脂基复合材料层合板在横向压缩及剪切荷载下展现明显塑性及非线性特点，在剪切方向上尤为显著。

此外，许多学者通过试验表明^[6-12]，与静态荷载下的力学性能相比，纤维增强树脂基复合材料在高速冲击或动态荷载下的力学性能有显著提升。Hosur等^[10]通过Hopkinson压杆对碳纤维/环氧复合材料层合板在3种不同的应变率下进行压缩试验，结果表明，与静态试验结果相比，复合材料层合板的强度和刚度均有大幅增加，且展示出了显著的率相关性。Grimes^[11]通过试验发现AS4/PEEK碳纤维/聚醚醚酮复合材料层合板在8 ${\rm{s}}^{-1}$应变率下的强度和断裂应变比静态荷载下相应数值分别约增加42%和25%；Daniel等^[12]对AS4/3501-6碳纤维/环氧树脂复合材料在不同应变率下的弹性模量和破坏强度进行测量，结果表明，由基体材料主导的横向及面内剪切模量和破坏强度均与应变率呈正相关。这是由于组成纤维增强树脂基复合材料的基体材料多为黏弹性材料，反映黏弹性材料力学性质的一些物理量，如材料强度、杨氏模量、剪切模量和泊松比等对应变率具有显著的敏感性。因此，为准确反映复合材料的物理特点，基体材料主导方向上的材料强度及弹性模量应考虑应变率效应。

由于组成成分力学性能的显著差异，复合材料存在复杂失效模式，包括纤维断裂、基体开裂和分层损伤。为判断复合材料单层板的失效，大部分破坏准则采用材料的应力状态达到某个限值来表征，仅由一个损伤判据表征宏观失效，如Tsai-Hill准则^[13]、Hoffman准则^[14]和Tsai-Wu准则^[15]等。针对纤维增强树脂基复合材料损伤模式复杂的特点，Hashin^[16-17]用4个损伤判断准则区分纤维拉伸、纤维压缩、基体拉伸和基体压缩失效4种不同的破坏模式，但并未对引起材料破坏的机制进行探讨。

Puck等^[18-19]基于Mohr-Coulomb断裂理论有效应力的概念提出了Puck准则，合理解释了复合材料在横向压缩和面内剪切作用下的基体失效机制。该准则的难点在于确定基体失效断裂面的角度，该角度可通过杨凤祥等^[20]提出的分区反抛物线插值法进行求解。与传统求解方法相比，该方法具有更高的计算精度和效率。Davila^[21]和Pinho等^[22-26]基于Puck准则通过寻找材料断裂面确定损伤主轴的方法^[18-19]和纤维扭结理论^[27]提出LaRC系列准则，将单层板的失效模式区分为基体失效、纤维扭结/劈裂失效和纤维拉伸失效。通过计算扭结带(Kink band)内部应力水平判断是否发生扭结带内部的纤维扭结/劈裂失效。LaRC系列准则既可以确定基体破坏时的断裂面损伤主轴，又能确定纤维扭结破坏时的断裂面损伤主轴，与Puck准则相比更能解释纤维增强树脂基复合材料的损伤机制。

为准确描述复合材料的非线性或塑性力学行为，同时考虑损伤对力学属性的影响，大部分学者基于塑性理论和损伤力学建立了不同的本构模型。Giannadakis等^[28]提出了包含微损伤影响的非线性黏弹性和黏塑性模型，较好地预测拉伸荷载下的非线性剪切应力-应变曲线。陈静芬^[29-30]提出了二维弹塑性渐进损伤本构模型，采用Sun和Chen^[31]提出的正交各向异性塑性模型考虑复合材料平面内横向和剪切方向的塑性变形，预测结果与试验结果吻合良好。薛康等^[32]建立了三维弹塑性损伤模型，模型考虑了材料的拉压不对称性，采用LaRC05准则判断材料损伤起始，并对IM7/8552碳纤维/环氧树脂不同偏轴角度单向板压缩试验进行数值模拟验证了模型的有效性。但上述模型未考虑复合材料的应变率效应。Chen等^[33]建立了率相关的二维弹塑性渐进损伤本构模型，塑性模型采用了Thiruppukuzhi和Sun^[34]提出的黏塑性屈服准则，该屈服准则考虑了应变率效应和塑性效应，但由于所采用的损伤初始准则为Hashin准则，对复合材料基体失效损伤萌生物理机制描述不够准确。

针对现有模型存在的问题，本研究建立了同时考虑复合材料非线性力学响应、应变率效应和损伤累积导致材料属性退化的三维率相关弹塑性损伤本构模型。分别采用最大应力准则判断复合材料层内的纤维拉伸初始损伤，LaRC05准则判断纤维压缩扭结/劈裂初始损伤，Puck基体失效准则判断基体初始损伤，通过改进塑性力学模型表征材料动态荷载下的非线性力学行为。为准确预测复合材料在动态荷载下的弹塑性力学响应，引入了率相关放大系数对准静态下的塑性强化函数进行修正。采用Bazant等^[35]提出的“断裂带模型”对已开发的本构模型软化段进行规则化，以减轻有限元分析结果的网格敏感性。采用杨凤祥等^[20]提出的分区反抛物线插值法对基体损伤初始断裂面角度及纤维扭结/劈裂平面角度进行求解。编写了包含显式积分数值算法的用户自定义材料子程序VUMAT，并嵌于有限元程序ABAQUS v6.14中，可用于预测横向和剪切方向展现显著塑性效应及应变率效应的正交各向异性纤维增强树脂基复合材料层合板的破坏行为，并预测其破坏强度。

1.   三维率相关弹塑性损伤本构模型

模型包括描述未损伤与损伤材料的应力-应变关系、包含应变率效应的塑性模型和描述损伤初始及演化的损伤模型。其中，损伤模型包括：(1)损伤初始判断准则判断损伤是否萌生和在加载条件下是否进一步演化；(2)损伤演化法则描述损伤的大小。

1.1   应力-应变关系

当材料发生损伤后，承受荷载的截面有效面积由于微裂缝的产生和发展而逐渐减小，从而引起材料属性(如弹性模量和泊松比)的退化。在连续介质损伤力学中，通过在材料刚度张量中引入损伤变量d来表征损伤。在外荷载不变的单轴加载条件下，含损伤材料的名义应力σ与有效应力$\bar {\sigma }$之间的关系如下：

式中：P为沿材料截面法线方向的外力；${A}_{0}$为材料截面的初始面积；${A}_{\rm{eff}}$为材料截面的有效面积。式(1)表明，当材料损伤发生时，外力由材料的有效面积抵抗。

假设材料具有显著的塑性效应，满足小应变假设，则总应变张量${\boldsymbol{\varepsilon }}$可以表示为

式中：${{\boldsymbol{\varepsilon }}}^{\rm{e}}$为弹性应变张量；${{\boldsymbol{\varepsilon }}}^{\rm{p}}$为塑性应变张量。

损伤与未损伤的应力-应变关系表示为

式中，符号“∶”为对两个张量进行缩并运算；${\boldsymbol{\sigma }}$为名义应力张量；${\bar{\boldsymbol{\sigma }}}$为有效应力张量；${\boldsymbol{C}}\left(d\right)$为含损伤单向复合材料层合板的四阶刚度张量；${\boldsymbol{C}}$为未损伤单向复合材料层合板的四阶刚度张量。在有限元应用中，常采用Voigt形式表示，本文采用如下形式：

式中：

式中：${d}_{1}$、${d}_{2}$、${d}_{3}$分别为纤维方向、平面内垂直于纤维方向和平面外法线方向描述纤维损伤、基体损伤和分层损伤的损伤变量；${d}_{12}$、${d}_{13}$、${d}_{23}$为3个剪切方向的损伤变量；${E}_{11}^{0}$、${E}_{22}^{0}$、${E}_{33}^{0}$、${G}_{12}^{0}$、${G}_{13}^{0}$、${G}_{23}^{0}$、${v}_{12}$、${v}_{21}$、${v}_{13}$、${v}_{31}$、${v}_{23}$、${v}_{32}$为未损伤复合材料单层的弹性模量和泊松比。

针对纤维增强树脂基复合材料横向和剪切方向弹性模量和强度随应变率增加而增加的特点，本文根据Daniel等^[12]的研究对不同应变率下纤维增强复合材料的横向和剪切方向的强度和弹性模量进行修正：

其中：${E}_{ij}^{*}$、${F}_{ij}^{*}$为未损伤单向复合材料层合板的率相关弹性模量和率相关强度；${\dot{\varepsilon }}_{0}$为准静态下的应变率，${\dot{\varepsilon }}_{0}=1\times {10}^{-4}{{\rm{s}}}^{-1}$；$\dot{\varepsilon }$为当前状态下的应变率；${E}_{ij}^{0}$、${F}_{ij}^{0}$为未损伤单向复合材料层合板准静态下的弹性模量和强度；${m}_{\rm{e}}$、${m}_{\rm{s}}$为应变率对弹性模量及材料强度的影响系数，通过试验数据拟合确定。为保持刚度矩阵的对称性，泊松比需要根据应变率的影响进行修正：

本研究与基体相关的材料属性，即弹性模量$E_{22}$、$E_{33}$，剪切模量$G_{12}$、$G_{13}$及材料强度$Y_{\rm{T}}$、$Y_{\rm{C}}$和$S_{12}$，按式(5)考虑应变率效应；纤维方向的力学行为不考虑应变率效应，相应的弹性模量和强度不进行率相关修正，即${E}_{11}^{*}={E}_{11}^{0}$，${v}_{12}^{*}={v}_{12}^{}$，${v}_{13}^{*}={v}_{13}^{}$，${X}_{\rm{T}}^{*}={X}_{\rm{T}}^{}$，${X}_{\rm{C}}^{*}={X}_{\rm{C}}^{}$。对相关模量进行率相关修正后的应力-应变关系表示为

式中：${{\boldsymbol{C}}}^{{*}}\left(d\right)$为率相关修正后的含损伤单向复合材料层合板的刚度张量；${{\boldsymbol{C}}}^{{*}}$为率相关修正后的未损伤单向复合材料层合板的刚度张量。

本文复合材料塑性效应通过塑性模型考虑，模型包括屈服准则、流动法则和强化法则。

1.2   塑性模型

在含损伤的材料中，外力是由材料的有效面积抵抗的。因此，可以合理地认为塑性变形发生在材料的未损伤区域。塑性模型采用有效应力$\bar {\sigma }$表示，同样地，等效塑性应变${\tilde {\varepsilon }}^{\rm{p}}$、等效应力$\tilde {\sigma }$也是基于有效应力空间的概念建立的。

Sun和Chen^[31]提出的塑性屈服准则能有效模拟复合材料层合板准静态下的各向异性塑性行为，本研究采用该模型的等效形式：

式中：$F^{\rm{p}}\left(\bar{\boldsymbol{\sigma}}\right)$是塑性势函数，由各有效应力分量表示；$\tilde {\sigma }\left({\tilde {\varepsilon }}^{\rm{p}}\right)$为各向同性强化函数；同时定义$\tilde {\sigma }$为材料等效应力，其值由当前有效应力状态的塑性势函数计算得到，即$\tilde {\sigma }= F^{\rm{p}}\left(\bar{\boldsymbol{\sigma}}\right)$ 。本研究将Sun和Chen^[31]提出的塑性势函数拓展到三维状态下的表达式形式如下：

在纤维方向，复合材料单层并不展现显著的塑性行为，模型未予考虑；而加载时在横向和剪切方向塑性变形的不同发展水平则通过材料参数${a}_{66}$描述。

假设复合材料塑性流动满足相关联流动法则，则塑性应变增量分量可表示为

其中：${\Delta }\lambda$是整个塑性加载历史的非负塑性流动因子；梯度矢量$\dfrac{\partial F}{\partial {\bar {\sigma }}_{ij}}$为塑性梯度，描述了塑性应变增量分量${\rm{d}}{\varepsilon }_{ij}^{\rm{p}}$的方向。

联立式(8)~(10)，可得：

同样地，假设等效塑性应变符合相关联强化法则：

式中，${h}^{\rm{p}}$描述了等效塑性应变演化的方向，由式(8)和式(12)可得${h}^{\rm{p}}=1$。

采用相关联强化法则可保证单位体积塑性功率相等：

将式(9)与式(11)代入式(13)可得：

Sun和Chen^[31]针对准静态加载情况下提出的强化函数的等效形式如下^[30-31]：

式中：材料参数${a}_{66}$为通过试算法使不同偏轴角度层合板试件的等效应力$\tilde {\sigma }$和等效塑性应变${\tilde {\varepsilon }}^{\rm{p}}$试验强化曲线与纤维铺设角度为90°的相应试验强化曲线重合在一起得到的最优值；材料参数β和${n}_{\rm{p}}$为通过对单向复合材料层合板试件的偏轴试验应力-应变曲线获得加载方向的应力与塑性应变关系后进行对数线性回归分析获得^{[31, 36]}。

若直接采用该强化函数形式，在模拟不同加载速率条件下的复合材料力学行为时将不能得到准确的模拟结果。因此，本文参考Daniel等^[12]提出的弹性模量强化对数函数形式(见式(5))对该塑性强化函数进行了修正，引入了率相关放大系数，修正后的强化函数如下：

式中：${\xi }_{\rm{s}}$为率相关放大系数；${m}_{\rm{p}}$为应变率对材料塑性行为的影响系数，为了简单考虑，本研究令${m}_{\rm{p}}={m}_{\rm{e}}$。

1.3   损伤模型

1.3.1   损伤初始与发展判断准则

为预测单层内各种损伤模式的损伤初始及评估加载过程中的有效应力状态，采用应力危险系数表征各种模式的损伤初始：

式中：${FI}_{I}$为反映不同破坏模式的加载函数；${r}_{I}$为损伤阈值，初始值为0。当${r}_{I}=1$时，该式定义了材料的损伤初始，当${r}_{I}>1$时，该式定义了材料损伤发展。ft、fc和m分别代表复合材料的不同破坏模式，即纤维拉伸破坏、纤维压缩破坏和基体破坏。

本文采用三维Puck准则^[18]中的基体失效准则、最大应力准则和LaRC05准则^[26]中纤维压缩失效准则分别判断基体开裂、纤维拉伸和纤维压缩的损伤初始。

(1)基体破坏模式

纤维增强复合材料单向层合板在横向应力及面内剪切应力的作用下基体损伤初始断裂面会发生在平行于纤维的方向，基体后续的损伤失效行为将沿着这个断裂面扩展。基于Mohr-Coulomb断裂理论，Puck等^[18]采用基体损伤初始断裂面局部坐标系上的应力分量建立了描述复合材料单向板基体损伤初始的判断准则：

图1给出了作用在基体断裂面上的各个应力分量，其中1-2-3为单向板的材料坐标系，l-n-t为基体断裂面的局部坐标系，两个坐标系的第一方向重合；$\theta$为基体初始断裂面角度，其取值区间为[−90°，90°]；${\bar {\sigma }}_{\rm{n}}$和${\varepsilon }_{\rm{n}}$分别为基体断裂面上的法向应力和法向应变；${\bar {\tau }}_{\rm{l}}^{}$和${\gamma }_{\rm{l}}^{}$分别为基体断裂面平行于纤维方向的剪切应力和剪切应变；${\bar {\tau }}_{{\rm{t}}}$和${\gamma }_{{\rm{t}}}$分别为基体断裂面上垂直纤维方向的剪切应力和剪切应变。断裂面上应力、应变分量与基体初始断裂面角度θ间存在如下关系：

图  1  复合材料单向板基体断裂面上应力分量

Figure  1.  Stress components on matrix fracture surface of unidirectional composite

下载: 全尺寸图片 幻灯片

式中：${S}_{\rm{L}}$、${S}_{\rm{T}}$为基体断裂面上沿纤维方向和横向的剪切强度；${Y}_{\rm{T}}^{{*}}$为根据式(5)进行率相关修正后的材料横向拉伸强度；${\eta }_{\rm{L}}^{}$、${\eta }_{\rm{T}}^{}$为摩擦系数；${S}_{\rm{L}}$、${S}_{\rm{T}}$、${\eta }_{\rm{L}}^{}$和${\eta }_{\rm{T}}^{}$可按下式求得：

式中：${\alpha }_{0}$为材料在横向压缩破坏时的基体断裂面角度，对于碳纤维增强树脂基复合材料，已有文献^[19]表明${\alpha }_{0}$的取值区间为[51°，55°]，本文取53°；${Y}_{\rm{C}}^{*}$和${S}_{12}^{{*}}$为根据式(5)进行率相关修正后的材料横向压缩强度和面内剪切强度。

(2)纤维拉伸破坏(${\bar {\sigma }}_{11}\geqslant 0$)：

式中，${X}_{\rm{T}}$为单向板纤维方向拉伸强度。

(3)纤维压缩破坏(${\bar {\sigma }}_{11}<0)$：

Schultheisz和Waas^[37]通过试验发现，当受到纤维方向的压缩荷载时，纤维增强树脂基复合材料的局部纤维会发生偏折变形，并称之为“纤维扭结(Fiber kinking)”。Pinho等^[26]根据以上现象将纤维压缩破坏分为纤维扭结(Fiber kinking)和纤维劈裂(Fiber splitting)两种形式，并基于Mohr-Coulomb断裂理论提出了纤维受压损伤初始判据判断这2种纤维受压损伤形式的初始：

式中，${\bar {\tau }}_{23}^{\varphi }$、${\bar {\tau }}_{12}^{\varphi }$、${\bar {\sigma }}_{22}^{\varphi }$为发生纤维扭结/劈裂时纤维扭结方向的应力分量。当材料受到的纤维方向压缩荷载较小时$\left(\left|{\bar {\sigma }}_{11}\right| < \dfrac{{X}_{\rm{C}}}{2} \right)$，纤维的失效由剪切应力主导，此时纤维发生劈裂损伤；当纤维方向压缩荷载足够大时$\left( \left|{\bar {\sigma }}_{11}\right|\geqslant \dfrac{{X}_{\rm{C}}}{2}\right)$，纤维发生扭结损伤。

为求得纤维扭结/劈裂方向的应力分量${\bar {\tau }}_{23}^{\varphi }$、${\bar {\tau }}_{12}^{\varphi }$、${\bar {\sigma }}_{22}^{\varphi }$，需要对材料坐标系O123进行两次坐标变换，首先绕材料坐标系1轴逆时针旋转$\varPsi$角得到纤维扭结/劈裂所在平面的坐标系O${1}^{\varPsi }{2}^{\varPsi }{3}^{\varPsi }$，再绕该坐标系的${3}^{\varPsi }$逆时针旋转φ角得到纤维扭结/劈裂方向，如图2所示。因此，应力分量需要从材料坐标系上的${\bar {\sigma }}_{ij}$转换为纤维扭结/劈裂所在平面的坐标系O${1}^{\varPsi }{2}^{\varPsi }{3}^{\varPsi }$的应力分量${\bar {\sigma }}_{ij}^{\varPsi }$，再转换至纤维扭结/劈裂方向坐标系上的应力分量${\bar{\sigma }}_{ij}^{\varphi }$。应力分量${\bar {\tau }}_{23}^{\varPsi }$、${\bar {\tau }}_{12}^{\varPsi }$、${\bar {\sigma }}_{22}^{\varPsi }$通过如下公式转换：

图  2  纤维扭结/劈裂物理模型

Figure  2.  Physical model for fiber kinking/splitting formation

下载: 全尺寸图片 幻灯片

式中，$\varPsi$为发生纤维扭结/劈裂平面与材料坐标系2轴间夹角，取值区间为[0°，180°]。而纤维扭结/劈裂方向坐标系上的应力分量${\bar {\sigma }}_{ij}^{\varphi }$与应力分量${\bar {\sigma }}_{ij}^{\varPsi }$有如下关系：

式中，$\varphi$为纤维偏折角，该值不能通过试验测量得到，Pinho等^[26]认为该值可看作与材料剪切强度及纵向压缩强度有关的材料参数，计算如下：

式中：${\rm{s}}{\rm{ign}}$为符号函数；${\varphi }_{\rm{C}}$为纯压状态下的纤维偏折角。

1.3.2   损伤演化

当损伤阈值${r}_{I}>1$时，材料损伤发展。随着有效应力的增加，材料损伤将进一步演化，宏观表现为材料的力学属性退化及承载力下降。采用线性渐进折减方案描述材料的损伤演化：

式中：${d}_{I}$为损伤状态变量，${d}_{{\rm{f}}{\rm{t}}}$、${d}_{{\rm{f}}{\rm{c}}}$和${d}_{{\rm{m}}}$分别为纤维拉伸、纤维压缩及基体开裂损伤变量，当材料未发生损伤时，${d}_{I}=0$；当材料完全失效时，${d}_{I}=1$，为防止有限元计算过程中由于刚度折减至0导致的材料刚度矩阵产生奇异，设定${d}^{*}=$0.99；${\varepsilon }_{{\rm{e}}{\rm{q}}}^{I}$为等效应变，${\varepsilon }_{{\rm{e}}{\rm{q}},0}^{I}$为发生初始损伤时的等效应变，${\varepsilon }_{{\rm{e}}{\rm{q}},{\rm{f}}}^{I}$为完全失效时的等效应变。

本模型中材料在纤维方向上的力学行为是线弹性的，在横向及剪切方向上用塑性模型描述材料的非线性力学行为，因此在发生损伤前，材料的力学行为分为线性和非线性响应；在损伤演化阶段均采用线性渐进折减方案。在线性本构关系中，应变只包含弹性应变；非线性本构关系，应变分为弹性应变和塑性应变两部分，假定材料卸载和再加载以卸载时的弹性模量按线弹性关系进行，如图3所示。图中，${\sigma }_{{\rm{e}}{\rm{q}},0}^{I}(I={\rm{f}}{\rm{t}},{\rm{f}}{\rm{c}},{\rm{m}})$为各失效模式发生初始损伤时的等效应力。

图  3  横向和剪切方向材料损伤本构关系

Figure  3.  Material damage constitutive relationships in the transverse and shear directions

d_I (I=ft,fc,m)—Damage variables corresponding to tensile and compressive failure in the fiber direction and matrix cracking, respectively; ${\sigma }^I_{{\rm{e}}{\rm{q}},0}$ and ${\varepsilon }^I_{{\rm{e}}{\rm{q}},0}$(I=ft,fc,m)—Equivalent stresses and equivalent strains when damage initiation criteria corresponding to tensile and compressive failure in the fiber direction and matrix cracking (see Eq. (17)) are met; ${\varepsilon }^I_{{\rm{e}}{\rm{q}},{\rm{f}}}$(I=ft,fc,m)—Equivalent strains at complete tensile and compressive failure in the fiber direction and complete matrix failure

下载: 全尺寸图片 幻灯片

无论是线性还是非线性本构关系，损伤演化过程皆由材料应变引起，当材料的应变能释放率等于其断裂韧性时，意味着该材料发生完全失效。不同尺寸的有限单元采用相同的材料损伤本构关系，将导致有限元计算结果具有显著的网格敏感性。本文采用Bazant等^[35]提出的“断裂带模型”，通过引入与单元尺寸相关的单元特征长度降低采用含渐进损伤本构关系进行有限元计算的数值分析结果的网格敏感性。材料完全失效时的等效应变可由下式得出：

其中：$G_{I, c}$为各损伤模式下的材料断裂韧性，一般通过试验测得；${L}_{\rm{C}}$为单元特征长度，其计算公式如下：

式中，V为单元体积。

各损伤模式在任意时刻、发生初始损伤时和完全失效时的等效应变，发生初始损伤时的等效应力及对应的材料断裂韧性按下列各式计算：

(1)纤维拉伸破坏：

式中，${G}_{{\rm{f}}{\rm{t}},{\rm{c}}}$为纤维拉伸断裂韧性。

(2)纤维压缩破坏：

式中，${G}_{{\rm{f}}{\rm{c}},{\rm{c}}}$为纤维压缩断裂韧性。

(3)基体破坏模式：

式中：${G}_{{\rm{I}}{\rm{C}}}$、${G}_{{\rm{I}}{\rm{I}}{\rm{C}}}$、${G}_{{\rm{I}}{\rm{I}}{\rm{I}}{\rm{C}}}$分别为基体损伤张开型裂纹(I型)、滑移型裂纹(II型)和撕开型裂纹(III型)对应的断裂韧性；${G}_{{\rm{m}}}$为混合模式下的等效基体断裂韧性。

此时，式(4)中损伤变量由下式计算：

式中，${\theta }_{{\rm{f}}{\rm{p}}}$为基体损伤初始后基体断裂面角度，可由1.3.3节求解方法计算得到。

1.3.3   断裂角求解方法

基体损伤准则及纤维压缩损伤准则的加载函数${F}_{{\rm{m}}}$和${F}_{{\rm{f}}{\rm{c}}}$是潜在断裂面角度$\theta$和$\varPsi$的一元函数。在不同的应力状态下，存在$\theta$和$\varPsi$的值使得对应的加载函数${F}_{{\rm{m}}}$和${F}_{{\rm{f}}{\rm{c}}}$达到最大值，因此需要求解令加载函数${F}_{{\rm{m}}}$和${F}_{{\rm{f}}{\rm{c}}}$达到最大值时的潜在断裂面角度$\theta$和$\varPsi$。

为求解Puck失效准则中基体损伤初始断裂面角度，杨凤祥等^[20]提出了分区反抛物线插值法(Selective range inverse parabolic interpolation algorithm, SRIPI)，与逐角度遍历法和分区黄金分割法相比，该方法具有更高的计算精度和效率。基体潜在断裂面角度$\theta$的求解问题与纤维扭结/劈裂潜在平面角度$\varPsi$的求解问题相同，本文采用分区反抛物线插值法对潜在断裂面角度$\theta$和$\varPsi$进行求解。现以求解基体潜在断裂面角度$\theta$为例，阐述该算法。该算法以步长为10°将基体潜在断裂面角度$\theta$取值区间[−90°，90°]等分为18个子区间，分别计算18个间隔点上函数${F}_{{\rm{m}}}$的数值，并比较每相邻的3个间隔点上函数${F}_{{\rm{m}}}$的值。如果中间间隔点上数值高于两侧间隔点上的相应值，则可以确定此区间内必定存在函数${F}_{{\rm{m}}}$的一个局部极大值。再采用反抛物线插值法求解每个存在局部最大值的区间内的相应局部最大值，最后对比所有的局部极大值可以得到全局最大值，此时对应的基体断裂面角度即为基体潜在断裂面初始断裂面角度$\theta$。

图4为分区反抛物线插值法示意图，存在局部最大值的区间内3个点的坐标分别为${{A}}\left( {\theta }_{1} ,{F}_{{\rm{m}}}\left({\theta }_{1}^{}\right) \right)$、${{B}}\left({\theta }_{2} ,{F}_{{\rm{m}}}\left({\theta }_{2}^{}\right)\right)$和${{C}}\left({\theta }_{3}, {F}_{{\rm{m}}}\left({\theta }_{3}^{}\right) \right)$，反抛物线插值法构造函数为

图  4  分区反抛物线插值法(SRIPI)示意图

Figure  4.  Illustration of selective range inverse parabolic interpolation algorithm (SRIPI)

下载: 全尺寸图片 幻灯片

对应的基体潜在断裂面角度$\theta$由下式得出：

对于加载过程中的每一个应力状态都存在一个基体潜在断裂面角度$\theta$使得${F}_{{\rm{m}}}$($\theta$)数值达到最大值，当${F}_{{\rm{m}}}$($\theta$)<1时，基体损伤尚未发生；当${F}_{{\rm{m}}}$($\theta$)=1时，基体损伤初始发生，对应的基体潜在断裂面角度$\theta$即被确定为基体损伤初始断裂面角度${\theta }_{{\rm{f}}{\rm{p}}}$。

纤维扭结/劈裂平面角度$\varPsi$按相同方法求解，将式(33)和式(34)中的${FI}_{{\rm{m}}}$和$\theta$分别替换为${FI}_{{\rm{f}}{\rm{c}}}$和$\varPsi$。

2.   数值实现方法

已开发的模型及其算法通过编写用户自定义材料子程序VUMAT嵌入到有限元显式计算程序ABAQUS/Explicit中，适用于对复合材料层合板进行渐进失效分析。

2.1   数值计算程序

在应变增量${{\Delta }\boldsymbol{\varepsilon} }_{n+1}$驱动下，将${t}_{{\rm{l}}{\rm{o}}{\rm{a}}{\rm{d}}}$时长内的加载过程离散为一系列的微小增量步，荷载以增量形式施加，将第$n$增量步的计算结果$\left\{{\boldsymbol{\varepsilon}}_{\left(n\right)},{\boldsymbol{\varepsilon}}_{\left(n\right)}^{\rm{e}},\right.$$\left. {\boldsymbol{\varepsilon} }_{\left(n\right)}^{\rm{p}},{\tilde {\varepsilon }}_{\left(n\right)}^{\rm{p}},{\boldsymbol{\sigma} }_{\left(n\right)},{\bar{\boldsymbol{\sigma}}}_{\left(n\right)},{\tilde {\sigma }}_{\left(n\right)},{d}_{{{I}},\left(n\right)},{\theta }_{\left(n\right)} ,{\varPsi }_{\left(n\right)}\right\}$作为第$n$+1增量步的初始条件，求解出第$n$+1增量步更新后变量值$\left\{{\boldsymbol{\varepsilon} }_{\left(n+1\right)},{\boldsymbol{\varepsilon} }_{\left(n+1\right)}^{\rm{e}},{\boldsymbol{\varepsilon} }_{\left(n+1\right)}^{\rm{p}},{\tilde {\varepsilon }}_{\left(n+1\right)}^{\rm{p}},{\boldsymbol{\sigma} }_{\left(n+1\right)}, {{\bar {\boldsymbol{\sigma }}}_{\left(n+1\right)},\tilde {\sigma }}_{\left(n+1\right)},{{d}_{{{I}},\left(n+1\right)}},\right.$$\theta _{(n+1)}^{},\left.{\varPsi }_{(n+1)}\right\}$，更新后的应力以及相关的状态变量通过用户子程序VUMAT传递至下一增量步作为初始条件，具体步骤如下：

(1)初始条件：

(2)计算应变率，更新材料参数：

式中：${u}_{1}$为沿试件长度方向施加的位移；${l}_{0}$为试件在加载方向上的原始长度；假设整个加载过程中应变率保持不变。

(3)弹性预测步：

其中，${\boldsymbol{\varepsilon} }_{(n+1)}^{{\rm{p}},{\rm{t}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{l}}}$、${\boldsymbol{\varepsilon} }_{(n+1)}^{{\rm{e}},{\rm{t}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{l}}}$、${\tilde {\varepsilon }}_{(n+1)}^{{\rm{p}},{\rm{t}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{l}}}$、${\widetilde {\sigma }}_{(n+1)}^{{\rm{t}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{l}}}$、${\bar {\boldsymbol{\sigma} }}_{(n+1)}^{{\rm{t}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{l}}}$和${F}_{(n+1)}^{{\rm{t}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{l}}}$分别为塑性应变、弹性应变、等效塑性应变、等效应力、有效应力和塑性屈服准则的试算值。

(4)塑性修正步：

判断${F}_{(n+1)}^{{\rm{t}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{l}}}$：若${F}_{(n+1)}^{{\rm{t}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{l}}}<0$，表明材料未发生塑性变形，令：

若${F}_{(n+1)}^{{\rm{t}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{l}}}\geqslant 0$，则材料发生塑性变形，进入塑性修正步：利用向后欧拉法更新应力使当前增量步下的应力状态重新满足屈服条件${F}_{(n+1)}^{{\rm{t}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{l}}}\left({\bar{\boldsymbol{\sigma} }}_{(n+1)}^{},{\tilde {\varepsilon }}_{(n+1)}^{\rm{p}}\right)\leqslant$TOL，TOL为容许误差，本文取$1\times {10}^{-6}$。

(a)初始化向后欧拉法所需要的参数：

(b)计算塑性屈服准则：

(c)计算塑性流动因子增量$\delta \left({\Delta }\lambda \right)$，更新塑性流动因子：

(d)判断${F}_{(n+1)}^{(k+1)}\leqslant$TOL，若成立，进入步骤(5)，否则设$k=k+1$，返回到步骤(4b)。

(5)损伤判断步：

(a)采用分区反抛物线插值法计算基体潜在断裂面角度${\theta }_{(n+1)}$及纤维扭结/劈裂潜在平面角度${\varPsi }_{(n+1)}$；

(b)计算基体断裂面上的有效应力${\bar {\sigma }}_{\rm{n}}$、${\bar {\tau }}_{{\rm{t}}}$、${\bar {\tau }}_{\rm{l}}$及纤维扭结/劈裂平面上的有效应力${\bar {\sigma }}_{22}^{\varphi }$、${\bar {\tau }}_{12}^{\varphi }$、${\bar {\tau }}_{23}^{\varphi }$；

(c)检查判断准则${F}_{I,(n+1)}$($I={{\rm{f}}{\rm{t}},{\rm{f}}{\rm{c},\rm{m}}}$)，若${F}_{I,(n+1)} \geqslant 1$，跳至步骤(6)，且当${F}_{{\rm{m}},(n+1)}=1$时，该时刻对应的潜在基体断裂面角度${\theta }_{(n+1)}$将被保存作为随后损伤演化阶段的固定值${\theta }_{{\rm{f}}{\rm{p}}}$，当${F}_{{\rm{f}}{\rm{c}},(n+1)}^{}=1$时，该时刻对应的潜在断裂面角度${\varPsi }_{(n+1)}$将被保存作为随后损伤演化阶段的固定值${\varPsi }_{{\rm{f}}{\rm{p}}}$；否则跳至步骤(7)。

(6)损伤修正步：

(7)更新名义应力：

已编写的用户自定义子程序VUMAT流程图如图5所示。

图  5  用户自定义子程序VUMAT流程图

Figure  5.  Flow chart of the user subroutine VUMAT

下载: 全尺寸图片 幻灯片

2.2   数值模拟

为验证已开发模型的有效性，对已报道的IM7/8552碳纤维/环氧树脂单向复合材料层合板在准静态(10⁻⁴s⁻¹)、中应变率(1 s⁻¹)和高应变率(300~1000 s⁻¹)加载下的偏轴压缩试验^[38]进行了数值模拟，并预测其极限荷载。试件铺层方式分别为[15°]₅₄、[30°]₅₄、[45°]₅₄、[60°]₅₄、[75°]₅₄和[90°]₅₄，试件长25.0 mm，宽7.6 mm，厚7.6 mm，试件几何形状及尺寸如图6所示。

图  6  IM7/8552碳纤维/环氧树脂偏轴压缩试件的几何形状及尺寸

Figure  6.  Geometry and dimensions of IM7/8552 carbon fiber/epoxy resin off-axial compressive specimens

下载: 全尺寸图片 幻灯片

有限元模型使用C3D8R实体单元进行模拟，由于试件为单向板，假定各层界面为完全粘结，且为了防止单元数量过多导致模型计算量太大和单元尺寸长宽高比太大导致计算结果出现大的误差，设置网格尺寸为0.95 mm×0.95 mm×0.95 mm，如图7所示。有限元模型左端设置约束限制左端平面在加载方向上的位移，并允许其在平面内发生变形；右端加载由位移控制，施加在其面外的一个参考点上，采用Create constraint方法对参考点与右端加载平面建立Equation线性约束，保证加载在参考点上的位移等于施加在加载平面上的位移。数值分析采用的IM7/8552碳纤维/环氧树脂复合材料的弹性模量、破坏强度和泊松比等均取自本算例所在文献[38]，材料断裂韧性取自同种材料的相应数值^[39]，塑性模型参数β、${n}_{\rm{p}}$和${a}_{66}$的数值由本文采用Sun和Chen^[31]的方法根据文献[38]中报道的偏轴试验应力-应变曲线确定。材料属性及模型参数值详见表1。

图  7  IM7/8552碳纤维/环氧树脂复合材料单向层合板试件有限元模型

Figure  7.  Finite element model of IM7/8552 carbon fiber/epoxy unidirectional laminated specimen

下载: 全尺寸图片 幻灯片

表  1  IM7/8552碳纤维/环氧树脂复合材料单向板属性

Table  1.  Material properties of IM7/8552 carbon fiber/epoxy unidirectional composite laminate

${E}_{11} ^0$/GPa	${E}_{22}^0={E}_{33} ^0$/GPa	${X}_{\rm{T}}$/MPa	${X}_{\rm{C}}$/MPa	${Y}_{\rm{T}}$/MPa	${Y}_{\rm{C}}$/MPa	${S}_{12}$/MPa
154	9	1725	2650	76.4	288	89
${\nu }_{12}={\nu }_{13}$	${\nu }_{23}$	${G}_{12}^0={G}_{13} ^0$/GPa	${G}_{23} ^0$/GPa	${G}_{{\rm{I}}{\rm{C}}}$/(N·mm−1)	${G}_{{\rm{I}}{\rm{I}}{\rm{C}}}={G}_{{\rm{I}}{\rm{I}}{\rm{I}}{\rm{C}}}$/(N·mm−1)
0.33	0.43	5.6	3.98	0.26	1.002
${G}_{{\rm{f}}{\rm{t}},{\rm{c}}}$/(N·mm−1)	${G}_{{\rm{f}}{\rm{c}},{\rm{c}}}$/(N·mm−1)	${m}_{\rm{e}}$	${m}_{\rm{s}}$	${a}_{66}$	β	${n}_{\rm{p}}$
80	80	0.035	0.055	2	794.233	0.1973
Notes: ${E}_{11} ^0$, ${E}_{22} ^0$, ${E}_{33} ^0$—Quasi static elastic moduli in the fiber, transverse and through-thickness directions; ${X}_{\rm{T}}$,${X}_{\rm{C}}$—Tensile and compressive strengths in the fiber direction; ${Y}_{\rm{T}}$, ${Y}_{\rm{C}}$—Tensile and compressive strengths in the transverse direction; ${S}_{12}$—Shear strength in 1-2 plane; ${\nu }_{ij}$,${G}_{ij} ^0$ (ij=12,13,23)—Poisson’s ratios and shear moduli in the 1-2,1-3, 2-3 planes; ${G}_{{\rm{I}}{\rm{C}}}$, ${G}_{{\rm{I}}{\rm{I}}{\rm{C}}}$, ${G}_{{\rm{I}}{\rm{I}}{\rm{I}}{\rm{C}}}$—Fracture toughnesses for Mode I, II and III failure in the through-thickness direction, respectively; ${G}_{{\rm{f}}{\rm{t}},{\rm{c}}}$, ${G}_{{\rm{f}}{\rm{c}},{\rm{c}}}$—Tensile and compressive fracture toughnesses in the fiber direction; ${m}_{\rm{e}}$, ${m}_{\rm{s}}$—Strain-rate parameters related to elastic moduli and failure strengths; ${a}_{66}$—Material parameter of the plastic model; β, ${n}_{\rm{p}}$—Coefficients that fit the experimental hardening curve.

下载: 导出CSV 

| 显示表格

采用本文提出的三维弹塑性损伤本构模型和强化函数考虑或未考虑应变率效应的模拟方法对偏轴角分别为15°、30°、45°、60°、75°和90°的IM7/8552碳纤维/环氧树脂复合材料单向板试件在准静态、中、高应变速率下的偏轴压缩试验进行模拟，得到的预测应力-应变曲线如图8所示，其中Present model为强化函数考虑应变率效应的模拟方法，而Reference model为强化函数未考虑应变率效应的模拟方法。

图  8  不同偏轴角度IM7/8552碳纤维/环氧树脂复合材料单向层合板在应变率分别为10⁻⁴ s⁻¹(准静态)、1 s⁻¹(中应变率)和高应变率压缩荷载下的试验与预测应力-应变曲线

Figure  8.  Experimental and predicted stress-strain curves of off-axis IM7/8552 carbon fiber/epoxy unidirectional laminated specimens subjected to compressive loadings at strain rates of 10⁻⁴ s⁻¹ (quasi-static), 1 s⁻¹ (intermediate) and high strain rates

下载: 全尺寸图片 幻灯片

由计算结果可知，与强化函数未考虑应变率效应的三维弹塑性损伤本构模型相比，强化函数考虑应变率效应的三维弹塑性损伤本构模型对IM7/8552碳纤维/环氧树脂复合材料单向层合板在中、高应变率加载条件下偏轴压缩试验的数值模拟精度有了明显的提升，各应变率加载条件下在15°、30°、45°、60°、75°和90°偏轴方向上的预测结果与试验结果吻合良好。强化函数未考虑应变率效应的三维弹塑性损伤本构模型，其材料刚度及材料强度考虑了应变率效应，因此在动态荷载下的预测应力峰值较准静态荷载下的预测应力峰值有所提升，但与试验结果仍存在一定的误差，而强化函数考虑应变率效应的模拟方法在动态荷载下的预测应力峰值与试验结果更为接近。数值模拟结果表明，类似于在动态荷载作用下，材料刚度及材料强度有显著提升，纤维增强树脂基复合材料在动态荷载作用下的塑性力学行为与准静态荷载作用下塑性力学行为也有着显著的差异，需要对塑性模型中的强化函数进行率相关修正才能得到合理的数值模拟结果。综上，本文开发的率相关三维弹塑性损伤本构模型能较准确地反映复合材料在准静态，中应变率和高应变率下的力学响应，为预测具有显著应变率效应和塑性效应的复合材料极限荷载提供了一种有效的计算方法。

3.   结 论

(1)基于弹塑性力学和连续损伤力学，建立了同时考虑复合材料非线性力学响应、应变率效应和损伤累积导致材料属性退化的三维弹塑性损伤本构模型。推导了本构模型的应变驱动显式积分算法以更新应力和与解答相关的状态变量，开发了包含数值积分算法的用户自定义材料子程序VUMAT，并嵌于有限元程序ABAQUS v6.14中。模型被用于横向和剪切方向力学行为表现出显著塑性效应和应变率效应的IM7/8552碳纤维/环氧树脂复合材料层合板偏轴压缩试验的数值模拟，预测结果与试验结果吻合良好。开发的本构模型能够较准确地反映此类复合材料层合板的力学行为，为复合材料构件和结构设计提供一种有效的分析方法。

(2)在动态荷载作用下，纤维增强树脂基复合材料的塑性力学行为与准静态荷载作用下塑性力学行为有着显著的差异，需要对材料弹性模量、破坏强度和塑性模型中的强化函数进行考虑应变率效应的修正才能得到合理的数值模拟结果。与塑性强化函数未考虑应变率效应的本构模型相比，本文提出的含改进塑性强化函数的本构模型对预测纤维增强树脂基复合材料在中应变率和高应变率下塑性力学响应与承载力峰值具有更高的准确性。适用于对横向和剪切方向上表现出显著塑性效应和应变率效应的纤维增强树脂基复合材料的数值模拟。