基于Lattice Boltzmann方法的頁巖張性裂縫滲流特征研究

曲冠政，周德勝，彭 嬌，Randy Doyle Hazlett，Siwei Wu

(1.西安石油大學，陜西 西安 710065；2.The University of Tulsa McDougall School of Petroleum Engineering，Tulsa，74104，USA)

0 引 言

頁巖儲層的開發是油氣開發領域的熱點問題[1-5]，通常采用水平井分簇壓裂技術使壓裂裂縫溝通誘導裂縫和天然裂縫，實現流體從微納米級基質到井筒的有效滲流。裂縫滲流特征及空間分布決定了儲層的滲流特征，明確裂縫中流體滲流特征是研究頁巖儲層滲流的基礎。平行板模型忽略粗糙性影響，首先被用于描述裂縫結構[6-9]，隨后意識到粗糙性對滲流的重要影響，Lomize、李士斌等[10-15]從各個角度修正或描述裂縫的粗糙性，但獲取的均是裂縫整體滲流特征，對裂縫內滲流物理場變化及特征參數分布等的認識尚不完善。目前關于裂縫內滲流特征的描述，可見報道的只有鞠楊等[10]采用人工手段計算滲流速度變化，但其裂縫物理模型是基于W-M函數生成的，本質上并非三維裂縫結構。根據巖石的破裂機理，頁巖儲層裂縫可歸納為張性裂縫和剪切滑移裂縫[16]。因此，采用物理實驗和數值模擬技術相結合的方式，明確張性裂縫結構對流體滲流的影響，為研究單裂縫結構及裂縫網絡結構滲流特征提供參考。

1 張性裂縫結構的獲取

裂縫結構的獲取和描述是研究流體在裂縫中滲流的關鍵，裂縫的獲取方法主要有巴西劈裂實驗法和函數生成法。函數生成法中，以W-M函數的應用最為廣泛，其優勢在于構造二維裂縫，但無法構造三維裂縫結構。為保證研究的可參考性，擬采用巴西劈裂實驗獲取裂縫結構。采用巴西劈裂實驗將取自美國Barnett盆地的頁巖巖心劈裂成人工裂縫，并用三維輪廓儀收集裂縫面高度分布數據；在所獲裂縫面上隨機提取表觀尺寸為2.6 mm×1.5 mm的區域用于研究分析，并在后續研究中將其劃分為離散滲流單元體。

2 流體在張性裂縫流動的模擬方法

2.1 Lattice Boltzmann模型

Lattice Boltzmann方法將流體宏觀運動離散為微觀尺度上粒子的遷移和碰撞，時間和空間的離散通過粒子離散速度聯系；在保證微觀粒子運動特征參數密度、動量和能量等的統計結果與宏觀特征參數保持一致的前提下，通過微觀粒子間的碰撞和遷移運動來推演宏觀流體運動[17-21]。在單相流的Lattice Boltzmann模型中，碰撞項采用BKG模型，離散速度模型中Qian等[17]提出的DdQm模型(d為空間維數，m為離散速度數)被廣泛采用。因此，流體在三維裂縫中的流動采用LBM-BKG方法D3Q19模型進行模擬。

2.2 準確性驗證

采用Lattice Boltzmann方法專業軟件PowerFLOW進行流體滲流的模擬計算，該軟件是目前唯一一款基于Lattice Boltzmann方法開發的計算流體力學的商業軟件。陳耀松、Wang、Du和曲冠政等的研究成果證明了該專業軟件計算的準確性和可靠性[22-25]。

2.3 模擬參數

取水作為模擬流體，溫度為293.15 K，運動黏度為1.007×10-6m2/s，流體密度為998.203 kg/m3；出口壓力設為大氣壓；格子單位下，運動黏度通過調整松弛時間來確定，模擬中運動黏度取1/30，格子單位下流體密度取0.22，格子速度為1/3，由雷諾數計算公式，可得初始格子速度為1/30。流動模擬方面設置入口流速和出口壓力邊界條件。不考慮流體壓縮性的影響。模擬裂縫體表觀尺寸為2.6 mm×1.5 mm，模擬中固定裂縫開度為0.1 mm。采用笛卡爾正交網格，為確保模擬精度，裂縫開度方向上分辨率取10。將縱向和裂縫面上分別設為光滑邊界條件和標準反彈邊界條件。

