曲冠政, 彭 嬌, 趙 凱, 高 偉, RANDY Hazlett,周德勝 ,曲占慶, 姜海龍

(1.西安石油大學(xué)石油工程學(xué)院,陜西西安 710065; 2.中國(guó)石油長(zhǎng)慶油田分公司油氣工藝研究院,陜西西安 610065;3.University of Tulsa,Tulsa USA 74104; 4.中國(guó)石油大學(xué)(華東)石油工程學(xué)院,山東青島 266580)

在描述裂縫結(jié)構(gòu)時(shí)通常將其等效為平板模型,但實(shí)際裂縫滲透率與平板模型結(jié)果偏差較大。目前關(guān)于裂縫描述及其滲流規(guī)律的研究還處于探索階段[1-4];同時(shí),現(xiàn)場(chǎng)裂縫具有大尺度匹配、小尺度非匹配的特點(diǎn)[5-6]?；诖?以頁(yè)巖為例,統(tǒng)計(jì)并描述其高度分布特點(diǎn),建立裂縫結(jié)構(gòu)模型,生成滿(mǎn)足研究條件的非匹配裂縫結(jié)構(gòu),采用Lattice Boltzmann方法,結(jié)合正交試驗(yàn)和灰色關(guān)聯(lián)分析系統(tǒng)地研究流體在三維非匹配性裂縫中的滲流規(guī)律。

1 裂縫模型構(gòu)造方法

1.1 Werierstrass-Mandelbrot函數(shù)法

Werierstrass-Mandelbrot函數(shù)(W-M函數(shù))在描述二維自相似結(jié)構(gòu)中具有廣泛應(yīng)用[7-9]?；诹芽p剖面線(xiàn)具有多尺度、自仿射的特征,引入W-M函數(shù)描述二維裂縫剖面結(jié)構(gòu),基本形式為

(1)

式中,z(x)為x位置處的高度,mm;G為尺度系數(shù);D為自仿射分形維數(shù);γ為與裂縫面輪廓頻譜密度有關(guān)的量，一般取1.5;n1為描述W-M函數(shù)下限截止頻率,其值可由γn1=1/Ls給出;Ls為測(cè)量樣本長(zhǎng)度,mm。

粗糙度計(jì)算方法為

(2)

式中,σ為粗糙度,mm;nl為下限截止頻率，與樣本長(zhǎng)度有關(guān)，nl=1/L,mm-1;nh為上限截止頻率,由測(cè)量?jī)x器精度Lr決定,mm-1。

確定分形維數(shù)后,利用粗糙度計(jì)算公式求出尺度系數(shù),結(jié)合W-M函數(shù)即可獲得具有任意粗糙度和分形維數(shù)的二維粗糙曲線(xiàn),將剖面線(xiàn)縱向延伸獲取裂縫結(jié)構(gòu)[7-9]。

1.2 基于Gaussian分布的裂縫構(gòu)造方法

采用巴西劈裂試驗(yàn)將取自Barnett頁(yè)巖的巖心Cabana-ZG-1、DARADO-ZG-1、punta-piedra-ZG-1 和Dorado-Sur12劈裂成人工裂縫,掃描裂縫面數(shù)據(jù)并重構(gòu),如圖1所示。裂縫面高度Gaussian分布擬合如圖2所示。

圖1 頁(yè)巖裂縫面重構(gòu)圖Fig.1 Rough fracture surface reconstitution

圖2中藍(lán)色部分為裂縫面高度分布直方圖,紅色曲線(xiàn)為依據(jù)裂縫面高度分布數(shù)據(jù)計(jì)算出的Gaussian正態(tài)分布曲線(xiàn)。4塊頁(yè)巖裂縫面高度分布特征與計(jì)算出的Gaussian分布曲線(xiàn)間有出入,但總體上符合Gaussian分布規(guī)律。考慮到巖心尺寸和測(cè)量手段等對(duì)結(jié)果的干擾,結(jié)合裂縫面高度分布的實(shí)際情況,可認(rèn)為其分布符合Gaussian分布規(guī)律。

根據(jù)頁(yè)巖裂縫面高度分布特點(diǎn),采用Gaussian分布函數(shù)構(gòu)造裂縫結(jié)構(gòu)。裂縫結(jié)構(gòu)中,特征描述參數(shù)分為裂縫開(kāi)度和裂縫面粗糙性描述參數(shù)。生成裂縫面結(jié)構(gòu)后,裂縫開(kāi)度可任意設(shè)定。因此關(guān)鍵問(wèn)題是裂縫面的構(gòu)造。構(gòu)建裂縫面時(shí),需要考慮的主要因素為裂縫面高度概率分布及其在三維空間的展布規(guī)律[10-11]。分別采用Gaussian分布函數(shù)和自相關(guān)函數(shù)描述裂縫面高度概率分布和高度在三維空間的展布規(guī)律。裂縫結(jié)構(gòu)示意圖如圖3所示。

