卸荷煤體分數(shù)階滲透率模型參數(shù)探討及其在滲流模擬中的應(yīng)用

趙家巍，周宏偉，劉澤霖，謝森林，趙文慧，鐘江城，5，史艷楠，王毅穎

(1.中國礦業(yè)大學(北京)力學與建筑工程學院，北京 100083；2.中國礦業(yè)大學(北京)能源與礦業(yè)學院，北京 100083；3.煤炭資源與安全開采國家重點實驗室，北京 100083；4.中國礦業(yè)大學(北京)應(yīng)急管理與安全工程學院，北京 100083；5.臨沂大學 土木工程與建筑學院，山東 臨沂 276000；6.河北工程大學 機械與裝備工程學院，河北 邯鄲 056038)

目前煤炭開采深度已達1 500 m，并且仍在繼續(xù)延深，未來深部資源開采將成為常態(tài)，在深部賦存環(huán)境下，煤巖體的力學行為呈現(xiàn)出強烈的非線性和不可逆特征，深部煤巖體非線性力學響應(yīng)等許多問題難以采用傳統(tǒng)的巖石力學理論來解釋和描述。深部開采擾動造成煤巖體強烈卸荷，煤體內(nèi)原生裂隙快速擴展，產(chǎn)生嚴重的損傷破壞，引起滲透率急劇變化。目前開采卸荷煤巖體的滲透率演化問題逐漸引起眾多學者的關(guān)注，并從不同角度分別建立了由體積應(yīng)變表征和有效應(yīng)力表征的兩大類滲透率模型。

由體積應(yīng)變表征的滲透率模型以KOZENY-CARMAN方程為基礎(chǔ)，如謝和平等在綜合考慮支承壓力、孔隙壓力和瓦斯吸附膨脹耦合作用的基礎(chǔ)上，推導得出了煤體增透率表達式；程遠平等在建立了考慮有效應(yīng)力和瓦斯吸附解吸變形等因素的、以應(yīng)變?yōu)樽兞康拿后w卸荷損傷增透理論模型；榮騰龍等通過引入了內(nèi)膨脹應(yīng)變系數(shù)的概念，基于DRUCKER-PRAGER破壞準則，建立了考慮煤體損傷破裂的滲透率演化模型；ZHANG等在等效應(yīng)變的框架下建立了雙重孔隙介質(zhì)各向同性損傷本構(gòu)方程，采用標量損傷變量的概念，開發(fā)出滲透演化的動態(tài)模型。

由有效應(yīng)力表征的滲透率模型以SHI-DURUCAN模型為基礎(chǔ)，如楊天鴻等通過引入煤體變形過程中應(yīng)力、損傷與滲透性演化的耦合作用方程，建立了滲流-應(yīng)力滲透率模型，并據(jù)此分析了深部采動影響下瓦斯抽采過程中煤層滲透性的演化和抽采孔周圍瓦斯壓力的變化規(guī)律；薛熠等通過對峰前和峰后滲透率進行分段表達，在峰后階段的模型中加入突變系數(shù)，結(jié)合煤體有效偏應(yīng)力的變化，建立了峰后損傷煤體滲透率模型；張雷等基于CUI-BUSTIN模型(C-B模型)建立了平均有效應(yīng)力作用下的滲透率演化模型；CHEN等利用改進的Logistic增長函數(shù)，建立了滲透率從彈性狀態(tài)到破壞后狀態(tài)變化的模型，提出的滲透率模型是關(guān)于控制煤體壓縮的平均有效應(yīng)力和控制煤體破裂的有效偏應(yīng)力的函數(shù)。

另外，還有一類模型是從宏-細觀雙尺度關(guān)聯(lián)研究的角度建立的，此類模型將細觀力學與斷裂力學相結(jié)合，考慮裂紋擴展和各項異性特征，明確宏-細觀變量之間的關(guān)系，建立起裂隙損傷滲透率模型，如SHAO、MALEKI、陸銀龍等所做的工作。這類模型一般具有相對清晰的物理背景，能夠觸及到問題的本質(zhì)，但其表達式往往過于復雜而不方便用于數(shù)值模擬和工程應(yīng)用。

