李文彬 ，張書(shū)金 ，衡獻(xiàn)偉 ,3

（1.貴州省煤礦設(shè)計(jì)研究院有限公司，貴州 貴陽(yáng) 550025；2.貴州省礦山安全科學(xué)研究院有限公司，貴州 貴陽(yáng) 550025；3.貴州大學(xué) 礦業(yè)學(xué)院，貴州 貴陽(yáng) 550025）

煤層氣是一種以甲烷為主要成分的清潔高效能源，對(duì)其開(kāi)采不僅可以減少溫室氣體排放、接替常規(guī)天然氣能源，同時(shí)也可減少煤礦災(zāi)害的發(fā)生，保障礦井安全生產(chǎn)等[1-2]。目前的研究認(rèn)為[3-5]，煤體是一種多孔介質(zhì)材料，其內(nèi)部廣泛分布的裂隙結(jié)構(gòu)是氣體運(yùn)移的重要通道，對(duì)產(chǎn)氣量有顯著影響。煤體內(nèi)的裂隙分布呈現(xiàn)明顯的隨機(jī)性和無(wú)序性，現(xiàn)有模擬儲(chǔ)層裂隙的方法主要分為等效連續(xù)介質(zhì)模型[6]和離散裂隙網(wǎng)絡(luò)模型[7]，2 種方法各有優(yōu)缺點(diǎn)，在實(shí)際應(yīng)用時(shí)，等效連續(xù)介質(zhì)模型對(duì)裂隙的假設(shè)過(guò)于理想化，無(wú)法描述真實(shí)裂隙的細(xì)節(jié)，而離散裂隙網(wǎng)絡(luò)模型恰好彌補(bǔ)了這一不足，可以充分刻畫(huà)裂隙的真實(shí)細(xì)節(jié)和分布特征，因而被廣泛應(yīng)用。其中，采用Monte Carlo 方法生成離散裂隙網(wǎng)絡(luò)模型是目前最多的方法，許多學(xué)者利用該方法分析了煤層氣在裂隙內(nèi)的流動(dòng)規(guī)律。然而傳統(tǒng)離散裂隙網(wǎng)絡(luò)模型仍然是假設(shè)裂隙為平滑裂隙，對(duì)真實(shí)裂隙細(xì)節(jié)的描述仍然存在不足。而真實(shí)煤體內(nèi)的裂隙絕大多數(shù)是凹凸不平的，氣體在其內(nèi)的流動(dòng)也不同于平滑裂隙網(wǎng)絡(luò)模型。因此，采用一種新的粗糙離散裂隙網(wǎng)絡(luò)建模方法建立裂隙網(wǎng)絡(luò)幾何模型，該方法相比Monte Carlo 方法可充分考慮裂隙粗糙特征，因而更加符合實(shí)際情況。為此，基于建立的裂隙網(wǎng)絡(luò)模型，結(jié)合數(shù)值模擬分析煤層氣在裂隙內(nèi)的流動(dòng)規(guī)律，探討不同裂隙結(jié)構(gòu)特征參數(shù)對(duì)煤層氣流動(dòng)特性的影響。

1 粗糙離散裂隙網(wǎng)絡(luò)模型

實(shí)際煤體內(nèi)的裂隙往往是彎曲的，具有一定的粗糙特性，且裂隙相互交叉形成裂隙網(wǎng)絡(luò)，組合形式多為網(wǎng)狀和不規(guī)則網(wǎng)狀，滲透率相對(duì)較好。為了更真實(shí)地建立符合實(shí)際工程情況的粗糙離散裂隙網(wǎng)絡(luò)模型，WU 等[8]和XIA 等[9]等提出了基于像元概率分解算法建立粗糙離散裂隙網(wǎng)絡(luò)模型的方法，根據(jù)其理論，在像元空間中，任意1 條裂隙都可以由點(diǎn)集合的形式組成，即：

式中：GF為裂隙點(diǎn)集合系統(tǒng)；Di為裂隙點(diǎn)的成長(zhǎng)方向；Mi為成長(zhǎng)概率；ti為成長(zhǎng)步長(zhǎng)。

