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

黃健豐，吳 璋，張振振，2，李永強(qiáng)，張 浩

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

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

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

1 工程概況

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

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

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

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

2.1 模型構(gòu)建流程

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

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

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

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

以研究區(qū)域地表模型構(gòu)建為例，首先通過三維GIS 軟件收集研究區(qū)域數(shù)字化地形圖，提取地表矢量化等高線數(shù)據(jù)，將等高線數(shù)據(jù)導(dǎo)入Rhino3D工作空間離散提取點源三維坐標(biāo)，地表等高線點源如圖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 采用插值蒙皮法構(gòu)造曲面，不改變點源的坐標(biāo)，其中子集T1可形成均勻曲面，由于子集T2點源不均造成曲面存在尖銳區(qū)域，不利于后續(xù)數(shù)值計算分析，插值法構(gòu)建曲面如圖4。

圖4 差值法構(gòu)建曲面Fig.4 Difference method to construct a surface

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

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

從地表某端點沿著z 軸負(fù)方向建立體模型，結(jié)合工程地質(zhì)資料及現(xiàn)場鉆孔勘查數(shù)據(jù)對地質(zhì)體進(jìn)行切割形成地層界面，通過切割、組合構(gòu)建的研究區(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)裂隙則位于工作面上覆巖層中，并隨著工作面的推進(jìn)不斷擴(kuò)展，直至覆巖形成新的平衡結(jié)構(gòu)，其中導(dǎo)水?dāng)嗔褞Ц矌r存在較大自由空間，為礦區(qū)地下水資源流失提供通道。庫區(qū)水體長期作用下沿沉陷盆地邊緣拉張裂隙帶補(bǔ)給含水層，通過原生裂隙、與導(dǎo)水?dāng)嗔褞чg接溝通形成的水力通道，對水庫運營構(gòu)成重大安全威脅。因此，必須查明采煤沉陷區(qū)導(dǎo)水?dāng)嗔褞Х植继卣鳌?/p>

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

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

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

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

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

4.2 原巖應(yīng)力場反演情況

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

4.3 8712 工作面導(dǎo)水?dāng)嗔褞Оl(fā)育特征

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

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

圖8 初始應(yīng)力分布Fig.8 Initial stress field distribution

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

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

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

在模型上表面按照正常蓄水位（903 m）施加靜水壓力，庫區(qū)水體荷載分布如圖10。模型采用Mohr-Coulomb 破壞準(zhǔn)則，將模型四周及底部作為邊界條件，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。導(dǎo)水?dāng)嗔褞ё畲蟀l(fā)育高度為142 m，比蓄水前增加3.76 m。

5 結(jié) 語

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

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

3）庫區(qū)水體荷載作用下，綜放采煤沉陷區(qū)地表垂向位移增加334.7 mm，導(dǎo)水?dāng)嗔褞Оl(fā)育高度增加3.76 m。