連續加載應力下真實裂縫流場和滲透率演化規律數值研究

梁運培，張懷軍，王禮春，秦朝中，田鍵，陳強，史博文

（1.重慶大學煤礦災害動力學與控制全國重點實驗室，重慶 400044；2.重慶大學資源與安全學院，重慶 400044；3.天津大學表層地球系統科學研究院，天津 300072）

水力壓裂技術已被大規模運用到非常規油氣、深部地熱資源等地下能源開采中，且取得了顯著成效[1-3]。通過水力壓裂作業，原本致密的地層在高水壓下產生大量人工裂縫，最終形成主裂縫、次裂縫以及天然裂縫相互串通的復雜裂縫網絡體系[4-9]。由于尺寸匹配、顆粒懸浮等問題，壓裂液中的一些支撐劑無法進入次裂縫以及天然裂縫中進行有效支撐[10]，這類未被有效支撐的裂縫極易在地層應力擾動下出現不同程度的閉合，從而基本上失去裂縫的導流能力。因此，研究裂縫在應力擾動下的滲透率演化特征成了揭示裂縫導流能力損傷機制的前提。

國內外學者通過開展裂縫樣品應力（加載與卸載）實驗，測試了不同應力作用下的裂縫滲透率變化特點[11-15]，這些信息對分析物質在裂縫內部的運移行為至關重要。但采用實驗方法難以直接探究應力作用下裂縫內部流場、流動通道與速度場的演化規律，因此，越來越多的學者開始采用數值模擬方法揭示裂縫在應力作用下的導流能力變化機制。學者們通過數值生成的裂縫，發現應力擾動下裂縫的導流能力主要與裂縫的平均開度、表面粗糙度以及接觸面積等因素有關，并在此基礎上建立起表征裂縫滲透率的相關經驗模型[16-20]。然而，通過數值隨機生成的裂縫通常在裂縫幾何形態上滿足平均化特征，與人工裂縫和部分天然裂縫的幾何形態表現出的強非均質性存在顯著差別。因此，通過直接數值模擬研究非均質裂縫滲透率演化行為存在適用性受限問題。

本研究首先利用巴西劈裂法制備真實人工裂縫巖樣，提取裂縫表面形貌基本特征參數后，再通過數值方法拼接裂縫表面形成具有非均質裂縫特征信息的單條裂縫結構，隨后，研究裂縫在有效應力下的滲透率演化特征，分析裂縫流場、接觸面積、空間相關長度等參數的響應行為。最后，將研究模擬結果與傳統經典經驗模型進行驗證，證實了采用真實裂縫獲得的非均質裂縫與直接數值生成的裂縫在滲透率演化行為以及滲透率模型匹配關系上存在明顯差別。

1 單裂縫開度分布及非均質性

1.1 巖石裂縫提取

裂縫在頁巖氣藏和致密砂巖氣藏中普遍存在，對天然氣的采出發揮了重要作用。因此，選取頁巖和致密砂巖作為研究對象，同時選取花崗巖作為對比樣，分析3種典型巖石中裂縫的單相滲流行為。首先，對3 類巖石鉆取直徑為50 mm、長度為100 mm 的柱塞巖樣用于開展力學性能測試分析，獲取彈性模量和泊松比，為建立接觸力學模型提供基準參數。隨后，鉆取直徑為25 mm、長度為50 mm 的3 類巖石標準柱塞巖樣用于人工造縫，獲取真實裂縫的幾何形貌等信息。

采用MTS815 巖石力學試驗設備測量50 mm×100 mm樣品的多組巖石力學參數，平均結果如表1所示。測試時，通過軸向位移控制加載方式，頁巖和致密砂巖的加載速率為0.1 mm/min，考慮花崗巖強度更高、變形更慢，對花崗巖的加載速率設定為0.05 mm/min。

表1 不同種類巖石基本力學參數測試平均結果Table 1 Average results of basic mechanicalparameters of different kinds of rocks