3 張性裂縫流體滲流LBM模擬結果

3.1 滲流參數

在具體表征粗糙裂縫結構中流體滲流時的主要參數包括衡量慣性力與黏性力的雷諾數，反映壓力和慣性力的歐拉數、泊肅葉數(雷諾數的函數)等[15]。以上參數均是基于流體特征參數(密度、黏度、裂縫開度)及滲流場參數(壓力和速度)計算出來的流動分析參數。由于流體密度、黏度及裂縫開度等均為定值，參數的變化主要取決于壓力和速度的分布，因此，研究集中考察壓力和速度的變化。

3.2 可視化滲流場

3.2.1 整體滲流特征

以滲流速度為0.4 m/s為例，對流體在粗糙裂縫內的滲流進行數值模擬實驗研究。由數值模擬實驗可知：當滲流速度為0.4 m/s時，整體上裂縫中流體壓力沿滲流方向呈線性下降趨勢，粗糙性并不影響壓力的線性趨勢；同一滲流縱向截面上滲流速度呈非均勻分布，呈現兩端小中間大的趨勢；裂縫面處滲流速度基本為0；粗糙性造成不同截面滲流速度差別很大，以截面中心處最大速度的差距尤為明顯；沿滲流流線上各處速度值并不相等，因流線的迂曲性而有差異。

3.2.2 三維方向上速度分布

為觀察滲流速度分布，將速度場三維剖分。在垂直于x方向選取的10條速度切片，同一切片上速度呈非均勻分布，若沿y方向上裂縫面粗糙度較小則速度分布相對均勻；若局部粗糙度較大則明顯形成滲流高速區和低速區，其原理與多孔介質滲流中流體沿大孔道突進的原理一致；在垂直于y方向均勻選取6條速度切片，速度分布形態與垂直于x方向選取的速度切片上的速度分布情況基本一致，粗糙性造成速度分布的非均勻性；在垂直于z方向的速度切片上也較明顯展示了速度分布的非均質性。綜合分析認為：裂縫中滲流穩定后，盡管滲流場宏觀速度為平均速度，但由于粗糙性影響，速度分布依據空間位置的不同而相異；粗糙度大的區域，滲流速度較小；粗糙度較小的區域，滲流速度較大。研究局部粗糙度較大的裂縫面處的流體滲流發現：裂縫面的粗糙性造成了流體滲流路徑的迂曲性；在裂縫面處，速度矢量受粗糙性影響明顯，在部分裂縫面形態顯著變化的區域表現出滲流速度的雜亂無序性，產生回流區；回流區的產生增加了流體滲流的阻力。

3.2.3 滲流速度等值面分布

考察滲流速度為0.4 m/s時裂縫中滲流速度等值面的分布：在平均滲流速度附近(0.2～0.5 m/s)時，該區間流線間是相互平行的，只有部分粗糙性變化程度較劇烈區域由于出現回流導致滲流速度的非一致性；當滲流速度偏離平均速度程度較大時，粗糙性對滲流速度等值面分布的影響較大， 滲流速度小于0.1 m/s的區域基本對應于粗糙性最大的部分，而滲流速度大于0.6 m/s的區域基本對應于粗糙性變化較為平緩的區域。

3.3 截線上滲流特征

