冲击载荷下平板玻璃裂纹扩展和破坏形态的数值模拟

王木飞，李志强

（1. 太原理工大学机械与运载工程学院应用力学研究所, 山西 太原 030024；2. 山西省材料强度与结构冲击重点实验室, 山西 太原 030024；3. 力学国家级实验教学示范中心, 山西 太原 030024）

玻璃材料具有高强度、高透光性、低密度、耐高温等优异的性能，在民用建筑、运输工具、航空航天等各个行业得到了广泛的应用。在现实生产生活中，使用最普遍的玻璃类型是平板钠钙玻璃，研究其在强载荷下的动态响应具有重要的意义。在理论方面，Krauthammer 等[1]采用等效厚度和等效刚度将玻璃结构简化成等效单自由度质量-弹簧模型进行研究；Yuan[2]等将玻璃结构在空中爆炸载荷作用下的响应解耦为玻璃的弹性变形和PVB 的大塑性变形两个变形阶段。在实验方面，Holmquist 等[3]进行了较大应变率范围内无机玻璃的平板撞击实验，得到其状态参数；Daryadel 等[4]使用改进的分离式霍普金森压杆（split Hopkinson pressure bar，SHPB）研究了4 种不同的无机玻璃的动态压缩行为，得到其压缩强度、能量吸收值和最大失效应变。在数值模拟方面，裂纹扩展问题是计算力学领域长期存在的难点和热点。虽然有限元法在传统应力分析方面具有很大的优势，但是对材料失效断裂方面的模拟略显不足。这种不足源于经典连续力学理论控制方程中的偏微分方程，在裂纹尖端和裂纹面上会遇到偏微分方程空间导数不存在的情况，造成裂纹尖端应力的奇异性。为了解决上述数值求解上的困难，一方面，人们提出了带有失效准则的单元删除法[5]、节点分离法[6]和扩展有限元法[7]等多种有限元数值处理手段；另一方面，为了避免模拟多条裂纹带来的数值不稳定，人们采用诸如近场动力学（peridynamic，PD）[8]、离散元法（discrete element method，DEM）[9]、光滑粒子流体动力学（smoothed particle hydrodynamics，SPH）法[10]等无网格方法。例如：Sun 等[11]、Nie 等[12]利用SHPB 对玻璃材料开展高应变率动态加载实验，并利用断裂损伤力学方法，考虑损伤演化的本构关系，结合有限元方法（finite element method，FEM），对实验进行了数值模拟；Zang 等[13]采用SPH 与FEM 相耦合的方法模拟玻璃碎片的飞溅过程；Deng 等[14]采用多物质ALE 方法对中空玻璃幕墙的抗爆特性进行了大规模计算模拟；Bobaru 等[15]根据Bless 等[16]的实验方案和结果，使用PD 对多层钠钙玻璃在高速冲击载荷下的损伤进行了系统研究；Ren 等[17]在用键基PD 模拟动态脆性失效分析中应用了三维不连续伽辽金有限元方法，对冲击载荷下框支玻璃的动态裂纹扩展进行了数值模拟。

研究人员对玻璃材料的理论计算、实验分析和数值模拟的处理手段丰富，取得了很多优异的成果。但是，玻璃材料在冲击荷载下的细观裂纹数值模拟研究尚不全面，鲜有文献将不同算法的模拟结果进行对比研究。本工作将采用单元删除法、不连续伽辽金近场动力学（discontinuous Galerkin peridynamic，DG-PD）法和无网格粒子近场动力学（meshless peridynamic，M-PD）法，以常见的平板玻璃板为例，研究其在子弹冲击载荷作用下的动态响应。通过定量定性分析，对比单元删除法、DG-PD 法和M-PD 法在模拟裂纹扩展方面的特点，并结合相关文献对结果的有效性进行验证。

1 单元删除法在平板玻璃中的应用

单元删除法是指模型中的单元满足预先设定的失效准则时，程序会自动删除相应的单元和节点，当一系列单元和节点被删除而相互贯通后就形成了可视的裂纹路径。其中，最常用的方式是添加最大主应力、最大主应变或最大拉应力等失效准则，在加载过程中若满足任何一种预置的失效准则，则材料发生破坏；另一种方式是采用有限元软件中自带失效准则的材料模型，将失效准则嵌入相应的材料本构，无需添加额外的失效准则。单元删除法在一定程度上适合模拟玻璃材料的损伤破坏，删除的单元可被看作产生裂纹。单元删除法虽然操作简单，计算效率较高，但是无法保证系统的质量、能量和动量守恒，使得体系的运动控制方程不成立，从而很难准确地反映真实的材料特性。为了模拟出裂纹的细观特征，应设置合适的网格密度和材料失效阈值，以确保网格在合适的时候才被删除，从而提高模拟结果的精确度。

