孔隙尺度各向異性與孔隙分布非均質(zhì)性對多孔介質(zhì)滲透率的影響機(jī)理

李滔，李閩，荊雪琪，肖文聯(lián)，崔慶武

（1.西南石油大學(xué)油氣藏地質(zhì)及開發(fā)工程國家重點(diǎn)實(shí)驗(yàn)室，成都 610500；2.中石化新星（北京）新能源開發(fā)有限公司四川分公司，成都 610500；3.中石化中原石油工程有限公司鉆井一公司，河南濮陽 457000）

0 引言

多孔介質(zhì)的固有滲透率（文中簡稱滲透率）由其孔隙結(jié)構(gòu)決定[1-2]。在油氣開采領(lǐng)域，儲集層滲透率對儲集層評價(jià)和開發(fā)有著重要影響[3]。學(xué)者采用解析法[4-5]、實(shí)驗(yàn)法[6-10]和模擬法[11-14]，建立了多種滲透率與孔隙結(jié)構(gòu)參數(shù)的關(guān)系式，多孔介質(zhì)滲透率主要與孔隙度（或連通孔隙度）、孔徑、迂曲度、比表面積和固相顆粒大小等有關(guān)?，F(xiàn)有滲透率模型存在較大差異，且只適用于各向同性、相對均質(zhì)的常規(guī)儲集層[15]。與常規(guī)儲集層相比，非常規(guī)儲集層（致密砂巖、頁巖等）巖石的孔隙結(jié)構(gòu)具有以下特點(diǎn)：①微納米級孔隙發(fā)育[16]，孔隙形態(tài)復(fù)雜[17-20]；②孔隙尺度各向異性強(qiáng)[21-23]；③孔隙分布非均質(zhì)性強(qiáng)[17,24-25]?？紫冻叨雀飨虍愋耘c孔隙分布非均質(zhì)性均顯著影響多孔介質(zhì)滲透率[22-23,26-27]，是非常規(guī)儲集層不同于常規(guī)儲集層的重要特征。Gu等[28]通過中子散射提出孔隙尺度各向異性的計(jì)算模型，但模型中包含了需測量的中子散射強(qiáng)度，極大限制了該模型的可應(yīng)用性。隨后，Wang等[26]基于形態(tài)學(xué)原理提出了具有較強(qiáng)適用性的多孔介質(zhì)各向異性和孔隙分布非均質(zhì)性計(jì)算模型。

目前對多孔介質(zhì)各向異性和孔隙分布非均質(zhì)性的定量研究相對缺乏，孔隙尺度各向異性與孔隙分布非均質(zhì)性影響多孔介質(zhì)滲透率的機(jī)理認(rèn)識不清[21]，基于微CT圖像重構(gòu)三維數(shù)字巖心定量評價(jià)儲集層孔隙尺度各向異性和孔隙分布非均質(zhì)性的研究也未見文獻(xiàn)報(bào)道。隨著成像技術(shù)的發(fā)展，利用物理實(shí)驗(yàn)方法可直接觀測多孔介質(zhì)的微納米級孔隙及其連通關(guān)系。采用四參數(shù)隨機(jī)生成（QSGS）算法能夠較高效地構(gòu)建孔隙尺度各向異性、孔隙分布非均質(zhì)性強(qiáng)的多孔介質(zhì)，且能反映天然多孔介質(zhì)絕大部分的形態(tài)復(fù)雜性[29]。同時，采用格子-玻爾茲曼方法（LBM），可直觀、高效地處理流體內(nèi)部及流-固間的相互作用[30]，實(shí)現(xiàn)多孔介質(zhì)孔隙尺度的模擬[26,31]。這些均為非常規(guī)儲集層滲透率影響機(jī)理的研究提供了條件。本文結(jié)合微CT掃描實(shí)驗(yàn)、QSGS算法和LBM，進(jìn)一步揭示孔隙尺度各向異性與孔隙分布非均質(zhì)性影響多孔介質(zhì)滲透率的機(jī)理。

1 致密砂巖微觀孔隙結(jié)構(gòu)性質(zhì)

1.1 三維數(shù)字巖心的建立

