复合材料由于高比强度高、高比刚度和优秀的抗疲劳性等性能^[1-3]，被广泛应用于航空航天、船舶、建筑和体育等军民两用领域^[4-5]。对于工程常用的复合材料层合板结构，铺层纤维的存在使得其纵向和横向力学性能得到显著增强，但层间性能较弱。这使得结构在受到一定冲击或结构的连接边界受到挤压作用时，极易产生层间分层。分层损伤的发生虽然不会导致结构突然的破坏，但结构刚度和强度的大幅度下降是致命的^[6-8]。因此，分层问题一直以来受到复合材料领域研究人员的广泛关注。在对复合材料结构进行设计和强度的分析过程中，需要对复合材料结构的分层损伤起始与扩展过程开展深入的研究^[9-11]。

数值模拟技术具有简便、快捷和成本低的特点，是复合材料结构设计领域的重要方法。现阶段针对单向层合板复合材料的分层扩展模拟已经得到广泛的研究。但是，对于工程中应用更为广泛，分层扩展机理更为复杂的多向层合板的分层扩展模拟研究还远远不足^[12]。预测复合材料单向层合板的分层失效时，现有的分层扩展准则通常选用恒定的层间断裂韧性作为基本性能参数^[13]，通过对比裂尖应变能释放率与层间断裂韧性的大小来判断分层是否扩展。而多向层合板失效机理复杂^{[14, 15]}，分层界面张开后裂尖后面出现大面积的纤维桥接。随着分层的进一步扩展，断裂韧性逐渐增加并最终达到一个稳定值，表现出显著的R曲线特征^[16-18]。现有分层扩展准则没有考虑这一特征，不适用于含纤维桥接多向层合板分层扩展行为的模拟。龚愉等^{[13, 16]}对现有分层扩展准则进行改进，提出将R曲线引入到准则中，并据此建立了考虑纤维桥接影响的分层扩展准则。但是改进的准则由于位置依赖性，可以用于一维分层的模拟而不适用于二维或三维分层扩展情况。为解决该问题，作者提出了考虑纤维桥接的新型三线性内聚力本构。该新型本构在运用过程中涉及内聚力模型关键参数的确定。

笔者旨在所提新型三线性内聚力本构的基础上，采用二维有限元模型对三线性本构模型中的关键参数开展数值研究，分析关键参数对数值模拟结果的影响，确定一套适用于多向层合板分层扩展数值模拟的新型本构内聚力模型界面参数取值方案，最后采用该取值方案确定的内聚力模型界面参数建立有限元模型，对复合材料多向层合板的分层扩展行为进行预测。

1 新型三线性本构内聚力模型简介

纤维桥接是导致复合材料层合板分层扩展阻力增加，产生R曲线效应的关键因素。准确表征纤维桥接区的纤维在整个分层失效过程中的应力应变关系对分层扩展行为的准确模拟非常重要。作者所提出的考虑纤维桥接影响的新型三线性本构内聚力模型如图 1所示。该本构关系从微观失效机理角度出发建立，充分体现了分层扩展的复杂失效过程。此外，具有参数简单、物理意义明确且不依赖于位置信息等优点。该新型内聚力模型满足以下两点假设：

1) 由2个分别代表基体和纤维失效的双线性本构关系组合而成，如图 1所示，其中紫色线段代表基体的分层失效，紫色线段所组成的三角形面积等于分层扩展的起始断裂韧性G_ini；橙色线段代表桥接纤维失效，同时橙色段所组成的三角形面积等于由纤维桥接引起的断裂韧性G_bri；两断裂韧性值之和等于纤维桥接完全发展后分层扩展的稳态断裂韧性G_prop；

2) 在基体发生失效形成新裂纹的时刻，桥接区的纤维裸露并承担载荷，桥接应力达到其最大值，即δ_b既是基体失效双线性本构中的临界失效位移，同时也是纤维桥接双线性本构中的损伤初始点处对应的位移。

新型三线性本构关系如式(1)~(3)所示。

$ \sigma = (1 - d)D_{ij}^0, $

(1)

