高伟 臧孟炎，2†

(1.华南理工大学 机械与汽车工程学院，广东 广州 510640;2.湖南大学 汽车车身先进设计制造国家重点实验室，湖南 长沙 410082)

前挡玻璃是汽车重要零部件之一，通常使用由上下两层玻璃和树脂材料(PVB)压制而成的夹层玻璃.夹层玻璃的力学特性，尤其是破坏特性在交通事故中对乘客的安全有很大影响，因此，学者们对汽车夹层玻璃的破坏特性进行了大量研究.许俊等［1-2］使用试验方法，对汽车夹层玻璃的冲击破坏问题进行了深入研究.Dharani 等［3］使用非线性有限元方法分析了汽车前挡玻璃受到人头冲击后的动态响应.Flocker 等［4］研究了夹层玻璃中的锥形裂纹及其对夹层玻璃受力的影响.Xu 等［5］使用扩展有限元方法研究了夹层玻璃受人头模型冲击后的破裂形式.由于玻璃的脆性破坏特性，可以描述连续介质向非连续介质转化的离散元方法被认为是一种更为有效的仿真分析方法.Oda 等［6-7］将离散元法应用于汽车玻璃梁的冲击破坏问题的仿真研究，建立了汽车玻璃梁的二维离散元模型，并对一系列冲击破坏问题作了深入的研究，取得了许多有意义的研究成果.Zang等［8-9］在此基础上建立了三维离散单元模型，开发了离散元与有限元耦合方法，使用摩尔-库伦准则，仿真研究了夹层玻璃的冲击破坏过程.当两离散单元之间的作用力满足摩尔-库伦准则时，其间弹簧断开，弹簧连接力为0，这样影响计算的稳定性与精确性.

考虑裂口发生和发展的内聚力模型广泛运用于有限单元法中，用来模拟脆性材料的破坏问题［10-12］.文中基于内聚力模型的基本概念，提出了适用于三维球形离散单元的破坏模型.将该破坏模型在自主开发的离散元与有限元耦合计算程序CDFP［9，13］中实现，通过比较夹层玻璃梁冲击破坏过程的仿真分析与试验结果，确认了仿真分析方法的有效性.

1 离散单元的运动方程

三维球形离散单元方法的基本思想是将连续体离散成大小相等的刚性球，即球形离散单元.离散单元受到周围单元的作用和外力，其运动由牛顿定律确定［14-15］.例如，无阻尼作用的单元i 的运动方程为

式中，mi、Ii为单元i 的质量和惯性张量，t、g 为时间和重力加速度向量，ui、ωi为单元i 的位移向量和角速度向量，Ri、fik分别为单元i 受到的外力和单元k的作用力，rik为fik的力矩，Ki为外力偶，Ni为与单元i 相邻的离散单元数目.值得说明的是，式(2)适用于球形(圆形)离散单元，角速度分量相互不耦合.

2 离散单元之间的相互作用

离散单元之间的相互作用一般分为两类:连接型和接触型.连接型考虑单元间没有间隙且符合变形协调条件，主要用来处理连续介质力学问题.接触型是散体特有的作用形式，单元间的作用力为接触面上的接触力.本质上讲，脆性材料的破坏过程是连续介质向非连续介质的转化过程［15］.为了更好地模拟这一转化过程，文中在连接型和接触型之间增加了一种作用类型——内聚型.当离散单元满足破坏准则时，其作用类型由连接型转化成内聚型，内聚力由内聚模型求出.

2.1 离散单元之间的应力计算

离散单元之间的法向应力σn和切向应力σs按以下公式计算:

式中，Fn、Fs分别为两离散单元的法向和切向连接力，A 为相邻离散单元之间的等效作用面积.

当法向应力满足以下关系时，连接型转化成内聚型:

式中，σc为内聚强度.

2.2 内聚力模型

常用的内聚力模型包括指数型［10］、双线性［12，16］、单调递减内聚强度型［11，17-18］等，文中使用由Ortiz 等［19］提出的内聚力模型.该模型通过有效量将内聚力与单元之间的分开量联系起来.

单元之间的有效分开量δeff为

式中，δn、δs为法向和切向分开量，η 为切向位移因子.

最大有效分离量按下式计算:

由图1 和式(7)可以看出，最大有效分开量单调增加，并可用来判断加载、卸载.如果δeff=且ff≥0，则为加载，否则，为卸载.图1 中，δc为临界分开量.

图1 有效内聚力与有效分离量关系图Fig.1 Relationship between effective traction and separation

单元之间的有效应力σeff为

式(9)、(10)可以写成以下形式:

由式(7)可知，δeff≤恒成立，故式(11)可写为

当最大分开量δmaxeff大于临界分开量δc时，有效应力σeff等于0，即

单元间的内聚应力为

单元间的法向、切向内聚力分别为

内聚力模型主要受3 个参数影响:内聚强度σc、临界分开量δc、剪切应力因子η.一般情况下，内聚强度取材料的抗拉强度［20］.临界分开量为

式中，Gc为临界能量释放率.