①四川盆地某氣藏的4塊致密砂巖巖心（H6-1、H6-2、H106-1和H106-2）的鑄體薄片、掃描電鏡和X-衍射實(shí)驗(yàn)結(jié)果表明，其孔隙空間由晶間微孔隙、粒間孔隙、粒內(nèi)/粒間溶蝕微孔隙和粒間微縫等組成；成分以巖屑長石砂巖、長石巖屑砂巖為主，多為中粒、細(xì)—中粒砂巖，次為中—粗砂巖，少量為粉砂巖；膠結(jié)物主要為伊蒙混層、伊利石和黏土雜基。

②從4塊巖樣的端面鉆取近似小圓柱（實(shí)際物理尺寸見表1），并對其進(jìn)行微CT掃描（掃描測試電壓150 kV，電流90 μA，掃描時間120 min，圖像分辨率620 nm），掃描次序與巖樣軸向一致。采用專業(yè)軟件對巖樣二維CT圖像進(jìn)行中值過濾和二值化處理，準(zhǔn)確獲取巖石骨架和孔隙空間，隨后經(jīng)圖像切割完成微觀結(jié)構(gòu)成像。圖1展示了H6-1和H106-1巖樣單元體（800像素×800像素×800像素，對應(yīng)實(shí)際物理尺寸0.496 mm×0.496 mm×0.496 mm）二值化前后以及孔隙空間的三維圖像?？梢钥吹街旅苌皫r孔隙形態(tài)復(fù)雜且極不規(guī)則。

表1 柱體巖樣實(shí)際物理尺寸

③對單元體進(jìn)行數(shù)值重構(gòu)，建立致密砂巖的三維數(shù)字巖心，并計(jì)算得到4個單元體的孔隙度分別為2.60%、4.49%、4.51%和2.61%。

1.2 各向異性和非均質(zhì)性的定量評價(jià)

基于已建立的三維數(shù)字巖心，采用Wang等[26]提出的多孔介質(zhì)各向異性模型：

和孔隙分布非均質(zhì)性模型：

定量評價(jià)致密砂巖的孔隙尺度各向異性與孔隙分布非均質(zhì)性。

孔隙尺度各向異性因子被定義為多孔介質(zhì)中所有孔隙的寬度之和與所有孔隙高度之和的比值，且各向異性強(qiáng)弱與A值成正比（A=1時表示各向同性）；孔隙分布非均質(zhì)性因子表示單元塊孔隙度的相對標(biāo)準(zhǔn)偏差（取決于觀測尺度，即單元塊大?。?，非均質(zhì)性的強(qiáng)弱與H值成正比。

圖1 基于CT圖像獲取的致密砂巖巖心三維圖像

圖2 孔隙尺度各向異性與數(shù)字巖心立方體邊長的關(guān)系

圖3 孔隙分布非均質(zhì)性與數(shù)字巖心立方體邊長的關(guān)系

選取線條間隔為5個格子，單元塊長、寬、高均為25個格子，計(jì)算得到4塊致密砂巖的孔隙尺度各向異性和孔隙分布非均質(zhì)性與數(shù)字巖心立方體邊長的關(guān)系（見圖2、圖3）。圖中可以看出盡管4塊致密砂巖三維數(shù)字巖心的孔隙度存在一定差異，但孔隙尺度各向異性、孔隙分布非均質(zhì)性的變化卻較為一致。圖2顯示，當(dāng)立方體邊長達(dá)到300個像素點(diǎn)時，4塊致密砂巖的孔隙尺度各向異性趨于穩(wěn)定，x、y和z方向（以巖樣軸向?yàn)閤，水平面內(nèi)垂直于x的方向?yàn)閥，垂直于x、y構(gòu)成平面的方向?yàn)閦，建立空間直角坐標(biāo)系）的A值分別為1.38～1.55、0.70～0.77和0.91～1.02，巖樣軸向的各向異性因子最大（這里將各向異性因子最大的方向定義為主流方向），其余2個方向相差不大。圖3顯示，立方體邊長達(dá)到600個像素點(diǎn)，致密砂巖孔隙分布非均質(zhì)性趨于穩(wěn)定，H值范圍為3.29～3.74。

2 三維隨機(jī)多孔介質(zhì)的LBM模擬