1.1 有限元模型

使用TrueGrid 软件作为前处理器进行参数化模型的建立，整个模型结构采用kg-mm-ms 单位制。单层平板玻璃的尺寸为600 mm×600 mm×6 mm（长×宽×厚），模型四周采用30 mm 宽的框支结构固定边界。玻璃部件的面内网格尺寸为3 mm×3 mm，考虑到玻璃材料弯曲效应的影响，在玻璃厚度方向划分为3 层网格，冲击头采用常规子弹形状。整个结构的有限元模型及弹头放大图如图1 所示。

图1 有限元模型及弹头放大图Fig. 1 Finite element model and enlarged view of the bullet

为了研究单层平板玻璃在不同子弹冲击速度下的动态损伤行为，分别采用带有失效准则的单元删除法和PD 法进行数值模拟。玻璃属于微弹脆性材料，在发生破坏时产生的应变很小，为了能够直观地观察裂纹演化过程，在单元删除法、DG-PD 法和M-PD 法中分别采用最大主应变云图、等效塑性应变和损伤因子表示损伤程度。

1.2 材料模型和损伤准则

目前，模拟玻璃、陶瓷等脆性材料时使用最广泛的动态本构模型是Johnson 和Holmquist 共同提出的*MAT_JOHNSON_HOLMQUIST_CERAMICS（JH-2）模型，它能够很好地反映脆性材料在高应力、大变形下的动态破坏行为。该模型描述的损伤演化进程为：加载初期，材料表现为弹性性质；当应力达到预先设定的屈服极限时，材料开始出现损伤；当损伤达到一定程度时，材料出现破碎。采用JH-2 材料模型和单元删除法，同时施加最大主应变和最大主应力失效准则来模拟裂纹的扩展，最大主应变和最大主应力分别设为0.001 2 和65 MPa[18]。

在JH-2 材料模型中，陶瓷等类型材料的等效屈服应力σ*可以表示为

表1 玻璃材料JH-2 材料模型参数[6]Table 1 JH-2 material model parameters of glass material[6]

边框和弹头均选用钢材料，采用LS-DYNA 中的刚体材料模型描述，表2 给出了边框和弹头的刚体材料模型参数，其中：ν 为泊松比，E为弹性模量。

表2 边框和弹头刚体材料模型参数Table 2 Rigid body material model parameters of frame and bullet

1.3 裂纹扩展和破坏形态

为了探究平板玻璃在子弹冲击载荷下的裂纹扩展和破坏模式，结合单元删除法和JH-2 材料模型进行数值模拟分析，弹体的冲击速度分别为5、8、12 和15 m/s。为了直观展现平板玻璃裂纹萌生、扩展进程，图2 给出了15 m/s 冲击速度下迎撞面和背撞面的损伤演化过程。由图2 可以看出：单元删除法呈现出明显的十字形裂纹，冲击接触区域出现空洞，但受自身算法的限制，并未出现玻璃碎片飞溅；裂纹在冲击接触点萌生，并不断向前扩展，与唐剑[6]采用1/4 对称模型所得模拟结果不同的是，裂纹并非完全对称，且扩展的先后顺序和长短略有不同；背撞面在4 ms 时刻出现十字形裂纹，明显早于迎撞面。这可能是由于在玻璃厚度方向划分了多层网格，从而考虑了其弯曲效应，再叠加四周夹支的边界条件，产生了反射拉伸波，进而造成背撞面损伤加剧。

图2 迎撞面和背撞面的损伤演化Fig. 2 Damage evolution of the top face and bottom face