在滲流場中隨機選取滲流截線，分析截線上的壓力和速度分布情況。隨機在z方向上高度為1.43 mm處的滲流平面上沿滲流方向和垂直滲流方向上任意選取3條水平截線(沿滲流方向上3條隨機滲流截線y分別為0.55、0.15、0.10 mm；垂直于滲流方向上的3條隨機截線x分別為0.75、1.70、2.00 mm)，分析其壓力和速度分布情況可知：壓力沿滲流方向總體呈斜率固定的線性下降趨勢且截線的斜率保持一致，但壓力曲線稍有波動并非呈純粹的線性下降趨勢；垂直于滲流方向，壓力分布隨坐標變化呈無規則波動，該模擬中，沿垂向壓力波動在1 Pa以內，約占滲流方向壓力變化的5%；沿滲流方向和垂直于滲流方向上，同一水平截線上滲流速度差別很大；粗糙性導致裂縫面高度分布的非均勻性，流線依托于裂縫面形態，裂縫形態導致滲流流線的迂曲性。

3.4 張性裂縫離散單元滲流特征

為明確張性裂縫粗糙性結構中不同區域滲流特征參數分布，將裂縫面沿滲流方向均分為5個區域，中心位置分別為0.15、0.45、0.75、1.05、1.35 mm，并用垂直于xy平面尺寸為0.2 mm×0.3 mm的65個長方形結構覆蓋滲流場，考察離散單元內的滲流特征。圖1為離散單元內滲流速度與位置的關系。由圖1可知：圖1a、b中最大滲流速度分別為0.45、2.30 m/s，分別位于2.70、0.15 mm處，最小滲流速度分別為0.31、1.54 m/s，分別位于2.70、1.35 mm處，滲流速度最大值與最小值的差距分別為31%、33%；同一滲流區域中，沿滲流方向各點的速度相異；不同滲流區域，同一滲流截面上不同位置的速度亦不同。因此，由于粗糙性影響滲流速度因空間位置而異，在垂直于滲流剖面上同樣有壓差存在，前面關于截面壓力分布的研究也證明了該觀點；模擬滲流速度為0.4、2.0 m/s時，速度分布趨勢基本一致。因此，在分析中可選取某一滲流速度下的滲流特征，從差異辨別度的角度考慮選取模擬滲流速度為2.0 m/s進行滲流特征分析。

圖1 微觀滲流速度分布

圖2為滲流區域各滲流速度矢量值。由圖2可知：x方向滲流速度矢量值與滲流總速度矢量值基本接近，而y方向滲流速度矢量值最小，基本為0；z方向滲流速度矢量值在0附近波動，波動幅度較大，如中心位置y=0.75 mm處滲流速度分布圖。分析顯示：整個滲流場沿y方向壓力略有波動，但y方向上壓差的驅動作用不足以抵消粗糙性導致的滲流場迂曲性的阻力作用，因此，y方向上滲流速度值較小；模擬滲流方向沿x方向，模擬速度為2.0 m/s，同時壓力梯度沿x方向保持恒定值，流體壓能及動能完全能克服粗糙性導致的迂曲性阻力作用及裂縫面的摩擦阻力作用等；同時，由于滲流路徑沿z方向的迂曲性必然會造成z方向上具有一定的滲流速度分量，該分量值與驅動壓差、滲流速度及裂縫面形態等相關。

圖2 滲流速度矢量圖

常用的表征裂縫面粗糙性參數包括粗糙度、迂曲度、分形維數[26-30]，粗糙度和分形維數是從裂縫面高度分布的復雜程度考慮，迂曲度是從滲流路徑考慮。考慮到數據點采集量，采用粗糙度定性描述裂縫面高度分布[29]。

中心位置分別為0.15、0.45、0.75、1.05、1.35 mm，所對應的沿滲流方向各區域粗糙度分別為0.038、0.056、0.055、0.037、0.077 mm。5個橫向滲流區域滲流特征見圖3。為考察垂直于滲流方向上的滲流特征，縱向劃分為13個相互平行的滲流區域，每個區域包含5個滲流單元，滲流特征見圖3b、d、f。由圖3a、b可知：沿滲流方向各區域滲流速度各異，分布為1.90～2.10 m/s，按照滲流速度由大至小對5個區域進行排序，所對應區域的中心位置依次為0.15、1.05、0.45、0.75、1.35 mm，5個區域粗糙度依次為0.038、0.037、0.056、0.055、0.077 mm。滲流速度基本與粗糙度呈反相關關系，個別區域稍有出入。分析表明：離散單元滲流速度具有非一致性，相鄰區域間由于速度差異而存在剪切力，同時裂縫面附近的回流區域與粗糙度及滲流路徑變化率有關，回流區對滲流速度有阻力作用；垂直于滲流方向各截面滲流速度也不一致，在1.96～2.02 m/s內波動，波動范圍比沿滲流方向波動程度小(圖3b)。