根據(jù)像元概率分解算法的基本理論，借助MATLAB 軟件編制二維粗糙離散裂隙網(wǎng)絡(luò)生成程序，以實(shí)現(xiàn)像元空間范圍內(nèi)多條裂隙的隨機(jī)生成。由式（1）可知，如果確定了裂隙的成長(zhǎng)方向、成長(zhǎng)概率、成長(zhǎng)步長(zhǎng)，裂隙點(diǎn)集合系統(tǒng)也就確定了。為了更好地說(shuō)明其具體實(shí)現(xiàn)過(guò)程，假定裂隙的成長(zhǎng)概率M1={0.67, 0.33, 0, 0, 0, 0, 0, 0}，為二項(xiàng)隨機(jī)，成長(zhǎng)步長(zhǎng)ti=15，則1 條在像元空間范圍[100 100]中裂隙點(diǎn)集合，粗糙離散裂隙網(wǎng)絡(luò)生成過(guò)程如圖1。

圖1 粗糙離散裂隙網(wǎng)絡(luò)生成過(guò)程Fig.1 Generation process of rough discrete fracture network

當(dāng)需要生成具有同樣特征的150 條裂隙點(diǎn)集合時(shí)，原理與生成1 條裂隙一致，只需要在程序中設(shè)置裂隙總數(shù)，裂隙在像元空間中隨機(jī)分布，結(jié)果如圖1（b）。由于真實(shí)煤巖體中的裂隙相互交錯(cuò)，為生成更為符合真實(shí)情況的裂隙網(wǎng)絡(luò)，這里裂隙的分布特征參考文獻(xiàn)[10-11]獲得的統(tǒng)計(jì)規(guī)律，即裂隙長(zhǎng)度冪律符合分布特征，傾角分布符合Fisher分布特征。在圖1（b）中增加150 條裂隙，成長(zhǎng)概率M2={0, 0, 0, 0, 0, 0, 0.33, 0.67}，成長(zhǎng)步長(zhǎng)ti=15，因此，具有不同成長(zhǎng)概率，數(shù)量為300 條的二維像元空間中的裂隙網(wǎng)絡(luò)集合系統(tǒng)如圖1（c）。

2 粗糙離散裂隙網(wǎng)絡(luò)滲透率模型

煤體內(nèi)的裂隙相互交叉形成裂隙網(wǎng)絡(luò)，對(duì)于任意復(fù)雜的裂隙網(wǎng)絡(luò)均是由一條條單裂隙所組成。目前描述流體在單裂隙內(nèi)的流動(dòng)多采用經(jīng)典的平行板立方定律[12]，即

式中：q為流過(guò)裂隙的總流量，m3/s；α為裂隙開(kāi)度，m；l為裂隙長(zhǎng)度，m；θ1、θ2分別為裂隙方位角和傾角，（°）；μ為流體動(dòng)力黏度，μPa·s；Δp為壓差，Pa；L0為特征單元體長(zhǎng)度，m。

MIAO 等[11]借鑒多孔介質(zhì)模型中孔隙網(wǎng)絡(luò)滲透率模型的推導(dǎo)方法，在上述平行板立方定律的基礎(chǔ)上，結(jié)合分形結(jié)合推導(dǎo)得到了適用于隨機(jī)分布的裂隙網(wǎng)絡(luò)滲透率模型，即：

式中：Kf為裂隙網(wǎng)絡(luò)滲透率，m2；β為比例系數(shù)，取值范圍一般為10-3～10-1；Dl為裂隙網(wǎng)絡(luò)的分形維數(shù)；?f為裂隙率；lmax為最大裂隙長(zhǎng)度，m。

然而，上述模型是假設(shè)煤巖體內(nèi)的裂隙為平滑裂隙，沒(méi)有考慮裂隙的粗糙或迂曲度特征。針對(duì)上述平行板裂隙網(wǎng)絡(luò)模型存在的不足，LUO 等[12]改進(jìn)了上述模型，在上述模型中增加了裂隙的粗糙特征，即：

式中：DTf為裂隙迂曲度分形維數(shù)。

式（4）考慮了裂隙的粗糙特征，因而也更加符合實(shí)際情況。因此，在后續(xù)章節(jié)將主要結(jié)合該滲透率模型和粗糙離散裂隙網(wǎng)絡(luò)建模方法分析氣體在裂隙內(nèi)的流動(dòng)特性。

