劉亞群，周宏偉，李翼虎，易海洋，薛東杰

（1.中國礦業大學 力學與建筑工程學院，北京 100083；2.中鐵第五勘察設計院集團有限公司，北京 102600）

淺埋煤層開采突水潰砂的顆粒流模擬研究

劉亞群1，周宏偉1，李翼虎2，易海洋1，薛東杰1

（1.中國礦業大學 力學與建筑工程學院，北京 100083；2.中鐵第五勘察設計院集團有限公司，北京 102600）

為研究淺埋煤層開采下突水潰砂的覆松散層滲透破壞、水砂耦合流動及突涌規律，指導預防突水潰砂災害的發生提供依據，采用二維顆粒流軟件PFC2D分別對上覆松散層滲透變形破壞、水砂混合流運移過程進行數值模擬研究，建立了淺埋煤層開采突水潰砂的數值模擬模型。結果表明：上覆松散含水層的初始水壓力與裂隙通道的寬度對松散層的滲透變形破壞過程有較大影響。當含水層初始水壓力增加時，裂隙處流體的流速、顆粒間的接觸應力均明顯增加，且呈現出非線性增長趨勢；水砂混合流的突出時間以及突出速度與裂隙通道的傾角、粗糙度相關。同時根據水砂混合流的Y方向突出速度大小的波動變化規律，預測了2種可能存在的水砂突出類型；含水層的初始水壓力影響突水潰砂的發生及運移過程,由此提出了通過疏放松散層水來預防突水潰砂災害發生的方法。

淺埋；突水潰砂；PFC2D；滲透變形；裂隙通道

0 引 言

在淺埋煤層開采過程中經常會遇到上覆松散層含砂含水率較高的情況，這時由上覆巖層的冒落帶和裂隙帶形成的裂隙通道導通松散含水層，在一定條件下松散含水層發生滲透變形破壞，水砂混合流很容易通過采動過程中產生的裂隙突涌到煤層開采的工作面，此類災害常不僅淹沒礦井及開采設備，而且對井下作業人員生命安全造成嚴重威脅。

關于淺埋深煤層開采突水潰砂的研究，中科院地質所提出“強滲通道”理論[1]，認為突水關鍵在于是否具備突水通道；張世凱[2]提出了厚含水松散層薄基巖條件下開采工作面涌水量變化具有特殊性，確定其最大涌水量對工作面安全開采具有重要意義，通過對該條件下涌水規律及頂板垮落特點進行分析，提出了最大涌水量的預計方法；梁燕等[3]對第三系弱膠結砂巖地板突水潰砂進行室內試驗，得到了不同試樣突水的臨界水頭梯度并分析了突水機理；湯愛萍等[4]利用設計的試驗裝置對某礦井的突水潰砂的機理研究，確定了與突水潰砂相關的因素，并進一步提出了防治此災害的方法；張敏江等[5]通過試驗研究了3種弱膠結的砂巖突水、涌砂的階段性的特點；隋旺華，蔡光桃，董青紅等[6-8]通過模型試驗得出含水層的初始水頭和突砂口張開程度是控制礦井工作面突砂量的關鍵因素；此外，楊偉峰[9]也通過室內試驗探討了水砂混合流在運移過程中孔隙水壓力變化的規律。

綜上所述，大量學者基于理論溫習和室內模型試驗研究了突水潰砂災害過程，得到了一定的相關控制因素和水砂混合流動規律。鑒于控制研究成本和發展數值計算，因此文中通過二維顆粒流軟件PFC2D分別對上覆松散層發生滲透破壞、水砂混合流運移規律進行模擬研究，探究突水潰砂過程中松散層滲透變形破壞以及水砂混合流在裂隙通道內的運移規律，建立淺埋深薄基巖煤層開采突水潰砂的數值模擬模型。

1 流體與顆粒相互作用的理論基礎

1.1 顆粒運動方程

