基于固有型内聚力模型模拟双层夹胶玻璃冲击断裂行为

双层夹胶玻璃（Laminated Glass，LG）是一种常见的复合材料，广泛应用于汽车和建筑行业。在交通事故中，LG作为汽车风挡玻璃对于行人的保护起到至关重要的作用：由于其良好的吸能特性，可以减少行人头部与汽车挡风玻璃之间的冲击作用，降低对行人的伤害；大多数玻璃碎片黏附在PVB胶层上，大大降低了因玻璃碎片飞溅造成的二次伤害。近年来，通过数值仿真方法模拟LG在不同工况下的力学性能成为国内外研究的热点，其中LG在冲击荷载作用下的裂纹扩展行为是目前研究的热点之一。Camanho等[1]为研究复合加载状态下复合材料的层间分离行为，在不同材料实体单元层间加入内聚力单元，运用BK（Benzeggagh-Kenane）损伤演化准则在位移模式下模拟复合材料的分层行为，并分别在单一和复合准静态加载模式下进行了复合材料的分层试验以验证模拟结果的合理性。Su等[2]提出在ABAQUS有限元软件中对材料模型进行网格划分，在生成的实体单元间插入被赋予软化分离准则、损伤初始准则以及损伤演化准则的内聚力单元模拟裂纹的扩展行为，并进行了劈裂试验、缺口梁的扭转试验、气缸推出试验以及缺口梁的冲击试验，验证了模拟结果的准确性。Xu等[3]从市面上的商用有限元软件对PVB夹层玻璃裂纹模拟分析的结构关系出发，将模拟结果与经典的赫兹应力模型的理论结果进行对比，验证了模拟分析结果的准确性。Peng等[4]采用LS-DYNA有限元软件，通过改变LG模型中玻璃与胶层的连接方式及网格的大小，讨论两者对冲击荷载作用下LG动态力学响应模拟的影响。Chen等[5]将非固有型内聚力单元插入到整个玻璃实体单元间模拟LG在冲击荷载作用下的裂纹扩展行为，并通过固有型内聚力单元模拟玻璃板与PVB胶层间的黏附作用；Pelfrene等[6]采用固有型内聚力本构模型研究玻璃层与PVB层间的黏结以及失效作用，随后将固有型内聚力本构模型应用于硬冲击下的单片玻璃破坏研究。Gao等[7]采用ABAQUS有限元软件，在玻璃实体单元间插入固有型内聚力单元，PVB夹层用Mooney-Rivlin（MR）超弹性本构模型，通过对比LS-DYNA与ABAQUS有限元软件中超弹性材料本构的参数差异，得出ABAQUS有限元软件中的材料参数设置，并运用ABAQUS模拟LG在冲击荷载作用下的冲击破坏行为。目前，LG冲击破坏性能研究的数值计算方法主要分为4种，分别是离散元/有限元耦合法、单元删除法、扩展有限元（XFEM）、内聚力模型法。离散元/有限元耦合法需考虑离散元区域与有限元区域交界处的几何尺寸差异问题；单元删除法难以真实反映裂纹的张合现象，仿真所得裂纹与试验所得裂纹差别较大；XFEM得到的裂纹模式与试验结果相差较大；而内聚力模型法通过优化脆性材料网格，可获得与试验结果更接近的宏观冲击破坏特性以及实现裂纹的张合现象。因此，本研究采用内聚力模型法模拟LG在冲击荷载作用下的裂纹扩展行为。由于Gao等[7]已研究了材料参数对LG抗冲击性能的影响，但在一些材料参数上样本数量不足，无法得出比较有规律性的结论。因此本研究先优化“过渡型”网格画法，然后系统地研究和分析LG在冲击荷载作用下裂纹扩展规律，并细化玻璃厚度和内聚力单元参数K的样本数量，进一步研究两者对于LG裂纹扩展的影响。

1 理论基础

1.1 内聚力破坏理论