由以上分析可知，微CT圖像能夠較準(zhǔn)確地反映巖樣孔隙空間的三維形貌以及孔隙間的配置情況，但受分辨率限制，微CT掃描通常難以完整獲取非常規(guī)儲集層中的納米級孔隙，由此建立的數(shù)字巖心孔隙空間不連通。對此，可在致密砂巖孔隙尺度各向異性、孔隙分布非均質(zhì)性分析結(jié)果的基礎(chǔ)上，采用QSGS方法構(gòu)建三維各向異性、非均質(zhì)性多孔介質(zhì)模型，并運(yùn)用三維19個離散速度方向（D3Q19）的多弛豫時間格子-玻爾茲曼（MRT-LB）模型模擬其中流體的運(yùn)移，研究對滲透率的影響機(jī)理。

2.1 多孔介質(zhì)構(gòu)建

2.1.1 QSGS方法

QSGS方法構(gòu)建三維多孔介質(zhì)的具體過程[29]可分為3步：①設(shè)定構(gòu)造區(qū)域與隨機(jī)分布固相生長核（Cs），Cs的概率（Pcd）不能大于最終構(gòu)建三維多孔介質(zhì)的孔隙度（φ），即Pcd≤φ；利用隨機(jī)數(shù)發(fā)生器生成[0，1]間的隨機(jī)數(shù)，隨機(jī)數(shù)小于Pcd的節(jié)點(diǎn)被指定為固相生長核；②給定不同方向的生長速度（Pk，其中k代表固相生長方向），利用隨機(jī)數(shù)發(fā)生器生成[0，1]間的隨機(jī)數(shù)，當(dāng)Cs節(jié)點(diǎn)的隨機(jī)數(shù)小于Pk時，則沿k方向的鄰點(diǎn)生長；共設(shè)計(jì)26個生長方向，可分為側(cè)線方向（6個）、對角線方向（12個）和直徑線方向（8個）；當(dāng)各方向的生長速度相同時，可獲得各向同性多孔介質(zhì)；當(dāng)各方向的生長速度不同時，可獲得各向異性多孔介質(zhì)；③重復(fù)步驟②直到達(dá)到設(shè)定的多孔介質(zhì)孔隙度。

為更好地研究孔隙分布非均質(zhì)性對三維多孔介質(zhì)滲透性的影響，可采用2個相互獨(dú)立的生成過程來構(gòu)建多孔介質(zhì)[26]：①細(xì)結(jié)構(gòu)生成較小的固相，固相生長核概率為Pr，孔隙度為Rφ；②粗結(jié)構(gòu)生成較大的固相，固相生長核概率為Pc，且Pc遠(yuǎn)小于Pr；2個生成過程相結(jié)合，任一過程占據(jù)的節(jié)點(diǎn)均被設(shè)為固相，直至達(dá)到設(shè)定的孔隙度。

2.1.2 孤立孔隙刪除

QSGS方法構(gòu)建的三維多孔介質(zhì)中包括孤立孔隙（或稱死孔隙）和連通孔隙，連通孔隙是流體運(yùn)移的通道。對四連通區(qū)域標(biāo)記算法[32]進(jìn)行改進(jìn)可以去除三維多孔介質(zhì)中的孤立孔隙：①掃描三維多孔介質(zhì)中的孔隙和固相，將孔隙標(biāo)記為0；②檢查每個孔隙像素點(diǎn)東、南、西、北、上和下6個方向的鄰近孔隙像素點(diǎn)，如果鄰近的孔隙像素點(diǎn)均未被標(biāo)記，則給該像素點(diǎn)設(shè)定一個新的標(biāo)記數(shù)，新標(biāo)記數(shù)從1開始遞增，且不重復(fù)；如果只有一個鄰近孔隙像素點(diǎn)被標(biāo)記，則將孔隙像素點(diǎn)設(shè)定為與其相同的標(biāo)記數(shù)；如果兩個或以上的鄰近孔隙像素點(diǎn)已被標(biāo)記，則將孔隙像素點(diǎn)設(shè)定為已標(biāo)記鄰近孔隙像素點(diǎn)中最小的標(biāo)記數(shù)；③檢查三維多孔介質(zhì)中所有的孔隙像素點(diǎn)，將相互連通的相鄰孔隙像素點(diǎn)均使用其中最小的標(biāo)記數(shù)標(biāo)記；④檢查各標(biāo)記數(shù)所代表的孔隙空間是否同時與三維多孔介質(zhì)出入口端相連。

2.2 LBM模型及驗(yàn)證

