基于THM耦合凍結煤層溫度場演化規律研究*

劉 軍，衛彥昭，王 蔚，張露偉

(1.煤炭安全生產與清潔高效利用省部共建協同創新中心，河南 焦作 454003；2.河南理工大學 安全科學與工程學院，河南 焦作 454003；3.河南理工大學 應急管理學院，河南 焦作 454003)

0 引言

煤與瓦斯突出是制約煤礦安全生產的主要因素之一。由于地應力、瓦斯壓力及煤體自身力學性質等因素，使得石門揭煤[1]時煤與瓦斯突出十分危險且強度較大，盡管采取防突措施，瓦斯突出事故仍時有發生。文獻[2-3]研究發現低溫冷凍煤體可以抑制瓦斯解吸，降低煤體瓦斯壓力，從而降低煤與瓦斯突出危險性。

近年來，部分學者對石門揭煤冷凍防突開展研究：馮濤等[4]、謝雄剛等[5-7]對石門揭煤突出煤層的凍結溫度場進行實驗研究和數值模擬，采用低溫鹽水管道循環的制冷方式對煤體進行凍結，發現凍結區形成以凍結管中心為圓心的凍煤圓柱；葉青等[8]、任明龍[9]、董若蔚[10]、周震等[11]通過分析液氮凍結對煤體的力學性質影響，發現凍結后煤體的彈性模量、抗壓強度和抗拉強度得到顯著提高，可有效降低石門揭煤突出危險；王賓賓[12]發現在低溫作用下，含水率對煤體強度有顯著影響，煤體內水凍結成冰，可有效提高煤體抗拉強度；陶云奇[13-14]、黃旭超[15]基于溫度對瓦斯含量的影響，建立含瓦斯煤流-固-熱三場耦合數學模型，發現在瓦斯抽采過程中，考慮溫度變化是必要的。

現有研究均為基于常溫瓦斯抽采或注熱促抽技術建立的三場耦合方程。煤是多孔介質，溫度場改變會引起瓦斯賦存狀態改變或煤體變形，瓦斯賦存狀態改變會引起瓦斯在煤層中的熱對流，進而導致瓦斯流動、煤層孔隙率和滲透率的改變。同時，煤層變形會引起溫度場及滲流場的改變，含瓦斯煤低溫冷凍過程中溫度場演化規律的熱流固(THM)耦合需要進一步研究。因此，本文通過建立應力場、滲流場、多孔介質傳熱場耦合模型，利用COMSOL軟件模擬煤層冷凍過程中的溫度演化規律。

1 THM耦合方程構建

煤體凍結過程是應力場、滲流場、溫度場相互耦合作用的結果，首先需要構建熱流固THM耦合方程。

1.1 應力場方程

煤體凍結過程中，煤基質應變由地應力、瓦斯吸附及低溫凍結3者相互耦合作用。由于水結冰后體積膨脹9%，模型選取地點山西晉城玉溪煤礦煤層原始含水率為1.14%，含水率低。為簡化模型，水分凍結后體積膨脹暫未考慮。應力場方程如式(1)所示[16]：

(1)

1.2 滲流場方程

瓦斯在煤層中的流動符合達西定律，如式(2)所示：

(2)

式中：u為煤層瓦斯流動速度，m/s；κ為煤層滲透率，mD；μ為瓦斯動力黏度，Pa·s；P為煤層瓦斯壓力，MPa。

煤巖體孔隙率θ如式(3)所示：

(3)

式中：θ0為煤體初始滲透率，%，本文取6.96 %；εP為吸附應變；e為煤體體積形變。

煤巖體滲透率κ如式(4)所示：

(4)

式中：κ0為煤層初始滲透率，mD，本文取3.81×10-7mD。

連續性方程如式(5)所示：

(5)

式中：ρ為煤層中瓦斯壓力為P時的密度，kg/m3；W為單位質量煤的瓦斯質量，kg/kg；t為時間，s。

氣體狀態方程如式(6)所示：

(6)

式中：ρn為煤層中瓦斯壓力為1個大氣壓(Pn)時的密度，0.717 kg/m3；z為煤層中瓦斯可壓縮系數，取1。

單位質量含瓦斯煤的瓦斯質量如式(7)所示：

(7)

通過式(2)～(7)聯立得到煤體滲流場方程，如式(8)所示：

(8)

式中：Ks為煤層骨架模量，Pa；m為滑脫因子，Pa，取7.6×105Pa。

1.3 溫度場方程