圖3中裂縫結(jié)構(gòu)由上下兩個(gè)裂縫面組成,裂縫面在局部可能會(huì)接觸。兩裂縫面可分別視為在平均高度平面h0+和h0-波動(dòng),兩個(gè)平面與xy平面平行。面Sp+和Sp-可表示為

(3)

裂縫開(kāi)度定義為

b(x)=z+(x)-z-(x) .

(4)

式中,b(x)為裂縫開(kāi)度,mm;z+(x)為上裂縫面x處高度,mm;z-(x)為下裂縫面x處高度,mm。

圖2 裂縫面高度Gaussian分布擬合圖Fig.2 Fitted Gaussian distribution of fractures surface height

圖3 裂縫結(jié)構(gòu)示意圖Fig.3 Sketch map of fracture structure

裂縫面的波動(dòng)隨其所在位置不同而差異很大。由于波動(dòng)是基于平均面,因此裂縫面波動(dòng)平均值為0,即

(5)

裂縫面間平均距離定義為

(6)

式中,bm為裂縫面間平均距離,mm。

采用Gaussian分布函數(shù)描述其高度分布概率,表達(dá)式為

(7)

(8)

式中,P(H)和p(h)分別為概率函數(shù)和概率密度函數(shù);h為裂縫高度,mm。

自相關(guān)函數(shù)是描述信號(hào)間關(guān)聯(lián)程度的函數(shù)。對(duì)于三維裂縫面,可表示為

(9)

式中,Ch+(x,y)和Ch-(x,y)為自相關(guān)函數(shù),兩者相等,mm;h+(x,y)和h-(x,y)為裂縫面高度分布,mm;u和v為滯后,mm;l為周期,mm。

相比于常規(guī)數(shù)域中求解問(wèn)題的局限性,在傅里葉域求解函數(shù)具有明顯優(yōu)勢(shì);自相關(guān)函數(shù)在傅里葉域中對(duì)應(yīng)于功率譜密度函數(shù),將信號(hào)采用不同時(shí)間域或空間域正弦信號(hào)的組合。正弦信號(hào)具有不同的波長(zhǎng)、振幅和相位,振幅的平方為功率,振幅與波數(shù)的關(guān)系曲線(xiàn)為功率譜,功率譜的正交化即為功率譜密度。功率密度函數(shù)譜矩能有效描述裂縫面粗糙度。功率密度函數(shù)譜矩[12]定義為

(10)

其中

k0=2π/λ0.

式中,mn為n階矩;k0為對(duì)應(yīng)于剖面線(xiàn)長(zhǎng)度λ0的波數(shù);S(k)為功率譜密度函數(shù),即自相關(guān)函數(shù)的Fourier變換。

確定Gaussian概率分布函數(shù)和自相關(guān)函數(shù)(功率譜密度函數(shù))后即可建立裂縫結(jié)構(gòu)的基本模型。功率譜密度函數(shù)與波數(shù)間關(guān)系[12]表達(dá)式為

S(k)=Ck2D-7.

(11)

考慮粗糙度沿裂縫面兩方向上的分布規(guī)律,引入各向異性系數(shù)b/a度量其非均質(zhì)性。當(dāng)b/a<1時(shí),非均質(zhì)性垂直于x方向;當(dāng)b/a=1時(shí),粗糙度為均質(zhì)分布;當(dāng)b/a>1時(shí),非均質(zhì)性平行于x方向。具體表達(dá)式為

S(k)=C[(kx/a)2+(ky/b)2]D-3.5.

(12)

采用Gaussian概率分布函數(shù)和自相關(guān)函數(shù)分別描述粗糙度分布概率、裂縫面粗糙度間的關(guān)聯(lián)性,并引入分形維數(shù)、非均質(zhì)性系數(shù)描述粗糙度在裂縫面的展布情況。以上函數(shù)及相關(guān)參數(shù)已足以描述拉張型粗糙裂縫結(jié)構(gòu),但研究表明:實(shí)際的粗糙裂縫兩面結(jié)構(gòu)間存在一定的非匹配性,傅里葉域中,長(zhǎng)波長(zhǎng)情況下兩裂縫面是匹配的;而當(dāng)波長(zhǎng)小于某一閥值時(shí),裂縫結(jié)構(gòu)表現(xiàn)為非匹配性。Brown[5]首次提出非匹配長(zhǎng)度的概念,Glover等[6]在Brown研究的基礎(chǔ)上發(fā)展了非匹配性表征方面的理論:在傅里葉域中,引入臨界波長(zhǎng)λc衡量其匹配情況。當(dāng)λ<λc/2時(shí),采用兩個(gè)獨(dú)立的隨機(jī)數(shù)據(jù)源(R1、R2)生成裂縫面,在此波長(zhǎng)范圍內(nèi)兩裂縫面完全非匹配;當(dāng)λ>λc/2時(shí),引入R1和R2的線(xiàn)性組合R3形成隨機(jī)數(shù)據(jù)源,表達(dá)式為

