基于Rhino3D-FLAC3D 耦合技術的庫區(qū)下伏綜放采煤沉陷區(qū)穩(wěn)定性分析

黃健豐，吳 璋，張振振，2，李永強，張 浩

（1.中煤科工集團西安研究院有限公司，陜西 西安710077；2.煤炭科學研究總院，北京100013）

隨著我國國民經(jīng)濟的持續(xù)發(fā)展，區(qū)域性交通工程、能源工程、水利水電工程等選址規(guī)劃常與采煤沉陷區(qū)產(chǎn)生疊壓現(xiàn)象，尤其部分采煤沉陷區(qū)工程地質(zhì)條件復雜、地表移動監(jiān)測數(shù)據(jù)缺失，使外部荷載作用下采煤沉陷區(qū)穩(wěn)定性評價存在巨大困難，對建（構）筑物的建設及運營構成重大安全隱患[1-3]。國內(nèi)外相關學者對采煤沉陷區(qū)穩(wěn)定性開展了大量研究工作，取得了系列研究成果。李東陽等[4-5]采用砌筑法建立了淺部采煤沉陷區(qū)三維相似實驗模型，研究了地面荷載作用下淺部采煤沉陷區(qū)覆巖破壞演變過程；楊利民等[6]提出了矩形截面采空區(qū)上方路基瞬時沉降預計公式，計算結果與現(xiàn)場監(jiān)測數(shù)據(jù)一致性良好；鄧喀中等[7-8]基于采空區(qū)覆巖移動規(guī)律建立空洞空隙預計模型，可為老采空區(qū)穩(wěn)定性評價提供參考；張志祥等[9-10]通過相似模擬實驗研究了交通荷載作用下采煤沉陷區(qū)覆巖“活化”規(guī)律；張向東等[11]基于彈性板理論建立了動靜荷載耦合作用下采煤沉陷區(qū)活化力學模型，研究了地震與采動耦合作用下采煤沉陷區(qū)覆巖變形特征；Seryakov[12]研究了采空區(qū)覆巖應力狀態(tài)，提出了一種圍巖破壞狀態(tài)估計方法；Rezaei[13]基于最小勢能原理研究了導水斷裂帶穩(wěn)定條件，推導了長壁開采地表最大下沉量預計公式。上述研究推動了采煤沉陷區(qū)在外加荷載作用下穩(wěn)定性研究的進展，但關于我國黃土高塬綜放采煤沉陷區(qū)在外部荷載作用下穩(wěn)定性研究還較少，不能完全滿足工程實踐需求。

針對上述問題，以紅巖河水庫下伏綜放采煤沉陷區(qū)為研究對象，聯(lián)合地表等高線數(shù)據(jù)及勘探地質(zhì)信息通過Rhino3D快速構建三維數(shù)值模型，將其導入FLAC3D數(shù)值模擬軟件，結合鉆孔沖洗液漏失量實測法，開展了綜放采煤沉陷區(qū)覆巖裂隙發(fā)育特征及庫區(qū)水體荷載作用下采煤沉陷區(qū)穩(wěn)定性模擬研究，以期為類似采煤沉陷區(qū)穩(wěn)定性評價及防治技術提供參考。

1 工程概況

紅巖河水庫位于彬州市西北部，是以城鎮(zhèn)居民生活、工業(yè)生產(chǎn)供水和攔沙為主要用途的Ⅲ等中型水利工程，能夠有效緩解居民生活供水及工業(yè)用水保證率低的現(xiàn)狀，為彬州市打造關中能源基地發(fā)揮巨大作用。勘查表明，庫區(qū)淹沒線與火石咀煤礦8712 綜放采煤沉陷區(qū)存在疊壓現(xiàn)象，采煤沉陷區(qū)覆巖變形對水庫蓄水及運營構成重大安全隱患，亟需研究水體荷載作用下采煤沉陷區(qū)穩(wěn)定性，研究區(qū)域地理位置圖如圖1。

圖1 研究區(qū)域地理位置圖Fig.1 Location map of the study area

火石咀煤礦采用傾斜長壁綜放開采技術，主采煤層為侏羅系延安組4-2煤層，4-2煤為近水平煤層，平均厚度為厚7.1 m，結構簡單。8712 工作面采用自然垮落法管理頂板，傾向長度為210 m，走向長度為2 250 m。地表起伏多變，溝壑縱橫，大部分被第四系、新近系松散層覆蓋，紅巖河溝內(nèi)側有華池組、洛河組巖層出露，井田內(nèi)地層由老至新依次為：侏羅系富縣組（J1f）、延安組（J2y）、直羅組（J2z）及安定組（J2a）；白堊系宜君組（K1y）、洛河組（K1l）、華池組（K2h）；新近系（N）；第四系（Q）。

2 三維數(shù)值模型可視化建模方法

2.1 模型構建流程