在内聚力破坏理论中，断裂发生是裂纹尖端或内聚力区域中物质间克服内聚作用逐渐分离的结果。如图1所示，在裂纹尖端附近的微小区域内，作用力达到材料承受极限时，孔穴开始生长，裂纹随之产生，材料逐渐软化[8]。因此，将裂纹尖端的微小区域视为内聚力区。在该区域内，材料受到裂纹尖端应力场作用产生微小的分离量，同时内聚力区表面受到张力作用。通常用表面张力和分离量有关的函数表征内聚力区域的力学性能，即内聚力本构模型。在采用内聚力模型模拟断裂问题时，事先无需设置预制裂纹，通过计算可得到裂纹的发生及裂纹尖端的位置、形状等，一定程度上削弱甚至解决了裂纹尖端应力奇异性[9]的问题。

图1 内聚力破坏过程
Fig.1 Cohesive failure process

2 有限元模型

图2 LG模型
Fig.2 Laminated glass model

LG由2层相同厚度的玻璃与1层PVB胶层构成，其中单层玻璃的厚度是2 mm，PVB胶层厚度是0.76 mm。考虑到模型对称性以及计算效率问题，模型取LG实际尺寸的1/4进行建模（见图2），LG水平面内尺寸为100 mm×40 mm，玻璃板支撑面内尺寸为40 mm×10 mm，厚度为4 mm，锤头由半球和圆柱体组合而成。有限元模型通过ABAQUS前处理软件建立，单位制采用mm-t-s。如图3所示，LG冲击破坏后其裂纹扩展路径呈现“中心塌陷，环径向裂纹交错并由中心向外辐射”的特点。由于固有型内聚力模型法将设置了初始损伤和演化准则的零厚度内聚力单元插入玻璃单元间，当LG冲击破坏后内聚力单元失效，玻璃单元间发生相对位移，即裂纹产生，所以通过该方法模拟LG裂纹扩展路径将依赖于玻璃网格的划分方式。为此，本研究通过对玻璃表面分区并局部布种实现“过渡网格”的划分，使其中心网格细化，从中心向外划分足够多的环向和径向的网格，以遵循玻璃裂纹扩展规律。同时考虑到对于冲击问题的数值仿真，网格的密度和质量对模拟结果有明显影响。经过多次调试，采用图4所示的网格划分方法，其中：玻璃中玻璃实体单元（C3D8）个数为35 152，内聚力单元（COH3D8）个数为95 676；锤头中实体单元（C3D8）个数为3 753；PVB胶层中实体单元（C3D8）个数为4 397；橡胶支撑中实体单元（C3D8）个数为756。上玻璃板、PVB胶层和下玻璃板在冲击方向（z方向）的网格大小分别为0.50、0.76、0.50 mm，LG被分割为9层。玻璃与胶层之间通过Tie约束实现连接。

图3 LG冲击破坏后的裂纹路径
Fig.3 Crack path of laminated glass after impact failure

图4 网格划分
Fig.4 Mesh generation

3 材料本构模型

3.1 内聚力单元本构模型

固有型内聚力本构模型最早应用于模拟复合夹层板间的黏附作用[10]，目前广泛应用于玻璃等脆性材料在高应变率、大变形、高冲击力下脆性断裂行为的仿真模拟。图5所示为双线性内聚力本构模型曲线。由图5可知，双线性内聚力本构模型模拟裂纹生成过程可分为两个阶段[11]：在弹性阶段（OA段），内聚力单元表面张力随等效分离量 δm线性增大，在该阶段内聚力单元呈现出线弹性的特性，其中K为罚刚度，表征内聚力单元的“软硬”程度；当 δm达到初始损伤分离量时，材料进入损伤状态，故称AC段为损伤演化阶段；当 δm达到完全失效分离量时，内聚力单元完全失效，形成裂纹。该内聚力模型在裂纹开始生成前就包含裂纹生成的初始损伤判定准则，故称其为固有型。常用的固有型内聚力单元本构模型主要有双线性型、梯形型以及指数型。根据内聚力单元间的法向、切向分离量是否合成为等效分离量可将内聚力单元失效模式分为单一模式和混合模式[12]：单一模式是内聚力单元表面的等效位移量等于法向分离量或切向分离量，混合模式是内聚力单元表面的等效位移量等于法向分离量和切向分离量的合成分离量。本研究采用混合模式的双线性内聚力本构模型研究LG在冲击作用下的裂纹扩展行为。

