李 峰,何子建,張 悅,沙芳菲,董新慧,劉 堯

(1.中國礦業(yè)大學(xué)(北京) 深部巖土力學(xué)與地下工程國家重點(diǎn)實(shí)驗(yàn)室，北京 100083; 2.中國礦業(yè)大學(xué)(北京) 應(yīng)急管理與安全工程學(xué)院，北京 100083; 3.中海油研究總院，北京 100028; 4.中國中材國際工程股份有限公司,北京 100102)

煤炭資源是我國的主要能源，煤礦事故災(zāi)害是制約煤礦安全生產(chǎn)的關(guān)鍵因素。大多數(shù)礦井屬于高瓦斯和煤與瓦斯突出礦井，一半以上的礦井具有瓦斯(煤塵)爆炸危險(xiǎn)性。在我國煤礦的災(zāi)害事故中，瓦斯、煤塵爆炸事故造成的傷亡人數(shù)占總傷亡人數(shù)的60%以上。瓦斯(煤塵)爆炸中礦井通風(fēng)構(gòu)筑物損毀嚴(yán)重，致使通風(fēng)系統(tǒng)遭到嚴(yán)重破壞而發(fā)生紊亂，導(dǎo)致災(zāi)情迅速擴(kuò)大波及井下其他區(qū)域，致使此類事故應(yīng)急處置與救援難度相當(dāng)大[1-2]。如果能保證關(guān)鍵的通風(fēng)構(gòu)筑物在事故發(fā)生時(shí)不被破壞，整個(gè)通風(fēng)系統(tǒng)依然可靠，實(shí)現(xiàn)通風(fēng)系統(tǒng)在災(zāi)變條件下的有效性、那么就可以為應(yīng)急處置與救援提供可能條件，就可以盡可能減少瓦斯(煤塵)爆炸事故造成的破壞。目前，國內(nèi)外學(xué)者在瓦斯爆炸沖擊波傳播規(guī)律、衰減規(guī)律、沖擊波破壞作用等方面進(jìn)行了理論與實(shí)驗(yàn)研究[3-7]。周杰[8]采用平面沖擊波方法，計(jì)算了泡沫材料有限元模型前部、內(nèi)部及后部監(jiān)控點(diǎn)的壓力波曲線，了解沖擊波與材料作用的反射、透射的特征。爆炸沖擊波傳播過程中在固定構(gòu)筑物處會(huì)出現(xiàn)多個(gè)波峰震蕩，在此種區(qū)域有著顯著的超壓與振蕩現(xiàn)象，對(duì)礦井通風(fēng)設(shè)施破壞效應(yīng)明顯，但固定構(gòu)筑物(如風(fēng)門、風(fēng)窗、密閉等)在爆炸沖擊波作用過程中具體的破壞原因、過程和機(jī)理尚不清楚。

筆者研究瓦斯爆炸沖擊載荷作用下矩形風(fēng)窗動(dòng)態(tài)響應(yīng)特性與破壞機(jī)理，對(duì)保護(hù)風(fēng)窗在瓦斯(煤塵)爆炸事故中不被破壞，實(shí)現(xiàn)通風(fēng)系統(tǒng)在災(zāi)變時(shí)期的有效性，為應(yīng)急處置與救援提供可能性，減少人員傷亡和損失具有重要的理論指導(dǎo)意義與實(shí)際價(jià)值。

1 矩形風(fēng)窗動(dòng)態(tài)響應(yīng)計(jì)算模型

1.1 矩形風(fēng)窗振動(dòng)控制方程

礦井矩形風(fēng)窗的尺寸大都符合薄板的尺寸要求(厚度與板面寬度的比值在1/5～1/80)，薄板模型可以抗彎、抗扭，也可以承擔(dān)平面內(nèi)的應(yīng)力，符合矩形風(fēng)窗的力學(xué)分析要求，以中性面建立坐標(biāo)系，如圖1所示。以薄板模型研究矩形風(fēng)窗在爆炸沖擊載荷作用下的動(dòng)態(tài)響應(yīng)特征與破壞機(jī)理。

圖1 矩形風(fēng)窗的坐標(biāo)系

在沖擊載荷q(x,y,t)作用下，根據(jù)瞬態(tài)平衡條件，得矩形風(fēng)窗受迫振動(dòng)的微分方程[9]:

(1)

其中，w(x,y,t)為撓度;D為風(fēng)窗的彎曲剛度;ρ為密度;t為時(shí)間。先求解式(1)的齊次方程，當(dāng)q(x,y,t)=0時(shí)，矩形風(fēng)窗自由振動(dòng)微分方程[10]為

(2)

假設(shè)矩形風(fēng)窗的振動(dòng)具有如下隨時(shí)間的諧振子變化規(guī)律[10]:

w(x,y,t)=W(x,y)eiωt

(3)

其中，ω為自振頻率;W(x,y)為矩形風(fēng)窗的撓度振型函數(shù)。把式(3)代入式(2)，變換得

(4)

其中，k4=ρhω2/D。由矩形風(fēng)窗的受力分析可得各物理量之間的關(guān)系為

若用W(x,y)表示各物理量，關(guān)系式為

(11)

(12)

式中，Mx，My為x,y方向的彎矩;Mxy為扭矩；Qx，Qy為x,y方向的橫向剪力;Vx，Vy為x,y方向的橫向總剪力;ν為泊松比。

1.2 矩形風(fēng)窗振動(dòng)Hamilton體系對(duì)偶方程

為了解耦控制方程中的物理量，需引進(jìn)Hamilton體系對(duì)偶方程[11-12]。由式(7)與式(12)可得

(13)

(14)

由(10)得

(15)

綜合上述可得

(16)

聯(lián)合式(5)，(12),(15)得

(17)

令向量v=[W,θ,-Vx,Mx]，可得Hamilton體系對(duì)偶方程:

v′=Hv

(18)

即

1.3 Hamilton體系對(duì)偶方程的求解

按照辛幾何方法[11]，利用分離變量法求解方程(18)，令

v=Y(y)X(x)

(19)

方程(18)有如下形式的解:

V(x)=ξ(x)ψ

(20)

將式(20)代入方程(18)中，得

(21)

因此

(22)

由此可得:μ與-μ均是哈密頓矩陣H的特征值，則方程(18)解的形式為

v=(a1eμx+b1e-μx)Y(y)

由此可得撓度振型函數(shù)W(x,y)為

W=(a1eμx+b1e-μx)W(y)

(23)

將(23)代入(2)得

(24)

方程(24)有如下形式解:W(y)=eλy，代入式(24)可得

(λ2+μ2)2=k4

(25)

因此，可得

λ12=±iα1，λ34=±α2

(26)

同理可得

μ12=±iβ1，μ34=±β2

(27)

由式(26)與(27)可得

(28)

由此可知:沿x,y方向的撓度振型分別為W(x)，W(y)，具體形式[13-15]可寫為

W(x)=a1eiβ1x+b1e-iβ1x+c1eβ2x+d1e-β2x=

A1cosβ1x+B1sinβ1x+C1coshβ2x+D1sinhβ2x

(29)

W(y)=a2eiα1y+b2e-iα1y+c2eα2y+d2e-α2y=

A2cosα1y+B2sinα1y+C2coshα2y+D2sinhα2y

(30)

2 矩形風(fēng)窗振動(dòng)的撓度振型函數(shù)

設(shè)矩形風(fēng)窗的尺寸為a×b，彈性模量為E。基于井下矩形風(fēng)窗的安裝特點(diǎn)，取風(fēng)窗四邊為固支邊界(C-C-C-C)。沿y方向有:W(x,0)=0，?W(x,0)/?y=0;W(x,b)=0，?W(x,b)/?y=0。

由式(30)可得

上式有解，則左側(cè)系數(shù)矩陣行列式為0，得出沿y方向的頻率方程為

(31)

令C2=1，聯(lián)合頻率方程解得

A2=-1，B2=k1α1/α2，C2=1，D2=-k1

因此

coshα2y-k1sinhα2y

同理

coshβ2x-k2sinhβ2x

因此方程(2)矩形風(fēng)窗振動(dòng)的撓度振型函數(shù)為

W(x,y)=C1C2W(x)W(y)=C1×C2×

(32)

3 矩形風(fēng)窗振動(dòng)的動(dòng)態(tài)響應(yīng)特征

3.1 控制方程求解