分別選取25 mm×50 mm 的頁巖、致密砂巖和花崗巖樣品各2塊，采用巴西劈裂法對樣品進行人工造縫（圖1）。造縫儀（型號AG-250kN IS）加載方式為軸向位移控制，試件的加載速率為0.1 m/min。為了使巖樣均勻受力，造縫時將巖樣水平放在凹槽中，造縫儀壓頭與巖樣表面保持接觸平整，使巖樣受壓形成張拉裂縫，隨后，采用Cronos Dual 三維光學掃描儀對人工裂縫進行掃描（設備分辨率為300 萬像素，掃描精度為0.02 mm），獲取裂縫表面的空間位置等形貌參數信息。

圖1 基于巴西劈裂法的人工裂縫樣品制作實物Fig.1 Photos of artificial fracture sample via Brazilian splitting method

將每個裂縫樣品上下2個表面的形貌數據，通過位置坐標進行拼接得到整合裂縫，整合裂縫的拼接需滿足2 個條件：①在垂直方向上，選擇上下2 個裂縫面的最高點和最低點，確保其在同一條豎直線上呈對準狀態；②在水平方向上，當上下2 個裂縫表面拼接后，2個裂縫面輪廓線能夠重合。通過插值方法對拼接后的裂縫表面進行長度x方向和寬度y方向網格單元劃分，網格單元尺寸為l1（長度）×lw（寬度）=0.16 mm×0.16 mm。在掃描提取裂縫時，由于出現了部分邊界損失，將這部分損失的邊界切除后，得到所有裂縫中尺寸最小的裂縫長度為46.88 mm，寬度為21.92 mm。為了方便計算和對比，將其余掃描裂縫的尺寸也處理為長度為46.88 mm，寬度為21.92 mm，裂縫拼接完成后，裂縫的開度可由下列公式得出：

式中：b為裂縫開度大小，單位mm；Z1、Z2分別為裂縫同一基準線上的上表面位置高程和下表面位置高程。在有效應力下，如果裂縫面發生相互接觸，則b=0。

1.2 裂縫開度分布特征

在開展數值研究前，對6條裂縫初始狀態下的裂縫開度進行分析，結果如圖2 和表2 所示。開度場云圖（圖2）表明，3 種樣品裂縫開度分布的空間相關性與巖石類型有關[16,21]。其中，花崗巖裂縫G1 和G2 開度分布呈現出左側開度小、右側開度大的特點；致密砂巖裂縫TS1和TS2開度整體呈現出與水平y方向為斜角45°分布特征；頁巖裂縫S1 和S2 開度則表現為水平y方向底部開度大。此外，裂縫開度大小概率密度分布統計顯示，花崗巖裂縫的開度分布范圍相對更寬，平均開度更大，表明裂縫開度分布相對不均勻，開度高低起伏較明顯；頁巖裂縫和致密砂巖裂縫開度范圍較為集中，表明開度大小分布相對集中。表2 給出了3 種巖石裂縫開度分布在x和y方向上的相關長度，表現為花崗巖裂縫開度空間分布相對集中，致密砂巖裂縫開度空間分布范圍相對更廣，頁巖裂縫開度空間分布則主要沿著裂縫y方向分布且y方向上裂縫開度的相關長度與寬度相等。

圖2 不同裂縫樣品開度場分布云圖及相應的裂縫開度概率密度分布Fig.2 Results of fracture aperture distribution and corresponding aperture probability density distribution

表2 不同裂縫樣品開度初始參數統計Table 2 Statistical analysis results of apertures of different fracture samples

2 裂縫接觸力學及滲流模型

2.1 接觸力學模型

學者們通常采用接觸力學模型數值研究裂縫在有效應力載荷下的變形規律，該模型遵循HOPKINS[22]、PYRAK-NOLTE等[23]、PETROVITCH等[24]、KLING等[25]提出的方法，且通過了大量實驗驗證。研究根據接觸力學模型分析裂縫應力敏感，在進行裂縫變形分析時，將裂縫面內的微凸起簡化為圓柱體，并假設在裂縫的上表面和下表面分別規定一塊基準平板，將有效應力均勻施加在該平板上，而在有效應力連續加載過程中，當裂縫滲透率為0時，停止分析，從而得到裂縫在不同有效應力下開度、粗糙度和相關長度等變化規律。根據研究方法，裂縫表面某一圓柱凸起i的變形計算公式為[20,24]：

