防水煤柱安全評價的流-固耦合分析

戰玉寶，劉煥芝，尤春安

( 1.山東科技大學資源與土木工程系，山東 泰安 271019；2.山東科技大學泰安校區，山東 泰安 271019 3. 山東科技大學土建學院，山東 青島 266510)

突水是煤礦防水煤柱失效的重要形式，引起的損失是巨大的。突水事故輕者會沖垮工作面、淹沒設備、增加礦井的排水負擔，嚴重者會造成人員傷亡、淹井等。全國已經發生許多起近松散層開采上部含水層造成的水砂突涌和淹井事故，造成重大經濟損失。防水煤柱失穩突水，不僅涉及到煤柱在水壓力作用下的力學性質與長期穩定性問題，還涉及到滲流場與應力場的耦合問題。[1]

含水松散層下采煤，是水體下采煤的重要方面。近年來，華東、華北、西北等地的許多煤礦，開采上限不斷提高，己經在不同礦區、不同富水程度的松散含水層下開展近松散層開采。對防水煤柱的合理留設進行了深入研究，[2-4]并取得一批重要理論成果。于進廣等[5]對中等含水層下留設防砂煤柱的安全開采機理進行了研究，廖學東、疏開生[6]對松散含水層下采煤合理煤巖柱高度的確定進行了探討，蔡榮[7]對水體下采煤巖層的移動與保護層理論及應用進行了研究，尹先尚[8]對防水煤柱留設對圍巖變形的影響進行了數值模擬方面的研究。中國礦業大學的狄乾生[9]和隋旺華[10]，結合“三下”采煤和開采覆巖移動的有限元分析，對工程地質模型進行了深入研究，建立了厚松散含水層下提高開采上限的工程地質模型。

巖土體地下滲流的流-固耦合問題的研究已經開展多年，但流-固耦合理論應用在防水煤柱安全評價方面的研究非常少，這方面的文獻能檢索到的幾乎沒有。本文就是應用有限差分軟件FLAC3D，對某礦煤礦四采二層煤防水安全煤柱進行滲流場與應力場的耦合分析。

1 計算分析方法與原理[11-13]

FLAC3D模擬多孔介質(如煤巖體)中流體流動時，流體的模擬獨立于結構計算。其主要通過孔隙水壓力的消散，引起巖體中位移的變化，這一過程包含兩種力學效果。第一、孔隙水壓力的變化，引起結構體中有效應力的變化；第二、孔隙水壓力的變化，又引起流體區域的變化。流體在孔隙介質中的流動依據Darcy定律，流-固耦合過程滿足Biot方程。

FLAC3D中，實現流-固耦合分析的幾個關鍵方程為[13]：

(1)平衡方程

對于小變形，流體質點平衡方程為：

(1)

式中：qi是滲流速度(m/s)；qv是被測體積的流體源強度(1/s)；ζ是單位體積孔隙介質的流體體積變化量。而

(2)

式中：M是Biot模量(N/m2)；p是孔隙壓力；α是Biot系數；ε是體積應變；T是溫度；β是考慮流體和顆粒熱膨脹系數(1/°C)。

(2)運動方程

流體的運動用Darcy定律來描述。對于均質、各向同性固體和流體密度是常數的情況，這個方程具有如下形式：

(3)

式中：k是介質的滲透系數(m2/pa·s)；ρf是流體密度(kg/m3)；gi(i=1，3)是重力加速度的三個分量(m/s2)。

(3)本構方程

流體的流動導致孔隙介質中，孔隙壓力p、飽和度s、體積應變e和溫度T的改變，孔隙流體方程為

(4)

體積應變的改變，引起流體孔隙壓力的變化，反過來，孔隙壓力的變化也會導致體積應變的發生。

孔隙介質本構方程的增量形式為

(5)

(4)相容方程

應變率和速度梯度之間的關系為

(6)

式中：υi是介質中某點的速度。

(5)邊界條件

在計算中，有四種類型的邊界條件，他們分別是：①給定孔隙水壓力；②給定邊界外法線方向流速分量；③透水邊界；④不透水邊界。

不透水邊界程序中默認，透水邊界采用下式給出：

(7)

式中：qn是邊界外法線方向流速分量；h是滲漏系數(m3/N·s)；p是邊界面處的孔隙水壓力；pe是滲流出口處的孔隙水壓力。

2 數值計算方案

2.1 工程概況

山東某礦礦井田西翼四采二層煤煤層露頭區，位于該礦井四采風井以西，煤層露頭東西走向長450m，走向92°，平均傾角25°。露頭區地面位于一條河的北岸，東部為該河所覆蓋。由于附近磚廠取土，村民挖砂，致使煤層露頭上覆的第四系表土覆蓋層受到不同程度的破壞。加上煤層淺部被非法開采、亂采和盜采，近年來四采區二層煤露頭區域地表明顯地下沉，形成連續的斑裂區，與地表形成直接的水力聯系。斑裂區位于該河北岸，是雨季洪水淹沒范圍，因此存在地表水泄漏威脅。

2.2 計算模型與參數

數值模擬時，模型取為340m×160m×1m，水平面內向右為X軸，垂直向上為Y軸，X、Y、Z軸成右手系。計算范圍取為-50m≤X≤290m，0m≤Y≤160m，-1m≤Z≤0m。對模型進行劃分，總共得到2106個單元，4388個節點。

