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

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

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

0 引 言

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

1 張性裂縫結(jié)構(gòu)的獲取

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

2 流體在張性裂縫流動(dòng)的模擬方法

2.1 Lattice Boltzmann模型

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

2.2 準(zhǔn)確性驗(yàn)證

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

2.3 模擬參數(shù)

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

3 張性裂縫流體滲流LBM模擬結(jié)果

3.1 滲流參數(shù)

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

3.2 可視化滲流場(chǎng)

3.2.1 整體滲流特征

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

3.2.2 三維方向上速度分布

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

3.2.3 滲流速度等值面分布

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

3.3 截線上滲流特征

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

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

為明確張性裂縫粗糙性結(jié)構(gòu)中不同區(qū)域滲流特征參數(shù)分布，將裂縫面沿滲流方向均分為5個(gè)區(qū)域，中心位置分別為0.15、0.45、0.75、1.05、1.35 mm，并用垂直于xy平面尺寸為0.2 mm×0.3 mm的65個(gè)長方形結(jié)構(gòu)覆蓋滲流場(chǎng)，考察離散單元內(nèi)的滲流特征。圖1為離散單元內(nèi)滲流速度與位置的關(guān)系。由圖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ū)域中，沿滲流方向各點(diǎn)的速度相異；不同滲流區(qū)域，同一滲流截面上不同位置的速度亦不同。因此，由于粗糙性影響滲流速度因空間位置而異，在垂直于滲流剖面上同樣有壓差存在，前面關(guān)于截面壓力分布的研究也證明了該觀點(diǎn)；模擬滲流速度為0.4、2.0 m/s時(shí)，速度分布趨勢(shì)基本一致。因此，在分析中可選取某一滲流速度下的滲流特征，從差異辨別度的角度考慮選取模擬滲流速度為2.0 m/s進(jìn)行滲流特征分析。

圖1 微觀滲流速度分布

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

圖2 滲流速度矢量圖

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

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

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

圖3 滲流場(chǎng)內(nèi)部區(qū)域滲流特征

4 結(jié) 論

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

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

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

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

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

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

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

[3] 張士誠，郭天魁，周彤，等.天然頁巖壓裂裂縫擴(kuò)展機(jī)理試驗(yàn)[J].石油學(xué)報(bào)，2014，35(3)：496-503.

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

[5] 程遠(yuǎn)方，李友志，時(shí)賢，等.頁巖氣體積壓裂縫網(wǎng)模型分析及應(yī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)的研制與應(yīng)用[J].特種油氣藏，2016，23(2)：136-139，158.

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

[9] 李小龍，肖雯，王凱，等.支撐劑在清潔壓裂液中的沉降規(guī)律[J].大慶石油地質(zhì)與開發(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] 張程賓，陳永平，施明恒，等.端面粗糙度的分形特征及其對(duì)微通道內(nèi)層流流動(dòng)的影響[J].物理學(xué)報(bào)，2009，58(10)：7050-7056.

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

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

[16] 鐘建華，劉圣鑫，馬寅生，等.頁巖宏觀破裂模式與微觀破裂機(jī)理[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方法的理論與應(yīng)用[M].北京：科學(xué)出版社，2009：98-100.

[19] 郭照立，鄭楚光.格子Boltzmann方法的原理及應(yīng)用[M].北京：科學(xué)出版社，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ì)算流體力學(xué)的新方向及其在工業(yè)上的應(yīng)用[J].中國科學(xué) E輯：技術(shù)科學(xué)，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ì)算[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ī)粗糙微通道內(nèi)流動(dòng)特性研究[J].物理學(xué)報(bào)，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.