與單松弛（SRT）模型相比，多松弛（MRT）模型的碰撞算子中引入了多個松弛時間，計(jì)算精度更高且穩(wěn)定性更好，也更適用于復(fù)雜多孔介質(zhì)滲流問題的研究。MRT-LB模型的演化方程為[30]：

式中M是19×19的轉(zhuǎn)換矩陣[33]，M-1是M的逆矩陣，可由高斯列主元消去法求得，S為對角矩陣：

式中非零矩陣元表示不同物理量（如：能量矩、能量平方矩、能量通量矩和應(yīng)力張量矩等）的弛豫速率，取值分別為[34]：

D3Q19模型中，平衡態(tài)分布函數(shù)可表示為[34]：

式中宏觀流體的密度、速度和運(yùn)動黏度可分別表示為：

模擬中流體流動均為層流（Re值較小），結(jié)合達(dá)西定律可推導(dǎo)得三維多孔介質(zhì)的滲透率計(jì)算公式[1]：

三維多孔介質(zhì)的迂曲度等價(jià)計(jì)算公式[35]為：

為驗(yàn)證以上LBM算法的準(zhǔn)確性，運(yùn)用D3Q19 MRT和D3Q19 SRT模型（τR取1）分別模擬了體心立方中的流體運(yùn)移。體心立方的孔隙度設(shè)為58.2%，邊長設(shè)為320 nm，網(wǎng)格數(shù)分別為16×16×16、32×32×32、48×48×48、64×64×64、80×80×80、96×96×96和100×100×100。采用（8）式計(jì)算滲透率，并將結(jié)果與體心立方的滲透率解析解[36]作對比（見圖4）可知，MRT模型模擬結(jié)果的準(zhǔn)確性高于SRT模型；由于體心立方中球體半徑與理論值始終存在一定的偏差[26]，隨著網(wǎng)格數(shù)的增加，模擬結(jié)果只能越來越接近解析解而不能完全一致，即使在16×16×16網(wǎng)格條件下，MRT模型模擬結(jié)果與解析解的相對偏差最大僅為4.09%，說明MRT模型可準(zhǔn)確模擬三維多孔介質(zhì)中的流體運(yùn)移。

圖4 LBM模擬滲透率與解析解結(jié)果對比

2.3 模擬結(jié)果

2.3.1 三維孔隙尺度各向異性多孔介質(zhì)的模擬

結(jié)合致密砂巖孔隙尺度各向異性的分析結(jié)果，同時為避免孔隙形態(tài)顯著變化對模擬結(jié)果的影響，將x與y、z方向設(shè)定為不同的固相生長速度（其中y和z方向相等），采用QSGS算法構(gòu)建三維孔隙尺度各向異性多孔介質(zhì)（基礎(chǔ)參數(shù)見表2）。為消除孔隙結(jié)構(gòu)隨機(jī)性對多孔介質(zhì)滲透率的影響，同一組參數(shù)均構(gòu)建10個隨機(jī)多孔介質(zhì)模型，共計(jì)180個。

表2 三維各向異性多孔介質(zhì)模型QSGS算法與MRT-LB模擬中的基礎(chǔ)參數(shù)

刪除其中孤立孔隙，將線條間隔設(shè)為4個格子，運(yùn)用（1）式計(jì)算多孔介質(zhì)的孔隙尺度各向異性因子，得到多孔介質(zhì)模型連通孔隙度、比表面積與孔隙尺度各向異性的關(guān)系（見圖5、圖6）。由圖可知，多孔介質(zhì)的各向異性因子介于0.98～1.62；各向異性因子由0.98增加至1.62時，多孔介質(zhì)比表面積增加0.33%～2.55%，連通孔隙度減小1.31%～2.00%；隨著Pcd的增加，多孔介質(zhì)中小孔隙增多，比表面積顯著增加（3種Pcd取值對應(yīng)的比表面積均值分別為4.1×105cm2/cm3，4.7×105cm2/cm3，5.3×105cm2/cm3）；連通孔隙度受Pcd的影響相對較小，與Pcd呈負(fù)相關(guān)。

圖5 各向異性模型比表面積與各向異性因子的關(guān)系

圖6 各向異性模型連通孔隙度與各向異性因子的關(guān)系

采用D3Q19 MRT-LB模型模擬三維各向異性多孔介質(zhì)中流體運(yùn)移（基礎(chǔ)參數(shù)見表2），模擬過程中的模擬偏差可表示為：