為了適應(yīng)采動煤巖體非線性力學響應(yīng)特征，在前述研究基礎(chǔ)上，ZHOU等構(gòu)建了分數(shù)階滲透率模型，體現(xiàn)了損傷效應(yīng)對煤體滲透性的影響作用，將經(jīng)典S-D模型的適用范圍拓展至峰后。筆者進一步對分數(shù)階滲透率模型所引入的2個參數(shù)進行深入探討，揭示其物理意義，并將分數(shù)階滲透率模型應(yīng)用于充填采動瓦斯流動規(guī)律的數(shù)值模擬中。

1 卸荷煤體應(yīng)力-滲透率模型分析

從力學本質(zhì)上講，應(yīng)力狀態(tài)改變是煤體滲透率變化的原因，因此，考慮應(yīng)力-滲流耦合作用的S-D模型更能從本質(zhì)上描述煤體滲透率演化規(guī)律，該模型的表達式為負指數(shù)關(guān)系式：

=exp[-3(-)]

(1)

式中，為滲透率;為有效應(yīng)力;下標0表示初始狀態(tài)，MPa；為割理壓縮系數(shù)，MPa，=1/-1/，其中為煤體裂隙的體積模量，MPa，為煤體的體積模量，MPa。

對式(1)兩邊取對數(shù)，并令=ln()，則式(1)可改寫為

=-3(-)

(2)

可以看出，以煤層氣抽采為背景建立的S-D滲透率模型，與Δ構(gòu)成線性關(guān)系，是可以回退的彈性模型。但該模型無法反映采動卸荷煤體的滲透率變化，根據(jù)現(xiàn)有滲透率研究成果，將卸荷滲透率演化趨勢繪制成如圖1所示的概化模型。

由圖1可知，在峰前階段，煤體滲透率隨有效應(yīng)力增加而緩慢降低，峰后階段產(chǎn)生損傷擴容，滲透性質(zhì)發(fā)生改變，在強烈卸荷作用下煤體滲透率急劇增加。

圖1 煤體卸荷滲透率演化概念模型[3]

為定量描述卸荷引起的滲透率變化，其中一個思路是通過引入透氣系數(shù)增大倍率和應(yīng)力敏感因子，得到滲流-應(yīng)力耦合模型，將透氣系數(shù)變換為滲透率后得到

(3)

其模型表征的滲透率變化曲線如圖2所示(以瓦斯壓力=5 MPa，有效應(yīng)力系數(shù)=0.8為例)，通過設(shè)定不同的應(yīng)力敏感系數(shù)和透氣性增大系數(shù)來適應(yīng)不同條件下開采擾動煤體應(yīng)力釋放引起的滲透性變化。

圖2 不同ξ和β下滲透率比隨應(yīng)力變化分布

煤體峰后滲透率強烈變化的實質(zhì)是卸荷導致煤體產(chǎn)生裂隙損傷，為此，有學者將損傷變量和突變系數(shù)嵌入S-D模型得到的滲透率模型為

(4)

其模型表征的滲透率變化曲線如圖3所示(初始有效應(yīng)力=40 MPa)，同樣能夠反映峰后滲透率的急劇變化。

圖3 不同Cf和λD下滲透率比隨應(yīng)力變化分布