图3 从定性角度展示了上述4 种不同冲击速度下平板玻璃的最终损伤模态。冲击速度为5 m/s时，弹头和玻璃接触处仅有淡淡的弹痕，未见裂纹萌生；冲击速度为8 m/s 时，出现明显的十字形裂纹，弹头未冲破玻璃而出现反弹现象；冲击速度为12 m/s 时，裂纹出现分叉现象，弹头微微嵌入玻璃内部；冲击速度为15 m/s 时，裂纹的对称性较好，弹头穿过玻璃，且边框处损伤较严重。由此可见，采用单元删除法模拟得到的玻璃裂纹主要以十字形裂纹形式扩展，在损失质量的前提下会造成裂纹的不稳定性。

图3 不同冲击速度下平板玻璃的最终损伤模态Fig. 3 Final damage of glass plates at different impact velocities

图4 从定量角度展示了不同冲击速度下损伤参数的对比。从玻璃和弹头接触区域损伤面积来看，5 m/s 速度下的玻璃几乎未受到破坏；8～12 m/s 速度下区域损伤面积的变化不大；而15 m/s 速度下玻璃在弹头的贯穿作用下，损伤面积明显增大。从玻璃最大裂纹扩展距离来看，5～8 m/s 速度下，最大裂纹扩展距离与速度大小正相关，8～15 m/s 速度下裂纹增长幅度很小。在冲击速度增大的过程中，玻璃结构中的应力和应变也随之增大，当达到预设的损伤阈值时，材料开始断裂。从子弹嵌入深度来看，在低于8 m/s 的冲击速度下，子弹将出现与速度正相关的反弹距离；随着冲击速度继续增大到12 m/s 时，弹头将嵌入玻璃内部；15 m/s 的冲击速度下，子弹完全贯穿玻璃。

图4 不同冲击速度下损伤参数的对比Fig. 4 Comparison of damage parameters at different impact velocities

综上所述，弹头的冲击速度对玻璃的裂纹扩展和破坏模态具有重要的影响；单元删除法可以在一定程度上模拟出玻璃结构的损伤形态，且计算效率较高，但是在捕捉裂纹分叉、贯通等细观方面的精确度不高。值得注意的是，模拟30、50 m/s 等中高速下单层平板玻璃的破坏时，网格畸变和单元负体积会造成计算终止，捕捉的裂纹形态也很不理想。由此可见，单元删除法在模拟中高速材料损伤破坏方面存在一定的缺陷。

2 DG-PD 法在平板玻璃中的应用

传统的有限单元离散是一种连续的网格划分，相邻单元之间共享节点和边界，阻止了单元之间的分离，从而无法进行断裂模拟。而不连续伽辽金有限单元中的相邻单元各自拥有独立的节点，在断裂发生前，它们的节点位置相互重合，断裂发生后，每个单元拥有各自的节点。而在使用有限单元离散PD 模型时，只有单元的高斯点之间具有键的连接关系，模拟断裂时，也只考虑高斯点之间的键是否断裂。图5 为二维有限单元断裂前后的节点分离示意图。虽然不连续伽辽金法比传统的有限元网格增加了节点的个数，但是相比经典连续介质力学模型，在PD 模型中引入不连续伽辽金单元有其天然的优越性。因为经典连续介质力学理论包含了位移的导数，而PD 理论不需要包含位移的导数，且不连续伽辽金单元离散在单元边界处允许位移的不连续。如果在经典连续介质力学模型中应用不连续伽辽金单元以模拟断裂问题，就需要考虑单元与单元边界处的通量约束，从而大幅增加计算的复杂度。

图5 二维有限单元断裂前后的节点分离示意图Fig. 5 Node separation before and after fracture of 2D finite elements

2.1 有限元模型和材料参数

采用与1.1 节中相同的有限元模型，采用DG-PD 方法研究平板玻璃的抗冲击性能，边框和弹头保持原有模型不变，唯一不同的是需要对平板玻璃进行节点分离。节点分离操作会造成节点增多，与单元删除法相比，计算量会明显地增大。LS-DYNA 提供了简单的节点分离操作，只需对目标部件中的节点集进行Detach 前处理即可。玻璃材料的DG-PD 模型的关键参数如表3 所示[15]，其中：Gt为临界断裂能量释放率；HSFAC 表示标准化的支持域大小，LS_DYNA 会自动调整HSFAC，确保材料点合适的邻域。这里，最重要的是断裂失效准则的选取，即断裂能量释放率，而临界能量释放率与临界伸长率直接相关，当伸长率超过临界值时，点对之间的键就不再有相互作用，并且不可恢复。

表3 玻璃材料DG-PD 材料模型参数Table 3 DG-PD material model parameters of glass material

