施 峰,汪華君,舒 才,王宏圖
(1.貴州工程應用技術學院 礦業工程學院,貴州 畢節 551700; 2.重慶大學 煤礦災害動力學與控制國家重點實驗室,重慶 400044)
煤與瓦斯突出嚴重威脅著我國煤礦的安全生產,而保護層開采作為迄今為止最經濟有效的防治煤與瓦斯突出的措施之一,至今仍是防治煤與瓦斯突出領域的研究熱點。國內外學者主要對保護層開采消突機理與消突效果展開研究[1-3],或在假定保護層開采其他條件不變的基礎上,以被保護層應變或應力改變量為考察指標,對保護效果受某一特定條件的影響進行單因素分析[4-6],而基于可綜合描述被保護層應力應變改變的被保護層應變曲線,進行保護效果的多保護效果指標、多影響因素的定量研究相對較少。
由于基于被保護層應變曲線形態的保護效果指標較多,并且被保護層的卸壓變形過程是受多因素影響的非線性過程[7],保護效果受多個保護層開采條件變化影響的定量研究難以開展。因此,將建立顯式的響應面模型(RSM)[8-10],代理保護效果與保護層開采條件間的不易顯式化的非線性關系;對上保護層開采進行變參數有限元迭代求解,確定響應面模型參數;對響應面模型求偏導,得到保護效果指標受參數影響的敏感度[11];通過對敏感度的量化研究,為保護層開采的安全設計提供理論依據。
南桐礦區內礦井主采煤層為K4,K6,目前開采水平深度約為650 m,平均傾角為45°。K6煤層位于K4煤層上方,層間距平均為35 m;K4煤層、K6煤層厚度平均分別為2.5和2.0 m。2層煤均具有突出危險性,礦井多選擇突出危險性較低的K6煤層進行上保護層開采,消除下部K4煤層突出危險性。
采用Ansys有限元軟件進行數值計算,為響應面模型的建立提供基本輸入參數。有限元計算采用圖1所示簡化的平面應變模型,單元類型為PLANE183,強度準則為摩爾-庫侖準則[3]。按試驗設計,進行多次不同幾何條件和邊界條件的數值計算。

圖1 上保護層開采保護效果敏感度分析的幾何模型Fig.1 Geometry model for FEA
以保護層平均開采條件的幾何參數確定變參數有限元模型的基本幾何模型,如圖1所示。有限元計算的幾何模型總體尺寸為600 m×400 m(長×高)。
基本幾何模型中,煤層與模型底邊夾角為45°;煤層間距35 m;保護層K6煤層開挖區域傾向長度固定為113 m,布置于K6煤層中央以減小邊界效應影響;在K6煤層工作面開挖位置設生死單元,并通過生死單元實現工作面開挖,通過在頂底板布置接觸單元實現頂底板接觸。
有限元模型整體施加自重應力,下部及左側為對稱邊界條件,上邊界施加垂直應力荷載補償研究區域的礦井實際垂直地應力,右側根據垂直地應力及側壓系數施加水平方向的斜坡應力荷載模擬水平地應力。
有限元模型采用的材料物理力學參數主要基于現場鉆孔得到的煤層柱狀圖和實測的巖石物理力學參數,如表1所示。

表1 數值模型各煤巖層的物理力學參數Table 1 Physical mechanic parameters for all stratum
保護層與被保護層間存在6層厚度及力學性質各不相同的巖層,具有組合介質結構的力學特性。在有限元模型中通過對材料參數進行厚度加權平均,簡化為單一等效復合巖層[3,13]。
保護層開采保護效果的影響因素較多,因此根據礦區地質條件變化的實際情況,選擇地應力參數:垂直地應力、側壓系數以及幾何參數:保護層厚度、復合巖層等效保護層間距、煤層傾角,共5個保護層開采參數作為試驗設計的輸入參數[4-6,12-14]。以上述5個保護層開采參數的平均值為中心,分別向正負方向展開成表2所示的5維輸入參數研究空間。
針對敏感度研究中保護層開采參數在平均值附近局部多參數變化的特點,響應面模型試驗設計方案采用可旋轉中心組合設計方案(CCD)[9]。該試驗方案試驗點空間分布參數α為2.378 4,每個因素含5個水平,總計43次試驗。其中中心點試驗1次、軸向試驗10次、對角線試驗32次。

