王彥海 時文峰 李清泉 陳波

(三峽大學(xué) 宜昌443000)

引言

在過去的幾十年里，在大量的煤炭資源被開采的同時還造成了大規(guī)模、大面積的采空區(qū)。 這些采空區(qū)極易造成上覆巖體的冒落、彎曲乃至斷裂，并且容易使圍巖的力學(xué)強度降低，故而導(dǎo)致采空區(qū)上方的建筑物地基承載力降低，嚴(yán)重的則會造成地表塌陷、大面積的農(nóng)田損毀、森林植被的破壞、水土流失加劇和土地沙漠化程度的增高等[1，2]，使得礦區(qū)居民的財產(chǎn)乃至生命安全受到嚴(yán)重威脅。 針對采空區(qū)地表塌陷與變形的研究問題，國內(nèi)外學(xué)者采用了包括概率積分法、巖石力學(xué)理論和數(shù)值模擬分析等技術(shù)和方法進行了一系列的研究[3-5]。 但由于前兩者方法的理論公式較為復(fù)雜，且計算過程較為冗長，在實際的施工中難以精準(zhǔn)測得各項參數(shù)，以及其他眾多的多變因素是此類方法無法考慮到的。 近年來，隨著力學(xué)機理研究的深入和計算機技術(shù)的發(fā)展，數(shù)值模擬方法在采空區(qū)沉降變形研究中得到了廣泛應(yīng)用。如陳贊成、柴靜靜、趙付萬、張德鵬等學(xué)者[6-9]就基于ANSYS 有限元分析軟件，對實際工程中的采空區(qū)進行仿真計算并取得相應(yīng)成果，對運用ANSYS 計算采空區(qū)的地表沉降做出貢獻和指導(dǎo)。但這些研究內(nèi)容僅為針對某一實際工程進行展開，所得結(jié)論不具備一定的普遍性，對于其他的工程而言也不具備一定的參考性。 又如孫超、容宇等學(xué)者[10-12]就基于ANSYS 對采空區(qū)地表沉降影響因素展開研究，得出了各因素與采空區(qū)地表沉降的影響關(guān)系。 但上述這些研究都有著考慮因素較少、且沒考慮到各因素變化下的相互影響作用、設(shè)計工況較少和對土層及采空區(qū)的建模較為簡略單一等缺點。

對此，本研究選取了數(shù)值模擬分析法，通過運用ANSYS 有限元軟件對采空區(qū)進行數(shù)值模擬分析，利用了多種土壤參數(shù)模擬出較為復(fù)雜的土壤環(huán)境，并模擬設(shè)計出共計29 種計算工況，對采空區(qū)地表變形進行計算分析，得出其地表變形與采深、采寬、采厚以及上覆巖性之間的規(guī)律。分別研究了單因素變化下與多因素變化下的各自規(guī)律并得出結(jié)論，且深入地對采空區(qū)地表變形影響因素進行分析，得出了各因素對地表沉降塌陷速度的影響，并且還得出各因素的影響主次關(guān)系。 因此對有限元數(shù)值模擬方法在采空區(qū)地表沉降研究中的應(yīng)用以及采空區(qū)地表的處理和礦區(qū)安全生產(chǎn)有一定指導(dǎo)意義。

1 采空區(qū)數(shù)值模擬分析方案設(shè)計