(13)

式中,R1、R2和R3為隨機(jī)數(shù)據(jù)源;λ為波長(zhǎng);β為匹配因子,代表最長(zhǎng)波長(zhǎng)時(shí)最大程度的匹配度;kc為非匹配波長(zhǎng)對(duì)應(yīng)的波數(shù)。

明確各參數(shù)后,以R1為生成上裂縫面的隨機(jī)數(shù)據(jù)源,R3為生成下裂縫面的隨機(jī)數(shù)據(jù)源。生成的兩裂縫面結(jié)構(gòu)分形維數(shù)、粗糙度、各向異性分布均一致,只有非匹配長(zhǎng)度存在差異。根據(jù)以上理論,依照研究需要可任意設(shè)置裂縫面分形維數(shù)、粗糙度、各向異性系數(shù)及非匹配長(zhǎng)度等參數(shù)。

考慮到計(jì)算資源問(wèn)題,生成裂縫結(jié)構(gòu)的表觀尺寸限定為3 mm×3 mm。根據(jù)研究需要共生成20組不同的裂縫面結(jié)構(gòu),部分粗糙裂縫面結(jié)構(gòu)如圖4所示。

圖4 粗糙裂縫面Fig.4 Rough surface morphology

圖5為生成的部分裂縫結(jié)構(gòu)及其開(kāi)度分布,裂縫面分形維數(shù)均為2.2,粗糙度0.1 mm;其中圖5(a)、(b)各向異性系數(shù)為1,圖5(c)、(d)各向異性系數(shù)為3。各裂縫結(jié)構(gòu)由于設(shè)定的裂縫結(jié)構(gòu)參數(shù)的差異性,導(dǎo)致所生成的裂縫平均開(kāi)度各異。

為直觀展現(xiàn)不同裂縫面間的差異,選取裂縫結(jié)構(gòu)某一剖面線(xiàn)分布定性考察其分布特點(diǎn),鞠楊等[2]和張程賓等[3]采用Weierstrass-Mandelbrot函數(shù)生成與圖6類(lèi)似的滿(mǎn)足粗糙度和分形維數(shù)的粗糙剖面線(xiàn)并將其橫向延伸,沿縱向影像得到匹配性裂縫結(jié)構(gòu);鄭宏等[4]則采取分形插值法生成滿(mǎn)足設(shè)定分形維數(shù)的粗糙裂縫結(jié)構(gòu),并截取某一剖面線(xiàn)并影像進(jìn)行粗糙裂縫滲流的研究,其在研究中也沒(méi)有考慮非匹配性的影響。

圖6為沿x=2 mm處截面y方向上各裂縫結(jié)構(gòu)剖面線(xiàn)分布。當(dāng)非匹配長(zhǎng)度為0.3 mm時(shí),沿y方向上裂縫上下面吻合很好,極個(gè)別點(diǎn)處裂縫開(kāi)度略有偏移;當(dāng)非匹配長(zhǎng)度為0.7 mm、可明顯見(jiàn)到部分區(qū)域裂縫面的非一致性,但從整體尺度上看兩裂縫結(jié)構(gòu)面基本保持一致,該情況下的裂縫開(kāi)度相比較與非匹配長(zhǎng)度0.3 mm情況下明顯增加,裂縫開(kāi)度的波動(dòng)幅度較之非匹配長(zhǎng)度0.3 mm情況增加;當(dāng)非匹配長(zhǎng)度為0.7 mm、裂縫面各向異性系數(shù)為3時(shí),相比于各向異性系數(shù)為1情況下,裂縫平均開(kāi)度略有增加,同時(shí)裂縫開(kāi)度波動(dòng)幅度明顯增加;而當(dāng)非匹配長(zhǎng)度為2.1 mm時(shí),裂縫面間差異明顯,裂縫間距較之非匹配長(zhǎng)度較小的情況明顯增加,波動(dòng)幅度亦明顯增加。>

圖5 構(gòu)造的裂縫結(jié)構(gòu)Fig.5 Constructed fracture structures

圖6 x=2 mm處裂縫結(jié)構(gòu)剖面線(xiàn)Fig.6 Fracture profiles along x=2 mm