當(dāng)模擬偏差小于10-6時，認(rèn)為模擬結(jié)果不再變化，結(jié)束迭代過程。

圖7為Pcd=0.002 5，A=1.55時各向異性多孔介質(zhì)的模擬結(jié)果。以各向異性多孔介質(zhì)的模擬結(jié)果為基礎(chǔ)，采用（8）式、（9）式分別計(jì)算三維各向異性多孔介質(zhì)模型的滲透率和迂曲度。分析滲透率、各向異性、迂曲度間的相互關(guān)系（見圖8—圖10）可知，各向異性多孔介質(zhì)模型的滲透率與各向異性因子呈正相關(guān)；滲透率與迂曲度負(fù)相關(guān)；迂曲度與孔隙尺度各向異性因子呈負(fù)相關(guān)；各向異性因子由0.98增加至1.62時，Pcd為0.002 5，0.005 0和0.010 0的多孔介質(zhì)滲透率分別增加80.20%，60.97%和57.38%，迂曲度分別減少17.72%，11.02%和6.06%。

圖7 三維各向異性多孔介質(zhì)模擬結(jié)果

圖8 各向異性模型滲透率與各向異性因子的關(guān)系

圖9 各向異性模型滲透率與迂曲度的關(guān)系

圖10 各向異性模型迂曲度與各向異性因子的關(guān)系

這里需要說明的是：①數(shù)字巖心模擬計(jì)算多孔介質(zhì)滲透率所得結(jié)果與單位網(wǎng)格長度密切相關(guān)，單位網(wǎng)格長度增加1個數(shù)量級時，多孔介質(zhì)的滲透率至少增加2個數(shù)量級；②非常規(guī)儲集層納米級孔隙發(fā)育，為更好研究孔隙的微觀結(jié)構(gòu)及其與滲透率的相關(guān)性，模擬中選取了較小的單位網(wǎng)格長度（10 nm），滲透率模擬計(jì)算結(jié)果與實(shí)際巖心的測量值差異較大，但滲透率與孔隙微觀結(jié)構(gòu)的相關(guān)性與單位網(wǎng)格長度的取值無關(guān)。

2.3.2 三維孔隙分布非均質(zhì)性多孔介質(zhì)的模擬

采用QSGS算法生成三維非均質(zhì)多孔介質(zhì)模型（基礎(chǔ)參數(shù)見表3）。這里x、y和z方向均設(shè)定相同的固相生長速度，忽略各向異性的影響，同時網(wǎng)格數(shù)和單位網(wǎng)格長度保持與三維各向異性多孔介質(zhì)模型一致。為消除孔隙結(jié)構(gòu)隨機(jī)性的影響，同一組參數(shù)下均構(gòu)建10個隨機(jī)多孔介質(zhì)模型，共計(jì)100個。構(gòu)建完成后刪除其中孤立孔隙，將單元塊的長、寬和高均設(shè)為16個格子，采用（2）式計(jì)算多孔介質(zhì)孔隙分布非均質(zhì)性因子，分析連通孔隙度、比表面積與孔隙分布非均質(zhì)性的關(guān)系（見圖11、圖12）可知，不同細(xì)結(jié)構(gòu)孔隙度多孔介質(zhì)的孔隙分布非均質(zhì)性因子的計(jì)算結(jié)果分布范圍有差距，當(dāng)Rφ=40%時，孔隙分布非均質(zhì)性因子相對較小，計(jì)算結(jié)果分布范圍為0.89～1.13；當(dāng)Rφ=60%時，非均質(zhì)性較強(qiáng)，計(jì)算結(jié)果分布范圍為1.03～1.42；多孔介質(zhì)的連通孔隙度受非均質(zhì)性的影響較弱，比表面積與非均質(zhì)性呈現(xiàn)明顯的負(fù)相關(guān)關(guān)系；當(dāng)Rφ分別取40%和60%時，非均質(zhì)性因子在其分布范圍內(nèi)增加，多孔介質(zhì)連通孔隙度變化幅度均小于4.57%，而比表面積變化幅度較大，分別減小了21.55%和17.76%。

表3 三維非均質(zhì)多孔介質(zhì)模型QSGS算法中的基礎(chǔ)參數(shù)

圖11 非均質(zhì)模型連通孔隙度與非均質(zhì)性因子的關(guān)系