2.2 裂纹扩展和破坏形态

采用DG-PD 法对平板玻璃的抗冲击性能进行数值模拟，采用等效塑性应变反映损伤状态。图6、图7、图8、图9 分别给出了弹体冲击速度为5、15、30、50 m/s 时不同时刻玻璃迎撞面的损伤演化。5 m/s 冲击速度下，玻璃出现轻微的十字形裂纹，弹头与玻璃接触的局部区域损伤严重；15 m/s 冲击速度下，十字形裂纹扩展到玻璃结构边界处，破坏程度进一步加大，有数道裂纹开始相互交叉和贯通；30 m/s 冲击速度下，玻璃呈现复杂的破坏模式，中心点损伤严重，产生大量的粒子碎片，边框处出现环状裂纹，玻璃内部多条裂纹分叉、贯通，细节捕捉效果较好；50 m/s 冲击速度下，中心点损伤区域进一步扩大，裂纹的数量明显增多，大量粒子开始飞溅。

图6 5 m/s 冲击速度下迎撞面的损伤演化Fig. 6 Damage evolution of the impact face at the impact velocity of 5 m/s

图7 15 m/s 冲击速度下迎撞面的损伤演化Fig. 7 Damage evolution of the top face at the impact velocity of 15 m/s

图8 30 m/s 冲击速度下迎撞面的损伤演化Fig. 8 Damage evolution of the top face at the impact velocity of 30 m/s

图9 50 m/s 冲击速度下迎撞面的损伤演化Fig. 9 Damage evolution of the impact face at the impact velocity of 50 m/s

综上所述，弹头冲击速度对平板玻璃的裂纹扩展和破坏模式有很大的影响，玻璃结构整体以十字形裂纹为主，随着速度的不断增大，裂纹扩展形式逐渐复杂化。此外，与单元删除法相比，DG-PD 法不仅保证了结构的质量守恒，还可以捕捉到玻璃碎片的飞溅，裂纹扩展形式更加多样。

图10 给出了不同冲击速度下中心点单元的速度变化曲线。从图10 可以看出，在5～15 m/s的低速冲击下，单元速度会在短时间内急速增大至与弹头相同的速度，并保持相对长的一段时间；在30、50 m/s 相对高速冲击下，单元速度在短时间内提升至42、59 m/s，略高于弹头的冲击速度，并在此速度下来回振荡。这可能是由于在高速冲击下单元和弹头不再接触而出现了分离，导致单元的进一步加速。此外，随着弹头冲击速度增大，单元速度增至最大速度所需的时间降低。由此可以解释1.2 节中在中高速冲击下单元删除法算法不稳定的原因，更加突出了DG-PD 数值算法的稳定性。

图10 采用DG-PD 法得到的不同冲击速度下中心点单元的速度变化Fig. 10 Velocity change of center element obtained by DG-PD method at different impact velocities

3 M-PD 法在平板玻璃中的应用

M-PD 法是基于离散粒子方法建立模型，并赋予这些粒子质量密度、速度、体积（V）等物料信息。粒子间以非局部力的方式在一定范围内相互作用，这些长程力中包含了材料所有的本构关系，图11 为M-PD 非局部相互作用示意图。通常情况下，采取施加临界伸长率S0来模拟玻璃结构的损伤破坏行为，并通过物质点已破坏的键与所有键的比值体现材料的损伤程度。此外，与DG-PD 法相比，M-PD 法无需进行单元节点分离操作。

图11 M-PD 非局部相互作用示意图Fig. 11 Schematic diagram of M-PD non-local interaction

3.1 粒子模型和材料参数

对1.1 节中的平板玻璃结构有限元模型进行“.k”文件的python 读写操作，将有限元模型中的每个单元节点转化为近场动力学中相应的离散粒子形式，包含坐标、体积和材料编号等物料信息，即(x,y,z,V, Block_ID)。与直接建立离散粒子模型相比，对有限元模型读写操作的转换更适合于复杂粒子模型的建立。边界条件的确定与上述两种方法略有区别，采用的是对玻璃部件边缘30 mm 处的粒子，以约束所有的自由度来代替框支结构，无网格粒子模型如图12 所示。

图12 M-PD 粒子模型及弹头放大图Fig. 12 M-PD particle model and enlarged view of the bullet

