基于Comsol 的礦井頂板水突出模擬研究

孫 凱，王 哲，呂 飛

（1.山東省煤炭泰山療養(yǎng)院，山東 泰安271000；2.兗煤萬福能源有限公司，山東 菏澤 274900）

0 引 言

礦井頂板突水問題一直是煤礦安全生產的重大隱患。采煤工作面頂板狀況對于煤礦安全生產具有重大意義。傳統(tǒng)的礦壓理論有：預成裂隙梁、鉸接巖梁、懸臂梁以及砌體梁模型等，運用這些理論研究采場頂板破壞也取得一定成果[1-3]。然而煤礦頂板破壞是一個極為復雜的過程，簡單的理論分析不能精確的分析這一過程。而這些理論無法深入研究頂板水突出時，流體的分布情況。研究表明在出現頂板水突出時，工作面會有大量水涌[4-5]。運用數值模擬的方法對這一過程進行模擬研究對于研究頂板突出問題具有一定的意義。前人在這方面做過大量的研究[6-9]。

Comsol Multiphysics 是一款耦合模擬仿真軟件，由于去具有開放界面，可結合MATLAB 對求解復雜多物理場模型，模擬出整個物理過程[10-11]。本文運用Comsol Multiphysic 中的PDE 模塊進行模型鑲嵌，結合MATLAB 編程，對復雜變量進行設定。模擬出礦井頂板突水過程中水分布情況已經整體變形。并將模擬結果和已做的實驗做對比，二者的位移及滲流規(guī)律基本一致。

1 耦合模型控制方程

礦井水突出是一個復雜的動力學過程。突水過程中，巖層斷裂。本文模型基于滲流—損傷力學。

1.1 應力場方程

在受力初始階段，彈性多孔介質孔隙流體壓力p及位移的本構方程為:

式中：G 為介質的剪切模量；v 為介質的泊松比；Fi和ui（i=x,y,z）分別為體力在位移i 方向的分量。

圖1 單軸應力狀態(tài)下的巖石的本構方程

設，，當介質的應力狀態(tài)滿足最大拉應力準則和摩爾庫倫準則時，其分別發(fā)生拉伸損傷和剪切損傷：

1.2 滲流場方程

由流體質量守恒方程和Darcy 定律得：

式中：εv為體積應變，I 為介質的孔隙率；βl為孔隙流體的體積模量（Pa）；t 為時間（s）；k 為連續(xù)介質固有的滲透率（m2）；μl為流體的動力粘性系數（Pa.s），ρl為流體密度（kg/m2）；g 為重力加速度（m/s2）。

2 頂板突水過程數值模擬

以上數學模型在時間和空間域都是非線性的，因此運用有限元軟件進行數值求解。本文應用運用MATLAB 強大的編程功能，聯(lián)合Comsol Multiphysics對以上模型進行求解。

2.1 幾何模型及參數

設計的幾何模型尺寸為：160m（長）×190m（高），如圖2 所示：

圖2 幾何模型

表1 模型巖石力學參數

表2 模型耦合參數

實驗過程中整個巖層左邊、右邊及底面均受固定約束。對頂面進行加壓，分析在壓力變化過程中煤層頂板突水的情況。

2.2 數值模擬結果

模擬得出整個巖層加壓至8MPa 時，開始出現頂板斷裂，有大量水突出，整個模型流體速度場分布如圖3、圖4 所示。

圖3 水速度場分布

圖4 頂板水速度場曲線圖

由圖3 得，采空區(qū)出現頂板跨落時，工作面的水速度場分布集中，表明在此位置有大量水涌出。頂板水流速度分布如圖4 所示。其模擬結果和前人的研究基本吻合。

圖5、圖6 分別為2MPa 和8MPa 外力作用下整個模型的應變。

圖5 2MPa 下應變

圖6 8MPa 下應變

由圖5、圖6 可以看出，在工作面推進到一定程度，直接暴露于采空區(qū)的上覆巖層不再受到下部煤層的支撐作用，在重力作用下發(fā)生彎曲、沉降、離層。煤層頂板覆巖變形破壞存在一個漸進的演變過程，是伴隨工作面向前推進，采空區(qū)范圍的不斷擴大，應力的不斷調整而發(fā)展變化的。當外界壓力達到一定值時，出現頂板斷裂和局部彎曲下沉。

圖7 垂直位移分布

圖8 頂板最終垮塌情況

圖7 為頂板跨落時，整個模型豎直方向位移。由于受工作面回采的影響，工作面豎直方向位移變化較大，位移變化和理論研究基本一致。工作面前方巖層下沉梯度明顯大于工作面后方巖層下沉梯度，巖層下沉呈非對稱分布。

實際實驗結果如圖8 所示，模擬研究同實驗研究結果基本吻合。

3 結 論

本文根據流固耦合相關理論，分析得出煤層頂板突水流固耦合模型。運用Comsol Multiphysics 結合MATLAB 進行數值模擬，得出頂板突水時，煤層、巖層模型體積應變、豎直方向上的位移及水速度場分布。模型的模擬變形同實際實驗結果基本吻合，驗證了數學模型的正確性。