由于煤中含有瓦斯，在含瓦斯煤冷凍過程中，不僅有煤基質進行熱傳導，瓦斯也在進行熱對流，因此溫度場方程采用多孔介質傳熱方程，如式(9)所示：

(9)

式中：等效比熱容CPeff=θCP，CH4+(1-θ)CP，coal，kJ/(kg·K)；等效導熱系數Keff=θKCH4+(1-θ)Kcoal,W/(m·K)；CP，CH4為甲烷比熱容，kJ/(kg·K)；CP，coal為煤體比熱容，kJ/(kg·K)；KCH4為煤導熱系數，W/(m·K)；Kcoal為煤導熱系數，W/(m·K)；?T為傳熱過程煤體溫差，℃；μ為瓦斯對流擴散速度，m/s。

應變內能如式(10)所示：

(10)

本文模擬忽略瓦斯吸附釋放的熱量，所以QCH4=0。將式(1)、式(8)和式(9)耦合，得到煤層凍結過程中的THM耦合模型。

2 凍結過程滲透率演化數值模擬

該模型以山西晉城玉溪煤礦3號煤層為研究對象，模型寬度、高度為3 m×3 m。CaCl2冰點比水低，且溶液配料方便，有良好的傳熱性能，根據凍結溫度合理配置水溶液濃度。凍結溫度-30 ℃需CaCl2溶液濃度為28%，因此采用氯化鈣(CaCl2)溶液作為冷凍介質，凍結孔直徑94 mm。應力場邊界條件為模型上下左右邊界位移為0，上邊界考慮煤層受上覆巖層重力作用，載荷為12.5 MPa；滲流場邊界條件為煤層頂底板為巖層，0通量，工作面初始瓦斯壓力1 MPa；溫度場邊界條件為煤層初始溫度為28 ℃，凍結孔邊界溫度-30 ℃。模型物理參數見表1。

表1 模型物理參數Table 1 Physical parameters of model

采用COMSOL二維平面模型，將式(1)、式(8)和式(9)分別使用固體力學模塊、PDE模塊、多孔介質傳熱模塊耦合進行數值模擬，得到煤體滲透率隨凍結時間變化規律，如圖1所示。由圖1可知，從凍結孔到模型邊界，滲透率逐漸降低，最后趨于穩定。由于地應力作用，凍結孔周圍產生大量裂隙，瓦斯滲透率高。隨距凍結孔距離增加，煤層受鉆孔應力影響逐漸減小，滲透率逐漸減小，最后趨于穩定，但由于受凍結溫度作用，煤基質收縮，穩定值均大于煤層原始滲透率3.81×10-7mD。以凍結50 d為例，凍結孔邊界處滲透率為5.73×10-7mD，凍結孔在地應力及溫度應力作用下產生裂隙，滲透率相對最高。從凍結孔邊界處到距離凍結孔0.3 m，滲透率快速減小至4.54×10-7mD，隨距凍結孔距離增加，凍結孔周圍裂隙逐漸減少，滲透率變化快。由距離凍結孔0.3 m處至煤層邊界，滲透率緩慢減小，至煤層邊界處滲透率為3.92×10-7mD，隨距凍結孔距離增加，凍結作用逐漸減小，煤體形變減小，滲透率變化緩慢，但由于煤層凍結50 d，邊界處溫度低于煤層初始溫度28 ℃，煤基質收縮，滲透率大于煤層原始滲透率3.81×10-7mD。

圖1 不同凍結時間滲透率變化規律Fig.1 Change law of permeability under different freezing time

沿凍結孔0.3 m處作直線，如圖1所示，沿直線提取不同凍結時間滲透率，得到凍結孔0.3 m處滲透率隨凍結時間變化，如圖2所示。由圖2可知，煤層滲透率隨凍結時間增加，先快速增大后緩慢增大。凍結10 d，滲透率由3.81×10-7mD快速增加至4.32×10-7mD，煤層初始溫度為28 ℃，凍結溫度-30 ℃，溫差較大，溫度場演化速度快，凍結作用強，煤基質快速收縮，滲透率快速增高；凍結10～50 d，滲透率緩慢增加至4.56×10-7mD，煤層經低溫凍結作用使其溫度降低，溫差逐漸減小，凍結作用逐漸減弱，煤基質收縮緩慢，滲透率緩慢增高。

圖2 距凍結孔0.3 m處滲透率隨凍結時間變化Fig.2 Change of permeability at 0.3 m away from freezing hole with time

