卸荷煤體分數階滲透率模型參數探討及其在滲流模擬中的應用

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

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

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

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

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

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

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

1 卸荷煤體應力-滲透率模型分析

從力學本質上講，應力狀態改變是煤體滲透率變化的原因，因此，考慮應力-滲流耦合作用的S-D模型更能從本質上描述煤體滲透率演化規律，該模型的表達式為負指數關系式：

=exp[-3(-)]

(1)

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

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

=-3(-)

(2)

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

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

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

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

(3)

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

圖2 不同ξ和β下滲透率比隨應力變化分布

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

(4)

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

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

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

2 雙/單參數分數階滲透率模型

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

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

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

(5)

式中，為與煤體割理壓縮系數相關聯的參數；為單位應力，MPa；為分數階階數；()為伽碼函數。

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

圖4 不同參數下分數階滲透率曲線

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

2.1 模型驗證與參數取值探討

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

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

圖5 分數階滲透率模型驗證結果

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

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

為獲得損傷變量的演化規律，需要構建損傷本構方程，在連續介質損傷力學中三維線彈性損傷本構方程可表示為

(6)

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

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

=2+(1-)

(7)

式中，為軸向應力;為圍壓;為軸向應變。

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

(8)

式中,為材料參數，由試驗確定；為廣義剪應變，其表達式為

(9)

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

常規三軸加載路徑下的廣義剪應變表達式為

(10)

式中，為徑向應變。

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

(11)

式中，為材料參數，由試驗確定。

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

表1 損傷參數擬合結果

圖6 理論損傷本構曲線和損傷演化曲線

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

2.2 單參數分數階滲透率模型

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

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

表2 模型參數計算結果

圖7 γ,ξc分布關系

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

(12)

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

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

3 分數階滲透率模型在數值模擬中的應用

3.1 煤體應力-損傷-滲流控制方程

將分數階滲透率模型應用于充填采動煤體瓦斯滲流規律數值模擬，借助文獻[13]建立的煤體應力-損傷-滲流控制方程展開研究，煤體流固耦合本構方程和流動方程分別為

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

(13)

(14)

3.2 充填采動煤體瓦斯滲流規律

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

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

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

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

圖9 煤體損傷演化規律

圖10 瓦斯越層流動規律

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

圖11 不同充填體彈性地基系數時上下層煤瓦斯滲流速度

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

因此，提高充填體彈性地基系數，從巖層控制方面來說，有利于降低采場整體擾動程度，而從瓦斯流動方面來說，提高充填體彈性地基系數，本層煤體瓦斯流動能力變差，增加了瓦斯壓力梯度，不利于瓦斯安全防治，但同時抑制住了鄰近層瓦斯流動，又有利于瓦斯安全防治。總而言之，對于多煤層開采的高瓦斯礦井或瓦斯突出礦井，實施充填開采會減小采煤工作面瓦斯濃度，降低瓦斯濃度超限風險，但需要注意隨之而來的本層煤瓦斯突出風險，必須加強本層煤采前預裂瓦斯抽采工作。

4 結 論

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

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

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