玻璃准静态破坏现象隐式和显式组合仿真研究

王 程，臧孟炎

(华南理工大学 机械与汽车工程学院，广州 510640)

据公安部网站统计，2020年上半年，我国汽车保有量已经达到2.7亿辆。风挡玻璃是重要的安全部件，由于生产过程中不可避免的质量缺陷而导致强度分布不均，汽车在行驶过程中，巨大载荷所引起的车架变形或者风压所形成的较大风载可能会导致玻璃的破坏，成为威胁司乘人员生命安全的一个重要因素[1-2]。所以，可靠高效的汽车风挡玻璃的准静态破坏性能评价对于保证汽车行驶的安全性至关重要。

关于汽车玻璃的准静态破坏问题，相关学者做了大量研究工作。玻璃的准静态破坏问题研究常用的试验方法是三点或四点弯曲方法，但由于玻璃存在质量缺陷，这两种方法无法准确评价玻璃的准静态破坏性能。因此，VITMAN等[3]最早提出了用于一定边长玻璃板的同轴双环加载试验，可准确评价玻璃的准静态破坏性能。RITTER等[4]使用同轴双环加载试验，评价了不同尺寸玻璃的准静态破坏性能。GULATI等[5]对新制的和老化的汽车玻璃进行同轴双环加载试验，研究玻璃老化对汽车玻璃破坏强度的影响。MÜLLER-BRAUN等[6]对边长为125 mm的平板玻璃进行了同轴双环试验，同时使用ANSYS隐式数值分析方法，研究玻璃的径向及切向破坏强度的分布。CASTORI等[7]对边长为400 mm的玻璃板做了不同支撑情况下的同轴双环试验，同时进行隐式有限元仿真分析等效应力分布，研究一种能够预测玻璃破坏强度的数学模型。综上所述，对于玻璃的准静态破坏问题研究，多数采用试验方法。即使有限元显式分析方法在风挡玻璃冲击破坏方面得到了广泛应用[8]，但玻璃准静态破坏的仿真分析尚无有效计算方法对应。主要原因在于准静态破坏问题现象时间过长，显式有限元方法难以实现，隐式有限元方法在破坏发生后的收敛性难以保证且计算效率极低。所以，如何准确且高效地对玻璃准静态加载破坏进行仿真分析，对提升汽车玻璃的开发设计水平具有重大意义。

基于上述考虑，提出玻璃准静态破坏仿真的隐式和显式组合仿真方法，充分发挥两种计算方法的优势。即在破坏前的加载阶段采用隐式有限元方法以较大的时间步长分析玻璃的变形，在即将发生破坏时改用显式有限元方法分析破坏现象。以玻璃的同轴双环加载试验为研究对象，使用内聚力模型描述玻璃裂纹的发生和传播，比较组合仿真分析方法和单纯隐式分析方法的计算精度和计算效率，为汽车风挡玻璃准静态破坏性能评价的有限元仿真应用提供指导和帮助。

1 内聚力模型

内聚力模型法是一种基于断裂力学，预测脆性材料的动态断裂的有效方法，由于可以描述裂纹的扩展和考虑破坏发生的能量耗散，近年来广泛应用于汽车风挡玻璃等脆性材料的冲击破坏仿真研究工作中[8-10]。内聚力模型首先由DUGDALE[11]和BARENBLATT[12]提出，BARENBLATT将内聚力模型应用于脆性材料的断裂中。HILLERBORG等[13]首次将内聚力模型应用到有限元计算中，模拟了脆性材料的断裂过程。

双线性型固有内聚力模型由MI等[14]提出，其法向和切向的本构曲线如图1所示。其牵引-分离法则假定：在外载荷的作用下，界面上作用牵引力随着裂纹面间相对位移的增大而线性增加。当内聚区域的牵引力达到最大牵引力之后，内聚区的裂纹面相对位移会持续增加，但界面上的作用牵引力开始线性下降，此时材料刚度下降，发生软化，裂纹开始萌生及扩展。当牵引力降为0时，裂纹面相对位移达到最大值，此时材料已经发生破坏。

图1 双线性型内聚力模型的法向及切向本构曲线

双线性型固有内聚力模型的本构方程为