根據(jù)在采空區(qū)地表變形問題的研究中發(fā)現(xiàn)，影響采空區(qū)穩(wěn)定的因素包括礦層的埋深條件、覆巖的力學(xué)性質(zhì)、采空區(qū)的幾何尺寸等。 其中對采空區(qū)地表變形影響尤為顯著的有：采掘方式、頂板管理方法、回采率的高低等采空區(qū)工程活動[13]。 在工程的實踐中，運用數(shù)值方法研究一個已經(jīng)探明且確定的采空區(qū)的地表變形，其數(shù)值計算模型的建立將會變得非常簡單。 但遺憾的是，在工程中對于采空區(qū)的幾何尺寸不但往往無法探明清楚，而且對于其回采率以及諸多的頂板管理方法也無法精確掌握。 故而想要詳細(xì)地了解采空區(qū)的地下挖開情況是比較困難的[14]。 為了能夠便捷地探究采空區(qū)地表變形的影響因素，本文研究簡化了采空區(qū)的實際模型。 將影響采空區(qū)地表變形的因素簡化為4 個，并均給出4 個不同的水平。 對于多因素多等級的實驗工況，其實驗工作量巨大，且各個因素之間還可能會有交互作用，相互產(chǎn)生影響。 在研究多因素多等級的實驗設(shè)計方法中，正交實驗設(shè)計是一種高效、快速且經(jīng)濟的多因素實驗設(shè)計方法。 故本研究選取了正交實踐設(shè)計，運用正交實驗表設(shè)計出四因素四水平的正交工況共計16 個。 分別計算出各工況下采空區(qū)地表沉降數(shù)值，并得出各因素對采空區(qū)地表沉降的影響趨勢函數(shù)曲線。 為驗證正交實驗設(shè)計運用在本研究的準(zhǔn)確性，又設(shè)計出在其余三因素等級不變的情況下，某一單因素等級變化的單因素實驗工況共計13 個。 并以單因素實驗工況的計算結(jié)果為基準(zhǔn)驗證正交實驗工況計算結(jié)果的準(zhǔn)確性。 以此完成各個因素對采空區(qū)地表變形的研究。

1.1 采空區(qū)的模型建立

為尋求各因素對采空區(qū)地表沉降的影響結(jié)果，本研究利用了模擬工況法，選取了采深、采厚、采寬和上覆巖性作為采空區(qū)地表沉降的主要的影響因素并將其劃分為4 個等級，將其進行匯總見表1。

表1 釆空區(qū)地表變形影響因素和等級Tab.1 Factors and grades of surface deformation

其中上覆巖性1 級對應(yīng)材料折減系數(shù)取1.0，上覆巖性2 級對應(yīng)材料折減系數(shù)取0.95，上覆巖性3 級對應(yīng)材料折減系數(shù)取0.9，上覆巖性4 級對應(yīng)材料折減系數(shù)取0.85。

1.采空區(qū)地表變形影響分析正交工況組合

為探究各影響因素對于采空區(qū)地表沉降的影響，本研究運用正交實驗設(shè)計研究，將4 個主要影響因素的4 個等級按照正交模型表組合，得出正交工況組合見表2。

表2 多因素正交工況L16(44)Tab.2 Multi-factor orthogonal condition table L16(44)

2.采空區(qū)地表變形影響單因素工況組合

為論證多因素正交實驗設(shè)計的準(zhǔn)確性，根據(jù)表1 中采空區(qū)地表變形因素和等級，將影響因素包括采深、采寬、采厚和上覆巖性4 個因素按4個等級組合起來，在其余因素不變的情況下，研究單因素對采空區(qū)地表變形影響規(guī)律。 單因素工況組合見表3。

表3 單因素工況組合Tab.3 Single factor working condition combination

為研究采深對采空區(qū)地表變形的影響規(guī)律，工況 1、2、3、4 保持了采厚，采寬以及上覆巖性不變，將采深逐步增加等級。 為研究采厚對采空區(qū)地表變形的影響規(guī)律，工況 1、5、6、7 保持采深，采寬以及上覆巖性不變，將采厚逐步增加等級。 為研究采寬對采空區(qū)地表變形的影響規(guī)律，工況 1、8、9、10 保持采深，采厚以及上覆巖性不變，將采寬逐步增加等級。 為研究上覆巖性對采空區(qū)地表變形的影響規(guī)律，工況1、11、12、13 保持采深，采厚以及采寬不變，將上覆巖性逐步增加等級。

2 ANSYS模型建立

