动荷载作用下预测玻璃板破坏强度和破坏时间的一种改进模型——Ⅰ.理论

柳锦春，于润清，唐德利

(1.解放军理工大学 国防工程学院，南京 210007;2.爆炸冲击防灾减灾国家重点实验室，南京 210007)

随着城市建设的高速发展，以及建筑设计理念的不断提升，玻璃门窗和玻璃幕墙等成为现代建筑重要的组成部分。但玻璃是脆性材料，在爆炸等强动载作用下易于破碎，且产生高速飞散的碎片，易对人和物造成伤害[1-2]。准确确定玻璃板的破坏应力，对玻璃门窗和玻璃幕墙的抗爆设计有着重要意义。目前有很多计算玻璃抗爆能力的方法，影响最广的是Beason的理论[3]，他基于Weibull分布进行复杂的理论推导，将任意形式荷载等效为恒定荷载，由此得出玻璃的最大抗力。由他的理论发展出美国和加拿大的玻璃抗爆设计规程。

但玻璃板的破坏是由板中的微小裂纹引起的[4]，裂纹的分布和裂纹的扩展都会影响玻璃板的破坏应力，在爆炸冲击荷载作用下更是如此，而且裂纹的分布和扩展又有着很大的随机性。本文首先对常用玻璃破坏模型进行了分析和讨论，指出现有玻璃破坏模型的适用性和不足。在借鉴已有模型[5]的基础上，通过蒙特卡洛法确定玻璃板表面的随机裂纹，运用有限元软件计算玻璃板的应力分布，再综合考虑裂纹的动态扩展和玻璃表面的预应力，确定玻璃的破坏准则，最终确定出玻璃板的破坏强度和破坏时间。基于此思路，本文提出了一种可以预测动荷载作用下玻璃板破坏应力和破坏时间的改进模型。

1 常用玻璃破坏模型介绍及讨论

玻璃的破坏应力有很大的离散性，一般采用Weibull分布描述破坏应力的分布[6-7]，形式如下：

(1)

其中:Pf是失效概率，σ0是破坏应力的均值，σm是最小破坏应力。A、a分别是有效体积(面积)和参考体积(面积)，m是参数。式(1)被称为Weibull三参数分布。

当σm=0时，式(1)简化为式(2)，式(2)被称为Weibull双参数分布。由于Weibull双参数分布也可以较好的描述其分布且形式更加简单，因此，很多预测玻璃破坏应力的模型[3,8-12]的形式是基于Weibull双参数分布的。

(2)

但是，这些模型存在以下不足点，一是材料参数较难准确确定，这是因为影响这两个参数的因素较多，且参数值容易受到试验误差的影响；二是Weibull双参数分布的两个参数表达的物理意义并不十分明确[5]；三是Weibull分布只有当裂纹的数量服从指数分布时才成立[13-14]；四是由Weibull分布发展的模型一般是用于描述玻璃板在准静态情况下破坏应力的分布，并不适用于分析玻璃板承受如爆炸荷载等动荷载作用的情况。

由于玻璃板的破坏是由裂纹引起的，因此可以由裂纹的开裂来计算破坏应力。Nurhuda等[5]根据此思想提出了一种计算玻璃板破坏应力的模型。他将玻璃板中裂纹简化为二维Griffith裂纹，不考虑裂纹的交叉，通过蒙特卡洛法模拟玻璃板裂纹的几何随机性，结合有限元分析软件，对玻璃板的破坏应力做出预测。玻璃板的破坏是由某一关键裂纹控制的，因此，该模型主要是确定出关键裂纹的位置和长度，由此估算出破坏应力。具体思路如下[5]：

(1) 随机确定裂纹的数目、长度、位置。

裂纹的数目服从Poisson分布

f(n)=[(ρA)n/n!]exp(-ρA)

(3)

式中：ρ是裂纹密度，n是裂纹数目。

设定裂纹的最大长度为278 μm，裂纹长度假定服从特定函数分布，即均匀分布，正态分布，Weibull分布(左偏分布)，对数正态分布(右偏分布)中的一种。

裂纹的位置是均匀随机分布。

(2) 建立关键裂纹确定标准，找出关键裂纹。

Griffith裂纹长度为aj，与破坏应力σj有如下关系式：

(4)

其中:KI是应力强度因子，Y是裂纹的形状参数。