由圖3c、d可知：橫向截面壓力總體呈線性均勻分布；相對于縱向壓力分布，橫向壓力梯度很小；沿滲流方向壓差總體呈線性下降趨勢，但相鄰考察點間壓差呈波動趨勢，波動趨勢與局部粗糙性有關；總之，2個方向壓力總體均呈線性分布，滲流阻力與滲流路徑迂曲程度有關，因此，壓力梯度隨粗糙度的增加而增加。由圖3e、f可知：x方向滲流速度比滲流總速度略小，橫向截面上的壓力梯度造成橫向上具有較小波動；y、z方向滲流速度分量基本在0附近；圖3f中，沿滲流方向上，x方向滲流速度依然占據滲流速度絕對比例，各剖面y方向滲流速度基本在0附近，由于滲流路徑的迂曲性造成z方向速度分量值為-0.50～0.50 m/s。

圖3 滲流場內部區域滲流特征

4 結 論

(1) 整體上壓力沿滲流方向呈線性分布，但受粗糙性影響壓力下降趨勢具有小幅波動性；受裂縫面粗糙性影響，在垂直于滲流方向上也存在呈線性分布的壓降，但相比于滲流方向的壓降很小，可忽略。

(2) 沿滲流方向任意截線上，壓力均呈波動狀線性分布且線性分布規律保持一致；垂直于滲流方向截線上，壓力分布受裂縫面粗糙性影響具有隨機性。

(3) 沿滲流方向和垂直于滲流方向上切片各處的滲流速度不同；同一滲流流線上不同位置處滲流速度略有差別；受裂縫面粗糙性變化影響，部分區域會出現速度的雜亂無序性分布，產生回流區；沿滲流方向和垂直于滲流方向上，任選一條平行于xy平面的截線，其各處速度值具有隨機性，該分布特點與粗糙性造成的流線分布特點有關。

(4) 當速度值在平均速度附近時，滲流速度等值面間基本上相互平行；當滲流速度偏離平均速度值較大時，速度等值面分布與裂縫面形態密切相關，滲流速度較低的等值面基本分布于裂縫面形態較粗糙位置，滲流速度較高的等值面基本分布于裂縫面形態分布平緩的區域。

(5) 滲流單元體中的滲流速度各不相同。由于粗糙性造成滲流路徑迂曲性的差異，滲流單元體中的滲流總速度并不相等。滲流路徑迂曲性造成三維方向上均有滲流速度，但沿滲流方向速度分量占總速度分量的絕對比重。

[1] 蔣裕強，董大忠，漆麟，等.頁巖氣儲層的基本特征及其評價[J].天然氣工業，2010，30(10)：7-12.

[2] 侯冰，陳勉，李志猛，等.頁巖儲集層水力裂縫網絡擴展規模評價方法[J].石油勘探與開發，2014，41(6)：763-768.

[3] 張士誠，郭天魁，周彤，等.天然頁巖壓裂裂縫擴展機理試驗[J].石油學報，2014，35(3)：496-503.

[4] 吳奇，胥云，張守良，等.非常規油氣藏體積改造技術核心理論與優化設計關鍵[J].石油學報，2014，35(4)：706-714.

[5] 程遠方，李友志，時賢，等.頁巖氣體積壓裂縫網模型分析及應用[J].天然氣工業，2013，33(9)：53-59.

[6] WITHERSPOON P A，WANG J S Y，IWAI K，et al.Validity of cubic law for fluid flow in a deformable rock fracture[J].Water Resources Research，1979：16(6)：1016-1024.

[7] 溫慶志，劉欣佳，黃波，等.水力壓裂可視裂縫模擬系統的研制與應用[J].特種油氣藏，2016，23(2)：136-139，158.