2 滲流模擬方法

研究流體滲流問(wèn)題的方法歸納起來(lái)包括物理實(shí)驗(yàn)法、理論研究方法及數(shù)值模擬方法。在處理微尺度問(wèn)題上,物理實(shí)驗(yàn)方法受實(shí)驗(yàn)成本、儀器精度、測(cè)量手段及實(shí)驗(yàn)條件等影響很大;理論研究方法通常采取各種簡(jiǎn)化條件,研究條件與流體滲流條件有較大出入;數(shù)值模擬方法可根據(jù)研究需要設(shè)置物理模型,并結(jié)合相應(yīng)滲流理論對(duì)研究目標(biāo)重點(diǎn)分析,能夠克服物理實(shí)驗(yàn)方法的缺陷,因此在研究目標(biāo)較為復(fù)雜且對(duì)實(shí)驗(yàn)要求較高時(shí),數(shù)值模擬方法是較好的選擇?？紤]到研究目標(biāo)的特征尺度及裂縫面的復(fù)雜性,采用Lattice Boltzmann方法模擬流體在裂縫中的滲流[13]。

2.1 Lattice Boltzmann模型

Lattice Boltzmann方法基于分子熱運(yùn)動(dòng)和統(tǒng)計(jì)理論,將無(wú)規(guī)則的分子熱運(yùn)動(dòng)簡(jiǎn)化為格子規(guī)則化的碰撞和遷移運(yùn)動(dòng)。Lattice Boltzmann方法中最關(guān)鍵的因素是格子運(yùn)動(dòng)規(guī)律和平衡態(tài)分布函數(shù)。在單相流模擬方面,LBM-BGK模型被廣泛接受;三維流體運(yùn)移模擬中,以Qian等[14]提出的D3Q19模型應(yīng)用最為廣泛。該模型演化方程、平衡態(tài)分布函數(shù)及粒子權(quán)重系數(shù)表達(dá)式分別為

(14)

(15)

(16)

式中,cs為格子聲速,模型中為(1/3)0.5;ei為i方向格子速度,i=0,1,…,18;feq為D3Q19模型平衡態(tài)分布函數(shù);fi為粒子密度分布函數(shù);r為粒子空間位置;ωi為i方向上的粒子權(quán)重系數(shù);τl為松弛時(shí)間;δ為時(shí)間步長(zhǎng)。

流體宏觀密度和速度分別為

(17)

(18)

邊界條件設(shè)置上,根據(jù)在裂縫中流體滲流的實(shí)際情況,將縱向上邊界設(shè)為光滑邊界條件;考慮到裂縫面的粗糙性,在裂縫面上采用標(biāo)準(zhǔn)反彈邊界。

根據(jù)量綱等效原則,裂縫滲透率和格子滲透率關(guān)系式為

(19)

式中,kf為裂縫滲透率,μm2;k為格子滲透率;r為格子分辨率。

采用Lattice Boltzmann方法專(zhuān)業(yè)軟件PowerFLOW進(jìn)行相關(guān)模擬計(jì)算。陳耀松等[15]、Wang等[16]、Du等[17]的研究成果證明了軟件的準(zhǔn)確性和可靠性。采用平板模型和圓管模型驗(yàn)證計(jì)算方法的可靠性,理論值與模擬值吻合很好。

2.2 計(jì)算網(wǎng)格劃分

空間網(wǎng)格采用笛卡爾正交網(wǎng)格,PowerFLOW軟件中是以分辨率的形式表示網(wǎng)格結(jié)構(gòu)。經(jīng)計(jì)算當(dāng)格子分辨率大于10以后,模擬結(jié)果基本保持不變。因此選取格子分辨率為10。

2.3 流體運(yùn)移模擬

首先固定分形維數(shù)、粗糙度和各向異性系數(shù),建立不同非匹配長(zhǎng)度的裂縫結(jié)構(gòu)模型,進(jìn)行流體滲流模擬;然后,從所建立的非匹配性裂縫模型中,以一定間隔選取其中4個(gè)裂縫結(jié)構(gòu),對(duì)比兩面完全吻合和非匹配情況下的流體滲流特征;此外,綜合考察分形維數(shù)、粗糙度和各向異性系數(shù)等裂縫特征參數(shù)對(duì)流體滲流的影響。

采用水作為模擬流體,溫度20 ℃,密度0.998×103kg/m3,運(yùn)動(dòng)黏度1.006×10-6m2/s,采取入口速度,出口壓力作為邊界條件。入口速度設(shè)為1.0 m/s,出口壓力采用標(biāo)準(zhǔn)大氣壓,即1.013 25×105Pa。為保證模擬結(jié)果的收斂,模擬時(shí)間設(shè)為20 000時(shí)間步長(zhǎng)。為保證模擬中流動(dòng)為達(dá)西滲流,模擬中恒定雷諾數(shù)為10。