圖12 非均質(zhì)模型比表面積與非均質(zhì)性因子的關(guān)系

運(yùn)用D3Q19 MRT-LB模型模擬非均質(zhì)多孔介質(zhì)中流體運(yùn)移（基礎(chǔ)參數(shù)見表2）。以非均質(zhì)模型的模擬結(jié)果為基礎(chǔ)，運(yùn)用（8）式和（9）式分別計(jì)算三維非均質(zhì)性多孔介質(zhì)的滲透率和迂曲度，分析滲透率、非均質(zhì)性、比表面積、迂曲度間的相互關(guān)系（見圖13—圖16）可知，滲透率與孔隙分布非均質(zhì)性因子、細(xì)結(jié)構(gòu)孔隙度呈正相關(guān)[26]；Rφ越大，多孔介質(zhì)中較大的孔隙越多，滲透率越大；當(dāng)Rφ分別取40%和60%時，非均質(zhì)性因子在其分布范圍內(nèi)增加，多孔介質(zhì)滲透率分別增加136.01%和344.53%；滲透率與比表面積無明顯相關(guān)性；滲透率與迂曲度呈明顯的負(fù)相關(guān)；迂曲度與非均質(zhì)性因子無明顯的相關(guān)性。

3 分析與討論

圖13 非均質(zhì)模型滲透率與非均質(zhì)性因子的關(guān)系

圖14 非均質(zhì)性模型滲透率與比表面積的關(guān)系

圖15 非均質(zhì)性模型滲透率與迂曲度的關(guān)系

三維隨機(jī)多孔介質(zhì)模擬中，三維各向異性多孔介質(zhì)中固相顆粒量級均為單位網(wǎng)格大小[37]，而非均質(zhì)性多孔介質(zhì)中，粗結(jié)構(gòu)固相生長核概率減小，固相顆粒的大小逐漸增加。在此基礎(chǔ)上，主要研究了滲透率、比表面積、迂曲度、連通孔隙度等物性參數(shù)與孔隙尺度各向異性、孔隙分布非均質(zhì)性的關(guān)系。

這里進(jìn)一步研究孔隙尺度各向異性、孔隙分布非均質(zhì)性對隨機(jī)多孔介質(zhì)孔隙大小和分布的影響，以及各向異性、非均質(zhì)性對復(fù)雜多孔介質(zhì)滲透率的影響機(jī)理。采用最大球法獲取多孔介質(zhì)的孔隙直徑、不同孔隙占孔隙空間的體積比例。最大球法獲取三維多孔介質(zhì)孔隙特征參數(shù)的步驟包括：①確定最大球半徑的相關(guān)概率；②尋找最大球；③刪除冗余球（包括單個及多個包含關(guān)系）等。即使多孔介質(zhì)中孔徑大小和分布均相同，孔隙間的配置關(guān)系也會導(dǎo)致滲透率的改變[1]，為消除這一影響，同一組基礎(chǔ)參數(shù)均構(gòu)建了10個隨機(jī)多孔介質(zhì)模型。

圖16 非均質(zhì)性模型迂曲度與非均質(zhì)性因子的關(guān)系

多孔介質(zhì)中流體的運(yùn)移主要發(fā)生在較大的連通孔隙中，因此采用最大球法抽取時不考慮只包含單個節(jié)點(diǎn)的孔隙。圖17、圖18為各向異性（Pcd=0.002 5）、非均質(zhì)性（Rφ=60%）多孔介質(zhì)模型的孔徑分布，橫坐標(biāo)為孔隙直徑（以網(wǎng)格數(shù)量表示），縱坐標(biāo)為同一組基礎(chǔ)參數(shù)下10個隨機(jī)多孔介質(zhì)模型孔徑分布中對應(yīng)孔隙體積所占比重的均值。分析圖17可知，不同生長速度（對應(yīng)各向異性）條件下多孔介質(zhì)的孔徑分布形態(tài)十分接近，相對偏差為1.78%～8.64%；圖18顯示出不同粗結(jié)構(gòu)固相生長核概率（對應(yīng)非均質(zhì)性）條件下多孔介質(zhì)的孔徑分布形態(tài)也十分接近，相對偏差為1.82%～6.09%；相對偏差大小與各向異性或非均質(zhì)性強(qiáng)弱無關(guān)。

圖17 各向異性模型的孔徑分布