采用摩爾-庫侖材料模擬巖石與煤層，計算參數根據該礦提供的巖石的物理力學性質選定。具體的參數見表1。

2.3 初始地應力場

根據煤礦提供的資料，初始地應力在垂直方向按巖體自重考慮，水平地應力根據勘測單位提供的側壓力系數為0.6考慮。

3 數值計算結果分析

3.1 不考慮流-固耦合時計算結果分析

此種情況下，只考慮采空區頂板以及防水煤柱受自重和上覆巖層的作用，分別對20m、24m、28m、32m不同長度的防水煤柱進行數值模擬，得到了開挖后采空區頂板的塑性區分布以及防水煤柱的塑性區分布(圖1)。

圖1 未考慮流-固耦合時不同長度的安全防水煤柱的塑性區

圖2 流-固耦合時不同長度的安全防水煤柱的塑性區分布

(1)采空區頂板上端20m范圍內塑性區較大，有向上發展的趨勢；在采空區下端30m范圍出現較大塑性區；中間塑性區較小，比較平緩；整個塑性區呈現為兩端高、中間低的“馬鞍”形狀。

(2)當煤柱大于24m時，煤柱上方采空區下端的塑性區比留設20m煤柱時的塑性區小，隨著煤柱長度的增加，塑性區分布開始趨于穩定。

(3)安全煤柱長度從20m一直到32m的情況下，整個煤柱均未全部進入塑性狀態。

3.2 考慮流-固耦合的計算結果分析

當上部采空區與地表形成水力連接時，考慮采空區頂板以及安全煤柱受自重和上覆巖層的作用，同時考慮水壓力的作用。分別對20m、24m、28m、32m不同長度的防水煤柱進行數值模擬，得到了開挖后采空區頂板的塑性區分布以及防水煤柱的塑性分布(圖2)

分析圖2可以得出：

(1)上部采空區頂板上端20m范圍內，塑性區有繼續向上發展的趨勢，在采空區下端30m范圍，也出現更大的塑性區，中間塑性區較小，比較平緩；整個塑性區呈現為兩端高中間低的“馬鞍”形狀。

(2)當煤柱留設20m和24m時，煤柱全部進入塑性狀態，并且水能夠透過煤柱及其上方的裂隙帶，形成水力通道，會出現水潰入下部采空區的現象，對安全生產構成威脅。

(3)當煤柱留設28m、32m時，煤柱未全部進入塑性狀態，雖然非塑性區范圍不大，但是可以通過采取治理措施加以控制，保證生產安全。

(4)考慮流-固耦合后，底板出現了不同程度的塑性區。

4 結論

對比不考慮流-固耦合與考慮流-固耦合作用兩種情況下防水煤柱的數值分析結果，可以得出如下結論：

(1)考慮流-固耦合時與不考慮耦合時，塑性區都呈現為兩端高、中間低的“馬鞍”形狀，但考慮耦合時塑性區的范圍比較大。

(2)考慮流-固耦合作用時，能很好的模擬底板的受力狀態與變形破壞規律，與不考慮流-固耦合相比，煤柱留設需加長。

(3)采用FLAC3D，可以很好地對防水煤柱進行流-固耦合數值模擬，模擬結果與實際情況更加接近。

[1] 張金才, 張玉卓, 劉天泉. 巖體滲流與煤層底板突水[M]. 北京: 地質出版社, 2001.

[2] 許傳峰, 王思棟, 紀憲江，等. 柴里煤礦留12米煤柱的開采上限研究[J] .河北建筑科技學院學報,2001, 18(2):68-70.

[3] 隋旺華. 煤層開采防水煤巖柱工程地質預計[J]. 江蘇煤炭, 1994,(2):21-25.

[4] 馬福祥, 隋旺華. 山東橫河煤礦分層開采防水煤柱留設[J].中國地質災害與防治學報, 1994, 5(增):15-17.

[5] 于進廣, 陳風祥,等. 中等含水層下留設防砂煤柱的安全開采機理[J].淮南工業學院學報, 2001, (2): 15-17.

[6] 廖學東, 疏開生. 松散含水層下采煤合理煤巖柱高度的確定[J] .淮南礦業學院學報, 1998.8(2):29-33.

[7] 蔡 榮. 水體下采煤巖層的移動與保護層理論及應用[D].徐州: 中國礦業大學, 2003.

[8] 尹先尚. 陷落柱防水煤柱留設對圍巖變形影響的數值模擬[J].煤炭學報, 2006, 4(2):179-182.

[9] 狄乾生，隋旺華，黃山民. 開采巖層移動工程地質研究[M]. 北京：中國建筑工業出版社, 1992.

[10] 隋旺華, 董青紅,等. 工程地質模型在防水煤巖柱研究中的應用[J].中國礦業大學學報, 1999，28(4):417-420.

[11] FLAC-3D,Fast Lagrangian Analysis of Continua in 3Dimensions, Version 3.0,User’s Manual, Itasca Consulting Group, Inc. USA.

[12] 劉 波，韓彥輝. FLAC原理、實例與應用指南[M] . 北京：人民交通出版社, 2005.

[13] 李廷春，李術才，陳衛忠，等. 廈門海底隧道的流固耦合分析[J] .巖土工程學報, 2004, 26(3):397-401.