笔者已将该破坏模型在自主开发的离散元有限元耦合计算程序CDFP 基础上编程实现，以进行脆性材料的破坏仿真分析.

3 夹层玻璃梁冲击破坏试验仿真

为考察汽车玻璃冲击破坏过程中裂纹发生和扩展的机理，笔者等进行了夹层玻璃梁冲击破坏试验［21］.由于普通汽车夹层玻璃厚度太小，难以捕捉到上下两片玻璃的破坏过程，因此该试验定制了厚度约为普通汽车玻璃4 倍的试验片(两块玻璃各厚10 mm，PVB 厚4 mm、长度为200 mm、宽度为10 mm).

图2 为夹层玻璃梁冲击破坏过程的照片［21］.由于高速相机容量的限制，进行了两次重复试验，图2是两次重复试验的结果.图2(a)中照片的时间间隔为20 μs，可以清楚地看到冲击侧玻璃裂纹的发生及传播.图2(b)中照片的时间间隔是100 μs，它记录了上下层玻璃裂口发生及传播的全过程.由于夹层玻璃试验样片一致性较差，导致图2 中的2 次重复试验结果裂纹形式差异较大，但是冲击发生后首先在冲击侧玻璃下表面产生向上传播的裂口、然后再发生下层玻璃自下向上的破坏，破坏趋势是一致的.仿真分析的目的在于利用开发的三维离散元内聚力破坏模型获得与试验现象一致的破坏趋势.

图2 夹层玻璃梁冲击破坏裂纹的发生与传播的试验结果［21］Fig.2 Occurrence and propagation of cracks of laminated glass beam obtained by experiment［21］

使用离散元与有限元耦合方法对以上夹层玻璃梁的冲击破坏过程进行仿真.计算模型如图3 所示.其中，夹层玻璃受到4 个支撑体支撑，上支撑体的上表面和下支撑体的下表面受到全约束.冲击块质量为1 kg.初始时，冲击块位于夹层玻璃梁中部正上方，且刚好与梁的上表面接触，并以3.13 m/s 的速度垂直冲击夹层玻璃梁.由于冲击速度较低，模型中所有材料均假设为线弹性材料，各材料物性参数见表1.剪切应力因子为0.1.

图3 夹层玻璃梁模型(单位:mm)Fig.3 Model of laminated glass beam (Unit:mm)

表1 材料物性参数Table 1 Material parameters

夹层玻璃梁计算模型也由3 层构成:上、下层是厚度为10 mm 的玻璃片;中间层是厚度为4 mm 的PVB 膜.上下层玻璃片离散成5 000 个半径为1 mm的离散单元;中间层PVB 划分为1 000 个六面体有限单元.上、下层玻璃片和中间层PVB 膜之间的粘着关系使用离散元与有限元耦合方法处理［9］.冲击块离散成54 个六面体单元，每个支撑体离散成60 个六面体单元，对于冲击块及支撑体与玻璃片之间的接触问题，使用离散元与有限元的接触方法处理［13］.

图4 夹层玻璃梁冲击破坏裂纹的发生与传播的仿真结果Fig.4 Occurence and propagation of cracks of laminated glass beam obtained by simulation

夹层玻璃梁冲击破坏过程的仿真结果如图4 所示.12 μs 时，在冲击侧玻璃片下侧中部首先出现裂纹(见图4(a)).图中裂纹的表述方法是:两个离散单元之间发生破坏时，将这两个离散单元用深黑色表示.这可认为是，冲击压缩波到达PVB 与冲击侧玻璃的结合面时，由于中间膜的密度和刚度过小，冲击压缩波被大部分反射而成为拉伸波，拉伸应力超过玻璃的破坏强度，因而在此发生了破坏［21］.之后，裂纹极快地向上扩展并迅速贯穿上侧玻璃(见图4(b)).由于穿过上侧玻璃的冲击波在PVB 中传播速度相对缓慢.104 μs 时，在下侧玻璃片下表面中部才出现裂纹(见图4(c)).随后，裂纹极快地向上扩展并迅速贯穿下侧玻璃(见图4(d)).

碰撞过程中冲击块与夹层玻璃梁之间的冲击力时间历程仿真结果如图5 所示.冲击开始后，冲击力f迅速增大并达到最大值，随后迅速下降并出现振荡，这一过程伴随着冲击点附近PVB 两侧玻璃的初始裂口发生和传播.在200 μs 后的相当一段时间里，冲击力为0，这是由于夹层玻璃破坏后，向下运动的速度比冲击体大，造成二者的分离.在930 μs 时冲击块与夹层玻璃再次接触，接触力再次出现，图4(e)为该时刻夹层玻璃梁的破坏情况.

图5 冲击块受到的冲击力的仿真结果Fig.5 Impact force of impact body obtained by simulation

对比试验结果图2，可以发现:尽管上、下层玻璃片的初始裂口发生时间与试验有差异，但裂纹产生的地点、上、下层玻璃的破坏顺序及裂口扩展方向均与试验结果基本一致，说明了内聚破坏模型用于汽车玻璃破坏机理研究的有效性.差异产生的原因主要在于仿真模型与试验样片的不一致性，如试验片PVB 两侧玻璃的加工面不可避免地存在毛细裂纹，而仿真模型目前无法考虑其影响.

