王彥海 時(shí)文峰 李清泉 陳波

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

引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 ANSYS模型建立

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

ANSYS 計(jì)算模型在水平方向取800m，垂直方向取200m，一共分為五層：地表覆蓋土層為黃土層，厚度為20m； 第二層為鈣質(zhì)泥巖層，厚度為35m； 第三層為砂質(zhì)泥巖，厚度為20m； 第四層為中粒砂巖，厚度為20m； 最后一層為粗砂巖，厚度為105m。 采空區(qū)的大小和位置由上述工況確定。 由于模型較大，為使計(jì)算結(jié)果更加精確，對(duì)采空區(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 計(jì)算結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

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

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

圖2 正交實(shí)驗(yàn)下各因素各等級(jí)沉降值Fig.2 Settlement values of various factors under orthogonal experiment

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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