3 煤層氣流動(dòng)數(shù)值模擬

由于實(shí)際煤層氣的流動(dòng)是一個(gè)極其復(fù)雜的過(guò)程，涉及的影響因素眾多，研究?jī)H分析裂隙結(jié)構(gòu)特征參數(shù)的變化對(duì)煤層滲透性的影響。因此，研究重點(diǎn)在于關(guān)注不同裂隙結(jié)構(gòu)特征對(duì)滲流的影響，故不考慮其它因素的影響，即假設(shè)氣體為單相流，遵循達(dá)西定律，忽略應(yīng)力、溫度等因素的變化對(duì)滲流的影響。

3.1 邊界條件和網(wǎng)格劃分及氣體流動(dòng)控制方程

為采用像元概率分解算法建立的粗糙離散裂隙網(wǎng)絡(luò)幾何模型，然后將建立的裂隙網(wǎng)絡(luò)模型導(dǎo)入到COMSOL 多物理場(chǎng)耦合軟件中進(jìn)行煤層氣流動(dòng)數(shù)值計(jì)算，煤層氣流動(dòng)數(shù)值模型如圖2。模型尺寸為100 mm×100 mm，在模型的左側(cè)設(shè)置壓力入口，入口壓力為1 000 Pa，模型右側(cè)為出口壓力，壓力為0，將模型下邊界固定，上邊界為不滲流邊界。然后對(duì)模型進(jìn)行網(wǎng)格劃分，這里采用自由三角形網(wǎng)格對(duì)其進(jìn)行劃分，對(duì)特別關(guān)注的裂隙區(qū)域進(jìn)行網(wǎng)格加密處理。三角形單元數(shù)量為57 617 個(gè)，邊單元數(shù)量為373 個(gè)，最小單元質(zhì)量0.393 1 mm，平均單元質(zhì)量0.824 6 mm。

圖2 煤層氣流動(dòng)數(shù)值模型Fig.2 Numerical modeling of coalbed methane flow

數(shù)值計(jì)算過(guò)程中，煤層氣的流動(dòng)控制方程采用COMSOL 內(nèi)置的Navier-Stokes 方程[13]，即：

式中：V為氣體流速，m/s； ?為哈密頓算子。

3.2 模擬結(jié)果

設(shè)定邊界條件及相關(guān)參數(shù)后，通過(guò)達(dá)西穩(wěn)態(tài)流求解控制方程（式（6）），可得到氣體壓力分布云圖以及氣體在裂隙內(nèi)的流速分布云圖，氣體壓力分布云圖如圖3，氣體在裂隙網(wǎng)絡(luò)中的流速分布云圖如圖4。

圖3 氣體壓力分布云圖Fig.3 Cloud diagram of gas pressure distribution

圖4 氣體在裂隙網(wǎng)絡(luò)中的流速分布云圖Fig.4 Cloud image of gas velocity distribution in fracture network

由圖3 可知，由于煤體內(nèi)裂隙的存在，使得壓力分布呈現(xiàn)非均勻分布特征。由圖4 可知，氣體的流動(dòng)主要發(fā)生在煤體內(nèi)相互連通的裂隙網(wǎng)絡(luò)內(nèi)，裂隙的分布特征影響著氣體的流動(dòng)特性。

4 參數(shù)分析及討論

由式（4）可知，裂隙網(wǎng)絡(luò)滲透率主要與裂隙率、分形維數(shù)、迂曲度分形維數(shù)、裂隙傾角、最大裂隙長(zhǎng)度以及比例系數(shù)等參數(shù)有關(guān)。因此，為了進(jìn)一步分析裂隙特征參數(shù)對(duì)煤層滲透性的影響規(guī)律，采用像元概率分解算法建立具有不同結(jié)構(gòu)特征的粗糙離散裂隙網(wǎng)絡(luò)模型，結(jié)合數(shù)模模擬方法分析不同裂隙特征參數(shù)對(duì)煤層滲透性的影響，裂隙特征參數(shù)見(jiàn)表1。