图5 双线性内聚力本构模型
Fig.5 Bilinear constitutive model of cohesion

由图5可知，双线性内聚力本构模型涉及到4个参数：等效分离量 δm、初始损伤分离量、完全失效分离量以及不同分离量对应的应力值 τi(i=1,2,3)。可通过以下公式计算上述4个参数。

（1）等效分离量 δm的计算公式如下

式中：δ1、 δ2、 δ3分别为Mode I、II和III所对应的分离量， δshear 为切向分离量，运算符 〈〉表示如下

为了表征玻璃的脆性，采用最大名义应力准则作为材料的损伤起始准则，即

式中： T1为Mode I所对应的材料抗拉强度， T2、 T3分别为Mode II和Mode III下对应的材料抗剪强度。假定Mode I、II和III下的K值均相同，弹性阶段的拉应力和剪应力 τ1、 τ2、 τ3的计算公式如下

式中分别为Mode I、II和III下对应的初始损伤分离量。

（2）分别将（1）式、（4）式和（5）式代入（3）式，可整理出材料初始损伤分离量的表达式

（3）本研究通过BK准则[13]表征材料的损伤演化过程，在该准则下假定Mode II、III下的能量释放率相等，即由BK准则的表达式可得完全失效分离量

式中分别表示Mode I、II和III下所对应的能量释放率[14]分别表示混合度和材料参数分别为Mode II和Mode III下材料完全失效分离量。

3.2 系统本构模型

在系统本构模型中将玻璃实体和钢性材质的锤头均视为线弹性体。而PVB胶层是一种高分子聚合物[15]，具有很好的柔韧性，其应力-应变曲线与橡胶材料类似，故PVB本构模型采用与橡胶材料相似的超弹性材料。PVB胶层和橡胶支撑均采用ABAQUS中的Mooney-Rivlin本构模型[16]。内聚力单元参数K=500 GPa/mm，T1=T2=T3=60 MPa，本研究中K值、玻璃刚度T以及能量释放率G的选取参考Chen等[5]的参数设置。材料物性参数如表1所示。

表1 各材料物性参数
Table 1 Physical parameters of each material

Note: ν is the Poission’s ratio, D1, C10 and C01 are material parameters.

4 结果分析

4.1 仿真工况

为了更好地了解冲击荷载下LG裂纹的产生和扩展规律，分别研究了LG玻璃厚度和固有型内聚力单元参数K值对裂纹扩展过程以及下玻璃板中心位移的影响。在不考虑PVB胶层和玻璃板间黏附作用失效、锤头质量为144 g以及冲击速度为6.78 m/s的情况下，利用ABAQUS有限元软件分别模拟了4组不同K值和5组不同玻璃厚度的LG冲击破坏响应，并分析其结果。分别在以下两种情况下讨论：（1）当确定玻璃厚度为2 mm时，变化K值；（2）当K值为500 GPa/mm时，变化玻璃厚度。LG标准件模型（简称LG标件）参数：K值为500 GPa/mm，玻璃厚度为2 mm，PVB胶层厚度为0.76 mm。具体工况如表2所示。

表2 模拟工况
Table 2 The simulated cases

4.2 裂纹模态分析