3 凍結孔布置方式凍結效果考察

凍結孔布置方式分為單孔凍結、三花眼凍結、四花眼凍結等形式，分別對其凍結效果進行研究分析。

3.1 單孔凍結效果

煤層冷凍有效影響半徑是指在冷凍過程中，通過冷量從凍結孔向四周傳遞，使煤層降溫至0 ℃以下。

應用COMSOL軟件建立單孔凍結物理模型，寬度、高度為3 m×3 m，在模型中心處布置凍結孔，直徑94 mm，如圖3所示。

圖3 單孔凍結物理模型Fig.3 Physical model of single-hole freezing

煤層單孔凍結效果如圖4所示，隨凍結時間增加，煤層凍結區域逐漸增大。凍結初期，煤層凍結區域快速增加，隨凍結過程進行，凍結區域緩慢增大。

圖4 單孔凍結溫度場分布Fig.4 Temperature field distribution of single-hole freezing

對煤層凍結孔兩側連線(圖3所示)凍結過程進行數值模擬，以連線左端點為起始點，得到不同凍結時間煤層溫度隨距離變化，如圖5所示。由圖5可知，隨凍結時間增加，煤層有效凍結區域逐漸增大，距凍結孔越遠，煤層溫度越高。以凍結10 d為例，凍結孔邊界至距凍結孔0.5 m處，距離凍結孔遠，凍結作用較小，煤基質溫度演化緩慢，滲透率變化緩慢，瓦斯對流傳熱緩慢，因此溫度由28 ℃降低至20.71 ℃；凍結孔邊界至距凍結孔0.12 m處，凍結作用大，煤基質熱傳導快，滲透率快速增大，瓦斯對流傳熱快，溫度快速減少至有效凍結溫度0 ℃。

圖5 不同凍結時間溫度隨距離變化Fig.5 Change of temperature with distance under different time

在圖5中沿凍結溫度孔0 ℃作1條直線，得到有效凍結半徑隨凍結時間變化如圖6所示。由圖6可知，有效凍結半徑隨凍結時間先快速增加，后緩慢增加。凍結0～10 d，有效凍結半徑快速增加至0.12 m，這是由于煤層初始溫度為28 ℃，凍結溫度-30 ℃，溫差較大，煤基質傳熱速度快，凍結作用強，煤基質快速收縮，滲透率快速增高，瓦斯對流速度快，瓦斯傳熱快；凍結10～60 d，有效凍結半徑緩慢增加至0.26 m，由于此階段煤層經低溫凍結作用導致溫度降低，溫差逐漸減小，煤基質傳熱緩慢，凍結作用逐漸減小，煤基質收縮緩慢，滲透率緩慢增高，瓦斯對流速度減緩，瓦斯在煤層中傳熱緩慢。有效凍結半徑與時間呈冪函數關系，其關系表達式為：

圖6 有效凍結半徑隨凍結時間變化Fig.6 Change of effective freezing radius with time

r=0.056 5t0.336(R2=0.999 51)

式中：r為有效凍結半徑，m；t為凍結時間，d。

3.2 三花眼(正三邊形)布置凍結效果

應用COMSOL軟件，建立三花眼(正三邊形)凍結物理模型，布置凍結孔，直徑94 mm，設定孔間距1 m，進行凍結效果分析。

煤層三眼花凍結物理模型如圖7所示，其凍結效果如圖8所示，隨凍結時間增加，煤層凍結區域逐漸增大。凍結孔之間相互影響，凍結25 d，孔間煤體全部凍結，即凍結半徑0.5 m，與單孔凍結25 d凍結半徑0.17 m相比，凍結半徑增大1.9倍。

圖7 三花眼凍結物理模型Fig.7 Physical model of three-flowered eye freezing

圖8 三花眼凍結溫度場分布Fig.8 Temperature field distribution of three-flowered eye freezing

對煤層凍結孔連線(圖7)凍結過程進行數值模擬，以連線左端點為起始點，得到三花眼布置孔間連線溫度變化，如圖9所示。由圖9可知，隨凍結時間增加，有效凍結范圍增大。凍結0～10 d，有效凍結半徑快速增加至0.14 m，孔間連線中心處溫度快速降至12.09 ℃，由于煤層初始溫度為28 ℃，凍結溫度-30 ℃，溫差較大，煤基質傳熱速度快，在凍結孔相互作用下，煤基質快速收縮，滲透率快速增高，瓦斯對流速度快，瓦斯在煤基質傳熱快；凍結10～25 d，有效凍結半徑緩慢增加至0.5 m，孔間連線中心處溫度緩慢降至-1.04 ℃。由于此階段凍結孔相互作用，溫差逐漸減小，煤基質傳熱緩慢，煤層傳熱緩慢，煤基質收縮緩慢，滲透率緩慢增高，瓦斯對流速度減緩，瓦斯在煤層中傳熱緩慢，直至煤體完全達到有效凍結溫度。

