熱－流－固耦合巖體三軸壓縮實驗數值模擬

2020-06-08 09:46:34楊長德毛金峰李金波張海東

煤礦安全 2020年5期

楊長德，王鵬，2，毛金峰，李金波，張海東

（1.新疆工程學院礦業工程與地質學院，新疆烏魯木齊830091；2.中國礦業大學（北京）能源與礦業學院，北京100083；

3.山西工程技術學院土木與建筑工程系,山西陽泉045000）

地下煤巖體在滲流力、地應力、開采擾動及溫度等各種不同組合的作用下，其內部結構與內部的物質含量會發生相應的改變，進而會導致其力學性能和滲透特性發生改變，這對煤炭的安全開采或是對礦井的穩定及維護都會產生重大的影響[1-2]。目前國內外學者對多場耦合下的煤巖體進行了大量的實驗研究和數值模擬研究，并取得了較為豐富的研究成果。在多場耦合實驗研究方面，謝和平等[3]采用煤巖采動卸壓增透試驗平臺研究了不同開采方式擾動下煤巖應力場-裂隙場-滲流場的差異性行為，對指導煤與瓦斯共采提供了一定的理論依據；江宗斌[4]等進行軸壓循環加卸載條件下的巖石應力-滲流耦合試驗，分析了巖石加卸載過程中的滲透性變化以及力學行為變化規律；楊金保[5]等開展了單裂隙花崗巖不同圍壓加、卸載和不同水力梯度作用下的多場耦合試驗，得出了應力歷史是影響裂隙巖體滲透特性的因素之一；陳天宇[6]等對含氣頁巖試件進行圍壓循環加卸載，滲透壓固定，無軸壓施加方式，得到了圍壓敏感性和各向異性對含氣頁巖變形以及滲透特性有重大影響；曹亞軍[7]等開展了不同圍壓和滲壓作用下滲流-應力耦合三軸流變試驗。在多場耦合數值模擬方面，趙延林[8]等建立了裂隙巖體溫度場-滲流場-應力場耦合的雙重介質模型，認為在一定時域下，存在起主導控制作用的雙場耦合系統規律；張樹光[9]等通過對裂隙巖體流-熱耦合傳熱過程進行三維數值模擬分析，得出了巖體滲透系數的變化對其內部溫度場的分布影響很大；劉泉聲[10]等通過開展多場耦合作用下巖體裂隙擴展演化規律的研究，得出THM 耦合及巖體變形、失穩全過程的數值模擬算法；康永水等[11]、馮梅梅[12]等、高娟[13]等對裂隙巖體凍融損傷特性以及多場耦合作用過程中巖體的應力-應變特征、裂隙發育規律進行了研究；王軍祥等[14-15]建立數值模型，利用反演的參數對隧道圍巖應力場、滲流場、損傷場分布規律及襯砌結構的受力特征進行了分析。通過以上研究成果可以看出，對煤巖體在多場耦合條件下的研究逐漸在形成一種常態。

隨著煤炭深部開采的日益增多，于“三高一擾動”條件下的深部煤巖體的力學作用機理和破壞失穩規律已不再適應于淺部煤巖體的發展力學規律[16-18]。為此，對溫度場、應力場、滲流場三場耦合條件下的三軸壓縮實驗進行數值模擬研究，進而獲得多場耦合條件下煤巖體的力學行為發展規律和破壞機理。

1 三場耦合下的數學模型

由于煤巖體自身具有不均質性，各向異性以及不可再生性，因此，實驗室中的每1 個煤巖體試件內部成分、均質分布、結構組成以及基本力學參數都不可能完全相同，即實驗室中的每1 個煤巖體試件都是不同的，并無復制性。目前體現煤巖體試件可復制性的研究手段有參數反演進行三維重構或者采用3D 打印技術，和CT 掃描進行三維重構[19]，再或者進行數值模擬研究。采用數值模擬研究建立滲流場-應力場-溫度場耦合的數學模型。

建立的三場耦合模型基于如下假設：巖體屬于均質連續介質體；滲流規律服從Darcy 定律；水的汽化不考慮，巖體全過程處于飽水狀態；熱質在固相、液相介質中傳遞方式以對流傳熱為主，各介質的比熱容及熱傳導系數不隨溫度和壓力而變化。

1）流體滲流控制方程。

式中：k 為滲透系數；μ 為動力黏度；▽p 為出入口速度差形成的滲透差；ρ1為流體密度；Qm為整個滲透過程流體流量。

2）溫度控制方程。

式中：▽T 為巖體傳熱過程形成的溫差；kf為巖體導熱系數；Cp為比熱容；u 為流體流動速度；ρ2為巖體密度；Q 為巖體本身縮傳熱量；Qted為流體吸收熱量。

3）巖體應力控制方程。

式中：S1、S2分別為巖體軸向壓力作用面面積和圍壓作用面面積；Fv1、Fv2分別為軸壓和圍壓；R 為巖體半徑；σz為軸向應力。

4）巖體溫度-應力-滲流耦合方程。

式中：T－ext為巖體點空間溫度；▽Fv為巖體耦合過程中作用力差；S 為對應力作用面積。

聯立以上控制方程，得到巖體三場耦合下的數學模型，結合相應條件，進行模型的求解。