图6对比了试验与仿真所得的冲击力随时间变化曲线，二者的趋势基本保持一致。由试验得到的冲击力-时间曲线可知：冲击开始时，冲击力迅速增大，并在冲击时间为0.081 ms时达到峰值583.9 N；随后冲击力随着冲击时间的增大而逐渐减小并趋于稳定。而仿真得到的冲击力-时间曲线中，冲击力峰值为653.5 N，对应的冲击时间为0.09 ms。将仿真得到的冲击力峰值与试验值进行对比发现，其相对误差为11.92%；同时冲击力峰值出现时刻二者也较为吻合。

图7给出了试验与仿真所得到的LG上表面的最终裂纹对比。试验和仿真所得的LG裂纹均包括局部粉碎裂纹、环向裂纹和径向裂纹，两者较为吻合；其中粉碎形裂纹均位于冲击点附近的圆形区域，试验、仿真中该区域直径与玻璃板宽度的比值分别为0.152和0.160，两者的相对误差为5.3%。

图6 试验与仿真所得冲击力曲线对比
Fig.6 Comparison of impact force curve obtained from test and simulation

图7 试验与仿真所得LG裂纹对比
Fig.7 Comparison diagram of laminated glass crack obtained by test and simulation

图8 冲击荷载作用下LG裂纹扩展的试验结果
Fig.8 Experimental result of LG crack growth under impact load

图8给出了LG标件在冲击荷载作用下的裂纹扩展过程以及试验结果。LG在冲击荷载作用30μs时，受到锤头冲击作用的LG中心发生凹陷，出现多条环径向裂纹，并产生许多细小的玻璃颗粒。随后径向裂纹继续向外扩展，并伴随少量相近径向裂纹间的环向裂纹产生。 330μs后，大量环向裂纹产生并伴随着径向裂纹向外延伸，继续产生更大半径的环向裂纹，裂纹最终延伸到玻璃板边缘，与试验获得的LG裂纹扩展模式较为吻合。

图9所示为LG破坏机理及其仿真结果。LG刚受到冲击荷载作用时，玻璃单元两端受弯矩作用，上、下玻璃板均上部受压下部受拉，但上玻璃板受压区域大于受拉区域，而下玻璃板受拉区域大于受压区域，其仿真结果模态未发生变化，如图9（a）所示；随着冲击过程的进行，下玻璃板受拉部分首先发生断裂，随即整个下玻璃板完全断开，仿真结果中下玻璃板单元间发生明显的相对位移而产生裂纹，如图9（b）所示；随后PVB胶层受力拉伸，上玻璃板受拉区域拉力增大并发生断裂，仿真结果中上玻璃板单元开始发生微小的相对位移，如图9（c）所示；最终上玻璃板完全断裂，仿真结果中上玻璃板发生明显的相对位移，LG在此处完全开裂，如图9（d）所示。单层玻璃的破坏机理与LG单层玻璃板的破坏机理相同，但由于PVB胶层的粘连作用使LG破坏后的玻璃碎片粘连在一起，不会发生飞溅，从而起到保护作用。从仿真结果可知，PVB胶层受到冲击作用后，其上下受压，中心受拉变形，由于未考虑PVB胶层的破坏，未对其设置初始损伤及其演化准则，因此PVB胶层没有发生断裂，呈现出类橡胶的超弹性材料特性。

图9 LG破坏机理图（左列）及其仿真结果（右列）
Fig.9 The failure mechanism diagram of LG (left column) and its simulation results (right column)

通过上述分析，比较试验和仿真所得到冲击力-时间曲线、LG裂纹形态以及破坏机理可知，仿真结果和试验结果较吻合。因此，采用基于固有型内聚力单元模拟LG在冲击荷载作用下裂纹扩展行为的方法是可行的。

4.3 罚刚度K对LG裂纹扩展的影响

图10所示为LG标件在不同K值下的裂纹模态。由图10可知，内聚力单元材料参数K值会影响裂纹扩展的模态。随着K值的增加，距LG中心较远的外围环径向裂纹更加明显，裂纹扩展的范围更大，但对于裂纹生成过程和数量影响不大。