式中：Tn和Tt分别表示裂纹内聚区域的法向牵引力和切向牵引力，对应的破坏应力为σmax和τmax；δn和δt分别为裂纹内聚区域两裂面间的法向及切向分离量；和分别为裂纹开始萌生时内聚区域两裂面间的法向及切向分离量；和分别为裂纹内聚区域两裂面间的法向及切向最终分离量。断裂能为双线性型固有内聚力模型的本构曲线与分离量坐标轴包围的三角形面积的大小，即：

而在实际的材料破坏中，材料的裂纹不是单一类型的破坏，而是某种混合模式。因此，为描述混合模式下的材料破坏，BENZEGGAGH和KENANE等[15]提出了一种开裂准则，称为BK开裂准则。此时，内聚区域裂纹的总断裂能为：

式中：β为破坏模式的混合度，一般为0.5；XMU为材料常数。这种材料开裂准则由于简单易用，所以广泛地用于复合材料的仿真破坏中。目前，有限元商用分析软件ABAQUS和LS-DYNA的隐式和显式计算方法中，均已嵌入表征玻璃脆性破坏的内聚力模型。

2 组合仿真分析方法

仔细研究夹层玻璃的准静态破坏问题，发现可以将这个过程分解为加载阶段和破坏阶段。其中加载阶段往往需要10～102 s甚至更长的现象时间，而破坏阶段的现象时间不足1 s。这是因为作为脆性材料的玻璃，发生破坏时裂纹的扩展速度极快[16]。如果单纯使用隐式求解器分析夹层玻璃的准静态破坏过程，则由于玻璃破坏的现象时间极短，极易导致玻璃破坏分析过程中出现收敛性问题。即便使用其他方法使其收敛，计算步长也会大幅度缩短，收敛性迭代计算效率将严重下降。反之，若单纯使用显式求解器分析，虽然无需考虑收敛性问题，但玻璃准静态加载的过长现象时间对应显式分析方法微秒级的时间步长，导致需要计算的步长数过多，发生数据溢出无法计算。所以，如果使用隐式分析和显式分析相结合的方法，即加载阶段使用隐式分析方法，而破坏阶段使用显式分析方法，就能充分发挥两个方法的优势：加载阶段使用隐式有限元方法可以有较大的时间步长，破坏阶段使用显式有限元方法避免了复杂的收敛性计算。

基于ABAQUS软件的玻璃准静态破坏现象的组合仿真分析方法如下：

（1）在ABAQUS的UVARM子程序中，定义玻璃内聚力单元等效拉应力破坏强度，且满足破坏强度后立即终止计算。通过编译得到新的ABAQUS执行文件。

（2）对嵌入内聚单元的玻璃准静态破坏仿真模型M1进行隐式有限元计算。当某个内聚单元满足破坏强度，UVARM子程序命令计算立即终止，得到夹层玻璃破坏开始时间。

（3）以玻璃破坏开始时间为计算结束时间，定义新的玻璃准静态破坏有限元模型M2，重新进行隐式计算。

（4）将M2计算最终时刻得到的玻璃节点位移和单元应力作为初始条件导入到显式计算仿真模型M3中，分析玻璃的破坏阶段。

3 数值仿真与评价

3.1 仿真模型

基于玻璃的同轴双环加载试验[7]，考虑结构和载荷的对称性，建立边长为400 mm，厚度为8 mm的1/4同轴双环加载试验仿真模型。其中厚度8 mm玻璃采用4层全积分六面体单元，网格分布为辐射状网格，玻璃单元数量为24 256个，如图2所示。在玻璃网格公共界面间插入零厚度的内聚力单元，使用内聚力单元模拟玻璃的破坏。加载环和支撑环的半径分别为75 mm和150 mm，它们的截面半径为2.5 mm。支撑环定义固定约束，给加载环定义0.023 mm/s的向下恒定速度，得到玻璃准静态加载模型M1，如图3所示。

图2 玻璃网格划分

图3 厚度8 mm玻璃同轴双环加载模型

玻璃选择线弹性材料模型。加载环和支撑环定义为刚体，材料参数均参照普通钢材。材料物性见表1[7]。

表1 玻璃、刚体环的材料物性

玻璃的内聚力模型定义为双线性型固有内聚力模型，破坏强度取99 MPa[7]， I型破坏的断裂能为10 N/m，II型、III型破坏的断裂能为50 N/m[8]。内聚力单元的破坏准则使用BK准则，XMU系数为2，罚刚度取500 GPa/mm[8]。