從唯象角度來看，原有的S-D模型表征的滲透率增大倍數(shù)有限，通過加入具有調(diào)節(jié)倍率功能的應(yīng)力敏感因子或突變系數(shù)恰是釋放了這種限度，從而得以在有關(guān)采動煤體瓦斯?jié)B流的數(shù)值模擬和工程實踐中應(yīng)用。但是，從細觀機理層面來看，現(xiàn)有應(yīng)力-滲透率模型中的參數(shù)或替代參數(shù)所代表的煤體隔理壓縮系數(shù)被設(shè)定為恒定值，即被視為不隨力學過程變化的靜態(tài)參數(shù)，與煤體損傷擴容、裂隙擴展過程中其割理壓縮系數(shù)不斷變化的事實不符，這也是S-D模型的適應(yīng)性被局限在彈性階段的原因。根據(jù)當前滲透率模型的發(fā)展趨勢，考慮損傷變量的卸荷滲透率模型構(gòu)建思路是可行的，但需要深入探討新增參數(shù)與損傷變量之間的內(nèi)在關(guān)系。

2 雙/單參數(shù)分數(shù)階滲透率模型

開采擾動煤體力學性質(zhì)呈現(xiàn)出強烈的非線性特征，對此，在前述模型基礎(chǔ)上進一步構(gòu)建應(yīng)力-損傷作用下的非線性滲透率模型。對于采動煤體，更多關(guān)注的是煤體峰后卸荷變形階段的滲透率變化，為簡化分析，假設(shè)損傷以應(yīng)力峰值為起點，峰前與應(yīng)力改變量Δ保持線性關(guān)系，峰后兩者呈非線性關(guān)系。

構(gòu)建非線性模型的方法有很多，而蘊含記憶效應(yīng)或歷史依賴性的分數(shù)階導數(shù)模型具備獨特優(yōu)勢，通過借鑒非達西滲流分數(shù)階建模思想，文獻[13]建立了滿足上述條件的雙參數(shù)分數(shù)階滲透率模型：

=exp{[-(-)]1+(2+)}

(5)

式中，為與煤體割理壓縮系數(shù)相關(guān)聯(lián)的參數(shù)；為單位應(yīng)力，MPa；為分數(shù)階階數(shù)；()為伽碼函數(shù)。

該模型揭示了采動卸荷煤體的滲透率變化與應(yīng)力之間的關(guān)系。分別以=05和=002為算例，借助Matlab計算程序可得出一組不同參數(shù)取值條件下滲透率的變化曲線(圖4)。

圖4 不同參數(shù)下分數(shù)階滲透率曲線

由圖4可以看出，分數(shù)階滲透率模型所描述的滲透率變化依賴于分數(shù)階導數(shù)的階數(shù)和割理壓縮系數(shù)相關(guān)的，并且對2個參數(shù)的敏感度較高。特別地，當=0時，分數(shù)階滲透率模型退化為經(jīng)典的S-D模型，表明S-D模型是分數(shù)階滲透率模型的一種特殊情形，=1時，分數(shù)階滲透率模型表現(xiàn)為S-D模型的冪函數(shù)形式，展現(xiàn)峰后滲透率的強非線性變化特征。

2.1 模型驗證與參數(shù)取值探討

為驗證所建立的卸荷煤體分數(shù)階滲透率模型的有效性，利用課題組已進行的不同圍壓下煤巖試樣的常規(guī)三軸壓縮滲流試驗的結(jié)果，進行擬合驗證。試驗煤塊取自平煤12礦己-31030工作面(埋深992～1 118 m)，加工成尺寸為50 mm×100 mm的標準圓柱形煤樣，利用MTS815巖石力學試驗系統(tǒng)進行測試。

采用非線性最小二乘法，對分數(shù)階滲透率模型和滲透試驗數(shù)據(jù)進行擬合分析，確定相關(guān)參數(shù)，滲透率擬合結(jié)果如圖5所示。

圖5 分數(shù)階滲透率模型驗證結(jié)果

對于峰前階段，由于假設(shè)煤體未產(chǎn)生損傷，=0，式(11)退回S-D模型，因此峰前采用S-D模型進行擬合，不再贅述。由圖5可見，分數(shù)階滲透率模型在峰后階段的擬合效果非常好，3個圍壓條件下的確定系數(shù)均大于0.99，說明所建立的分數(shù)階滲透率模型能夠適應(yīng)于卸荷煤體滲透率強烈非線性變化特征。