淺埋深煤層的上覆松散層可以看作是由砂粒和水組成的兩相體，因此對于淺埋深煤層開采突水潰砂的研究，要分別確定顆粒、流體的控制方程。顆粒流研究的離散單元法是基于這樣一個思想：假定單個時間步內，顆粒的運動只對直接相鄰的顆粒產生影響，而不會傳播給其他不相鄰的顆粒。因此顆粒的運動方程可以描述為

（1）

（2）

1.2 流體控制方程

根據連續體的質量守恒方程和流體動量守恒方程

（3）

（4）

（5）

結合方程（3）（4）（5），得到不可壓縮流體的局部平均的納維斯托克斯方程

（6）

1.3 流體與顆粒耦合流動

淺埋深煤層開采突水潰砂的顆粒流研究，關鍵是確定流體與顆粒耦合流動過程中的相互作用力，包括：拖拽力、浮力、虛擬質量力、巴塞特力、壓力梯度等。考慮文中研究內容，這里只選取拖拽力、浮力作為流體對顆粒的主要作用力。在PFC中流體與顆粒流動的耦合計算，實質上是一系列數據的交換，包括孔隙率、拖拽力、流體速度等。

在顆粒流程序PFC中假定，同一個流體單元內拖拽力、浮力平均的施加到每一個顆粒的中心，其中拖拽力的表達式如下

（7）

式（7）中拖拽力系數Cd，經驗系數χ分別為[10-11]

（8）

（9）

式（8）（9）中Rep為雷諾數，雷諾數的表達式如下

（10）

因此，在流體與顆粒耦合流動的計算過程中，流體對于顆粒施加的力為拖拽力與浮力的合力，可以表示為

（11）

2 上覆松散層滲透破壞的PFC2D模擬

2.1 顆粒流程序PFC2D

二維顆粒流軟件PFC2D是通過離散元法來模擬二維圓盤顆粒介質的運動及其相互作用，通過利用連續介質的方法求解復雜變形的真實問題。在PFC中，顆粒間的相互作用被視為一個動態過程，顆粒間的接觸力和位移是通過跟蹤單個顆粒的運動得到的。在計算過程中，采用時步算法在每個顆粒上反復使用運動方程，在每一個接觸上反復使用力-位移方程，并持續更新墻體的位置[12]。

由于PFC潛在的高效率、塊體可破裂、對模擬的位移大小沒有限制等優點，PFC逐漸成為模擬固體力學和顆粒流問題最有效的手段之一。

2.2 材料微觀物理力學參數的標定

利用PFC2D進行數值模擬需要通過校核過程選擇合適的微觀參數，主要是通過數值模擬的響應與實際材料的測試結果進行對比[13]。文中選取粗砂作為研究對象，并通過對粗砂的級配和微觀參數的標定，來使其符合粗砂的宏觀物理性質。通過虛擬雙軸試驗來匹配參數，由于砂土材料的特殊性質，采用循環加載的方式，在模型到峰值強度80%后，重新施加荷載。圖1為循環加載下的軸向偏應力隨軸向應變變化的關系曲線，從圖中可以看出材料明顯的塑性特征，近似符合實際的粗砂偏應力與軸向應變曲線。

圖1 偏應力隨軸向應變變化的關系曲線Fig.1 Curve of deviatoric stress change vs.axial strain

通過上述虛擬雙軸試驗，文中獲得了接近粗砂實際宏觀物理力學性質的顆粒比重、粒徑范圍、摩擦系數等一系列微觀參數值以及PFC2D模擬中粗砂的級配曲線，具體微觀參數數值與顆粒級配曲線分別見表1，如圖2所示。

圖2 PFC中顆粒的級配曲線Fig.2 Particle percentage curve in PFC

表1 PFC模擬中選用顆粒的微觀參數

2.3 PFC2D數值模擬模型