圖18 非均質(zhì)性模型的孔徑分布

多孔介質(zhì)的孔隙結(jié)構(gòu)參數(shù)與孔隙尺度各向異性、孔隙分布非均質(zhì)性的關(guān)系（見表4）可歸納為：①三維各向異性多孔介質(zhì)模型中，Pcd相同，其連通孔隙度、孔徑、比表面積和固相顆粒大小等幾乎不受各向異性影響或無明顯相關(guān)性，但與迂曲度呈明顯負(fù)相關(guān)。②三維非均質(zhì)性多孔介質(zhì)模型中，Rφ相同，其連通孔隙度、孔徑受非均質(zhì)性的影響較小，其值波動幅度不大；迂曲度與非均質(zhì)性雖無明顯相關(guān)性，但其值波動幅度明顯；多孔介質(zhì)非均質(zhì)性與固相顆粒大小有關(guān)，進(jìn)而影響比表面積，比表面積與非均質(zhì)性呈明顯負(fù)相關(guān)。

表4 孔隙結(jié)構(gòu)參數(shù)與各向異性、非均質(zhì)性的相關(guān)性

多孔介質(zhì)的孔隙尺度各向異性顯著影響迂曲度，孔隙分布非均質(zhì)性顯著影響比表面積，進(jìn)而顯著影響多孔介質(zhì)中流體的流動動態(tài)：①主流方向各向異性因子越大，多孔介質(zhì)中孔隙長軸沿主流方向的取向性更強(qiáng)（孔隙長軸更傾向于朝向主流方向）[21]，流動路徑更趨平直，迂曲度顯著減小，流體流過多孔介質(zhì)所消耗的能量減少，多孔介質(zhì)的滲流能力增強(qiáng)，迂曲度與各向異性的強(qiáng)相關(guān)是各向異性影響滲透率的根本原因。②孔隙分布非均質(zhì)性增強(qiáng)，多孔介質(zhì)中固相顆粒尺寸越大，孔隙聚集現(xiàn)象更加明顯[26]，孔隙間發(fā)生相互重疊的概率增加，比表面積將顯著減小；通常多孔介質(zhì)孔隙度相近，而固相顆粒越大或比表面積越小時，流體流過多孔介質(zhì)所受到的壁面摩擦阻力越小[38]，多孔介質(zhì)的滲透率將越大。③雖圖14顯示三維非均質(zhì)性多孔介質(zhì)模型的滲透率與比表面積無明顯相關(guān)性，但總的趨勢表現(xiàn)為滲透率隨非均質(zhì)性的增加而增加；這是因?yàn)橛厍纫彩怯绊懚嗫捉橘|(zhì)滲透率的重要因素[39]，迂曲度與多孔介質(zhì)中固相顆粒的排列方向和位置（尤其是其中較大的固相顆粒）有關(guān)[2]，三維非均質(zhì)性多孔介質(zhì)模型中，隨著非均質(zhì)性的增強(qiáng)，生成的固相顆粒間差異越大，分布位置的隨機(jī)性可能導(dǎo)致迂曲度隨機(jī)變化，可增大也可減小，導(dǎo)致迂曲度與非均質(zhì)性無明顯相關(guān)性；孔隙分布非均質(zhì)性對多孔介質(zhì)滲透率的影響是比表面積和迂曲度共同作用的結(jié)果。④采用迂曲度與比表面積的乘積分析其與非均質(zhì)性的關(guān)系（見圖19）顯示，其與非均質(zhì)性呈明顯負(fù)相關(guān)，說明迂曲度與比表面積共同作用，是非均質(zhì)性影響多孔介質(zhì)滲透率的內(nèi)在原因。

圖19 非均質(zhì)性模型迂曲度、比表面積乘積與非均質(zhì)性因子的關(guān)系

對模擬結(jié)果進(jìn)行數(shù)值擬合，得到三維各向異性、非均質(zhì)性多孔介質(zhì)模型（未細(xì)分Rφ）的滲透率與迂曲度滿足冪函數(shù)關(guān)系（見表5）。

系數(shù)a和b與多孔介質(zhì)孔隙結(jié)構(gòu)有關(guān)，且b值與多孔介質(zhì)的孔徑呈負(fù)相關(guān)，與比表面積呈正相關(guān)。擬合結(jié)果顯示相關(guān)系數(shù)平方值均接近70%，最高超過80%，精度較高，可用于巖心滲透率的近似估算。