3 模擬結(jié)果分析

3.1 非匹配長(zhǎng)度影響

如表1所示,共進(jìn)行20組非匹配長(zhǎng)度影響的模擬,模擬中裂縫表觀尺寸3 mm×3 mm,分形維數(shù)2.2,粗糙度0.1 mm,各向異性系數(shù)均為1。表1中的裂縫開(kāi)度為控制兩裂縫面不接觸所需的最小平均開(kāi)度。隨非匹配長(zhǎng)度增加,維持裂縫不接觸所需的平均開(kāi)度總體呈增加趨勢(shì)。

表1 非匹配裂縫模擬數(shù)據(jù)

圖7為裂縫滲透率與非匹配長(zhǎng)度的關(guān)系曲線(xiàn)。由圖7可得:裂縫滲透率隨非匹配長(zhǎng)度增加呈迅速增長(zhǎng)并在后期趨于穩(wěn)定的趨勢(shì)。對(duì)比理想平板模型滲透率,在相同開(kāi)度情況下裂縫滲透率的變化趨勢(shì)基本與理想平板模型一致;粗糙性的影響導(dǎo)致裂縫滲透率偏離理想平板模型結(jié)果且偏差較大。

圖7 裂縫滲透率與非匹配長(zhǎng)度關(guān)系Fig.7 Relationship between fracture permeability and mismatch length

考慮粗糙性影響的裂縫滲透率計(jì)算公式[18-20]概括起來(lái)為

(20)

式中,τ為迂曲度;A和B為系數(shù)。

式(20)中采用粗糙度和迂曲度修正裂縫粗糙性的影響,但研究中并未考慮非匹配性的影響。參考裂縫面匹配的裂縫滲透率計(jì)算模型,修正非匹配性的影響公式為

(21)

式中,k′為非匹配性裂縫滲透率,μm2;α為裂縫面非匹配性綜合修正系數(shù)。

非匹配性修正系數(shù)反映了非匹配性對(duì)裂縫滲透率的影響程度,該值越大則表明非匹配裂縫滲透率越大。非匹配長(zhǎng)度與非匹配性綜合修正系數(shù)關(guān)系如圖8所示。受粗糙性影響,實(shí)際裂縫滲透率約為理想平板模型滲透率的20%～50%;以裂縫結(jié)構(gòu)表觀尺寸3 mm為界,當(dāng)非匹配長(zhǎng)度小于3 mm時(shí),非匹配性綜合修正系數(shù)隨非匹配長(zhǎng)度變化有較大幅度的波動(dòng),當(dāng)非匹配長(zhǎng)度大于3 mm后,非匹配性綜合修正系數(shù)基本趨于穩(wěn)定。

圖8 非匹配長(zhǎng)度與非匹配性綜合修正系數(shù)關(guān)系Fig.8 Relationship between mismatch length and mismatch comprehensive coefficient

3.2 非匹配性裂縫與匹配性裂縫結(jié)構(gòu)對(duì)比

固定分形維數(shù)2.2、粗糙度0.05 mm、各向異性系數(shù)1.0,考察非匹配長(zhǎng)度0.3、1.1、1.9和2.7 mm情況下非匹配性裂縫結(jié)構(gòu)及采用非匹配性裂縫結(jié)構(gòu)兩裂縫面分別組成匹配性裂縫結(jié)構(gòu)等3種情況下的裂縫滲流能力,結(jié)果如圖9所示。兩種匹配性裂縫結(jié)構(gòu)的滲流能力基本保持一致,而非匹配性裂縫的滲透率明顯比匹配性裂縫滲透率低。

如圖10所示,裂縫結(jié)構(gòu)可視為由眾多平板模型串聯(lián)而成。流體在裂縫中的滲流阻力,一方面來(lái)自于眾多平板模型基本單元體結(jié)構(gòu)的滲流阻力,另一方面來(lái)自于相鄰平板模型基本單元體結(jié)構(gòu)間的銜接滲流阻力。銜接區(qū)域形態(tài)的突然改變更容易出現(xiàn)如圖11 中裂縫面改變率較大區(qū)域形成的回流區(qū),回流區(qū)的出現(xiàn)極大地增加了流體的滲流阻力。

對(duì)于非匹配性裂縫結(jié)構(gòu),同樣將其等效為平板模型結(jié)構(gòu)分析其滲流阻力:除具有平板模型基本單元體和基本單元體間的銜接滲流阻力外,還包括非匹配區(qū)域的滲流阻力;由于結(jié)構(gòu)間的非匹配性,在非匹配區(qū)域裂縫結(jié)構(gòu)形態(tài)變化率較大,造成非匹配性裂縫滲透率低于匹配性裂縫滲透率。