文中在PFC2D中建立了數值模擬模型，首先建立四組墻單元，使其構成長寬為0.2 m×0.15 m矩形結構，在墻單元圍成的矩形內按粗砂級配生成2 806個顆粒，并使其在重力下達到穩定狀態，以此模擬淺埋煤層開采的上覆松散層。此外，在此矩形空間內創建了300個大小為10 mm×10 mm的流體單元，以此模擬承壓含水層。（其中水的密度為1 000 kg/m3，粘度為1.0×10-3Pa·s）

這里只討論粗砂在水壓力作用下產生突水潰砂的過程，重點研究在突水潰砂初期，不同裂縫寬度、不同初始水壓下砂體顆粒之間的接觸應力分布與水流的流速分布。文中在模型的底部設置不同寬度的3組裂縫（見表2），定義了5組不同的初始水壓力（見表3），根據選取的底部裂縫和初始水壓力的個數，相互組合建立了15組數值模擬模型。

表2 底部裂縫的寬度值

表3 初始水壓力的值

2.4 結果分析

根據模擬結果可以看出隨著潰砂的發生，模型內部顆粒間最大接觸應力隨著砂粒穩定流出后逐漸趨于平穩，流體在底部裂隙處的最大流速在一定范圍內上下波動。下面分別取部分不同裂隙寬度、不同初始水壓下模型的接觸應力與流速分布圖，結果如圖3，4，5所示。

圖3 裂隙為20 mm時，不同水壓力（a）～（c）0.1 MPa，0.3 MPa和0.5 MPa下顆粒接觸應力與流體流速分布Fig.3 Width is 20 mm，distribution of particles contact stress and fluid velocity under different water pressure（a）～（c）0.1，0.3 and 0.5 MPa，respectively

圖4 裂隙為25 mm時，不同水壓力（a）～（c）0.1 MPa，0.3 MPa和0.5 MPa下顆粒接觸應力與流體流速分布Fig.4 Width is 25 mm，distribution of particles contact stress and fluid velocity under different water pressure（a）～（c）0.1，0.3 and 0.5 MPa，respectively

圖5 裂隙為30 mm時，不同水壓力（a）～（c）0.1，0.3和0.5 MPa下顆粒接觸應力與流體流速分布Fig.5 Width is 30 mm，distribution of particles contact stress and fluid velocity under different water pressure

從圖5可以看出，相同裂隙寬度的模型內，顆粒接觸應力隨著水壓力的增大而逐漸增大，接觸應力成枝狀分布形成不穩定的力鏈。模型底部的接觸應力伴隨著潰砂現象的發生，分布網絡的形態沒有較大的變化。隨著潰砂開始，顆粒間接觸應力逐漸趨于平穩，底部裂隙處的流速成上下波動狀態。模型在潰砂發生后25萬步，不同裂隙寬度下粒間接觸應力的均值與裂隙處流體流速的均值如圖6，圖7所示。

圖6 不同裂隙下初始水壓力與顆粒間接觸應力曲線Fig.6 Initial water pressure under different fractures between particles and the contact stress curves

圖7 不同裂隙下初始水壓力與裂隙處流體流速曲線Fig.7 Initial water pressure under different fractures and fracture fluid velocity curve

通過簡化實際復雜的模型，在PFC2D內建立簡化數值模擬模型進行研究，得出了上覆松散層在不同初始水壓的力、不同的導水裂縫寬度下滲透變形破壞時，砂體內部接觸應力以及流體流速的規律。從圖6，7可知，松散含水層的初始水壓力由0.1 MPa增加到0.5 MPa時，不同裂隙寬度下顆粒間接觸應力均增加5倍左右，同時裂隙處流體流速均增加2.5倍左右，且相同水壓力條件下裂隙寬度對于顆粒間接觸應力與裂隙處流速的甚微；顆粒間接觸應力與裂隙處流體流速隨含水層初始水壓力增加非線性的增加，且不同裂隙寬度下的增加趨勢一致。