FLAC3D是Itasca 公司基于連續(xù)介質(zhì)理論及顯式有限差分法開發(fā)的三維數(shù)值模擬軟件，廣泛應用于巖土工程、礦業(yè)工程、地質(zhì)工程等領域，尤其適用于復雜應力條件下巖土體變形破壞機理、地質(zhì)災害孕災機制及安全風險評估等研究，然而該軟件存在前期建模能力較弱、靈活性差等問題，難以滿足構建大型復雜地質(zhì)體三維數(shù)值模型的需求。針對上述問題，眾多學者[14-18]對三維數(shù)值模型構建方法展開研究，取得系列研究成果。但多依托數(shù)個建模軟件，操作復雜，且通過計算機語言編寫數(shù)據(jù)轉化接口，對使用者有較高編程水平要求。將矢量化地形圖、鉆孔地質(zhì)信息及采掘工作面輪廓導入Rhino3D工作空間，通過布爾運算對已有信息進行處理、融合，構造相應的地質(zhì)體和結構面，三維地質(zhì)體建模流程如圖2。

圖2 三維地質(zhì)體建模流程Fig.2 3D geological body modeling process

2.2 三維數(shù)值模型構建

從不連續(xù)的原始數(shù)據(jù)點源重新建立連續(xù)的地質(zhì)面，通常采用插值法或近似擬合法。插值法通過給定數(shù)據(jù)點剛性地形成數(shù)據(jù)曲面，不改變數(shù)據(jù)點的位置坐標，在原始點源數(shù)據(jù)較少時外推精確度較差；近似擬合法是通過某種數(shù)學方法構建最大近似平滑表面，這種方法會改變部分數(shù)據(jù)點的坐標，不依賴于點源數(shù)量，但難以達到高精度標準。因此，采用2 種方法相結合構建復雜地質(zhì)體數(shù)值模型。

以研究區(qū)域地表模型構建為例，首先通過三維GIS 軟件收集研究區(qū)域數(shù)字化地形圖，提取地表矢量化等高線數(shù)據(jù)，將等高線數(shù)據(jù)導入Rhino3D工作空間離散提取點源三維坐標，地表等高線點源如圖3。

圖3 地表等高線點源Fig.3 Point sources of contour lines

提取點源初始數(shù)據(jù)集合為T=｛Rs，（s=0，1，2，…，m-1）｝，分布均勻子集為T1=｛Ki，（i=0，1，2，…，r）｝，分布不均勻子集為T2=｛Kj，（j=0，1，2，…，m-1-r）｝。集合T 采用插值蒙皮法構造曲面，不改變點源的坐標，其中子集T1可形成均勻曲面，由于子集T2點源不均造成曲面存在尖銳區(qū)域，不利于后續(xù)數(shù)值計算分析，插值法構建曲面如圖4。

圖4 差值法構建曲面Fig.4 Difference method to construct a surface

Rhino3DNURBS 技術通過曲線起點、終點、控制點及編輯點完成曲線的編輯，曲線的形狀主要由曲線的控制點和編輯點決定，通過Rhino3DNURBS 技術利用給定的精度以及控制點源數(shù)據(jù)可實現(xiàn)地表插值-逼近擬合完整的NURBS 擬合地表曲面如圖5。

圖5 NURBS 擬合地表曲面Fig.5 NURBS fitting surface surface

從地表某端點沿著z 軸負方向建立體模型，結合工程地質(zhì)資料及現(xiàn)場鉆孔勘查數(shù)據(jù)對地質(zhì)體進行切割形成地層界面，通過切割、組合構建的研究區(qū)域三維數(shù)值模型如圖6。x 軸方向1 640 m，y 軸方向910 m，z 軸方向430～662 m，8712 開采工作面走向長度為600 m，傾向長度為210 m，采厚為7.1 m，六面體網(wǎng)格共686 504 個，節(jié)點共713 460 個。

圖6 研究區(qū)域三維數(shù)值模型Fig.6 3D numerical model of study area

3 綜放采煤沉陷區(qū)覆巖裂隙發(fā)育特征

基于開采損害理論，采煤沉陷區(qū)覆巖裂隙主要分為邊緣裂隙和動態(tài)裂隙。邊緣裂隙一般位于沉陷盆地的邊緣區(qū)域，以“帶狀”、“橢圓狀”的形態(tài)分布在沉陷盆地的邊緣；動態(tài)裂隙則位于工作面上覆巖層中，并隨著工作面的推進不斷擴展，直至覆巖形成新的平衡結構，其中導水斷裂帶覆巖存在較大自由空間，為礦區(qū)地下水資源流失提供通道。庫區(qū)水體長期作用下沿沉陷盆地邊緣拉張裂隙帶補給含水層，通過原生裂隙、與導水斷裂帶間接溝通形成的水力通道，對水庫運營構成重大安全威脅。因此，必須查明采煤沉陷區(qū)導水斷裂帶分布特征。