式中：Wi為粗糙圓柱凸起i的總變形；wii為圓柱凸起i受作用力后的法向變形；wij是圓柱凸起i由鄰近凸起j引起的變形。

在圓柱i施加的力為fi，則圓柱i的法向變形wii為：

式中：θ為角度，單位（°）；I1為圓柱i在法向應力作用下的法向變形函數，單位m；r0為圓柱的半徑，單位m；r為圓柱i中心的距離，單位m。

其中：

式中：fi為法向應力，單位N；E為彈性模量，單位GPa；v為泊松比。

wij可以用下式表達：

式中：r1和r2是圓柱j沿r方向到圓柱i中心的距離，單位m；θ1和θ2為圓弧和圓柱體j之間的交點；I2為圓柱i由臨近凸起j在擠壓作用下的橫向變形函數，單位m。

其中：

除了2個裂縫面的變形外，凸起長度在法向應力fi的作用下同樣被壓縮，則：

被壓縮后，上下裂縫面距離D和凸起變形Wi之和需等于凸起變形后的長度，則：

式中：D為上下裂縫面距離，單位m。

2.2 單相滲流模型

單相滲流模型將裂縫內的流動空間劃分為多個長方體，如圖3a所示，每一個長方體內的流量遵循質量守恒定律。在不可壓縮假設下，對每個開度單元采用體積守恒定律，通過Hagen–Poiseuille 方程計算裂縫內流體壓力分布[26-27]（圖3）：

圖3 裂縫開度網格示意圖Fig.3 Schematic diagram of fracture aperture grid

由于考慮裂縫表面之間的起伏，加入修正因子，將式（9）修正為：

式中：l′ij為i到j的傾斜距離，單位m；li和lj分別為i、j網格長度的一半，單位m。

對每個網格i上的水力傳導系數gi通過局部立方定律得到：

式中：gi為對每個網格i上的水力傳導系數，單位m3/Pa；bi為i點的開度；Ai為i的橫截面積，單位m2；μ為黏度，單位Pa·s，研究中為單相水的流動，黏度設置為1×10-3Pas；lw為網格寬度，單位m。

式（11）為確定每個開度網格上的水力傳導系數，然而流體總是從一個網格i流進另一個網格j，因此，不能只用一個網格的水力傳導系數計算壓力，需要對其做調和平均：

將計算得到的gij代入公式（10），同時將裂縫入口壓力設置為10 Pa，出口壓力設置為1 Pa，通過求解得到裂縫中流體的壓力分布情況，將壓力分布帶入式（10）得出裂縫內流體流量的分布，并根據壓力分布計算裂縫入口流量Q。

通過立方定律[28]計算出水力開度bh：

式中：bh為水力開度，單位m；Q為裂縫入口流量，單位m3/s；?p為入口與出口之間的壓力差，單位Pa；w為裂縫寬度，單位m；L為裂縫長度，單位m。

根據平板模型，裂縫滲透率可以用下式表示[19,28]：

式中：k為裂縫滲透率，單位m2。

3 結果與討論

3.1 不同有效應力下滲透率變化行為

由于同類裂縫樣品在不同有效應力下流場和滲透率演化模擬結果趨勢一致，因此，針對3 種裂縫樣品各選取其中1塊作為代表樣品，分析裂縫內流體流場隨有效應力變化的響應特征，所有流動方向為從左往右。模擬結果表明（圖4）：由于裂縫開度的強非均質性，在有效應力作用下花崗巖和致密砂巖流場演化方向與流體流動方向呈現為垂直交叉特點，而頁巖裂縫的流場變化方向與流體流動方向平行，這種由于裂縫開度非均質性和沿流動方向上相關長度的影響對滲透率的控制作用，將在進一步研究中討論。

圖4 有效應力加載過程中不同裂縫內流場演化結果Fig.4 Results of flow field distribution under loading effective stress