圖9 裂縫結(jié)構(gòu)滲透率對(duì)比Fig.9 Permeability comparison among different fracture structures

圖10 裂縫等效平板模型 Fig.10 Equivalent parallel plate model of rough fracture

圖11 裂縫滲流速度矢量Fig.11 Velocity vector in rough fracture

3.3 非匹配性裂縫滲透率影響因素敏感性分析

正交試驗(yàn)設(shè)計(jì)可定性確定各因素對(duì)考察目標(biāo)的影響程度,最顯著特點(diǎn)是以具有代表性的有限個(gè)方案反映大量方案中所包含的本質(zhì)規(guī)律和矛盾主次[36]。非匹配性裂縫結(jié)構(gòu)中,裂縫開(kāi)度是影響滲透率的最主要因素,其他影響因素有非匹配長(zhǎng)度、粗糙度、分形維數(shù)和各向異性系數(shù)。其中由于裂縫開(kāi)度與非匹配長(zhǎng)度存在正相關(guān)關(guān)系,裂縫開(kāi)度的波動(dòng)受迫于非匹配長(zhǎng)度的變化。因此,進(jìn)行正交試驗(yàn)設(shè)計(jì)時(shí)兩個(gè)因素中只選取非匹配長(zhǎng)度,即考慮非匹配長(zhǎng)度、粗糙度、分形維數(shù)和各向異性系數(shù),采用四因素三水平正交試驗(yàn)設(shè)計(jì)確定模擬方案如表2所示,表3為計(jì)算的裂縫滲透率及其基本參數(shù)數(shù)據(jù)。

表2 正交試驗(yàn)設(shè)計(jì)因素

表3 Lattice Boltzmann方法模擬滲透率Table 3 Permeability value based on LatticeBoltzmann simulation

為全面分析各因素對(duì)裂縫滲透率的影響程度,在進(jìn)行相關(guān)灰色關(guān)聯(lián)分析計(jì)算時(shí)考慮裂縫開(kāi)度的影響?；疑P(guān)聯(lián)分析方法對(duì)樣本容量和分布的規(guī)律性均沒(méi)有要求,并且分析結(jié)果與定性分析結(jié)果相吻合,不會(huì)出現(xiàn)量化分析結(jié)果與定性分析結(jié)果不符的情況。因此以上試驗(yàn)設(shè)計(jì)方案和灰色關(guān)聯(lián)分析所采用數(shù)據(jù)方案不會(huì)影響分析結(jié)果。具體步驟如下[21]。

將裂縫滲透率視為目標(biāo)數(shù)列:X0={X0(s)|s=1,2,3,…,9},非匹配長(zhǎng)度、粗糙度、分形維數(shù)、各向異性系數(shù)、裂縫開(kāi)度視為因素?cái)?shù)列Xi={Xi(s)|s=1,2,3,…,9},其中s為因素取值個(gè)數(shù),i為因素?cái)?shù)列個(gè)數(shù)。為消除量綱影響,將因素?cái)?shù)列和目標(biāo)數(shù)列進(jìn)行均值化無(wú)量綱處理,經(jīng)過(guò)均值化方法處理的各指標(biāo)數(shù)據(jù)構(gòu)成的協(xié)方差矩陣既可以反映原始數(shù)據(jù)中各指標(biāo)變異程度上的差異,也包含各指標(biāo)相互影響程度差異的信息。數(shù)據(jù)無(wú)量綱化處理后,計(jì)算因素?cái)?shù)列Xi(s)與目標(biāo)數(shù)列X0(s)的關(guān)聯(lián)系數(shù),表達(dá)式為

(22)

式中,ρ為分辨系數(shù),取0.5;Δi(s)=|X0(s)-Xi(s)|為第某時(shí)刻目標(biāo)數(shù)列與因素?cái)?shù)列的絕對(duì)值差。

根據(jù)式(22)可求得因素?cái)?shù)列Xi與目標(biāo)數(shù)列X0的關(guān)聯(lián)系數(shù)序列ξi={ξi(s)|s=1,2,3,…,9},利用平均值法可得各因素的關(guān)聯(lián)度γi為

(23)

