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

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

0 引 言

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

1 張性裂縫結構的獲取

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

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

2.1 Lattice Boltzmann模型

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

2.2 準確性驗證

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

2.3 模擬參數(shù)

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

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

3.1 滲流參數(shù)

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

3.2 可視化滲流場

3.2.1 整體滲流特征

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

3.2.2 三維方向上速度分布

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

3.2.3 滲流速度等值面分布

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

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)，分析其壓力和速度分布情況可知：壓力沿滲流方向總體呈斜率固定的線性下降趨勢且截線的斜率保持一致，但壓力曲線稍有波動并非呈純粹的線性下降趨勢；垂直于滲流方向，壓力分布隨坐標變化呈無規(guī)則波動，該模擬中，沿垂向壓力波動在1 Pa以內，約占滲流方向壓力變化的5%；沿滲流方向和垂直于滲流方向上，同一水平截線上滲流速度差別很大；粗糙性導致裂縫面高度分布的非均勻性，流線依托于裂縫面形態(tài)，裂縫形態(tài)導致滲流流線的迂曲性。

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

為明確張性裂縫粗糙性結構中不同區(qū)域滲流特征參數(shù)分布，將裂縫面沿滲流方向均分為5個區(qū)域，中心位置分別為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%；同一滲流區(qū)域中，沿滲流方向各點的速度相異；不同滲流區(qū)域，同一滲流截面上不同位置的速度亦不同。因此，由于粗糙性影響滲流速度因空間位置而異，在垂直于滲流剖面上同樣有壓差存在，前面關于截面壓力分布的研究也證明了該觀點；模擬滲流速度為0.4、2.0 m/s時，速度分布趨勢基本一致。因此，在分析中可選取某一滲流速度下的滲流特征，從差異辨別度的角度考慮選取模擬滲流速度為2.0 m/s進行滲流特征分析。

圖1 微觀滲流速度分布

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

圖2 滲流速度矢量圖

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

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

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

圖3 滲流場內部區(qū)域滲流特征

4 結 論

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

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

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

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

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

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

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

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

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

[5] 程遠方，李友志，時賢，等.頁巖氣體積壓裂縫網(wǎng)模型分析及應用[J].天然氣工業(yè)，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] 溫慶志，劉欣佳，黃波，等.水力壓裂可視裂縫模擬系統(tǒng)的研制與應用[J].特種油氣藏，2016，23(2)：136-139，158.

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

[9] 李小龍，肖雯，王凱，等.支撐劑在清潔壓裂液中的沉降規(guī)律[J].大慶石油地質與開發(fā)，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] 李士斌，陳波濤，張海軍，等.清水壓裂自支撐裂縫面閉合殘留寬度數(shù)值模擬[J].石油學報，2010，31(4)：680-683.

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

[16] 鐘建華，劉圣鑫，馬寅生，等.頁巖宏觀破裂模式與微觀破裂機理[J].石油勘探與開發(fā)，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] 陳耀松，單肖文，陳滬東，等.計算流體力學的新方向及其在工業(yè)上的應用[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].石油勘探與開發(fā)，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.