表1 裂隙特征參數(shù)Table 1 Parameter values of fracture characteristics

4.1 裂隙率和分形維數(shù)

裂隙率和分形維數(shù)是表征裂隙復(fù)雜程度的重要參數(shù)，其中裂隙率為裂隙的總面積與測(cè)量煤巖樣品的比值，其值越大表明裂隙數(shù)量越多，越容易形成更多的流動(dòng)通道。同樣，分形維數(shù)也是表征裂隙復(fù)雜程度常用的定量表征參數(shù)，是裂隙占有度量空間的多少，裂隙分形維數(shù)越大，裂隙分布越密集，越容易形成更多的滲流通道。為具體分析裂隙率和分形維數(shù)的變化對(duì)煤層滲透性的影響，固定其它裂隙參數(shù)不變（其中每個(gè)案例中裂隙的傾角為45°、迂曲度分形維數(shù)1.015、比例系數(shù)為0.014 以及最大裂隙長(zhǎng)度為17 mm），采用像元概率分解算法生成具有不同裂隙結(jié)構(gòu)參數(shù)的粗糙離散裂隙網(wǎng)絡(luò)幾何模型，裂隙率和分形維數(shù)的變化范圍分別為0.1～0.25、1.55～1.85。

裂隙率與裂隙網(wǎng)絡(luò)模型出口流速的關(guān)系如圖5，分形維數(shù)與裂隙網(wǎng)絡(luò)模型出口流速的關(guān)系如圖6。

圖5 裂隙率與裂隙網(wǎng)絡(luò)模型出口流速的關(guān)系Fig.5 Relationship between fracture rate and outlet velocity of fracture network model

圖6 分形維數(shù)與裂隙網(wǎng)絡(luò)模型出口流速的關(guān)系Fig.6 Relationship between fractal dimension and outlet velocity of fracture network model

從圖5 可以看出：隨著裂隙率的增大，模型出口流速也隨之增大，通過(guò)對(duì)模擬數(shù)據(jù)擬合發(fā)現(xiàn)，二者基本呈現(xiàn)較好的冪律增長(zhǎng)關(guān)系。

由圖6 可知：裂隙網(wǎng)絡(luò)分形維數(shù)的增加使得模型出口速度也隨之增大。對(duì)其數(shù)據(jù)進(jìn)行擬合發(fā)現(xiàn)，分形維數(shù)與出口流速之間同樣展現(xiàn)了較好的冪律增長(zhǎng)關(guān)系。

4.2 迂曲度分形維數(shù)和最大裂隙長(zhǎng)度

分析裂隙迂曲度分形維數(shù)的變化對(duì)煤層滲透性的影響所需參數(shù)見(jiàn)表1。固定其它裂隙參數(shù)不變，迂曲度分形維數(shù)和最大裂隙長(zhǎng)度的變化范圍分別為1.0～1.02 和11～26 mm。迂曲度分形維數(shù)與裂隙網(wǎng)絡(luò)模型出口流速的關(guān)系如圖7，最大裂隙長(zhǎng)度與裂隙網(wǎng)絡(luò)模型出口流速的關(guān)系如圖8。

圖7 迂曲度分形維數(shù)與裂隙網(wǎng)絡(luò)模型出口流速的關(guān)系Fig.7 Relationship between tortuosity fractal dimension and outlet velocity of fracture network model

圖8 最大裂隙長(zhǎng)度與裂隙網(wǎng)絡(luò)模型出口流速的關(guān)系Fig.8 Relationship between maximum fracture length and outlet velocity of fracture network model

由圖7 可知：迂曲度分形維數(shù)與流速之間呈現(xiàn)了負(fù)相關(guān)關(guān)系，即隨著迂曲度分形維數(shù)的增大，裂隙網(wǎng)絡(luò)模型出口流速逐漸減小，這是因?yàn)橛厍确中尉S數(shù)表征的是裂隙的彎曲程度，迂曲度分形維數(shù)值越大，說(shuō)明裂隙越彎曲（或粗糙），則氣體在裂隙內(nèi)的流動(dòng)阻力也就越大，因此，出口流速逐漸減小。