經(jīng)計(jì)算,各因素的關(guān)聯(lián)系數(shù)為0.58～0.82;裂縫開(kāi)度的影響程度最大,其關(guān)聯(lián)系數(shù)為0.813 3,是影響非匹配性裂縫滲透率的主要因素,但并不占絕對(duì)主導(dǎo)地位;其他依次為粗糙度0.743 7,分形維數(shù)0.671 4,各向異性系數(shù)0.660 4和非匹配長(zhǎng)度0.586 8;非匹配長(zhǎng)度、粗糙度、分形維數(shù)和各向異性系數(shù)對(duì)非匹配性裂縫滲透率的影響均不可忽略;粗糙度,即裂縫面高度分布偏離程度為影響非匹配性裂縫滲透率的次主要因素;分形維數(shù),即粗糙度在裂縫面分布規(guī)律的影響位于第三位;最后為非匹配長(zhǎng)度和各向異性系數(shù),兩者對(duì)非匹配性裂縫滲透率的影響程度相近。

在目前的壓裂開(kāi)發(fā)設(shè)計(jì)中,通常只考慮裂縫開(kāi)度與儲(chǔ)層油氣滲流條件的匹配性。而灰色關(guān)聯(lián)分析結(jié)果顯示:裂縫面粗糙度及粗糙度在裂縫面的整體分布規(guī)律(分形維數(shù))及橫縱向分布規(guī)律(各向異性系數(shù))、非匹配長(zhǎng)度均對(duì)裂縫滲流具有不可忽略的影響。相比于兩面完全吻合的裂縫結(jié)構(gòu),在裂縫平均開(kāi)度相同情況下,匹配性裂縫的滲透率比非匹配性裂縫的略大;隨粗糙度和分形維數(shù)增加,裂縫滲流阻力明顯增加,滲透率降低。目前的壓裂設(shè)計(jì)中,研究主要集中于所形成裂縫結(jié)構(gòu)的縫長(zhǎng)和縫寬等宏觀尺度參數(shù),而對(duì)于裂縫結(jié)構(gòu)面參數(shù)對(duì)滲流能力的研究較少。因此在以后的研究中,建議加強(qiáng)儲(chǔ)層巖石物性特征(泊松比、彈性模量等)及壓裂泵注程序等對(duì)裂縫面粗糙性影響方面的研究。在保證裂縫開(kāi)度的情況下,盡可能避免使裂縫結(jié)構(gòu)產(chǎn)生非匹配性,同時(shí)降低裂縫面的粗糙度和分形維數(shù)。

3.4 各向異性系數(shù)的影響

固定分形維數(shù)2.2、粗糙度0.10 mm情況下,考察各向異性系數(shù)與粗糙裂縫滲透率的關(guān)系,結(jié)果如圖12所示。隨各向異性系數(shù)增加,滲透率呈先上升后下降趨勢(shì);在裂縫面粗糙度分布各向同性時(shí),模擬裂縫滲透率值并不是最大;實(shí)例中,當(dāng)各向異性系數(shù)為0.7時(shí),模擬滲透率值最大,為402.45 μm2。

圖12 各向異性系數(shù)與裂縫滲透率關(guān)系Fig.12 Relationship between anisotropy coefficient and permeability

4 結(jié) 論

(1)采用統(tǒng)計(jì)學(xué)公式Gaussian分布函數(shù)描述裂縫面結(jié)構(gòu)分布特點(diǎn),考慮分形維數(shù)、非匹配長(zhǎng)度、粗糙度、各向異性等因素,建立更符合實(shí)際的真三維裂縫結(jié)構(gòu)模型。實(shí)例中,裂縫滲透率僅為理想平板模型滲透率的20%～50%,裂縫面的粗糙性對(duì)流體滲流具有較大的阻力作用。

(2)采用上、下面結(jié)構(gòu)影像得到的匹配性裂縫滲透率值基本保持一致,非匹配性裂縫滲透率比匹配性裂縫滲透率值略小。將裂縫結(jié)構(gòu)視為平板模型結(jié)構(gòu)的串聯(lián)組合,分析匹配性裂縫和非匹配性裂縫滲流阻力的差別。非匹配性裂縫銜接區(qū)域形態(tài)的突然改變更容易出現(xiàn)回流區(qū)域,回流區(qū)域增加了流體的滲流阻力。

(3)裂縫開(kāi)度和粗糙度是影響非匹配性裂縫滲透率的最主要因素;其他依次為分形維數(shù)、各向異性系數(shù)和非匹配長(zhǎng)度。粗糙度分布規(guī)律同樣對(duì)非匹配性裂縫滲透率具有重要作用。

(4)實(shí)例中,隨各向異性系數(shù)增加,粗糙裂縫滲透率呈先上升后下降趨勢(shì),且當(dāng)各向異性系數(shù)取0.7時(shí),粗糙裂縫滲透率取最大值。

:

[1] 速寶玉,詹美禮,趙堅(jiān).仿天然巖石裂隙滲流的實(shí)驗(yàn)研究[J].巖土工程學(xué)報(bào),1995,17(5):19-24.