非齊次方程(1)解的形式可寫成

(33)

將式(33)代入式(1)中可得

q(x,y,t)

由于

因此

兩邊同時(shí)乘以Wm(x,y)，根據(jù)撓度振型函數(shù)的正交性[16]可得

并在整個(gè)平面區(qū)域上積分，得

(34)

因此

根據(jù)杜哈梅積分可得

假設(shè)爆炸沖擊載荷作用在矩形風(fēng)窗上的壓力可看作瞬時(shí)均布載荷，具有如下形式:

q(x,y,t)=q0δ(t-t1)

因此

由此可得沖擊載荷作用下矩形風(fēng)窗振動(dòng)控制方程(1)的解為

(35)

3.2 矩形風(fēng)窗動(dòng)態(tài)響應(yīng)特征

3.2.1矩形風(fēng)窗振動(dòng)參數(shù)

取鋼制矩形風(fēng)窗的參數(shù)為:h=0.06 m，ρ=2 800 kg/m3，E=72 GPa，ν=0.3，a×b=3 m×3.6 m。聯(lián)立式(28),(31)，運(yùn)用牛頓迭代法求解四邊固支條件下(C-C-C-C)矩形風(fēng)窗的主振型函數(shù)，其前10階振動(dòng)參數(shù)計(jì)算值見表1。

表1 四邊固支(C-C-C-C)矩形風(fēng)窗主振型函數(shù)前10階振動(dòng)參數(shù)

3.2.2主振型撓度的動(dòng)態(tài)分布特性

由式(34)與式(35)計(jì)算可得:沖擊載荷下四邊固支(C-C-C-C)矩形風(fēng)窗前10階振動(dòng)參數(shù)值，見表2。由表2可得:1階、5階、6階振動(dòng)函數(shù)為主振型函數(shù)，其振動(dòng)的固有頻率分別為308，975，1 309 Hz;系數(shù)(An/q0)比較可得:5階系數(shù)約為1階系數(shù)的9.4%，6階系數(shù)約為1階系數(shù)的4.9%，因此取前10階主振型完全滿足分析要求。3種主振型的撓度分布，如圖2所示,圖2中Am為常數(shù)系數(shù)，矩形風(fēng)窗撓度的主要分布范圍分別為風(fēng)窗中心區(qū)域、兩條短邊附近與兩條長邊附近區(qū)域。

表2 沖擊載荷下四邊固支(C-C-C-C)矩形風(fēng)窗前10階振動(dòng)參數(shù)值

圖2 主振型撓度動(dòng)態(tài)分布

圖3 主振型應(yīng)力動(dòng)態(tài)分布

圖4 主振型的最大主應(yīng)力與應(yīng)變

3.2.3應(yīng)力與應(yīng)變動(dòng)態(tài)分布特性

(1)應(yīng)力、應(yīng)變動(dòng)態(tài)分布。1階、5階、6階振動(dòng)函數(shù)為主振型，其應(yīng)力動(dòng)態(tài)分布如圖3所示,圖3中An為常數(shù)系數(shù)。1階主振型x方向、y方向與xy扭應(yīng)力集中區(qū)域分別為:長邊中點(diǎn)區(qū)域與風(fēng)窗中心區(qū)域、短邊中點(diǎn)區(qū)域與風(fēng)窗中心區(qū)域、矩形邊角附近區(qū)域;5階主振型x方向、y方向與xy扭應(yīng)力集中區(qū)域分別為:四邊中點(diǎn)及其附近區(qū)域與風(fēng)窗四角附近區(qū)域、短邊中點(diǎn)及其附近區(qū)域、矩形邊角附近區(qū)域與其中一條長邊中點(diǎn)附近區(qū)域;6階主振型x方向、y方向與xy扭應(yīng)力集中區(qū)域分別為:短邊中點(diǎn)及其附近區(qū)域、四邊中點(diǎn)及其附近區(qū)域與風(fēng)窗四角附近區(qū)域、矩形風(fēng)窗邊角附近區(qū)域與一條短邊中點(diǎn)附近區(qū)域;1階、5階、6階主振型的x方向、y方向應(yīng)力在矩形風(fēng)窗中心區(qū)域均存在較明顯的應(yīng)力集中現(xiàn)象。