由圖8 可知：隨著裂隙長(zhǎng)度的逐漸增大，模型出口流速也隨之增大。這是因?yàn)榱严蹲畲箝L(zhǎng)度的增加意味著裂隙之間更易形成網(wǎng)絡(luò)，使得裂隙的連通性整體性增強(qiáng)，進(jìn)而影響著煤層滲透性。

4.3 比例系數(shù)和裂隙傾角

同樣，為分析比例系數(shù)和裂隙傾角對(duì)煤層滲透性的影響，固定裂隙長(zhǎng)度、迂曲度分形維數(shù)等參數(shù)不變，比例系數(shù)變化范圍為0.006～0.08，傾角變化范圍為10°～85°。比例系數(shù)與裂隙網(wǎng)絡(luò)模型出口流速的關(guān)系如圖9，裂隙傾角與裂隙網(wǎng)絡(luò)模型出口流速的關(guān)系如圖10。

圖9 比例系數(shù)與裂隙網(wǎng)絡(luò)模型出口流速的關(guān)系Fig.9 Relationship between proportional coefficient and outlet velocity of fracture network model

圖10 裂隙傾角與裂隙網(wǎng)絡(luò)模型出口流速的關(guān)系Fig.10 Relationship between fracture inclination angle and outlet velocity of fracture network model

由圖9 可知：隨著比例系數(shù)的增加，裂隙網(wǎng)絡(luò)模型出口流速逐漸增大。其原因在于，當(dāng)裂隙長(zhǎng)度不變時(shí)，比例系數(shù)與裂隙的開(kāi)度成正相關(guān)關(guān)系，比例系數(shù)的增加意味著裂隙開(kāi)度的增大，因此氣體在裂隙內(nèi)的流動(dòng)更為容易，導(dǎo)致裂隙網(wǎng)絡(luò)模型出口流速增大。此外，由圖中數(shù)據(jù)的變化可知，比例系數(shù)的增大對(duì)裂隙網(wǎng)絡(luò)流速的影響最為顯著。

由圖10 可知：隨著傾角的逐漸增大，裂隙網(wǎng)絡(luò)模型出口流速逐漸減小，這是由于較小裂隙傾角將更有利于氣體流動(dòng)，而傾角越大，氣體流動(dòng)阻力越大，可見(jiàn)煤層中裂隙傾角分布方向決定裂隙網(wǎng)絡(luò)的滲流方向。

綜上，為了更好地觀察裂隙特征參數(shù)與流速之間的關(guān)系，對(duì)模擬結(jié)果進(jìn)行非線性擬合，擬合結(jié)果表明，裂隙率、分形維數(shù)、最大裂隙長(zhǎng)度以及比例系數(shù)與裂隙網(wǎng)絡(luò)出口流速之間呈現(xiàn)較好的冪律增長(zhǎng)關(guān)系，而迂曲度分形維數(shù)、裂隙傾角與裂隙網(wǎng)絡(luò)出口流速之間呈現(xiàn)冪律降低的關(guān)系。

5 結(jié) 語(yǔ)

1）詳細(xì)介紹了粗糙裂隙離散裂隙網(wǎng)絡(luò)的建模理論及其實(shí)現(xiàn)過(guò)程，并基于該方法建立了更符合工程實(shí)際的粗糙離散裂隙網(wǎng)絡(luò)模型。

2）由于煤體內(nèi)裂隙網(wǎng)絡(luò)的存在，氣體壓力分布呈現(xiàn)非均勻分布特征，氣體的流動(dòng)主要發(fā)生在相互連通的裂隙網(wǎng)絡(luò)中，煤層滲透率主要為裂隙網(wǎng)絡(luò)所控制。

3）裂隙率、分形維數(shù)、最大裂隙長(zhǎng)度以及比例系數(shù)與裂隙網(wǎng)絡(luò)出口流速之間呈現(xiàn)較好的冪律增長(zhǎng)關(guān)系，而迂曲度分形維數(shù)、裂隙傾角與裂隙網(wǎng)絡(luò)出口流速之間呈現(xiàn)冪律降低的關(guān)系；在這些參數(shù)中，比例系數(shù)、裂隙傾角和最大裂隙長(zhǎng)度對(duì)煤巖滲透率的影響最為顯著。