表5 三維多孔介質(zhì)的滲透率與迂曲度關(guān)系式

綜上所述，為更好地描述各向異性、非均質(zhì)性對流體流動的影響，表征非常規(guī)儲集層的滲透率模型必需包含孔隙尺度各向異性、孔隙分布非均質(zhì)性這2項(xiàng)反映其微觀孔隙結(jié)構(gòu)性質(zhì)的參數(shù)[15]。

4 結(jié)論

致密砂巖孔隙形態(tài)復(fù)雜，孔隙尺度各向異性、孔隙分布非均質(zhì)性顯著，其中主流方向各向異性因子最大，其他方向相近。

各向異性影響多孔介質(zhì)中孔隙長軸的取向性及流體流動路徑，沿各向異性因子大的方向迂曲度小、流體流動消耗能量小，迂曲度與各向異性的強(qiáng)相關(guān)是各向異性影響滲透率的根本原因。

非均質(zhì)性對滲透率的影響表現(xiàn)為迂曲度與比表面積的共同作用，比表面積、迂曲度的乘積與非均質(zhì)性呈明顯負(fù)相關(guān)，孔隙分布非均質(zhì)性越強(qiáng)，乘積值越小，滲透率越大。

復(fù)雜多孔介質(zhì)的滲透率與迂曲度滿足乘冪關(guān)系式，擬合精度較高，可用于巖心滲透率的近似估算。

符號注釋：

a——擬合系數(shù)，10-3μm2；A——各向異性因子，無因次；b——擬合系數(shù)，無因次；cs——聲速，m/s；Cs——固相生長核；d——固相顆粒大小，m；D——固相生長速度，f；Dx、Dy、Dz——x，y，z方向的固相生長速度，f；ei——格子速度，m/s；fi（x，t），fj——離散分布函數(shù)，kg/m3；feq,j——平衡態(tài)分布函數(shù)，kg/m3；g——孔隙或子單元塊的序號；H——孔隙分布非均質(zhì)性因子，無因次；hg——第g個孔隙的高度，m；i，j——離散速度方向；k——QSGS算法中固相生長方向；K0——滲透率，10-3μm2；M——轉(zhuǎn)換矩陣；M-1——M的逆矩陣；n——多孔介質(zhì)劃分的子單元塊數(shù)量，塊；nx，ny——一條沿x或y方向穿過多孔介質(zhì)的直線，其單位長度所遇到的平均孔隙數(shù)，個；N——多孔介質(zhì)中孔隙總數(shù)量，個；o，p，q——多孔介質(zhì)中x，y，z方向的坐標(biāo)；Pc——粗結(jié)構(gòu)固相生長核概率，f；Pcd——固相生長核概率，f；Pk——不同方向的生長速度，f；Pr——細(xì)結(jié)構(gòu)固相生長核概率，f；r——孔隙半徑包含格子數(shù)，個；Re——雷諾數(shù)，無因次；se，sm，sq，sv，sε，sπ——能量矩、三階矩、能量通量矩、應(yīng)力張量矩、能量的平方矩、四階矩的弛豫速率，無因次；S——比表面積，cm2/cm3；S——對角矩陣；t——時間，s；u——流體速度，m/s；u（o，p，q）——多孔介質(zhì)（o，p，q）處流體速度，m/s；u（o，p，q，t）——t時刻多孔介質(zhì)（o，p，q）處流體速度，m/s；ux（o，p，q）——多孔介質(zhì)（o，p，q）處x方向的流體速度，m/s；——x方向平均流體速度，m/s；wg——第g個孔隙的寬度，m；wi——權(quán)系數(shù)，無因次；x，y，z——空間直角坐標(biāo)系的3個方向；δt——時間步長，s；δE——模擬結(jié)果相對偏差，無因次；δL——單位網(wǎng)格長度，m；?p——壓力梯度，Pa/m；μ——動力黏度，Pa·s；ν——運(yùn)動黏度，m2/s；ρ——流體密度，kg/m3；τ——迂曲度，無因次；τR——無因次弛豫時間；φ——多孔介質(zhì)總孔隙度，%；φg——第g個子單元塊的孔隙度，%；φeff——連通孔隙度，%；φR——細(xì)結(jié)構(gòu)孔隙度，%。