圖9 三花眼布置孔間連線溫度變化Fig.9 Temperature change of connection between holes in three-flowered eye arrangement

同一凍結時間，距凍結孔越遠，溫度變化越緩慢。以凍結25 d為例分析，凍結孔邊界至孔間連線中心處，隨距離增大，溫度由-30 ℃增大至-1.04 ℃，溫升速率先增大后減小，由于凍結孔邊界在地應力及溫度應力作用下產生裂隙，滲透率高，瓦斯對流速度快，降溫效果好，隨距離增加，裂隙減少，瓦斯對流速度減慢，溫度變化速度減慢。

3.3 四花眼(正四邊形)布置凍結效果

建立四花眼(正四邊形)凍結物理模型，布置凍結孔，直徑94 mm，設定孔間距1 m，對角線間距1.41 m，進行凍結效果考察，如圖10所示。

圖10 四花眼凍結物理模型Fig.10 Physical model of four-flowered eye freezing

煤層四花眼凍結效果如圖11所示，隨凍結時間增加，煤層凍結區域逐漸增大。隨著凍結時間增加，凍結孔之間相互影響，凍結25 d，孔間煤體全部凍結，即凍結半徑0.7 m，與單孔凍結25 d凍結半徑0.17 m相比，凍結半徑增大3.2倍.與三花眼凍結半徑0.5 m相比，凍結半徑增大0.4倍。

圖11 四花眼凍結溫度場分布Fig.11 Temperature field distribution of four-flowered eye freezing

對煤層凍結孔對角線(圖10所示)凍結過程進行數值模擬，以連線左端點為起始點，得到四花眼布置孔間連線溫度變化，如圖12所示。由圖12可知，隨凍結時間增加，有效凍結范圍隨之增大。凍結0～10 d，有效凍結半徑快速增加至0.1 m，孔間連線中心處位置溫度快速降低至15.85 ℃，由于煤層溫差大，煤基質傳熱速度快，在凍結孔相互作用下，滲透率快速增高，瓦斯在煤基質傳熱快；凍結10～25 d，有效凍結半徑緩慢增加至0.7 m，孔間對角線中心處位置溫度緩慢降至-0.81 ℃。此階段溫差逐漸減小，煤基質傳熱緩慢，煤層傳熱緩慢，煤基質收縮緩慢，滲透率緩慢增高，瓦斯對流速度減緩，瓦斯在煤層中傳熱緩慢，直至煤體完全達到有效凍結溫度。

圖12 四花眼布置對角線溫度變化Fig.12 Diagonal temperature change of four-flowered eye arrangement

同一凍結時間，距凍結孔越遠，溫度變化越緩慢。以凍結25 d為例，凍結孔邊界至孔對角線中心處，隨距離增大，溫度由-30 ℃增大至-0.81 ℃，溫升速率先增大后減小，由于凍結孔邊界在地應力及溫度應力作用下產生裂隙，滲透率高，瓦斯對流速度快，降溫效果好；隨距離增加，裂隙減少，瓦斯對流速度減慢，溫度變化速度減慢。

4 結論

1)煤層冷凍過程中，滲透率呈先快速增大后緩慢增大趨勢，有效凍結半徑隨時間增加，先快速增加后緩慢增加，2者呈冪函數關系為r=0.056 5t0.336。

2)三花眼布置凍結孔，25 d凍結完成，有效凍結半徑增大至0.5 m，與單孔凍結相比，凍結半徑增大1.9倍。孔間連線中心處溫度降低至-1.04 ℃。

3)四花眼布置凍結孔時，25 d凍結完成，有效凍結半徑增大至0.7 m，與單孔凍結相比，凍結半徑增大3.2倍；與三花眼凍結相比，凍結半徑增大0.4倍;孔間對角線中心處處溫度降低至-0.81 ℃。

4)多孔布置凍結孔時，隨凍結時間增加，有效凍結范圍增大；同一凍結時間，距凍結孔越遠，溫度變化越緩慢，對凍結煤層溫度場演化規律研究有指導意義。