研究分析了3 個裂縫樣品的滲透率隨有效應力變化的模擬結果（圖5）：其中，花崗巖裂縫樣品G1 的初始滲透率為5.64×103μm2，模擬結束后滲透率降至6.58×10-5μm2，降幅為8 個數量級；致密砂巖裂縫樣品TS1 初始滲透率為1.51×103μm2，模擬結束后滲透率降至3.39×10-3μm2，降幅為6 個數量級；頁巖裂縫樣品S1 的初始滲透率為1.27×103μm2，模擬結束后滲透率降至3.55×101μm2，降幅在2 個數量級。結合圖2可以看出：花崗巖和致密砂巖裂縫樣品的滲透率隨有效應力的增加下降程度更為明顯，而頁巖裂縫的滲透率變化則相對較小，主要是頁巖S1 裂縫下邊緣開度較大，在有效應力下一直難以壓實，這部分未被壓實的開度對整個裂縫的導流能力貢獻大，所以滲透率變化較緩慢。

圖5 有效應力下不同巖石裂縫滲透率演化結果Fig.5 Evolution of fracture permeability with increasing loading effective stress

3.2 有效應力作用下裂縫特征參數響應

3.2.1 裂縫開度

研究分析了有效應力變化下3 種裂縫樣品的裂縫開度演化過程（圖6）：裂縫開度演化與裂縫開度的非均質性有關，花崗巖G1 在有效應力不斷增加過程中，裂縫右側開度降幅有限，而裂縫左側更容易發生閉合；致密砂巖TS1 在有效應力作用下，裂縫內整體開度變化相對均勻；頁巖S1的開度變化則表現為從y方向上的一側往另一側逐漸壓合的過渡行為。結果表明，有效應力加載過程中，裂縫開度變化行為與圖4中的流場分布演化呈現出良好的對應關系。

圖6 不同有效應力下裂縫開度演化結果Fig.6 Evolution of fracture aperture distribution under different effective stress

3.2.2 空間相關長度

根據裂縫在x和y方向上的實際縫長和縫寬大小，將裂縫開度分布與x和y方向上的空間相關長度做歸一化處理，得到歸一化后裂縫開度的空間相關長度介于0~1，其中，0 代表沒有空間相關長度，1 代表x或y方向上的裂縫開度分布貫穿整個裂縫。分析在初始狀態下，只有頁巖S1 裂縫開度的相關長度在y方向上貫穿裂縫，且所有裂縫的相關長度在x和y方向上呈現顯著差異，說明3 類裂縫開度表現出各向異性[29]。當施加的有效應力小于0.5 MPa 時，所有裂縫的空間相關長度變化較??；隨著有效應力超過0.5 MPa，所有裂縫的空間相關長度開始隨有效應力的增加呈下降趨勢，并最終趨于平緩。分析結果表明（圖7），致密砂巖裂縫的相關長度和頁巖裂縫沿y方向的相關長度降幅明顯，而花崗巖裂縫的空間相關長度則相對緩慢。由表1可知，花崗巖的彈性模量為26 GPa，高于致密砂巖和頁巖，且花崗巖裂縫開度分布不均勻，因而有效應力加載過程中，花崗巖裂縫整體上不容易被壓縮，其裂縫右側還有大面積的開度未被壓縮（圖6），因此，花崗巖的相關長度減小緩慢。

圖7 裂縫相關長度隨有效應力變化的演化結果Fig.7 Changes of fracture correlation length with increasing effective stress

3.2.3 接觸面積

分析有效應力加載過程中裂縫巖樣接觸面積的變化特征，結果表明（圖8）：隨著有效應力的增加，所有裂縫接觸面積占比逐漸增大且呈現近似線性關系，但不同裂縫接觸面積增長速率差異越發明顯。在相同的有效應力下，由于花崗巖的強度高和開度分布不均勻，花崗巖裂縫的接觸面積遠小于致密砂巖和頁巖裂縫的接觸面積，而致密砂巖TS1裂縫開度的相關長度較小且巖樣的強度更低，因此，在有效應力下致密砂巖裂縫接觸面積增長的速率最快。