從模擬結果分析可得，當淺埋煤層開采過程中發生突水潰砂，相同的底部裂隙寬度下水砂突涌到開采面的速度隨含水層水壓力的增加而增大，水砂混合物快速突涌到工作面，形成較大的瞬間動力，對于巷道支護系統及液壓支架的沖擊更猛烈，更容易造成壓架事故。同時，水砂突涌速度的增加也提高了應對和處理突水潰砂事故的難度。

3 水砂混合流在裂隙通道內運移的PFC2D模擬

3.1 裂隙通道的設計及模擬方案

當上覆松散含水層發生滲透破壞以后，水砂混合流通過由冒落帶和裂隙帶形成的裂隙通道，突涌到采煤工作面，裂隙通道的形態對于水砂混合物的流速、運移規律有直接影響。因此對水砂混合流在裂隙通道內的運移規律進行研究，是探究煤層開采突水潰砂規律的重要內容之一。考慮到實際的地質情況下，裂隙通道形態復雜多樣，無論進行室內試驗還是數值模擬，都難以完全復制裂隙形態。因此文中選取理想單裂隙作為研究對象，選取影響突水潰砂的幾個主要因素進行研究，即裂隙的張開度、粗糙度以及傾角，進行數值模擬研究。其中裂隙通道的粗糙度用粗糙度系數JRC（Joint Roughness Coefficient）表示，裂隙通道的傾角以水平方向為基準表示，具體模型設計見表4.

表4 裂隙通道類型組合

首先，在模型上部裂隙通道內填充砂體并使其在自身重力下達到平衡。然后，刪除隔墻單元并施加流場，流場的網格沿著上下2部分裂隙通道分布。由于標準輪廓曲線是位圖格式的，而PFC2D中墻單元是矢量形式的，所以需要進行矢量化的處理。文中選用軟件Illustrator對其進行矢量化處理，后導入AutoCAD來達到曲線簡化并提取曲線坐標的目的。通過PFC2D中的Fish語言和提取的曲線坐標，編寫出模型裂隙通道的形態，圖8為張開度為20 mm，JRC值為11，夾角為60°的模擬裂隙通道。

圖8 模擬裂隙通道示意圖Fig.8 Schematic of fracture channel

3.2 結果分析

文中主要監測模型下部裂隙通道初始位置為（0，0.2）處、初始位置為（0，0.15）處以及初始位置為（0，0.1）3個監測點處的顆粒突涌出裂隙瞬間的Y方向流速與運算時步的關系（圖中黑線、紅線、藍線分別代表通過入口、中部、出口處的顆粒）。為了方便，把粗糙度與夾角簡寫，例如：粗糙度1，夾角45°的裂隙通道簡寫成（1，45）以此類推。監測結果如圖9所示，

圖9 3個監測顆粒Y方向的速度與運算時步的曲線Fig.9 Graph of the speed of three monitoring particles Y direction and operation steps

顆粒的Y方向流速與運算時步的關系監測說明，隨著裂隙粗糙度的增加，同一位置顆粒突涌出裂隙通道的時間也隨之增加，顆粒Y方向速度隨之而減小；隨著裂隙角度的增加相同位置顆粒突涌出裂隙通道的時間變短，且同一位置顆粒Y方向的速度隨之增加。

此外，當裂隙的傾角和粗糙度發生變化時，速度大小的波動頻率明顯發生改變。裂隙的粗糙度越大、傾角越小時，顆粒Y方向的速度的波動頻率越大；而裂隙的粗糙度越小、傾角越大時，則相反。由此，可以預測出水砂混合流在裂隙通道內可能存在2種流動形態。

1）快速突出型，即水砂混合物受裂隙影響較小，水砂混合流以較快的速度、較短的時間突涌到工作面，這種情況對于巷道支護的破壞嚴重，典型的情況為粗糙度為1，夾角90°的裂隙通道（圖12（b））；