[8] 溫慶志，胡藍霄，翟恒立，等.滑溜水壓裂裂縫內砂堤形成規律[J].特種油氣藏，2013，20(3)：137-139，158.

[9] 李小龍，肖雯，王凱，等.支撐劑在清潔壓裂液中的沉降規律[J].大慶石油地質與開發，2015，34(2)：95-98.

[10] LOMIZE G M.Water flow through jointed rock [M].Moscow：Gesenergoizdat，1951：122.

[11] STEPHEN R Brown.Transport of fluid and electric current through a single fracture[J].Journal of Geophysical Research，1989，94(B7)：9429-9438.

[12] ERDINC Eker，SERHAT Akin.Lattice Boltzmann simulation of fluid flow in synthetic fractures[J].Transport in Porous Media，2006，65(3)：363-384.

[13] 張程賓，陳永平，施明恒，等.端面粗糙度的分形特征及其對微通道內層流流動的影響[J].物理學報，2009，58(10)：7050-7056.

[14] 李士斌，陳波濤，張海軍，等.清水壓裂自支撐裂縫面閉合殘留寬度數值模擬[J].石油學報，2010，31(4)：680-683.

[15] 鞠楊，張欽剛，楊永明，等.巖體粗糙單裂隙流體滲流機制的實驗研究[J].中國科學 E輯：技術科學，2013，43(10)：1144-1154.

[16] 鐘建華，劉圣鑫，馬寅生，等.頁巖宏觀破裂模式與微觀破裂機理[J].石油勘探與開發，2015，42(2)：242-249.

[17] QIAN Y H.Lattice BGK models for Navier-Stokes Equation [J].Europhysics Letters，1992，17(6)：479-484.

[18] 何雅玲，王勇，李慶，等 .格子Boltzmann方法的理論與應用[M].北京：科學出版社，2009：98-100.

[19] 郭照立，鄭楚光.格子Boltzmann方法的原理及應用[M].北京：科學出版社，2009：121-123.

[20] ZHAO Liguo，Zhao T S.Lattice boltzmann model for incompressible flows through porous media[J].Physical review E，2002，66(3)：162-173.

[21] YOSHINO M，INAMURO T.Lattice boltzmann simulation for flow and heat/mass transfer problems in a three-dimensional porous structure[J].International journal for numberical methods in fluids，2003，43(2)：183-198.

[22] 陳耀松，單肖文，陳滬東，等.計算流體力學的新方向及其在工業上的應用[J].中國科學 E輯：技術科學，2007，37(9)：1107-1116.

[23] WANG Yiwei，WANG Yang，AN Yiran，et al.LBM calculation of high-speed train aerodynamics[J].Science China Technological Sciences，2008，38(11)：1795-1804.

[24] DU Tezhuan，LI Xiangqun，ZHANG Xialiang，et al.Lattice Boltzmann Method used for the aircraft characteristics computation at high angle of attack[J].Science China Technological Sciences，2010，53(8)：2068-2073.

[25] 曲冠政，曲占慶，HAZLETT R D，等.頁巖拉張型微裂縫幾何特征描述及滲透率計算[J].石油勘探與開發，2016，43(1)：115-120.

[26] BELEM T，HOMAND-ETIENNE F.Quantitative parameters for rock joint surface roughness[J].Rock Mechanics and Rock Engineering，2000，33(4)：217-242.

[27] ZHANG Zhenyu，JAN Nemcik，QIAO Qiuqiu，et al.A model for water flow through rock fractures based on friction factor [J].Rock Mechanics and Rock Engineering，2015，48(2)：559-571.

[28] 閆寒，張文明，胡開明，等.隨機粗糙微通道內流動特性研究[J].物理學報，2013，62(17)：174701-174706.

[29] LAMAS L N.An experimental and analytical study of the roughness of granite joints[J].Hydraulic Engineering，2013，33(4)：240-248.

[30] MAKURAT A，GUITERREZ M.Fracture flow and fracture cross flow experiments[C].SPE36732，1996：1-9.