Nurhuda取KI=KIC为裂纹破坏临界条件，即

(5)

其中:λj是任意处应力与板中最大应力比，σmax是板中最大应力。

当最大裂纹处于板中最大主应力处，玻璃板有最小破坏应力σmin，即

(6)

由式(5)和式(6)，可得

(7a)

(7b)

其中:rj是任意处裂纹长度与裂纹最大长度之比。

在求出每条裂纹的Sj后，取Sj值最小的裂纹为关键裂纹。在确定关键裂纹后，回代入式(4)可以求出破坏应力。

Nurhuda 等提出的模型虽然对裂纹的随机性做了充分的研究，但是没有考虑裂纹的动态扩展，此外，Nurhuda假设玻璃板最大应力在板中央点，而在动荷载作用下，板中央的应力并不一定是最大的，所以Nurhuda的模型不能很好地分析玻璃板承受动荷载的情况。

2 改进的玻璃破坏模型

基于上述分析，我们对Nurhuda提出的模型进行改进。在这里，将玻璃板视为二维板，不考虑裂纹的交叉及相互耦合，主要考虑以下四点来模拟玻璃的破坏。

(1) 考虑玻璃板在荷载作用下的几何非线性。文献表明[15]，薄板在承受较大均布荷载时表现出明显的非线性，若用线性薄板理论计算可能会高估薄板的挠度和应力值。

(2) 考虑表面预压应力对玻璃板破坏的影响。

(3) 考虑裂纹的动态扩展。

(4) 考虑裂纹的几何随机性(用蒙特卡洛法模拟)。

该改进模型采用ABAQUS软件计算玻璃板的应力。在计算时不考虑裂纹的影响，仍将玻璃板设为均值材料，提取出玻璃板每一点的应力时程曲线组成应力矩阵，作为原始分析数据；在玻璃板的各点运用蒙特卡洛法随机“虚加”裂纹，根据破坏准则对应力矩阵进行分析，找出最先破坏的点，并将此点此刻的应力作为破坏应力，此刻作为破坏时间，模型结构示意图如图1所示。

图1 模型结构示意图

模型的裂纹随机性通过编制随机程序产生，下面主要介绍Griffith的动态扩展、玻璃表面预应力和破坏准则的确定。

2.1 Griffith裂纹的动态扩展

图2 裂纹开裂速度与应力强度因子的关系

裂纹的扩展是有一定的条件的[16]，当KIKIC时，裂纹迅速扩展，此时可认为裂纹破坏。

在一定环境下，开裂速度可以表示为

(8)

其中:n，v0是与环境和材料有关的参数。

裂纹在扩展时可分为三个阶段，如图2，随着KI的增大，n的值是变化的，但是，由于在裂纹扩展的第二阶段和第三阶段的开裂速度较大，所以用第一阶段的n值代替第二三阶段的n值对估算的最终破坏时间影响并不大。

2.2 玻璃表面的预压应力

常用的单层玻璃有浮法玻璃，半钢化玻璃和钢化玻璃，不同性质的玻璃由于加工工艺的不同导致在玻璃的横截面上存在预压力[17]，如图3。一般而言，玻璃表面有预压应力，内部则是预拉应力。不同性质的玻璃板表面的预压力不同[18]，见表1。

表1 不同种类玻璃表面预压力值[18]

图3 钢化玻璃截面抛物线应力分布[17]

当玻璃板在受力时，表面的拉应力首先要抵消预压应力，由于玻璃板内部的预压应力小于表面的预压应力，因此，一旦表面被破坏，玻璃板就会迅速破坏。所以，可以只考虑表面预压应力的影响，忽略玻璃板内部预压应力的大小。对于不同的性质的玻璃板，本文取预压应力的最小值。

2.3 模型的破坏准则

图 4 裂纹的方向与裂纹处的应力

本文模型的破坏准则是指在各点应力时程曲线组成应力矩阵中确定出最先破坏的裂纹所遵循的规则。

如图4，将玻璃视为线弹性材料，设定好弹性模量和泊松比，用ABAQUS有限元软件计算σ11、σ12、τ12的瞬间大小，再用蒙特卡洛法随机模拟出裂纹的方向θ，由此就可以计算出裂纹处法向应力大小。

(9)

若忽略剪应力的影响，则得