4 结语

文中根据内聚力模型的基本概念，提出了适用于三维球形离散单元的破坏模型，并将该模型在自主开发的离散元与有限元耦合计算程序CDFP 中实现，仿真分析夹层玻璃梁的冲击破坏过程，并将仿真结果与试验结果进行比较，发现裂纹产生的位置、先后顺序及传播路径均与试验结果一致，说明该模型可以用于汽车夹层玻璃的破坏机理研究.

［1］许骏，刘博涵，葛东云，等.低速冲击下的聚乙烯醇缩丁醛夹层风挡玻璃动态响应研究［J］.兵工学报，2010，31(增刊1):136-139.Xu Jun，Liu Bo-han，Ge Dong-yun，et al.Research on dynamic response of PVB laminated windshield subjected to low velocity impact ［J］.Acta Armamentarii，2010，31(Supp11):136-139.

［2］Xu J，Li Y，Liu B，et al.Experimental study on mechanical behavior of PVB laminated glass under quasistatic and dynamic loadings ［J］.Composite Part B:Engineering，2011，42(2):302-308.

［3］Dharani L R，Ji F S.Dynamic analysis of normal impact of occupant head on laminated glass ［R］.Detroit，Michigan:SAE，1998.

［4］Flocker F W，Dharani L R.Modelling fracture in laminated architectural glass subject to low velocity impact［J］.Journal of Materials Science，1997，32(10):2587-2594.

［5］Xu J，Li Y B，Chen X，et al.Characteristics of windshield cracking upon low-speed impact:numerical simulation based on the extended finite element method［J］.Computation Materials Science，2010，48(3):582-588.

［6］Oda J，Zang M Y，Mori T，et al.Simulation of dynamic fracture behavior of laminated glass by using discrete element method ［C］∥Trans 8th Calculation Dynamics Symp JSME.［S.l.］:［s.n.］，1995:429-430.

［7］Oda J，Zang M Y.Analysis of impact fracture behavior of laminated glass of bi-layer type using discrete element method［J］.Key Engineering Materials，1998(145/146/147/148/149):349-354.

［8］Zang M Y，Lei Z，Wang S F.Investigation of impact fracture behavior of automobile laminated glass by 3D discrete element method［J］.Computational Mechanics，2007，41(1):73-83.

［9］Lei Z，Zang M Y.An approach to combining 3D discrete and finite element methods based on penalty function method［J］.Computational Mechanics，2010，46(4):609-619.

［10］Xu X，Needleman A.Numerical simulations of fast crack growth in brittle solids ［J］.Journal of the Mechanics and Physics of Solids，1994，42(9):1397-1434.

［11］Camacho G T，Ortiz M.Computational modelling of impact damage in brittle materials［J］.International Journal of Solids and Structures，1996，33 (20/21/22):2899-2938.

［12］Zavattieri P，Espinosa H.Grain level model analysis of crack initiation and propagation in brittle materials［J］.Acta Materialia，2001，49(20):4291-4311.

［13］Zang M Y，Gao W，Lei Z.A contact algorithm for 3D discrete and finite element contact problems based on penalty function method ［J］.Computational Mechanics，2011，48(5):541-550.

［14］王端宜，赵熙.沥青混合料单轴压缩试验的离散元仿真［J］.华南理工大学学报:自然科学版，2009，37(7):37-41.Wang Duan-yi，Zhao Xi.Simulation of uniaxial compression test for asphalt mixture with discrete element method［J］.Journal of South China University of Technology:Natural Science Edition，2009，37(7):37-41.

［15］Liu K，Gao L.The application of discrete element method in solving three-dimensional impact dynamics problems［J］.Acta Mechniaca Solida Sinica，2003，16(3):256-261.

［16］Geubelle P H，Baylor J.Impact-induced delamination of laminated composites:a 2D simulation［J］.Composites Part B:Engineering，1998，29(5):589-602.

［17］Pandolfi A，Ortiz M.An efficient adaptive procedure for three-dimensional fragmentation simulations ［J］.Engineering with Computers，2002，18(2):148-159.

［18］Zhang Z，Paulino G H，Celes W.Extrinsic cohesive modelling of dynamic fracture and microbranching instability in brittle materials［J］.International Journal for Numerical Methods in Engineering，2007，72(8):893-923.

［19］Ortiz M，Pandolfi A.Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis［J］.International Journal for Numerical Methods in Engineering，1999，44(9):1267-1282.

［20］Papoulia K D，Sam C H，Vavasis S A.Time continuity in cohesive finite element modeling［J］.International Journal for Numerical Methods in Engineering，2003，58(5):679-701.

［21］臧孟炎，雷周，尾田十八.汽车玻璃的静力学特性和冲击破坏现象［J］.机械工程学报，2009，45(2):268-272.Zang Meng-yan，Lei Zhou，Oda Juhachi.Static characteristic and impact fracture behaviour of automobile glass［J］.Journal of Mechanical Engineering，2009，45(2):268-272.