2）波動突出型，即水砂混合物受裂隙影響，在運移過程中受裂隙通道的阻力影響，突出速度和突出時間明顯減弱，水砂混合流的瞬間動力相比快速突出型明顯減小，典型的情況為粗糙度為19，夾角45°的裂隙通道（圖12（e））。

4 結 論

文中利用二維離散元軟件PFC2D，通過2個方面對淺埋煤層開采過程中發生突水潰砂地質災害問題的規律進行數值模擬，對其滲透變形破壞以及水砂混合流在裂隙通道內的運移規律，可通過對相關量的監測來進行分析研究。文中通過對淺埋煤層突水潰砂的模擬，得到以下結論

1）上覆松散含水層的初始水壓力以及裂隙通道的寬度影響上覆松散層的滲透變形破壞。當含水層初始水壓力增加時，裂隙處流體流速、顆粒間的接觸應力均明顯增加，且呈現出非線性增長，一旦裂隙通道導通，就會產生突水潰砂災害，初始水壓力以及裂隙寬度的增加會使突水潰砂愈劇烈；

2）通過對水砂混合流在裂隙通道內的運移進行數值模擬，可發現裂隙通道的傾角以及粗糙度影響水砂混合流的突出速度以及突出時間，同時通過對顆粒Y方向速度的監測，預測了2種可能的水砂突出類型，即快速突出型、波動突出型；

3）通過對淺埋深煤層開采突水潰砂的模擬分析可知，含水層的初始水壓力影響突水潰砂的發生及運移。由此，提出可以通過疏放底層水來降低含水層初始水壓力及松散層的含水率，從而降低突水潰砂災害影響和預防突水潰砂災害發生的重要手段之一。

References

[1] 許學漢，王 杰．煤礦突水預測預報研究[M]．北京：地質出版社，1992．

XU Xue-han，WANG Jie.Research on coal mine water inrush prediction[M].Beijing:Geological Press，1992.

[2] 張世凱.厚含水松散層薄基巖開采涌水規律及預計[J].河北煤炭，1999（2）：1-2.

ZHANG Shi-kai.Thin bedrock mining thick watered loose bed of water gushing and expected[J].Hebei Coal，1999（2）：1-2.

[3] 梁 燕，譚周地，李廣杰．弱膠結砂層突水、涌砂模擬試驗研究[J].西安公路交通大學學報，1996，16（1）：19-22．

LIANG Yan，TAN Zhou-di，LI Guang-jie.Simulation test research on water and soil outbursts of weak binding soil[J].Journal of Xi’an Highway Communication University，1996，16（1）：19-22.

[4] 湯愛平，董 瑩，譚周地，等.振動作用下礦井突水涌砂機理的研究[J].地震工程與工程震動，1999，19（2）：132-135．

TANG Ai-ping，DONG Ying，TAN Zhou-di，et al.Mechanism of sandy-silt seepage-deformation in mine under vibration[J].Earthquake Engineering and Engineering Vibration，1999，19（2）：132-135.

[5] 張敏江，張麗萍，姜秀萍．弱膠結砂層突涌機理及預測研究[J]．金屬礦山，2002，16（10）：48-50．

ZHANG Min-jiang，ZHANG Li-ping，JIANG Xiu-ping.Study on the inrushing mechanism of weak cemented quicksand layer and its forecasting[J].Metal Mine，2002，16（10）：48-50.

[6] 隋旺華，蔡光桃，董青紅．近松散層采煤覆巖采動裂縫水砂突涌臨界水力坡度試驗[J].巖石力學與工程學報，2007，26（10）：2 084-2 091.

SUI Wang-hua，CAI Guang-tao，DONG Qing-hong.Experimental research on critical percolation gradient of quicksand across overburden fissures due to coal mining near unconsolidated soil layers[J].Chinese Journal of Rock Mechanics and Engineering，2007，26（10）：2 084-2 091.