表2 響應面模型輸入參數取值范圍Table 2 Input parameters of response surface model
將保護效果考察指標分為范圍相關、位置相關、卸壓程度相關3類。
范圍相關考察指標,用于表征膨脹變形曲線與保護準則的關系,包括:1)卸壓范圍,即被保護層垂直于層面的膨脹變形量大于0的范圍;2)充分卸壓保護范圍長度,參考《防突規定》,定義為被保護層垂直于層面的膨脹變形量大于3‰的區域傾向長度;3)上部卸壓角,即充分卸壓保護范圍上界與保護層開挖工作面回風巷的連線與保護層的夾角;4)下部卸壓角,即充分卸壓保護范圍下界與保護層開挖工作面機巷的連線與保護層的夾角。
卸壓程度相關考察指標,用于表征膨脹變形曲線的高度,包括:1)保護范圍平均膨脹變形量,即充分卸壓保護范圍內垂直層面的膨脹變形量的積分中值;2)卸壓范圍平均膨脹變形量,即卸壓范圍內垂直層面的膨脹變形量的積分中值;3)被保護層膨脹變形極值,即膨脹變形曲線上的最大值。
位置相關考察指標,用于表征膨脹變形曲線沿層面方向的形態,包括:1)保護范圍中心位置,即充分卸壓保護范圍中心到保護層開挖區域中心線的距離;2)卸壓范圍中心位置,即卸壓范圍中心到保護層開挖區域中心線的層面距離;3)膨脹變形極值位置,即被保護層膨脹應變極值位置到保護層開挖區域中線的距離。
對圖1中K4煤層膨脹變形進行抽樣,共得到42個如圖2所示被保護層沿層面的膨脹變形曲線。通過“凸”形頂部位置、底部截距的中心位置及曲線兩側的保護范圍上下邊界可確定膨脹變形的整體形態。

圖2 保護層開采平均開采條件下被保護層膨脹變形曲線Fig.2 Strain plot of protected seam of average condition
響應面功能函數是進行保護效果指標受保護層開采影響參數敏感度分析的基礎。本文研究的敏感度問題參數空間較小,敏感度指標變化程度不大、規律性較強,因此選擇經典多項式響應面功能函數[9-10]。
各個敏感度考察指標選擇不同的響應面多項式功能函數,從基本的線性多項式開始逐漸增加影響因素的相互作用項。最終確定以下考察指標采用完全二次多項式功能函數擬合:下部卸壓角、上部卸壓角、保護范圍長度、保護范圍中心位置、卸壓范圍長度、卸壓范圍中心位置、卸壓范圍內極值、極值位置;以下考察參數采用包含因素項及因素交互關系項的多項式功能函數擬合:保護范圍內平均卸壓程度、卸壓范圍內平均卸壓程度。二次多項式功能函數的一般形式如式(1)所示,其中因素平方項系數取0即退化為僅包含因素項及因素交互關系項的多項式功能函數。
在構造響應面功能函數時,為方便影響因素的對比分析及減小求解過程中的截斷誤差,將表2參數變化區間線性等比例映射到區間[-1,+1]。
映射處理的參數作如下規定:記A為層間距,m;B為煤層傾角,(°);C為保護層厚度,m;D為上邊界應力,MPa;E為側壓系數;因素兩兩相乘表示兩因素交互作用;coeff為因素對應的系數;intercept為常數項;則參數線性等比例映射后的二次多項式響應面功能函數可表示為:
coeff(D2)*D2+coeff(E2)*E2
(1)
表3為式(1)對應的5因素10考察指標的多項式響應面功能函數擬合結果。表中每1列為1個考察參數的功能函數,每一行為因素對應的系數。
表3中位置參數>0說明位置偏上山方向,<0說明位置偏下山方向。
上保護層開采條件在平均開采條件附近變化,因此這里研究平均開采條件,即A,B,C,D,E分別取0值時,以響應面功能函數偏導數表示的敏感度:
(2)
式中:s為敏感度;Y為保護效果指標的響應面功能函數;X0為CCD實驗設計的中心點。
從式(2)可知,試驗方案中心點處的偏導數分別為表3中前5行。表3中卸壓程度相關指標的敏感度數量級均較小,小于10-3,因此不作分析。