需要注意的是，模型中的是與煤體割理壓縮系數(shù)相關(guān)聯(lián)的參數(shù)，S-D模型之所以不能適應(yīng)峰后滲透率的演化，正是由于在峰前階段，煤的割理壓縮系數(shù)是由其基本性質(zhì)所確定的常數(shù)，而在峰后階段，煤體產(chǎn)生損傷擴容而對割理壓縮系數(shù)造成了影響。因此在確保滲透率模型有效的同時，還需要建立起和與損傷變量的關(guān)系。

為獲得損傷變量的演化規(guī)律，需要構(gòu)建損傷本構(gòu)方程，在連續(xù)介質(zhì)損傷力學中三維線彈性損傷本構(gòu)方程可表示為

(6)

式中，為主應(yīng)變;為切應(yīng)變;為彈性模量;為泊松比;為主應(yīng)力;為Kronecker符號。

由式(6)得到常規(guī)三軸加載路徑下巖石的損傷本構(gòu)方程

=2+(1-)

(7)

式中，為軸向應(yīng)力;為圍壓;為軸向應(yīng)變。

損傷變量的表達式采用MAZARS提出的混凝土材料在受壓條件下的損傷演化方程，該方程應(yīng)用簡便，并且符合本文以峰值點為損傷起點的假設(shè)，其表達式為

(8)

式中,為材料參數(shù)，由試驗確定；為廣義剪應(yīng)變，其表達式為

(9)

式中，為偏應(yīng)變張量；為應(yīng)變張量；tr()為跡函數(shù)；為Kronecker符號張量。

常規(guī)三軸加載路徑下的廣義剪應(yīng)變表達式為

(10)

式中，為徑向應(yīng)變。

有學者對Mazars方程進行修正，將適用對象擴展到巖石，更好地展示巖石從低圍壓到高圍壓的脆-延性轉(zhuǎn)化特征，提高對損傷閾值后巖石損傷演化過程的擬合精度，修正后的Mazars損傷演化方程為

(11)

式中，為材料參數(shù)，由試驗確定。

根據(jù)損傷本構(gòu)方程和損傷演化方程，對前述加卸載滲流試驗得出的應(yīng)力應(yīng)變曲線進行擬合，確定損傷演化參數(shù)(表1)，得到的理論損傷本構(gòu)曲線和損傷演化曲線如圖6所示。

表1 損傷參數(shù)擬合結(jié)果

圖6 理論損傷本構(gòu)曲線和損傷演化曲線

由表1可以看出，圍壓由10 MPa增加至15 MPa，損傷應(yīng)變閾值提高，圍壓由15 MPa增加至25 MPa，損傷應(yīng)變閾值不變，但損傷演化曲線明顯變緩，高圍壓增加了其延性。

2.2 單參數(shù)分數(shù)階滲透率模型

分數(shù)階滲透率模型對和這2個參數(shù)有著高度敏感性，2個參數(shù)共同決定著滲透率曲線的走勢，盡管取得了不錯的擬合精度，但通過試驗數(shù)據(jù)擬合得出的參數(shù)值并非唯一確定，即對于峰后某個卸荷點處的滲透率比值，對應(yīng)著多個和組合，而在該卸荷點處的損傷變量是唯一確定的，換句話說，和的組合關(guān)系受的控制，為此需借助損傷變量作為中間橋梁進一步探討和的取值關(guān)系。

根據(jù)損傷演化曲線和分數(shù)階滲透率曲線，可以確定峰后任意一點處的損傷變量值和該處的滲透率比值。在不唯一確定的情形中，其值隨的變化而變化，由于階數(shù)和損傷變量的取值范圍都是[0,1]，不妨假設(shè)=，則的值可以反算得出，見表2，兩參數(shù)關(guān)系的分布情況如圖7所示。