σθ=σ11cos2θ+σ22sin2θ

(10)

式(10)即为Weibull和Overend采用的法向应力值[7,17]。

由于玻璃表面存在预压应力，因此还要考虑预压应力f的影响，式(9)可写为式(11):

(11)

裂纹的法向应力σθ与裂纹的长度a是耦合的，在计算时要先进行解耦，本文假设裂纹在很小的时间增量步内应力不变，由式(4)、式(8)可得裂纹长度与应力之间的迭代关系式：

(12)

其中:aj是第j时刻的裂纹长度，σj-1是j-1时刻裂纹法向拉应力，Δt为时间增量步。

式(12)在迭代时，根据时间输入相应的应力值(σ11、σ12、τ12)，此应力值是由ABAQUS预先计算的，在计算应力时并不考虑裂纹，因此可以说，裂纹是“虚加”的。通过式(12)的迭代计算，既可以得出各点应力强度因子的时程曲线，认定KI=KIC时为破坏临界点，此时的荷载作用时间即为该点的破坏时间，取破坏时间最小的点为最先破坏点，最先破坏点对应的破坏时间和破坏应力视为玻璃板的破坏时间和破坏应力。

2.4 改进模型的整体计算思路

(1) 运用蒙特卡洛法在玻璃板的上下表面“布置”假定的随机裂纹。裂纹的数目、长度、位置的模拟参照Nurhuda模型，裂纹方向假设为随机均匀分布。

(2) 根据式(12)进行裂纹长度的迭代，根据破坏准则判断是否开裂。

(3) 比较所有裂纹的破坏时间，取破坏时间最短的裂纹为关键裂纹，将此时刻作为玻璃板的破坏时间，并输出此时刻关键裂纹处的应力作为玻璃板的破坏应力。

(4) 玻璃板的破坏应力具有概率性，为得到不同破坏应力对应的概率大小，重复(1)～(3)步，计算N块玻璃板的破坏时间和破坏应力，由此可求得累计概率(CPD)Pf：

(13)

按上述方法进行模拟需要确定两个未知参数和一个未知分布，即裂纹分布密度、最大裂纹长度和裂纹长度的分布。在计算时可以先假设裂纹长度分布，结合试验确定在该长度分布下最优的分布密度和最大长度，然后对比不同长度分布模拟的结果，取与试验结果最相近的裂纹长度分布和相应的裂纹分布密度和最大长度。

该改进模型可以同时计算出玻璃板的破坏时间和破坏应力。此外，该模型在计算玻璃板应力时，采用的是有限元方法，可以充分考虑玻璃板的受力特性，且不受玻璃板形状和边界条件的限制。

2.5 改进模型的验证

由于同时给出破坏应力和破坏时间的实验数据较少，下面分别对采用本文模型计算的破坏时间和破坏应力进行试验验证对比。

Gavanski等[19]研究了浮法玻璃板的破坏时间，通过试验统计了破坏时间的累计概率。本节选择合适的裂纹长度分布函数、裂纹分布密度和最大的裂纹长度，对玻璃板的破坏时间进行模拟，未知参量的确定见文献[20]。

取玻璃的弹性模量E=73 GPa，泊松比v=0.25，密度ρ=2 450 kg/m3，玻璃板四边简支，尺寸为1 m×1 m×0.006 m，加载速率为6.5 kPa/s。取长度分布函数为对数正态分布，均值0.272 5，方差1.616，裂纹最大长度230 μm，裂纹密度为6 m-2，共进行取15组模拟，每组20个，取模拟结果的平均值与试验值对比，如图5所示。由图5可知，模拟结果与实验结果有较好的一致性。

图5 模拟的结果与试验结果的对比

取Nurhuda的试验[5]计算玻璃板破坏应力的累计概率。加载速率取为10 kPa/s，玻璃板尺寸为2 m×0.4 m×0.008 m，玻璃板四边简支，改进模型的参数确定及取值同上，计算结果如图6。由图6可知，模拟结果与实验结果基本吻合。

图6 模拟的结果与试验结果的对比

3 结 论

