動荷載作用下預測玻璃板破壞強度和破壞時間的一種改進模型
——Ⅰ.理論

柳錦春，于潤清，唐德利

(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.