SU Baoyu, ZHAN Meili, ZHAO Jian. Study of fracture seepage in the imitative nature rock [J].Chinese Journal of Geotechnical Engineering,1995,17(5):19-24.

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

JU Yang, ZHANG Qingang, YANG Yongming, et al. An experimental investigation on the mechanism of fluid flow through single fracture of rock[J].Science China (Technological Science),2013,43(10):1144-1154.

[3] 張程賓,陳永平,施明恒,等. 表面粗糙度的分形特征及其對(duì)微通道內(nèi)層流流動(dòng)的影響[J]. 物理學(xué)報(bào),2009,58(10):7050-7056.

ZHANG Chengbin, CHEN Yongping, SHI Mingheng, et al. Fractal characteristics of surface roughness and its effect on laminar flow in microchannels[J]. Acta Physica Sinica, 2009,58(10):7050-7056.

[4] 樊火,鄭宏.基于MRT-LBM的分形裂隙網(wǎng)絡(luò)滲流數(shù)值模擬[J].中國(guó)科學(xué)(技術(shù)科學(xué)),2013,43(12):1338-1345.

FAN Huo, ZHENG Hong. MRT-LBM based numerical simulation of seepage flow through fractal fracture networks[J].Science China(Technological Science), 2013,43(12):1338-1345.

[5] BROWN S R. Simple mathematical model of a rough fracture[J]. Journal of Geophysical Research,1995,100(B4):5941-5952.

[6] GLOVER P W J, MATSUKI K, HAYASHI K, et.al. Synthetic rough fractures in rocks[J].Journal of Geophysical Research,1998,103(85):9609-9620.

[7] GANTI S, BHUSHAN B. Generalized fractal analysis and its applications to engineering surfaces[J].Wear,1995(180):17-34.

[8] MAJUMDAR A, TIEN C L. Fractal characterization and simulation of rough surfaces[J].Wear,1990(136):313-327.

[9] THOMA T R, ROSEN B G, AMINI N. Fractal characterization of the anisotropy of rough surfaces[J].Wear,1999(232):41-50.

[10] BROWN S R, SCHOLZ H. Closure of rock joints[J]. Journal of Geophysical Research,1986,9(B5):4939-4948.

[11] BROWN S R, SCHOLZ C H. Broad bandwidth study of the topography of natural rock surfaces[J].Journal of Geophysical Research,1985,90(12):575-582.

[12] POWER W L, TULLS T E. Euclidian and fractal models for the description of rock surface roughness[J].Journal of Geophysical Research,1991,96(B1):415-424.

[13] 張磊,姚軍,孫海,等.利用格子Boltzmann方法計(jì)算頁(yè)巖滲透率[J].中國(guó)石油大學(xué)學(xué)報(bào)(自然科學(xué)版),2014,38(1):87-91.

ZHANG Lei, YAO Jun, SUN Hai, et.al. Permeability calculation in shale using Lattice Boltzmann method[J]. Journal of China University of Petroleum(Edition of Natural Science),2014,38(1):87-91.

[14] QIAN Y H. Lattice BGK models for Navier-Stokes equation[J].Europhysics Letters,1992,17(6):479-484.

[15] 陳耀松,單肖文,陳滬東,等.計(jì)算流體力學(xué)的新方向及其在工業(yè)上的應(yīng)用[J].中國(guó)科學(xué)(技術(shù)科學(xué)),2007,37(9):1107-1116.

CHEN Yaosong, SHAN Xiaowen,CHEN Hudong, et al. New direction of computational fluid dynamics and its application on industry[J].Science China(Technological Sciences), 2007,37(9):1107-1116.

[16] 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.

[17] DU Tezhuan, LI Xiangqun, ZHANG Xialing, 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.

[18] LOUIS C. A study of groundwater flow in jointed rock and its influence on stability of rock masses[D].London:Imperial College of Science and Technology,1969.

[19] LOUIS. Rock hydraulics in rock mechanics[M].New York:Springer-Verlag,1974.

[20] XIAO Weimin, XIA Caichu, WANG Wei, et al. Combined effect of tortuosity and surface roughness on estimation of flow rate through a single rough joint[J].Journal of Geophysics and Engineering,2013,10(4):45015-45023.

[21] 梁濤,常毓文,郭曉飛,等.巴肯致密油藏單井產(chǎn)能參數(shù)影響程度排序[J].石油勘探與開(kāi)發(fā),2013,40(3):357-362.

LIANG Tao, CHANG Yuwen, GUO Xiaofei, et al. Influence factors of single wells productivity in the Bakken tight oil reservoir[J].Petroleum Exploration and Development,2013,40(3):357-362.