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

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

0 引言

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

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

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

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

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

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

表1 柱體巖樣實際物理尺寸

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1.1 QSGS方法

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

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

2.1.2 孤立孔隙刪除

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

2.2 LBM模型及驗證

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

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

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

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

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

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

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

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

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

2.3 模擬結(jié)果

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

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

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

刪除其中孤立孔隙，將線條間隔設為4個格子，運用（1）式計算多孔介質(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取值對應的比表面積均值分別為4.1×105cm2/cm3，4.7×105cm2/cm3，5.3×105cm2/cm3）；連通孔隙度受Pcd的影響相對較小，與Pcd呈負相關(guān)。

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

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

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

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

圖7為Pcd=0.002 5，A=1.55時各向異性多孔介質(zhì)的模擬結(jié)果。以各向異性多孔介質(zhì)的模擬結(jié)果為基礎，采用（8）式、（9）式分別計算三維各向異性多孔介質(zhì)模型的滲透率和迂曲度。分析滲透率、各向異性、迂曲度間的相互關(guān)系（見圖8—圖10）可知，各向異性多孔介質(zhì)模型的滲透率與各向異性因子呈正相關(guān)；滲透率與迂曲度負相關(guān)；迂曲度與孔隙尺度各向異性因子呈負相關(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ù)字巖心模擬計算多孔介質(zhì)滲透率所得結(jié)果與單位網(wǎng)格長度密切相關(guān)，單位網(wǎng)格長度增加1個數(shù)量級時，多孔介質(zhì)的滲透率至少增加2個數(shù)量級；②非常規(guī)儲集層納米級孔隙發(fā)育，為更好研究孔隙的微觀結(jié)構(gòu)及其與滲透率的相關(guān)性，模擬中選取了較小的單位網(wǎng)格長度（10 nm），滲透率模擬計算結(jié)果與實際巖心的測量值差異較大，但滲透率與孔隙微觀結(jié)構(gòu)的相關(guān)性與單位網(wǎng)格長度的取值無關(guān)。

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

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

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

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

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

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

3 分析與討論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

復雜多孔介質(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方向的坐標；Pc——粗結(jié)構(gòu)固相生長核概率，f；Pcd——固相生長核概率，f；Pk——不同方向的生長速度，f；Pr——細結(jié)構(gòu)固相生長核概率，f；r——孔隙半徑包含格子數(shù)，個；Re——雷諾數(shù)，無因次；se，sm，sq，sv，sε，sπ——能量矩、三階矩、能量通量矩、應力張量矩、能量的平方矩、四階矩的弛豫速率，無因次；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——空間直角坐標系的3個方向；δt——時間步長，s；δE——模擬結(jié)果相對偏差，無因次；δL——單位網(wǎng)格長度，m；?p——壓力梯度，Pa/m；μ——動力黏度，Pa·s；ν——運動黏度，m2/s；ρ——流體密度，kg/m3；τ——迂曲度，無因次；τR——無因次弛豫時間；φ——多孔介質(zhì)總孔隙度，%；φg——第g個子單元塊的孔隙度，%；φeff——連通孔隙度，%；φR——細結(jié)構(gòu)孔隙度，%。