加载环和支撑环与玻璃间的接触，定义为面与面接触（Surface to Surface Contact），摩擦因数为0.1。在破坏阶段显式分析模型中，添加玻璃单元自接触，类型为通用接触（General Contact），摩擦因数为0.9，已将玻璃裂纹间的自接触考虑在内。

3.2 组合仿真计算

首先以M1模型隐式计算确定玻璃的初始破坏发生时间，然后以M2模型第2次隐式计算得到玻璃破坏发生前的等效应力状态，如图4所示。

图4 破坏发生前的等效应力云图

读入隐式分析结果，建立破坏阶段显式分析仿真模型M3，进行仿真分析，最终得到玻璃的破坏模式、破坏载荷和破坏位移。

3.3 仿真结果分析及比较

3.3.1 破坏模式分析

厚度8 mm玻璃同轴双环加载仿真模型使用组合计算方法和纯隐式计算方法得到的玻璃破坏仿真结果与试验结果，如图5所示。由图5可知，组合仿真计算方法得到的玻璃破坏模式图5a与试验结果图5c呈现良好的一致性：加载环内部区域玻璃发生全面粉碎性破坏，加载环和支撑环之间的玻璃产生大量径向裂纹且加载环和支撑环附近有显著的周向裂纹，而加载环外部区域玻璃的径向裂纹明显减少。由图5b可知，纯隐式仿真所得到的玻璃破坏情况仅在加载环内部区域与试验结果一致，其他区域的玻璃没有观察到破坏现象。

将图3所示模型的玻璃厚度由8 mm变为6 mm（删除一层玻璃六面体单元），使用组合仿真方法得到的厚度6 mm玻璃双环试验仿真结果与试验结果如图6所示。由图6可知，组合仿真方法仍然得到了与对应试验一致性良好的分析结果。

图5 两种仿真方法得到的破坏现象与试验结果

图6 厚度6 mm玻璃双环试验仿真与试验结果

综上所述，就破坏模式来看，纯隐式计算方法显然不能有效预测复杂的玻璃准静态破坏问题，而组合分析方法得到了一致性的仿真计算结果。

3.3.2 破坏载荷与破坏位移分析

表2为厚度8 mm与6 mm玻璃同轴双环加载试验和仿真所得到的破坏位移和破坏载荷对比。由表2可知，无论是破坏位移还是破坏载荷，仿真与试验结果的误差都在10%以内，呈现良好的一致性，从定量方面说明了组合仿真分析方法的有效性。

表2 试验与仿真得到的破坏位移和破坏载荷

3.3.3 计算效率分析

前述仿真结果使用具有Intel(R) Xeon(R) Silver 4210 2.2 GHz处理器的计算服务器，8核并行计算。表3为厚度8 mm和厚度6 mm玻璃模型分别使用组合分析方法和单纯隐式分析方法的计算时间对比。尽管组合分析方法在加载阶段计算了两次。由表3可知，无论是厚度8 mm玻璃模型还是厚度6 mm玻璃模型，组合仿真方法所消耗的CPU时间要远小于单纯使用隐式计算的总计算时间。这是因为玻璃同轴双环加载试验的破坏现象时间极短，而隐式算法在分析玻璃的破坏现象时，需要大幅降低计算步长进行迭代计算，而且一个步长的迭代计算需要多次迭代才能收敛，导致计算时间延长，计算效率降低。

表3 组合仿真和纯隐式仿真的计算时间对比

4 结论

（1）针对玻璃准静态破坏仿真问题，提出了隐式和显式组合仿真分析方法。通过用户材料子程序确定夹层玻璃初始破坏发生时间，通过两次隐式计算和一次显式计算，利用隐式计算方法和显式计算方法分别在玻璃加载阶段和破坏阶段的优势，实现提高精度和计算效率的目的。

（2）以厚度为8 mm和厚度为6 mm的玻璃同轴双环准静态加载破坏试验为研究对象，分别使用组合仿真分析方法和纯隐式仿真方法进行了仿真分析。组合仿真分析方法的破坏模式与试验结果基本一致，破坏位移与破坏载荷与试验结果误差小于10%，计算效率显著高于纯隐式计算方法，充分说明了组合仿真方法在玻璃准静态破坏现象仿真分析方面的有效性。