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

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

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

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

1 矩形風窗動態(tài)響應計算模型

1.1 矩形風窗振動控制方程

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

圖1 矩形風窗的坐標系

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

(1)

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

(2)

假設矩形風窗的振動具有如下隨時間的諧振子變化規(guī)律[10]:

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

(3)

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

(4)

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

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

(11)

(12)

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

1.2 矩形風窗振動Hamilton體系對偶方程

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

(13)

(14)

由(10)得

(15)

綜合上述可得

(16)

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

(17)

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

v′=Hv

(18)

即

1.3 Hamilton體系對偶方程的求解

按照辛幾何方法[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 矩形風窗振動的撓度振型函數(shù)

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

由式(30)可得

上式有解，則左側系數(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)矩形風窗振動的撓度振型函數(shù)為

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

(32)

3 矩形風窗振動的動態(tài)響應特征

3.1 控制方程求解

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

(33)

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

q(x,y,t)

由于

因此

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

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

(34)

因此

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

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

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

因此

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

(35)

3.2 矩形風窗動態(tài)響應特征

3.2.1矩形風窗振動參數(shù)

取鋼制矩形風窗的參數(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)，運用牛頓迭代法求解四邊固支條件下(C-C-C-C)矩形風窗的主振型函數(shù)，其前10階振動參數(shù)計算值見表1。

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

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

由式(34)與式(35)計算可得:沖擊載荷下四邊固支(C-C-C-C)矩形風窗前10階振動參數(shù)值，見表2。由表2可得:1階、5階、6階振動函數(shù)為主振型函數(shù)，其振動的固有頻率分別為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ù)，矩形風窗撓度的主要分布范圍分別為風窗中心區(qū)域、兩條短邊附近與兩條長邊附近區(qū)域。

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

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

圖3 主振型應力動態(tài)分布

圖4 主振型的最大主應力與應變

3.2.3應力與應變動態(tài)分布特性

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

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

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

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

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

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

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

圖6 主振型最大剪應力

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

3.2.5矩形風窗動態(tài)破壞過程

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

4 矩形風窗動態(tài)破壞過程模擬

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

5 矩形風窗動態(tài)破壞機理探討

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

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

6 結 論

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

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

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