圖8 有效應力作用下接觸面積與滲透率變化結果Fig.8 Changes of contact area and fracture permeability under effect of effective stress

通過分析接觸面積占比與滲透率之間的關系，對比圖8a 和圖8b 可知：裂縫接觸面積的增加與裂縫滲透率的減小并未呈現出良好的正相關性，其中，頁巖裂縫有效應力加載后的最終接觸面積增加幅度更明顯，但其滲透率降幅卻相對最小。結合流場分布和裂縫開度演化可知，相比接觸面積本身大小，裂縫接觸面積在方向上的變化對滲透率的控制更加關鍵。因此，針對強非均質性開度分布的裂縫，僅僅利用裂縫接觸面積變化分析有效應力作用對裂縫滲透率的影響存在局限性，需要考察裂縫接觸面積在空間方向上的變化特性，進而綜合分析對滲透率的控制作用。針對這種現象，有學者提出了無效滲流區域概念，即裂縫網格中除了接觸的這部分區域變成無效流動區域外，還將裂縫網格中局部流量低于裂縫整體平均流量1%的區域也劃分為無效流動區域，這2 類區域對裂縫整體流動的貢獻忽略不計[16]，反之，對裂縫滲流有貢獻的區域定義為有效滲流區域。因此，建議在今后的工作中，采用有效滲流區域作為關鍵參數表征裂縫導流能力隨有效應力變化的影響。

3.3 研究案例與經驗公式匹配性

為了驗證傳統模型是否可以預測通過巴西劈裂法形成的非均質性裂縫的滲透率演化行為，采取常用的ZIMMERMAN 等[19]提出的經典經驗公式（簡稱ZB模型）進行分析，該模型表達如下：

ZB模型得到的計算結果與研究數值計算結果對比，結果表明（圖9）：隨著有效應力的增大，由于相對花崗巖G1 和頁巖S1 裂縫樣品，致密砂巖裂縫TS1 開度的非均質性程度更低且相關長度更小，能在有效應力下基本滿足均勻接觸，因而能夠與經典模型較好匹配，因此，除了致密砂巖TS1的數據擬合較好，其余通過經驗公式計算的結果與本研究數據結果有較大差異。另外，雖然經驗公式計算結構與研究實際結果存在較大誤差，但是可以看出兩者結果在相對低有效應力下的趨勢是一致的，說明在低有效應力作用下，裂縫開度的非均質性對滲透率演化的影響較小。因此，傳統的經驗公式在分析預測不滿足平均化特征的非均質裂縫應力敏感下的滲透率演化時，存在明顯的不適性。

圖9 研究數據與經驗公式計算本文數據結果對比Fig.9 Comparison between results calculated by our method and ZB model

4 結論

本研究中分別利用花崗巖、致密砂巖和頁巖樣品，采取巴西劈裂法制成真實裂縫巖樣，通過三維光學掃描儀獲取裂縫上下表面微觀結構特征參數后，數值拼接形成整合裂縫，并開展了連續應力加載過程中裂縫流場和滲透率演化行為研究，取得了以下主要結論和認識：

1）采用巴西劈裂法形成的裂縫具有較強的非均質性，且不同種類的巖樣在裂縫開度空間分布上表現出不同的非均質性，這種表現會影響裂縫在有效應力下的滲流演化規律。

2）裂縫開度、接觸面積和空間相關長度是影響裂縫滲透率的重要參數。針對均質裂縫，有效應力加載過程中，裂縫滲透率降低趨勢與該參數線性關系明顯，尤其是與接觸面積密切相關。隨著裂縫非均質性結構的加強，這種線性關系逐漸失效，表現為即使有些裂縫開度沒有閉合，但未能形成滲流通道，滲透率仍然降低。

3）傳統的經驗公式在研究裂縫應力敏感時，裂縫開度分布滿足平均化條件，因而擬合程度較好。但在擬合本文利用巴西劈裂法形成的真實裂縫時，擬合程度隨有效應力增加而偏差嚴重。因此，對于存在非均質性特性或裂縫開度并沒有整體滿足平均化尺度時，采用傳統經驗公式研究裂縫導流能力變化會存在顯著的不適應性。