由此可得:應(yīng)力集中的區(qū)域主要有:① 四邊中點(diǎn)及其附近區(qū)域;② 矩形風(fēng)窗中心區(qū)域;③ 矩形風(fēng)窗四角附近區(qū)域。由于應(yīng)變、彎矩動(dòng)態(tài)分布與應(yīng)力動(dòng)態(tài)分布相似，因此可根據(jù)應(yīng)力集中區(qū)域推斷出矩形風(fēng)窗的主要破壞區(qū)域。

(2)最大主應(yīng)力與應(yīng)變動(dòng)態(tài)分布。1階、5階、6階主振型最大主應(yīng)力與應(yīng)變動(dòng)態(tài)分布如圖4所示。由圖4可知:1階主振型最大主應(yīng)力與應(yīng)變集中區(qū)域分別為:四邊中點(diǎn)區(qū)域與矩形風(fēng)窗的中心區(qū)域、矩形風(fēng)窗中心區(qū)域;5階主振型最大主應(yīng)力與應(yīng)變集中區(qū)域分別為:四邊中點(diǎn)及其附近區(qū)域、長邊中點(diǎn)區(qū)域與短邊中點(diǎn)附近區(qū)域;6階主振型最大主應(yīng)力與應(yīng)變集中區(qū)域分別為:四邊中點(diǎn)及其附近區(qū)域、短邊中點(diǎn)區(qū)域與長邊中點(diǎn)附近區(qū)域;1階、5階、6階主振型的最大應(yīng)力與應(yīng)變?cè)诰匦物L(fēng)窗中心區(qū)域也均存在較明顯的應(yīng)力集中現(xiàn)象。

由此可得:最大主應(yīng)力、最大主應(yīng)變主要集中區(qū)域?yàn)榫匦物L(fēng)窗四邊中點(diǎn)及其附近區(qū)域與矩形風(fēng)窗的中心區(qū)域。最大主應(yīng)力/應(yīng)變?yōu)榈谝恢鲬?yīng)力/應(yīng)變，其集中區(qū)域與應(yīng)力、應(yīng)變、彎矩的集中區(qū)域保持一致，也可由此推斷出矩形風(fēng)窗的主要破壞區(qū)域。

3.2.4橫向總剪力動(dòng)態(tài)分布特性

沿矩形風(fēng)窗厚度方向的剪力與彎矩產(chǎn)生的剪力之和為橫向總剪力，其分布如圖5所示，結(jié)合圖3(a)與圖4(a)可知:橫向總剪力值與主應(yīng)力值相比較小(大約相差兩個(gè)數(shù)量級(jí))，對(duì)矩形風(fēng)窗動(dòng)態(tài)破壞的影響作用相對(duì)較小，分析時(shí)可以忽略。

圖5 1階主振型橫向總剪力

圖6 主振型最大剪應(yīng)力

圖7 四邊固支矩形風(fēng)窗(C-C-C-C)動(dòng)態(tài)破壞過程

3.2.5矩形風(fēng)窗動(dòng)態(tài)破壞過程

基于第三強(qiáng)度準(zhǔn)則，矩形風(fēng)窗中最大剪應(yīng)力動(dòng)態(tài)分布如圖6所示，由圖6可知:1階主振型的最大剪應(yīng)力主要分布在四邊中點(diǎn)區(qū)域與矩形中心區(qū)域;5階主振型的最大剪應(yīng)力主要分布在兩條短邊中點(diǎn)區(qū)域與沿矩形長軸中心區(qū)域;6階主振型最大剪應(yīng)力主要分布在兩條長邊中點(diǎn)區(qū)域與沿矩形短軸中心區(qū)域。由此分析可得破壞過程如下:沖擊載荷作用下四邊固支(C-C-C-C)矩形風(fēng)窗發(fā)生損傷破壞的位置位于四邊中點(diǎn)區(qū)域與矩形風(fēng)窗中心區(qū)域，主要沿著風(fēng)窗長軸的中心區(qū)域破壞，其次是沿著短軸的中心區(qū)域破壞。

4 矩形風(fēng)窗動(dòng)態(tài)破壞過程模擬