采用鉆孔沖洗液漏失量觀測法獲取導水斷裂帶發(fā)育高度，判別標準為：①鉆孔沖洗液漏失量急劇增加；②吸風現(xiàn)象明顯。在采煤沉陷區(qū)地表布設探測孔K-14，地面高程為910 m，孔底高程為468.7 m，通過統(tǒng)計分析鉆孔沖洗液消耗曲線圖（圖7），孔深為306 m 時漏失量由0.14 L/（m·s） 急劇上升為18.49 L/（m·s），吸風現(xiàn)象明顯，認為已經(jīng)到達導水斷裂帶頂端。確定8712 工作面導水斷裂帶高度為135.3 m，裂采比為19.05。

圖7 鉆孔沖洗液消耗曲線圖Fig.7 Consumption curve of drilling fluid

4 庫區(qū)下伏綜放采煤沉陷區(qū)穩(wěn)定性

4.1 邊界條件及力學參數(shù)選取

采用Mohr-Coulomb 彈塑性本構破壞準則，利用FLAC3D軟件新功能對數(shù)值模型表面進行分組，將模型East 和West 面作為x 軸方向的邊界條件，模型South 和North 作為y 軸方向的邊界條件，Bottom 面約束z 軸邊界條件，TOP 面為自由邊界。參考工程地質(zhì)勘查資料及相鄰井田巖石物理力學實驗數(shù)據(jù)，確定模擬的巖體物理力學參數(shù)見表1。

4.2 原巖應力場反演情況

初始應力場計算采用彈性本構準則，計算數(shù)值模型僅在自重作用下的應力狀態(tài)，得到研究區(qū)域的初始應力分布如圖8。沿工作面傾向中線y=505 m布置剖面，最大垂直應力為13.97 MPa，應力場呈“波浪形”層狀分布，是由于黃土塬凸起埋深增加造成的。

4.3 8712 工作面導水斷裂帶發(fā)育特征

為深入分析隨工作面推進覆巖變形破壞全過程，采用分步開挖的方式對8712 工作面開采過程進行模擬，賦予開挖部分為空本構模型，每次工作面推進40 m，共開挖10 步，累計開采400 m。沿工作面傾向中線y=505 m 布置剖面，8712 工作面塑性區(qū)演化分布如圖9。

表1 巖體物理力學參數(shù)Table 1 Rock physical and mechanical parameters

圖8 初始應力分布Fig.8 Initial stress field distribution

圖9 8712 工作面塑性區(qū)演化分布Fig.9 Evolution distribution of 8712 working face plastic zone

通過分析8712 工作面推進過程的塑性區(qū)演化分布，采動覆巖先以剪切破壞為主，導致頂板裂隙初步發(fā)育，逐步發(fā)展成為拉伸破壞產(chǎn)生垮落帶。隨著工作面的推進，拉伸破壞區(qū)的范圍逐漸增大，上部剪切破壞區(qū)域也逐漸擴大，工作面前方煤壁和開切眼處的圍巖主要發(fā)生了剪切破壞，覆巖塑性區(qū)分布整體成“馬鞍”型。塑性區(qū)發(fā)育高度為113.92～138.24 m，最大裂采比為19.47，與鉆孔沖洗液漏失量觀測法實測導水斷裂帶高度（135.3 m）基本一致，反映出三維數(shù)值模型、本構條件及巖體力學參數(shù)選取的合理性。

4.4 庫區(qū)水體荷載作用下采煤沉陷區(qū)穩(wěn)定性

在模型上表面按照正常蓄水位（903 m）施加靜水壓力，庫區(qū)水體荷載分布如圖10。模型采用Mohr-Coulomb 破壞準則，將模型四周及底部作為邊界條件，TOP 面為自由邊界。

圖10 庫區(qū)水體荷載分布Fig.10 Water load distribution in the reservoir area

庫區(qū)水體荷載對采煤沉陷區(qū)的覆巖位移及塑性區(qū)影響如圖11、圖12。

圖11 蓄水前后覆巖位移對比Fig.11 Comparison of overburden displacement before and after water storage

圖12 蓄水前后覆巖塑性區(qū)對比Fig.12 Comparison of plastic area of overburden before and after water storage

庫區(qū)水體荷載作用下，采煤沉陷區(qū)覆巖破壞范圍增大，破壞逐漸向上發(fā)展，使原本處于穩(wěn)定狀態(tài)的采煤沉陷區(qū)產(chǎn)生“活化”現(xiàn)象，沿工作面傾向中線y=505 m 布置剖面，水庫蓄水后，地表垂向位移最大值為686.4 mm，比蓄水前垂直位移增加334.7 mm。導水斷裂帶最大發(fā)育高度為142 m，比蓄水前增加3.76 m。

5 結 語

1）采用Rhino3D-FLAC3D耦合技術結合8712 綜放采煤沉陷區(qū)地形特征及地層巖性快速構建了研究區(qū)域三維數(shù)值模型。

2）8712 綜放采煤沉陷區(qū)覆巖導水斷裂帶整體呈“馬鞍”狀，最大發(fā)育高度為138.24 m，裂采比為19.47，與鉆孔沖洗液漏失量觀測法實測值（135.3 m）基本一致。

3）庫區(qū)水體荷載作用下，綜放采煤沉陷區(qū)地表垂向位移增加334.7 mm，導水斷裂帶發(fā)育高度增加3.76 m。