[7] 隋旺華，董青紅，蔡光桃，等．采掘潰砂機理與預防[M]．北京：地質出版社，2009．

SUI Wang-hua，DONG Qing-hong，CAI Guang-tao，et al.Sand mining collapse mechanism and prevention[M].Beijing：Geological Press，2009.

[8] 隋旺華，董青紅．近松散層開采孔隙水壓力變化及其對水砂突涌的前兆意義[J].巖石力學與工程學報，2008，27（9）：1 908-1 916．

SUI Wang-hua，DONG Qing-hong.Varition of pore water pressure and its precursor significance for quicksand disasters due to mining near unconsolidated formations[J].Chinese Journal of Rock Mechanics and Engineering，2008，27（9）：1 908-1 916.

[9] 楊偉峰.薄基巖采動破斷及其誘發水砂混合流運移特性研究[D].徐州：中國礦業大學，2009.

YANG Wei-feng.Overburden failure in thin bedrock and characteristics of mixed water and sand flow induced by mining[D].Xuzhou：China University of Mining and Technology，2009.

[10] Di Felice R D.The voidage function for fluid-particle interaction systems[J].Int.J.Multiphase Flow，1994,20（1）：153-159.

[11] Xu B H，Yu A B.Numerical simulation of the gas-solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics[J].Chemical Engineering Science，1997，52（16）：2 785-2 809.

[12] 陳 俊，張 東，黃曉明.離散元顆粒流軟件（PFC）在道路工程中的應用[M].北京：人民交通出版社，2015.

CHEN Jun，ZHANG Dong，HUANG Xiao-ming.Discrete element particle flow software（PFC）application in road engineering[M].Beijing:People’s Traffic Press，2015.

[13] Itasca Consulting Group，Inc.PFC3Duser’s manual[R].Minneapolis，USA：Itasca Consulting Group，Inc.2006.

Water-sand inrush simulation under shallow coal mining based on the particle flow code

LIU Ya-qun1，ZHOU Hong-wei1，LI Yi-hu2，YI Hai-yang1，XUE Dong-jie1

（1.SchoolofMechanicsandCivilEngineering，ChinaUniversityofMiningandTechnology，Beijing100083,China； 2.ChinaRailwayFifthSurveyandDesignInstituteGroupCo.，Ltd.，Beijing102600,China）

For the research of the seepage-deformation of aquifer,the flow regularity of mixed water-sand flow,the inrush law of mixed water-sand under the shallow coal seam mining,and providing a basis of guidance to prevent water-sand inrush disasters from happening,the particle flow code in two dimensions （PFC2D） is employed to simulate the seepage deformation failure of loose overburden and the characteristics of mixed water-sand flow,respectively.Establish the numerical simulation model of water-sand inrush under shallow coal seam.Results show that: The intial water pressure of overlying loose aquifer and the width of fracture channels affect the process of seepage-deformation failure.When the intial water pressure of aquifer increased the speed of fluid in the fracture and the contact stress between particles increase apparently and nonlinearly.The inrush time and speed of mixed water-sand flow have link with the rough and obliquity of fracture channels.Meanwhile two potential types of the mixed water-sand inrush were carried out according to the fluctuation of theYdirction inrush speed.Since the initial water pressure of aquifer affect the happen and flow process of the water-sand inrush,draining away bottom water may give an potential approach to preventing and solving the disaster of water-sand inrush.

shallow buried；water-sand inrush；PFC2D；seepage-deformation；fracture channels

2015-04-12 責任編輯：劉 潔

國家自然科學基金（51134018）;國家重點基礎研究發展計劃（2011CB201201）;中央高校基本科研業務經費（2010YL07）

劉亞群（1990-），男，河北秦皇島人，碩士研究生，E-mail：yqtt1909@yeah.net

10.13800/j.cnki.xakjdxxb.2015.0502

1672-9315（2015）05-0534-07

TD 82

A