由于采空區(qū)地表的實際沉降問題較為復(fù)雜，模型建立難度高，故本文研究將其簡化為平面應(yīng)變問題，進行彈塑性數(shù)值模擬分析[15]。 首先用到了solid45 單元建立土壤模型，為了模擬采空區(qū)周圍巖層的真實應(yīng)力狀態(tài)需用到ANSYS 單元生死功能，每次的模擬計算過程分為兩步進行：第一步先建完土壤模型后，將土壤模型和采空區(qū)模型分別賦予不同的單元號并均保留參數(shù)用于模擬開挖前巖層僅有自重作用時的應(yīng)力狀態(tài)，以獲取模擬計算區(qū)域的初始應(yīng)力狀態(tài)以及僅由重力而引起的初始位移； 第二步通過ANSYS 單元生死功能，將采空區(qū)部分的單元“殺死”，模擬礦體被采出的情況，并計算出此狀態(tài)下的應(yīng)力和變形。 最后通過ANSYS 后處理中的荷載步相減求出采空區(qū)地表僅由開采礦體引起的真實變形。

ANSYS 計算模型在水平方向取800m，垂直方向取200m，一共分為五層：地表覆蓋土層為黃土層，厚度為20m； 第二層為鈣質(zhì)泥巖層，厚度為35m； 第三層為砂質(zhì)泥巖，厚度為20m； 第四層為中粒砂巖，厚度為20m； 最后一層為粗砂巖，厚度為105m。 采空區(qū)的大小和位置由上述工況確定。 由于模型較大，為使計算結(jié)果更加精確，對采空區(qū)所在的單元采用了ANSYS 網(wǎng)格局部細(xì)化功能。 模型底部采用固定邊界施加全部約束，模型兩側(cè)采用滾軸邊界僅施加垂直于該面的約束，取其開采深度為100m、采寬為100m、采厚為2.5m 的工況，并將各材料賦予不同屬性并用不同顏色標(biāo)出的模型示意如圖1 所示。 各層的巖土力學(xué)參數(shù)見表4。

圖1 ANSYS 建模示意Fig.1 Schematic diagram of ANSYS modeling

表4 各巖土層及填充材料物理力學(xué)參數(shù)Tab.4 Physical and mechanical parameters of each rock layer and filling material

3 計算結(jié)果與分析

3.1 正交工況結(jié)果匯總分析

根據(jù)表2 采空區(qū)計算模型正交模型表L16(44)運用ANSYS 計算出正交模型下16 個工況的結(jié)果并將其匯總見表5。

表5 正交工況L16(44)計算結(jié)果Tab.5 Calculation results of orthogonal working conditionL16(44)

將正交工況4 個影響因素中每個因素的4 個等級的4 個工況采空區(qū)地表沉降的極差計算值進行匯總，地表沉降正交極差分析見表6。

表6 地表沉降正交極差分析Tab.6 Surface settlement orthogonal range analysis

通過其中極差一欄的數(shù)據(jù)大小比較，可得出如下結(jié)論：

(1)采深因素影響分析：采深的4 個等級其ANSYS 模擬計算得出的采空區(qū)地表最大沉降由0.99mm 降低至 0.42mm。 減小比例依次為：11.1%、44.3%、14.3%。 其沉降大小隨著采深的增加而明顯減小，其原因可能由于采空區(qū)頂板及其上部的部分巖層與整體分離，破碎成小塊巖塊而不規(guī)則地填充了采空區(qū)，且采空區(qū)上方裂縫帶上本身可以自成平衡壓力拱，采空區(qū)越深，其壓力拱效果則越顯著。

(2)采厚因素影響分析：采厚的4 個等級其ANSYS 模擬計算得出的采空區(qū)地表最大沉降由0.34mm 增加至 1.09mm。 增大比例為：50%、64.7%、29.7%。 其原因可能由于隨著采厚的增加會引起較大的塌落高度，隨之裂縫帶、彎曲帶的影響范圍則會上升，對地表沉降的影響將必然增加。

(3)采寬因素影響分析：采寬的4 個等級其ANSYS 模擬計算得出的采空區(qū)地表最大沉降由0.25mm 增加至 1.42mm。 增大比例為：16%、179%、75%。 其原因可能由于隨著開采寬度的增加，采空區(qū)工作面寬度超過其極限值時，控制層將會因拉應(yīng)力超過其抗拉強度而斷裂，會隨著下覆巖層的冒落而下沉，地表沉降會明顯增加。