玻璃粒子离散间距设为 Δx=3 mm，近场半径为δ=3.01Δx。根据弹性模量和近场半径[15]计算出玻璃材料的微模量c=4.19×1019N/m6，得到键的临界伸长率S0=7.17×10-5。玻璃材料M-PD 模型的关键参数如表3、表4 所示。

表4 玻璃材料M-PD 材料模型参数Table 4 M-PD model parameters of glass

3.2 裂纹扩展和破坏形态

采用M-PD 法对平板玻璃的抗冲击性能进行数值模拟，采用损伤因子D来反映损伤状态。图13展示了30 m/s 冲击速度下平板玻璃的损伤演化。可以看出，玻璃构件在弹头的冲击下出现裂纹萌生、扩展和贯通等现象，冲击点附近局部区域损伤最严重，并伴有径向裂纹和环状裂纹。0.4 ms 时，玻璃最先产生十字形裂纹；2.2 ms 时，出现大量径向裂纹，边框处开始有轻微的损伤；7.2 ms 时，裂纹不断向外扩展，损伤程度也逐渐加重，向外延伸的裂纹相互贯通形成环状裂纹，边框处也形成一道明显的环状裂纹；20 ms 时，裂纹形状几乎没有任何改变，只是损伤程度有所增强。

图13 30 m/s 冲击速度下迎爆面的损伤演化Fig. 13 Damage evolution images of the top face at the impact velocity of 30 m/s

3.3 近场域大小对平板玻璃动态响应的影响

近场域的大小是近场动力学的关键所在，对近场动力学与连续介质力学中的接触力模型做出了区分，从而避免了传统有限元控制方程中微分的奇异性。为了探究近场域大小对平板玻璃动态响应的影响，以冲击速度30 m/s 为例，分别模拟了常见的3 种粒子间距下平板玻璃的破坏情况，如图14 所示。就损伤情况来看，2 倍粒子间距下的玻璃损伤程度远远弱于3 倍和4 倍粒子间距，捕捉裂纹数量较少，而3 倍和4 倍粒子间距下的裂纹模态几乎相同；就计算时长来看，3 种工况的计算时长为1.5、2.7 和5.2 h，随着粒子间距的增大，所需的计算资源显著增加。综上所述，在保持计算精度的情况下，应该尽量降低计算成本，故3 倍粒子间距是相对优异的选择，与文献[8]中的观点一致。

图14 不同近场域大小下的平板玻璃损伤Fig. 14 Damage of glass plate under different horizon sizes

3.4 冲击速度对平板玻璃动态响应的影响

冲击速度对平板玻璃的损伤模式具有决定性的影响，是研究玻璃材料抗冲击性能最直观的体现。

图15 不同冲击速度下平板玻璃的最终损伤破坏Fig. 15 Final damage of glass plate at different impact velocities

4 结 论

采用3 种不同方法对平板玻璃在冲击载荷下的动态响应进行数值模拟，通过对比分析，得出以下结论。

(1) 单元删除法作为经典连续力学中模拟脆性材料的常用方法，可模拟玻璃在冲击载荷下的主要形貌，即十字形裂纹，但在捕捉裂纹分叉和贯通等方面存在不足；低速冲击下裂纹形貌的变化不大，而在中高速冲击过程中，会出现网格畸变和计算不收敛的情况；模拟过程中未见玻璃碎片飞溅，无法通过碎片飞溅情况评估其安全性能。

(2) DG-PD 法进一步丰富了裂纹的扩展形式，可以同时模拟出低、中、高速下玻璃的裂纹扩展模式，但商业软件中的可选参数有限；通过DG-PD 法可以获得明显的环状裂纹和径向裂纹，且裂纹具有很高的对称性，能捕捉到裂纹的分叉、贯通等现象；冲击点和边框处有大量玻璃碎片飞溅，通过监测碎片飞溅情况，可为结构防护提供参考。

(3) M-PD 法通过粒子的形式建立模型，可通过调节近场域、本构模型和断裂准则等可选参数充分发挥近场动力学的优势；M-PD 法能捕捉到径向裂纹和环状裂纹，裂纹的对称性较好，与DG-PD 法得到的模拟结果具有很高的一致性；对比单元删除法和DG-PD 法获得的损伤模态，并综合考虑计算资源，近场域取3 倍粒子间距较为合适；低速冲击下的裂纹形态单一，而高速冲击下的裂纹扩展形态复杂。