表2 模型參數(shù)計算結(jié)果

圖7 γ,ξc分布關(guān)系

從圖7可以看出，和這2個參數(shù)受圍壓變化影響較小，主要還是受損傷變量的影響，2者近似呈線性關(guān)系，由式(5)表達的分數(shù)階滲透率模型可進一步寫成

(12)

從而得到一個只含損傷變量的分數(shù)階滲透率模型，將原模型中的靜態(tài)參數(shù)修改為動態(tài)參數(shù)，更加符合卸荷煤體的實際情況。需要指出的是，對于不同性質(zhì)的煤巖材料，和的關(guān)系式應(yīng)根據(jù)試驗重新確定。

至此，通過深入探討分數(shù)階滲透率模型所引入?yún)?shù)和的相互關(guān)系，得到了損傷變量作為中間控制條件的和的取值關(guān)系，將煤體隔理壓縮系數(shù)的靜態(tài)替代參數(shù)修正為符合峰后變形規(guī)律并受損傷變量所控的動態(tài)參數(shù)，彌補了現(xiàn)有應(yīng)力-滲透率模型中所引入?yún)?shù)與損傷變量不關(guān)聯(lián)的缺陷，使其物理意義更加清晰。

3 分數(shù)階滲透率模型在數(shù)值模擬中的應(yīng)用

3.1 煤體應(yīng)力-損傷-滲流控制方程

將分數(shù)階滲透率模型應(yīng)用于充填采動煤體瓦斯?jié)B流規(guī)律數(shù)值模擬，借助文獻[13]建立的煤體應(yīng)力-損傷-滲流控制方程展開研究，煤體流固耦合本構(gòu)方程和流動方程分別為

[?+(+)??·]+?+=0

(13)

(14)

3.2 充填采動煤體瓦斯?jié)B流規(guī)律

以山西焦煤集團東曲煤礦28802充填開采工作面為工程背景，建立充填采動瓦斯越層流動概化數(shù)值模型，采用多物理場仿真軟件COMSOL Multiphysics中的平面應(yīng)變模塊(Solid Module)和偏微分方程模塊(PDE Module)進行耦合求解，將前述煤體應(yīng)力—損傷—滲流控制方程嵌入PDE模塊中，分數(shù)階滲透率模型由伽瑪函數(shù)編譯求解。計算域和邊界條件如圖8所示，充填工作面位于8號煤層，其下賦存9號煤層，煤巖體相關(guān)參數(shù)取值參考文獻[13]。將充填體假設(shè)為彈性地基支承，分別模擬彈性地基系數(shù)為10，50，90 MN/m三種條件下上下2層煤體瓦斯?jié)B流規(guī)律，重點分析充填體彈性地基系數(shù)對下鄰近層煤體瓦斯流動的影響規(guī)律。

圖8 充填采動瓦斯越層流動概化數(shù)值模型

采用瞬態(tài)求解器，時間步階中的求解時間設(shè)為0～300 min，時間步長為1 min，采用全耦合算法對應(yīng)力場和滲流場一次耦合求解。采用參數(shù)化掃描模擬隨采充填過程，掃描步長設(shè)定為5 m。計算得出的固體力學部分(煤體損傷)和滲流力學部分(瓦斯流動)的結(jié)果分別如圖9，10所示。

從圖9可知，隨彈性地基系數(shù)增大，充填體對頂板和底板的支承能力增強，采場整體擾動程度降低，本層煤體和下層煤體的損傷范圍明顯減小，充填體彈性地基系數(shù)由10 MN/m增加至90 MN/m時本層煤損傷范圍由20 m減小至10 m，下層煤損傷范圍由30 m減小至20 m。從圖10可知，隨彈性地基系數(shù)增大，瓦斯?jié)B流速度減慢，煤體內(nèi)瓦斯流出量減少。

圖9 煤體損傷演化規(guī)律

圖10 瓦斯越層流動規(guī)律