2 網格劃分和邊界條件

選取模型尺寸為25 mm(直徑）×100 mm（高），采用3 節點三棱柱體單元進行網格剖分，共剖分為94 個節點，156 個三角形單元。實體單元網格剖分如圖1。

圖1 實體單元網格剖分圖Fig.1 Mesh generation diagram of entity cells

上邊界由溫度梯度傳熱計算出熱-流邊界條件，下邊界初始溫度為75 ℃，外部溫度為25 ℃，四周為對流傳熱邊界。上邊界為滲流入口，流速0.009 m/s，下邊界為出水口，流速為0.005 m/s。左右兩邊界為封閉狀態，無水流速，整個過程選值分別為1×10-10、1×10-9m2、1×10-8m2的滲透系數進行模擬。右邊界為圍壓4、5、6 MPa 應力載荷邊界，上邊界為軸壓載荷和指定位移邊界，軸壓載荷邊界通過指定位移產生的應力積分函數定義，下邊界為固定約束邊界。計算中所用到的計算參數和數據見表1。

表1 計算模擬參數表Table 1 Calculation of simulation parameters table

3 數值模擬結果與分析

采用有限元軟件對數學模型進行求解[20]，滲流-應力-熱耦合下的巖體試件溫度等溫線分布、表面溫度分布、應力分布以及滲流速度如圖2～圖5。

圖2 巖體溫度等值線圖Fig.2 Temperature contour map of rock mass

圖3 巖體溫度表面圖Fig.3 Temperature surface of rock mass

圖4 應力分布圖Fig.4 Stress distribution diagram

圖5 滲流速度表面圖Fig.5 Surface diagram of seepage velocity

從圖2 和圖3 可以看出，在熱-應力-滲流三場耦合作用下，巖體熱傳遞并非是均勻傳遞的，溫度變化由巖體底部向上部非線性傳遞，因而巖體各局部位置溫度不一，但整體溫降趨勢比較顯著。同時由溫度最大值與最小值可以看出，由于滲流場的介入，整個多場控制下的耦合系統存在熱量散失，這主要是因為滲透液吸收了巖體部分熱量所致。

從圖4 可以看出，在熱-應力-滲流三場耦合控制系統下，巖體應力分布呈現端部應力集中，但其應力分布是呈對稱分布的。這也就表明了在多場耦合條件下，巖體試件實驗存在場耦合作用的顯著影響以及明顯的端部效應現象。

從圖5 可以看出，在多場耦合控制機制下，由于受到應力作用，巖體內部孔隙結構發生改變，即滲流通道發生改變，從而形成了連續均質巖體的非均質滲流，并且在巖體內部產生了流速閥值，且已遠遠超過了初始滲流速度。

3.1 滲流場分析

由不同圍壓梯度下巖體各指定位移節點下滲透速率表面圖如圖6。可知，在應力場和溫度場耦合作用下，原本均質巖體均勻流動轉變成非均勻流動，流體在巖體內流速閥值受圍壓影響而不同。圍壓越大，流體在巖體中的流速閥值越低。流體滲流速度增值按圍壓值升高梯度增值成正向比例增加。圍壓越高，巖體內流體流速閥值點數越多，而且在各圍壓梯度下的流體速度閥值均遠超過流體的初始速度，造成這現象的原因可能是三軸應力下的滲透巖體內部均勻結構發生了改變，巖體內部產生了新生裂隙，進而貫通了巖體內部的新增滲流通道，從而改變了流經巖體的滲透速率，形成了以下巖體滲流場的速度分布狀態。由此可以看出作用于巖體的應力場對流經巖體的滲流場產生的影響是比較明顯的。

圖6 指定位移節點下滲透速率表面圖Fig.6 Surface diagram of penetration rate under specified displacement node

3.2 溫度場分析

通過滲流場、溫度場、應力場三場耦合下的三軸壓縮實驗進行模擬，在模型中心位置每10 mm 截取1個點，共計11 個空間點，獲得了在不同圍壓，不同滲透系數下的溫度數據列，不同滲透系數下溫度變化圖如圖7，不同圍壓梯度下溫度變化圖如圖8。

圖7 不同滲透系數下溫度變化圖Fig.7 Temperature changes under different permeability coefficients

圖8 不同圍壓梯度下溫度變化圖Fig.8 Temperature changes under different confining pressure gradients

由圖7 和圖8 可以看出在同一圍壓不同滲透系數下，由于水的滲流作用，各工況下的巖體溫度均發生非線性降低，而且巖體的滲透系數與巖體溫降梯度成正向比例關系。另外，作用于巖體的圍壓越大，巖體溫降梯度越大。而在同一滲透系數不同圍壓梯度下巖體的降溫梯度是相差不大的，但其溫度降低趨勢均是一致的。盡管巖體在應力場和滲流場綜合作用下溫度產生了非線性降低，但巖體的最低溫度還是高于外界溫度，這表明巖體內部溫度場傳熱產生了一定的耦合作用。綜合上述分析表明，滲流場和應力場的耦合作用會影響巖體溫度場的分布，同時溫度場對滲流場和應力場也會起到一定的附和耦合作用。從以上分析可以得出滲流場對溫度場的影響是起主導作用的，而應力場對溫度場的影響作用并非十分顯著。