(4)上覆巖性因素影響分析：上覆巖性的4個等級其ANSYS 模擬計算得出的采空區(qū)地表最大沉降由0.64mm 增加至0.91mm。 增大比例為：-15%、29.6%、30%。 其原因可能由于隨著上覆巖性的折減越大時，對應(yīng)的巖性越為軟弱，在其他條件相同的情況下，堅硬的上覆巖層地表沉降值將小于軟弱的上覆巖層地表沉降值。

(5)根據(jù)表 6 可知采深、采厚、采寬、上覆巖性四者的極差數(shù)據(jù)有：1.18 ＞0.74 ＞0.57 ＞0.37。 則在本模擬環(huán)境下，影響地表位移的4 個因素中，按照其對地表沉降影響的主次排序依次為：采寬＞采厚＞采深＞上覆巖性。

(6)給出多因素變化下正交工況采空區(qū)地表沉降隨采深、采厚、采寬以及上覆巖性變化的趨勢如圖2 所示。

圖2 正交實驗下各因素各等級沉降值Fig.2 Settlement values of various factors under orthogonal experiment

綜上所述，多因素變化下正交實驗工況的計算結(jié)果完全符合單因素工況的各因素對采空區(qū)地表沉降的影響規(guī)律。 說明正交實驗工況計算數(shù)據(jù)較為準(zhǔn)確。

3.2 單因素工況計算結(jié)果匯總分析

通過ANSYS 軟件計算出共13 個單因素工況各工況的最大垂直移動、水平移動，見表7。 4個單因素影響下各自變化對地表沉降及水平移動的影響如圖3 所示。

表7 單因素工況計算結(jié)果匯總Tab.7 Summary of calculation results of single factor working conditions

圖3 各單因素變化的影響Fig.3 Effect on various factors

通過上述單因素變化下工況計算結(jié)果可得，其結(jié)果與正交實驗所得結(jié)果趨勢一致，表明正交實驗設(shè)計工況的可行性與準(zhǔn)確性。 且根據(jù)控制變量法下的單因素工況結(jié)論，還能得出以下結(jié)論：

(1)隨著采深的增加采空區(qū)地表沉降減少的速率越慢。

(2)隨著采厚的增加采空區(qū)地表沉降增加的速率越快。

(3)隨著采寬的增加采空區(qū)地表沉降增加的速率越快。

(4)由于折減系數(shù)差值較小，導(dǎo)致上覆巖性對于采空區(qū)地表沉降的影響很小，這里難以體現(xiàn)出其對地表沉降變化速率的影響。

4 結(jié)論

本文研究運用了有限元分析軟件ANSYS 對模擬設(shè)計出的正交實驗工況進行數(shù)值模擬分析，分別構(gòu)建了各工況下的采空區(qū)三維有限元計算模型，且單因素變化下工況組合和多因素變化下的正交工況組合所計算出的結(jié)果均表現(xiàn)一致。 對各工況采空區(qū)地表沉降以及水平移動的結(jié)果進行了比較分析，從而得出以下結(jié)論：

1.對于探究采空區(qū)的地表穩(wěn)定性而言，運用數(shù)值模擬的方法來計算其地表變形數(shù)值是可行的，只要能給出所需的計算參數(shù)，即可模擬出采空區(qū)地表的各項形變，以預(yù)先了解采空區(qū)地表變形大小，從而可以為后續(xù)工作提供必要的指導(dǎo)。

2.在本模擬環(huán)境下，各影響因素對采空區(qū)地表變形的影響趨勢為：采深對采空區(qū)地表沉降的數(shù)值呈負(fù)相關(guān)，且隨著采深的增加采空區(qū)地表沉降減少的速率越慢。 而采厚、采寬和上覆巖性均對采空區(qū)地表沉降的數(shù)值呈正相關(guān)，且隨著采厚與采寬的增加，采空區(qū)地表沉降增加的速率越快。

3.在本模擬環(huán)境下，通過對正交工況組合的計算結(jié)果，建立了極差分析表得出各影響因素對地表沉降影響的主次排序依次為：采寬＞采厚＞采深＞上覆巖性。