根據(jù)圖11顯示結(jié)果，充填體彈性地基系數(shù)由10 MN/m增加至90 MN/m時，本層煤瓦斯?jié)B流速度由0.13 m/s降低至0.045 mm/s，下層煤瓦斯?jié)B流速度由0.018 mm/s降低至0.003 1 mm/s。根據(jù)圖12顯示結(jié)果，充填參數(shù)變化時，本層煤采動影響范圍內(nèi)的瓦斯含量的衰減總量一致而衰減速度不同，這是由采動和落煤共同影響導致的，本層煤瓦斯衰減總量均由20 kg/m減少至4 kg/m，隨充填體彈性地基系數(shù)增加，瓦斯含量衰減速度減緩，臨近工作面的瓦斯壓力梯度大，從而增加了瓦斯突出風險；當充填體彈性地基系數(shù)為10 MN/m時，下層煤采動影響范圍內(nèi)的瓦斯含量衰減程度最大，由20 kg/m衰減至18.1 kg/m，當充填體彈性地基系數(shù)為90 MN/m時，下層煤采動影響范圍內(nèi)的瓦斯含量衰減程度最小，由20 kg/m衰減至19.2 kg/m。

圖11 不同充填體彈性地基系數(shù)時上下層煤瓦斯?jié)B流速度

圖12 不同充填體彈性地基系數(shù)時上下層煤瓦斯衰減曲線

因此，提高充填體彈性地基系數(shù)，從巖層控制方面來說，有利于降低采場整體擾動程度，而從瓦斯流動方面來說，提高充填體彈性地基系數(shù)，本層煤體瓦斯流動能力變差，增加了瓦斯壓力梯度，不利于瓦斯安全防治，但同時抑制住了鄰近層瓦斯流動，又有利于瓦斯安全防治?？偠灾瑢τ诙嗝簩娱_采的高瓦斯礦井或瓦斯突出礦井，實施充填開采會減小采煤工作面瓦斯?jié)舛龋档屯咚節(jié)舛瘸揎L險，但需要注意隨之而來的本層煤瓦斯突出風險，必須加強本層煤采前預(yù)裂瓦斯抽采工作。

4 結(jié) 論

(1)分析了經(jīng)典S-D模型不能適應(yīng)開采卸荷煤體的原因是割理壓縮系數(shù)不能隨應(yīng)力損傷改變，導致其表征的滲透率增大倍數(shù)有限，與采動卸荷煤體的實際不符，現(xiàn)有的滲透率模型通過加入具有調(diào)節(jié)倍率功能的應(yīng)力敏感因子或突變系數(shù)恰是釋放了這種限度而得以應(yīng)用。

(2)根據(jù)改進的Mazars損傷準則，獲得了煤體卸荷過程中損傷變量的演化規(guī)律，借此探討了分數(shù)階滲透率模型中2個參數(shù)的關(guān)系，進而得到了只含損傷變量的分數(shù)階滲透率模型，將模型中的靜態(tài)參數(shù)修正為動態(tài)參數(shù)，符合損傷擴容過程中煤體割理壓縮系數(shù)不斷變化的事實，彌補了現(xiàn)有應(yīng)力-滲透率模型中所引入的靜態(tài)參數(shù)與損傷變量不關(guān)聯(lián)的缺陷，使其物理意義更加清晰，并且可以比較方便地應(yīng)用于數(shù)值模擬和工程實踐。

(3)將分數(shù)階滲透率模型應(yīng)用于充填采動煤體瓦斯?jié)B流規(guī)律數(shù)值模擬，獲得了充填采動煤體瓦斯?jié)B流規(guī)律，提高充填體彈性地基系數(shù)，加劇了本層煤的瓦斯突出風險而減小了下層煤向本層煤涌入的采煤工作面瓦斯體積分數(shù)。因此，實施充填開采以后，存在利害兩面性，特別是對于瓦斯突出礦井，充填前需加強瓦斯突出風險的評估工作。