本文在借鉴已有模型[5]的基础上，提出了一种能够同时预测破坏时间和破坏应力的改进模型。改进模型通过蒙特卡洛法确定玻璃板表面的随机裂纹，运用有限元软件计算玻璃板的应力分布，通过随机裂纹的破坏计算出破坏应力和破坏时间。在运用有限元软件计算应力时，考虑了玻璃板的几何非线性；在计算裂纹的动态扩展时，考虑了Griffith裂纹的动态扩展，推导出裂纹开裂长度与破坏应力之间的迭代关系式；在计算裂纹开裂时，考虑了裂纹开裂的条件(KIKIC时,裂纹迅速扩展)和玻璃板的表面预应力的影响。此改进模型有三个需要确定的未知参量，即裂纹分布密度、最大裂纹长度和裂纹长度的分布，这些参量的确定见论文的第Ⅱ篇[20]。最后对此改进模型进行了验证分析，结果表明该改进模型可以较准确地同时计算预测玻璃板的破坏应力和破坏时间。由于此改进模型是运用有限元软件进行玻璃板的应力分析，因此其不受玻璃板形状和边界条件的限制，使用范围较广。

[1] Norville H S，Harvill N．Glass-related injuries in oklahoma city bombing[J]．ASCE：Journal of Performance of Constructed Facilities，1999，13(2)：50-56.

[2] Norville H S，Conrath E J．Considerations for blast resistant glazing design[J]．ASCE：Journal of Architectural Engineering．2001，7(3)：80-86.

[3] Beason. A failure prediction model for window glass[R], Institute for Disaster Research, Texas Tech University, USA, 1980.

[4] Brown W G. A load duration theory for glass design [R]. NRCC 12354, Canada, National Research Council of Canada, Division of Building Research, 1972.

[5] Nurhuda I, Lam N T K, Gad E F, et al. Estimation of strength in large annealed glass panels [J]. Int. J. Solids Struct. , 2010, 47(18-19): 2591-2599.

[6] Behr R A, Karson, Minor. Reliability analysis of window glass failure pressure data [J]. Struct. Safety, 1991, 11: 43-58.

[7] Weibull. A statistical distribution functions of wide applicability [J]. Appl. Mech., 1951, 18: 293-297.

[8] Beason W L, Morgan. Glass failure prediction model [J]. J. Struct. Eng., 1984, 110(2): 197-212.

[9] Sedlacek G, Blank K, Gusgen J. Glass in structural engineering[J]. J. Struct. Eng., 1995, 73(2): 17-22.

[10] Dalgliesh W A, Taylor. The strength and testing of window glass [J]. Can. J. Civ. Eng., 1990, 17: 752-762.

[11] Charles R J. Static fatigue of glass 1&2[J], J. Appl. Phys, 1958, 29(11): 1549 1560.

[12] Fischer-Cripps A C, Collins. Architectural glazing: Design standards and failure models [J]. Build. Environ, 1995, 30(1): 29-40.

[13] Todinov M T. Is Weibull distribution the correct model for predicting probability of failure initiated by non-interacting flaws[J]. Int. J. Solids Struct., 2009, 46:887-901.

[14] Todinov M T. Selecting designs with high resistance to overstress failure initiated by flaws [J]. Computational Materials Science, 2008, 42: 306-315.

[15] Mehmet Z A. Laminated glass plates: revealing of nonlinear behavior[J]. Computers and Structures, 2003, 81:2659-2671.

[16] Evans A G , Wiederhorn S M. Proof testing of ceramic materials-an analytical basis for failure prediction [J]. Int. J. Fract, 1974, 10(3): 379-392.

[17] Overend M, Parke G A R, Buhagiar D. Predicting failure in glass-a general crack growth model [J]. Journal of Structural Engineering, 2007, 8:1146-1155.

[18] Wurm J. Glass structures: design and construction of self-supporting skins [M]. Berlin: Springer, 2007.

[19] Gavanski E, Gregory A K, ASCE M. Glass breakage tests under fluctuating wind loads[J]. American Society of Civil Engineers, 2011, 3:34-41.

[20] 柳锦春, 于润清, 唐德利. 动荷载作用下预测玻璃板破坏强度和破坏时间的一种改进模型——Ⅱ.参数确定及验证[J]. 振动与冲击, 2014,33(13):82-85,116.

LIU Jin-chun, YU Run-qing, TANG De-li. Improved prediction model for failure strength and failure time of glass panels under dynamic loading —— Ⅱ. Parameter determination and verification [J]. Journal of Vibration and Shock, 2014,33(13):82-85,116.