$ \delta = \left\{ {\begin{array}{*{20}{c}} {0, }&{\delta \le {\delta _0}, }\\ {\left( {1 - \frac{{{K_{AB}}}}{{{K_0}}}} \right)\left( {1 - \frac{{{\delta _0}}}{\delta }} \right), }&{{\delta _0} < \delta \le {\delta _{\rm{b}}}, }\\ {\left( {1 - \frac{{{K_{BC}}}}{{{K_0}}}} \right)\left( {1 - \frac{{{\delta _{\rm{f}}}}}{\delta }} \right), }&{{\delta _{\rm{b}}} < \delta \le {\delta _{\rm{f}}}, }\\ {1, }&{\delta > {\delta _{\rm{f}}}, } \end{array}} \right. $

(2)

式中：d为全局界面状态损伤变量，D_ij⁰为起始刚度标量并可定义为

$ D_{i j}^{0}=\bar{\delta}_{i j} K_{0}。$

(3)

式中：参数K₀为罚刚度，等于表征基体分层失效的内聚力单元界面初始刚度K₁与表征纤维桥接的内聚力单元界面初始刚度K₂之和。K_AB和K_BC分别为图 1中所示直线AB和直线BC的斜率，其值可通过式(4)~(5)计算得到，该式中σ₀和σ_b分别为表征基体分层失效的双线性本构的界面强度和表征桥接纤维的双线性本构的桥接强度。δ₀和δ_f分别为基体开始产生损伤的临界损伤位移和纤维桥接充分发展后桥接纤维失效的临界失效位移。

$ {{K_{AB}} = \frac{{{\sigma _{\rm{b}}} - \left( {{\sigma _0} + {K_2}{\delta _2}} \right)}}{{{\delta _{\rm{b}}} - {\delta _0}}}, } $

(4)

$ {{K_{BC}} = \frac{{{\sigma _{\rm{b}}}}}{{{\delta _{\rm{b}}} - {\delta _{\rm{f}}}}}}。$

(5)

2 新型内聚力模型关键参数介绍 2.1 内聚力区的单元个数

裂纹尖端存在一个微小的内聚力区，如图 2所示。内聚力模型是用于表征内聚力区力学行为的界面张开位移关系曲线。完整的内聚力模型包含损伤起始、损伤演化、完全失效3个阶段。为了准确地描述分层扩展损伤失效的整个过程，内聚力区必须包含足够的内聚力单元个数。而内聚力单元本质上是相互分离的2个界面，如图 3所示。2个界面间的相对位移δ决定了内聚力单元在损伤失效过程中的不同阶段或状态。在张开型分层模式下，界面间相对位移在裂纹扩展方向是不同的，但在宽度方向分布的内聚力单元的界面间相对位移是相同的。因此，内聚力区单元个数上是指二维内聚力区所划分的单元数目，如图 2所示。

当内聚力区的网格划分尺寸过大时，内聚力区单元数目过少，将无法准确体现内聚力模型的3个阶段特征，从而造成分层扩展损伤失效过程预测的失真。当内聚力单元尺寸过小时，虽然能准确描述分层失效过程，但计算成本大大增加。为了划分适当的内聚力区网格尺寸，需要预估内聚力区的长度。在早期研究中，用于估算各向同性材料内聚力区长度的方程为

$ l_{\mathrm{cz}}=E \frac{G_{\mathrm{C}}}{\left(\sigma^{0}\right)^{2}}, $

(6)

式中：E是材料的杨氏模量，G_C是临界应变能释放率，σ⁰是界面强度。对于正交各向异性材料，研究普遍采用将E取为其横向模量。在预估了内聚力区长度之后，即可根据内聚力区需要划分的单元个数确定合适的单元尺寸。

2.2 界面初始刚度

内聚力单元本构关系中，界面初始刚度是用来确保内聚力单元在未损伤之前上下界面间保持刚性连接的非物理量。在理想情况下，界面初始刚度应该是无限大的，但在有限元分析中，有限大小的界面初始刚度值才能保证迭代收敛。关于合理的界面初始刚度的选取，目前可从两种角度考虑，一种是从材料的韧性和脆性角度，一种是考虑内聚力界面变形对整体结构刚度的影响^[19]。

从材料的韧性和脆性角度出发。在确保界面初始刚度的选取满足分层萌生时的界面相对位移小于等于分层扩展时的界面相对位移时。界面初始刚度越大，材料的塑性越显著，反之，材料的脆性性能越显著。其次，考虑到内聚力界面的变形也会贡献于整体结构的变形，内聚力界面在分层时的变形要足够小，以确保内聚力界面不会对整体结构的刚度损伤产生影响。目前，界面初始刚度的准确选取尚没有通用的方法，普遍的选取范围是在10¹² N/m³到10¹⁵ N/m³之间。

2.3 内聚力单元黏性系数

内聚力模型中采用黏结单元时，由于加载过程中存在材料的软化和刚度退化，计算过程中会出现严重的收敛性问题。研究者者通常采用本构方程的黏性规则化，该方法通过引入黏性，在足够小的时间增量下可使软化材料的雅可比矩阵始终保持为正，以改善收敛性，如方程(7)所示。

$ {{\dot D}_{\rm{v}}} = \frac{1}{\mu }\left( {D - {D_{\rm{v}}}} \right), $

(7)

式中：D_v为对应于黏性的界面状态损伤变量，D为不考虑黏性时的界面状态损伤变量，μ是黏性系数。ABAQUS^®内聚力单元模块中已经集成该方法，内聚力单元黏性系数可在单元属性编辑过程中直接输入。研究者普遍选取黏性系数在10^-6~10^-3之间，新型三线性本构关系黏性系数选取也以此为参考。

3 新型内聚力模型关键参数研究

从降低数值计算时间的角度出发，以二维有限元模型为例，对新型三线性本构内聚力模型中的3个关键界面参数，包括内聚力区的单元个数、表征基体失效的内聚力单元的界面初始刚度和内聚力单元黏性系数，进行有限元数值分析并讨论参数的影响规律。

采用双悬臂梁试验装置开展拉伸载荷下的Ⅰ型分层测试，提取分层扩展过程中的载荷位移曲线。其中，载荷为双悬臂梁加载点的支反力，位移为上下臂加载点的相对位移。试样由4层单向铺层层合板0₂//0₂构成，每层铺层厚度为1.5 mm，符号“//”表示预制分层位置。试样左侧的中部为预制分层位置，在制造过程中通过铺设Teflon薄膜来得到预制分层。如图 4所示，试样几何尺寸如下：长度100 mm，单臂厚度3 mm，预置裂纹长度30 mm。采用有限元软件ABAQUS^®进行有限元模拟，层板上下臂采用二维平面应变单元CPE4I建立，单臂长度方向上设置200个单元，厚度方向上设置2个单元。分层采用放置于层合板中间位置的一层薄的内聚力单元进行模拟，内聚力单元采用ABAQUS^©内聚力模块并结合UMAT子程序进行二次开发实现，单元类型为二维四节点内聚力单元COH2D4，长度方向根据所考察的内聚力区单元个数决定单元划分情况(在长度方向设置280个单元时，内聚力区包含3个单元)，厚度方向设置的单元数为1。为了更好的网格适用性，内聚力单元的上下面分别与上臂的下表面和下臂的上表面之间建立绑定约束。边界条件和载荷情况如图 4所示，表 1同时列出了所研究层合板的基本力学参数^[20]。表征基体失效的双线性本构的界面参数列于表 2中。参照此表设计了新型三线性本构内聚力模型界面参数的数值研究方案。

三线性本构内聚力模型中的断裂韧性起始值G_ini取值100 N/m，扩展断裂韧性的稳定值G_prop取为460 N/m，使得两者的平均值等于表 2中双线性内聚力模型中的断裂韧性。根据参考文献[21]，内聚力单元界面强度σ₀和桥接强度σ_b分别取为57 MPa和2 MPa。因此可以计算得到表征桥接纤维失效的双线性本构的界面初始刚度K₂，其值为5.7×10¹¹ N/m²。具体的参数研究方案如表 3所示。对每个界面参数选取3个典型值进行考察。

对层合板分层扩展过程中cohesive铺层中内聚力区单元个数N_e进行研究，所得结果如图 5所示。可以看出内聚力区单元个数对数值结果有着显著的影响。当单元个数较少，如N_e=3时，数值结果明显偏大。进一步增加单元个数，内聚力区单元个数从6增加到9时，数值预测的极限载荷逐渐降低，继续增加单元个数到12时，预测的载荷位移曲线没有变化。极限载荷及其对应的张开位移变化都在5%以内，趋于稳定值。考虑到较多的单元个数带来计算时间的增加，认为内聚力区单元个数取值为6是合适的。图 6给出了不同内聚力单元界面初始刚度K₁下的数值结果，可以看出界面初始刚度对数值结果影响并不显著。当K₁=10¹⁴ N/m³时，载荷位移曲线在载荷的下降阶段出现了一定程度的震荡，这主要是由单元内张力的寄生震荡造成的。考虑到较低的界面初始刚度并不能很好地表征未损伤内聚力单元与相邻层之间的刚性连接，认为K₁取值为10¹⁵N/m³是合理的，这也与一般文献中的做法相一致^[12]。

不同内聚力单元黏性系数下的数值模拟结果如图 7所示。可以看到，较大的黏性系数将导致极限载荷点附近位移载荷曲线的软化和明显偏高的峰值载荷。当μ≤10^-5时，黏性系数对模拟载荷位移曲线的影响可以忽略。考虑到过小的黏性系数将显著增加计算成本，μ取值为10^-5。通过以上数值研究，可以确定新型三线性本构关键界面参数的取值方案为：内聚力区单元个数N_e为6，表征基体失效的内聚力单元界面初始刚度K₁为10¹⁵ N/m³，内聚力单元黏性系数μ为10^-5。

4 参数方案的验证

采用第二节建立的界面参数取值方案建立有限元分析模型，对T700/QY9511多向层合板的分层扩展行为进行数值模拟，提取分层扩展过程中的载荷位移曲线。将模拟结果同试验结果进行对比，验证参数取值方案的有效性。层合板的铺层方案为：0₁₆//(+5/-5/0₆)_S，几何尺寸为：长度160 mm，宽度25 mm，单臂厚度2.08 mm，预置分层长度35 mm。T700/QY9511的材料属性为：E₁₁=130 GPa，E₂₂=E₃₃=10.4 GPa，G₁₂=G₁₃= 6.36 GPa, G₂₃=4 GPa，υ₁₂=υ₁₃=0.289，υ₂₃=0.3。此外基体QY9511的拉伸强度为100 MPa。该层合板的起始断裂韧性为100 J/m²，扩展断裂韧性的稳定值为890 J/m^{2 [12]}。

图 8所示为T700/QY9511多向层合板的有限元模型。按照参考文献[21]中的结论，表征基体分层失效双线性本构的界面强度σ₀取值为55 MPa。表征桥接纤维失效的双线性本构的桥接强度σ_b取值为0.4 MPa。有限元模型中，试样宽度和厚度方向分别设置4个和16个单元，长度方向上靠近预制裂尖处每个单元尺寸为0.5 mm，其余部分每个单元尺寸为2 mm。其中厚度方向上每个铺层均设置一个单元来更好的捕获层间应力大小。层合板采用Composite Layup模块设置每个单层的铺层角度和材料等属性。层合板采用C3D8三维实体单元进行建模，并以自上而下扫掠的方式生成模型网格。界面层采用单元类型为COH3D8，内聚力区的单元个数为6，每个单元尺寸为0.125 mm，厚度方向不再划分单元，内聚力单元黏性系数设为10^-5。同样采用自上而下扫掠的方式生成界面层的网格。边界条件和载荷情况为：右端固定，左端的上下悬臂分别施加位移载荷。为改善数值计算的收敛性，在Step模块中分两步施加位移载荷，首先施加一个较小的位移使得层合板的上下臂之间稍微张开，然后施加实际试验过程中位移载荷最大值。为方便从计算结果中提取位移和相应的支反力信息，在载荷施加位置通过建立参考点来定义集合，这样就可以在后处理过程中获得模拟的载荷位移曲线。

数值模拟方法与试验得到的载荷位移曲线均绘制于图 9中。两者对比，发现预测的载荷位移曲线在起始的线弹性阶段和分层扩展阶段均与试验曲线接近，两者吻合很好。预测的极限载荷为106.9 N，与试验结果101.8 N之间的相对误差为5%，预测结果的准确性验证了所建立新型三线性内聚力模型界面参数取值方案的有效性。

5 结论

采用数值方法研究了考虑纤维桥接影响的新型三线性本构内聚力模型的界面参数。对新模型中的关键界面参数，包括：内聚力区单元个数、表征基体失效的内聚力单元界面初始刚度和内聚力单元黏性系数的影响进行研究，获得了含纤维桥接多向层合板分层扩展行为模拟的内聚力模型界面参数取值方案。采用该取值方案对T700碳纤维增强复合材料多向层合板分层扩展行为进行模拟，模拟的载荷位移曲线与试验结果吻合较好，说明所建立的新型三线性本构内聚力模型参数选取方案是有效的。采用这种方法获得的参数取值适用于文中所研究复合材料的不同界面多向层合板I型分层扩展行为模拟。对于其他材料类型，该参数取值具有一定的参考意义。此外，本文中所建立的内聚力模型参数方案选定流程具有一定的通用性。