基于LS-DYNA軟件，采用PLASTIC_ KINEMATIC材料模型[17]，模擬分析矩形風(fēng)窗在沖擊載荷作用下的動(dòng)態(tài)破壞過程，如圖7所示，四邊固支矩形風(fēng)窗(C-C-C-C)動(dòng)態(tài)破壞過程如下:首先風(fēng)窗短邊中心出現(xiàn)裂紋發(fā)生損傷破壞、風(fēng)窗長邊中心發(fā)生損傷破壞，然后是風(fēng)窗中心區(qū)域發(fā)生破壞，破壞方向主要沿著長軸中心方向發(fā)展，其次會(huì)沿著短軸中心方向發(fā)展。最終主要的破壞區(qū)域?yàn)轱L(fēng)窗的四周、中心區(qū)域、沿長軸與短軸中心區(qū)域;數(shù)值模擬分析的結(jié)果與理論分析結(jié)果完全吻合。

5 矩形風(fēng)窗動(dòng)態(tài)破壞機(jī)理探討

沖擊載荷作用下四邊固支(C-C-C-C)矩形風(fēng)窗振動(dòng)以1階為主(主頻)，5階與6階次之。應(yīng)力、應(yīng)變、彎矩動(dòng)態(tài)分布相似，由圖2～4可知:3個(gè)主振型的主要集中區(qū)域有:四邊中點(diǎn)及其附近區(qū)域、矩形風(fēng)窗中心區(qū)域、矩形風(fēng)窗四角附近區(qū)域。由圖8可得:最大主應(yīng)力與主應(yīng)變主要集中區(qū)域?yàn)榫匦物L(fēng)窗四邊中點(diǎn)及其附近區(qū)域、矩形風(fēng)窗的中心區(qū)域。由圖6(a)可得:1階主振型四邊中點(diǎn)區(qū)域最大剪應(yīng)力明顯大于中心區(qū)域的最大剪應(yīng)力，因此首先發(fā)生破壞的位置應(yīng)是四邊中點(diǎn)區(qū)域，其次是矩形風(fēng)窗的中心區(qū)域;由圖6(b)可知，除了四邊中點(diǎn)區(qū)域的最大剪應(yīng)力集中以外，5階主振型沿矩形長軸中心區(qū)域也具有最大剪應(yīng)力集中現(xiàn)象，因此，破壞的位置應(yīng)是沿著長軸中心區(qū)域擴(kuò)展;同時(shí)由圖6(c)可知:6階主振型沿矩形短軸中心區(qū)域也具有最大剪應(yīng)力集中現(xiàn)象，因此破壞的位置沿著短軸中心區(qū)域也會(huì)有相對(duì)小范圍的擴(kuò)展。

結(jié)合圖7綜合分析可知:由1階主振型(主頻)動(dòng)態(tài)分布特征可以推斷出矩形風(fēng)窗的起始破壞位置為矩形四邊中心區(qū)域，由5階、6階主振型可以判斷破壞的發(fā)展方向與趨勢為沿著長軸或短軸中心區(qū)域擴(kuò)展，3個(gè)主振型綜合構(gòu)成了沖擊載荷下四邊固支(C-C-C-C)矩形風(fēng)窗振動(dòng)的主振型函數(shù)。

6 結(jié) 論

(1)基于薄板模型，建立了矩形風(fēng)窗自由振動(dòng)的Hamilton體系對(duì)偶方程;基于振型函數(shù)的正交性與杜哈梅積分，得出了沖擊載荷作用下矩形風(fēng)窗振動(dòng)的撓度動(dòng)態(tài)分布方程。

(2)基于牛頓迭代法，得出了四邊固支條件下(C-C-C-C)矩形風(fēng)窗的主振型函數(shù);得出了沖擊載荷作用下四邊固支矩形風(fēng)窗振動(dòng)主要由1階、5階、6階主振型構(gòu)成，其振動(dòng)的固有頻率分別為308，975，1 309 Hz。

(3)基于第三強(qiáng)度理論，得出了沖擊載荷作用下四邊固支矩形風(fēng)窗發(fā)生損傷破壞的起始位置是四邊中點(diǎn)區(qū)域，然后是矩形風(fēng)窗中心區(qū)域;發(fā)展趨勢是主要沿著矩形風(fēng)窗長軸的中心區(qū)域破壞，其次是沿著短軸的中心區(qū)域破壞。