表3 響應面功能函數Table 3 Performance functions of response surface
1)范圍相關保護效果考察指標的敏感度分析
圖3為范圍相關保護效果考察指標分析,圖3(a)下部卸壓角受層間距影響最敏感,敏感度為1.02;圖3(b)上部卸壓角受側壓系數影響最敏感,敏感度為1.70;圖3(c)中保護范圍長度受層間距影響敏感度最高,為-1.89;上下部卸壓角受層間距影響敏感度均超過1,而傾向保護范圍長度對層間距敏感度為負,說明對于該保護層開采工程實例,隨層間距增加,上下部卸壓角均增大而傾向保護范圍長度減小,表現出卸壓角與傾向保護范圍長度間的非線性關系。

圖3 范圍相關指標受保護層開采條件影響的敏感度Fig.3 Sensitivity of three protection index vs mining factors
圖4是不同開采條件敏感度的分組對比圖。組內對比顯示,保護效果指標受開采參數影響的敏感度存在類似的特征:傾向保護范圍長度>上部卸壓角>下部卸壓角;各組間柱狀圖整體高度對比顯示,保護效果指標受開采參數影響的敏感度總體為側壓系數>層間距>垂直應力>保護層厚度>煤層傾角。

圖4 上下部卸壓角及傾向保護范圍長度敏感度分析Fig.4 Sensitivity column chart of boundary relief angles and length of protected range
從圖5可知,卸壓范圍長度受層間距影響最敏感,敏感度為2.99,對其他4個開采參數敏感度均較小;傾向保護范圍長度與卸壓范圍長度受層間距影響的敏感度為異號,說明膨脹變形量準則從零逐漸提高到3‰,傾向保護范圍長度受層間距的敏感度必存在符號改變的趨勢轉折點。

圖5 傾向保護范圍長度、卸壓范圍長度敏感度分析Fig.5 Sensitivity column chart of length of protected range & swelling range of the protected seam
2)位置相關保護效果考察指標的敏感度分析
圖6中位置相關效果指標受5個保護層開采參數影響的敏感度均小于1,其中側壓系數影響的敏感度最高(0.57,0.66,0.7),受層間距影響的敏感度次之(0.34,0.46,0.57)。位置相關效果指標受開采參數變異影響的敏感度總體上小于范圍相關保護效果指標的對應值,說明膨脹變形曲線特征位置表示的被保護層卸壓分布的敏感度不及卸壓角和保護范圍長度。

圖6 保護層開采效果位置相關指標的敏感度分析Fig.6 Sensitivity analysis column chart of protection effect indicators related to position
1)響應面功能函數回歸效果檢驗
表4顯示各響應面功能函數的回歸效果。表4中響應面功能函數的擬合系數及修正擬合系數均達到0.91以上,說明所選10項保護效果指標的多項式響應面功能函數回歸效果均較好。
2)基于物理相似模擬的敏感度分析方法的驗證
為驗證利用基于數值模擬的保護效果響應面模型進行敏感度分析的正確性,進行工程條件類似的不同層間距的物理相似模擬[15],得到保護范圍長度受層間距變異影響的敏感度為-0.555,與本文敏感度結果接近。

表4 響應面功能函數回歸效果Table 4 Regression effect of RSM functions
1)保護開采參數敏感度分析顯示:范圍相關保護效果指標受保護層開采參數影響的敏感度最高,位置相關保護效果指標次之,卸壓程度相關保護效果指標的敏感度最小。
2)范圍相關效果指標中,上部卸壓角受側壓系數影響最敏感,敏感度為1.70;下部卸壓角受層間距影響最敏感,敏感度為1.02;傾向保護范圍長度受層間距影響最敏感,敏感度為-1.89。