图10 LG标件在不同K值下的裂纹模态
Fig.10 Crack modes of LG specimens at different K values

4.4 玻璃厚度对LG裂纹扩展的影响

图11给出了相同K值下不同玻璃厚度对LG裂纹扩展的影响。由图11可知：随着玻璃厚度的增加，玻璃表面环向及径向裂纹条数减少，其抗冲击性能明显提高；但玻璃厚度达到4 mm后，随着玻璃厚度的增加，各组裂纹个数差别不明显。因此，当玻璃厚度达到4 mm后再增加玻璃厚度，其抗冲击性能差异不明显，对玻璃抗冲击性能的提升不大。

4.5 下玻璃板的中心位移

图12显示在玻璃厚度为2 mm的情况下，不同K值对LG下玻璃板中心位移u的影响。随着K值的增加，LG下玻璃板的中心位移变化范围逐渐减小；但当K=1 250 GPa/mm时，LG下玻璃板的中心位移变大，与上述趋势不符。由此说明，K值影响玻璃板的刚度，K值越大玻璃板的中心位移越小；当K值大于1 000 GPa/mm以后，两者不成正相关关系。

图13显示在K值为500 GPa/mm的情况下，不同玻璃板厚度对LG下玻璃板中心位移的影响。随着玻璃厚度的增加，LG下玻璃板的中心位移逐渐减小；当玻璃厚度为3 mm时，LG下玻璃板中心开始发生回弹作用；当玻璃厚度为4 mm时，LG下玻璃板中心开始发生明显的回弹作用。因此，LG玻璃板越厚，回弹的位移越大，其刚度越大。

图11 不同玻璃厚度的LG裂纹模态（K=500 GPa/mm）
Fig.11 LG crack modes with different glass thicknesses (K=500 GPa/mm）

图12 不同K值的LG中心位移时程曲线
Fig.12 LG center displacement time-history curves with different penalty stiffnesses K

图13 不同玻璃厚度的LG中心位移时程曲线
Fig.13 LG center displacement time-history curves with different thicknesses of glass layer

5 结 论

在不考虑PVB胶层与玻璃板间的黏附作用的情况下，采用固有型内聚力模型，利用ABAQUS有限元软件模拟了LG在冲击荷载作用下的裂纹扩展过程，分析材料参数罚刚度K和玻璃厚度对LG响应的影响，得到以下结果。

（1）当LG玻璃厚度为2 mm、PVB胶层厚度为0.76 mm、K值为500 GPa/mm时，其裂纹扩展的过程是：当落锤接触玻璃表面时，玻璃表面迅速产生中心小范围的环向和径向裂纹并产生玻璃颗粒；随后径向裂纹首先向外扩展；环向裂纹在相邻径向裂纹间扩展，并随着时间的推移相互连接形成较大的环形裂纹。

（2）当玻璃厚度为2 mm、PVB胶层为0.76 m时，改变K值会影响LG裂纹扩展的模态以及下玻璃板的中心位移。随着K值的增大，LG裂纹扩展的范围增大，但裂纹数量及裂纹模态无明显变化；下玻璃板的中心位移减小，但当K值大于1 000 GPa/mm后，两者关系不明显。因此，K值对LG裂纹扩展过程影响不明显，但对于玻璃的刚度影响较大。

（3）当K值为500 GPa/mm、PVB胶层为0.76 m时，改变玻璃厚度也会影响LG裂纹扩展的模态和下玻璃板的中心位移。随着玻璃厚度增加，LG裂纹扩展的范围减小并且裂纹个数减少，同时下玻璃板的中心位移减小；当玻璃厚度为3 mm时，LG下玻璃板中心开始发生回弹作用；当玻璃厚度为4 mm时，LG下玻璃板中心开始发生明显的回弹作用。由此可知，玻璃厚度对LG裂纹扩展以及刚度有较为明显的影响。