多孔介質(zhì)孔隙結(jié)構(gòu)模型及其在滲流分析中的應用

胡黎明，林丹彤，張鵬偉，2，張興昊，郭豪豪，Jay N Meegoda，董曉強

(1.清華大學 水利水電工程系 水沙科學與水利水電工程國家重點實驗室，北京 100084；2.北京交通大學 土木建筑工程學院，北京 100044；3.新澤西理工學院 土木與環(huán)境工程系，美國新澤西 07102； 4.太原理工大學 土木工程學院，太原 030024)

隨著人類經(jīng)濟和社會的高速發(fā)展，目前人類環(huán)境面臨許多新的科學與技術(shù)問題：土壤地下水污染日趨加重，環(huán)境污染修復工作任重道遠；溫室效應加劇，二氧化碳地下封存成為緩解溫室效應的有效手段；全球能源需求不斷攀升，頁巖氣等非常規(guī)資源是世界能源開采的重要發(fā)展方向。上述問題與地質(zhì)巖土材料中的滲流密切相關(guān)。地質(zhì)巖土材料是自然界最常見的多孔介質(zhì)，多孔介質(zhì)多相滲流特性的研究具有廣泛的理論價值及應用前景。

地下水是水資源的重要組成部分，根據(jù)2020年《中國水資源公報》，全國供水量的15.4%來自地下水，其中淺層地下水占95.7%[1]，高效地修復被污染的地下水是環(huán)境巖土工程領(lǐng)域的重要課題。污染物在地下水中通常以溶質(zhì)形式存在，其在地下水中的運動過程涉及對流、機械彌散、分子擴散和吸附等多種作用[2]；在采用曝氣法等常用修復方法進行地下水原位修復的過程中，通常涉及土壤介質(zhì)中的氣液兩相運動，以及物質(zhì)擴散、溶解、揮發(fā)、吸附等傳質(zhì)問題[3]；近年來納米材料作為一種新興的環(huán)境修復材料也開始被應用到地下水的原位修復當中，以納米鐵為代表的納米修復材料在地下水中的有效粒徑通常處于膠體的尺寸范圍[4-6]。上述溶質(zhì)、氣液兩相流、膠體在地下水中的運移主要存在于毫米-微米尺度的孔隙中。

隨著經(jīng)濟社會的不斷發(fā)展，化石能源不斷燃燒導致大量的二氧化碳氣體被排入大氣，人類逐漸面臨著海平面上升、全球氣候變暖等全球性氣候問題，我國也提出了“碳中和”與“碳達峰”的目標[7]。目前，二氧化碳地質(zhì)封存是減少二氧化碳排放的重要手段之一，該方法是指將二氧化碳轉(zhuǎn)變?yōu)槌R界體并將其注入深部地層中進行長期儲存[8]。超臨界狀態(tài)二氧化碳具有液相較高的密度，同時又保持氣相較低的粘滯系數(shù)，在巖體孔隙中復雜的兩相滲流和溫度、力學和化學反應等形成了多場耦合作用。二氧化碳通常封存于地下巖層中，其孔隙尺度在微米量級。

作為一種非常規(guī)的化石能源，頁巖氣的開采已有200年的歷史[9]，是世界能源的新領(lǐng)域[10]。頁巖儲層的孔隙尺度處于納米尺度，頁巖氣、孔隙水以及劈裂液的多相流動以及氣壓變化導致儲層壓縮、滲透率降低產(chǎn)生了多場耦合復雜過程；頁巖儲層的高壓、高溫特性進一步使頁巖氣流動呈現(xiàn)出非達西流動的特征，準確預估頁巖氣產(chǎn)量、優(yōu)化開采設(shè)計成為該領(lǐng)域的研究難點[11]；頁巖氣儲層滲透率低，具有較好的密封特性，近年來有學者提出采用二氧化碳驅(qū)替頁巖氣，同時實現(xiàn)頁巖氣開采與二氧化碳地質(zhì)封存[12]。

上述問題涉及的孔隙流體運動、物質(zhì)運移等物理現(xiàn)象，與多孔介質(zhì)孔隙結(jié)構(gòu)以及多相滲流特性密切相關(guān)。傳統(tǒng)的多孔介質(zhì)滲流研究多集中于宏觀物理現(xiàn)象的唯象描述與概化模擬，對于其微觀機理與物理本質(zhì)的研究尚不完善?；谶B續(xù)介質(zhì)假設(shè)的宏觀模型忽略巖土介質(zhì)的微觀孔隙結(jié)構(gòu)特性，通過設(shè)定宏觀參數(shù)來概化巖土介質(zhì)復雜的微觀特性，無法反映微觀尺度下的多相滲流、阻滯、彌散、指進效應等現(xiàn)象。巖土介質(zhì)的滲流特性與其孔隙微觀結(jié)構(gòu)密切相關(guān)。近年來，多孔介質(zhì)孔隙結(jié)構(gòu)模型(Pore Structure Model)得到了發(fā)展，區(qū)別于傳統(tǒng)宏觀模型中對多孔介質(zhì)為理想連續(xù)體的假設(shè)，孔隙結(jié)構(gòu)模型將多孔介質(zhì)中的固相骨架以及固相骨架之間的孔隙空間加以區(qū)分，并基于多孔介質(zhì)中復雜的孔隙結(jié)構(gòu)對滲流問題進行分析和計算。

按照對孔隙結(jié)構(gòu)的描述及滲流計算方法由難到易，可以將目前常用的孔隙結(jié)構(gòu)模型分為孔隙重構(gòu)模型(Pore Structure Reconstruction Model)、孔隙網(wǎng)絡(luò)模型(Pore Network Model)和等效孔隙網(wǎng)絡(luò)模型(Equivalent Pore Network Model)三種?？紫吨貥?gòu)模型是對真實多孔介質(zhì)內(nèi)部的孔隙結(jié)構(gòu)進行重構(gòu)，并直接求解滲流控制方程；孔隙網(wǎng)絡(luò)模型是在孔隙重構(gòu)模型的基礎(chǔ)上對孔隙結(jié)構(gòu)中的孔隙和孔喉單元進行區(qū)分，從而對其中滲流計算進行簡化；等效孔隙網(wǎng)絡(luò)模型則是在孔隙網(wǎng)絡(luò)模型的基礎(chǔ)上進一步將孔隙排列在規(guī)則網(wǎng)格格點上，基于多孔介質(zhì)孔隙特征統(tǒng)計參數(shù)建立簡化的等效模型。本文將對孔隙重構(gòu)模型、孔隙網(wǎng)絡(luò)模型和等效孔隙網(wǎng)絡(luò)模型的特點、構(gòu)建方法以及在滲流分析方面的應用情況進行總結(jié)，梳理該領(lǐng)域亟待解決的問題和孔隙結(jié)構(gòu)模型未來發(fā)展前景，為相關(guān)研究的開展提供參考。

1 孔隙重構(gòu)模型

孔隙重構(gòu)模型是指借助數(shù)值模擬或計算機圖形學等手段，對實際巖土介質(zhì)所占據(jù)空間中的固相骨架和孔隙結(jié)構(gòu)的位置進行數(shù)字化，進而以孔隙結(jié)構(gòu)作為邊界，采用數(shù)值模擬的方法對滲流過程控制方程進行直接求解的模型。孔隙重構(gòu)模型是對實際巖土介質(zhì)孔隙結(jié)構(gòu)的直接重現(xiàn)，對孔隙結(jié)構(gòu)的刻畫較為精細，同時由于其求解過程是在連通的復雜孔隙區(qū)域內(nèi)進行直接的計算，計算代價通常較大，主要應用于較小的計算尺度。

1.1 孔隙重構(gòu)模型的建立

目前建立孔隙重構(gòu)模型的主要手段包括離散元重構(gòu)法和序列圖像重構(gòu)法兩種。離散元重構(gòu)法是基于顆粒材料的粒徑級配、孔隙率等基本參數(shù)對多孔介質(zhì)的空間結(jié)構(gòu)進行模擬，主要適用于粗粒土；系列圖像重構(gòu)法則是針對特定樣本進行掃描，基于一系列二維的數(shù)字圖像對其三維孔隙結(jié)構(gòu)進行重構(gòu)的方法，不僅適用于粗粒土，也適用于巖石。

1.1.1離散元重構(gòu)法

離散元重構(gòu)法是指通過離散元軟件對一系列顆粒進行堆積，以模擬一定粒徑級配、孔隙率、阻尼和摩擦角等參數(shù)條件下粗顆粒巖土介質(zhì)的孔隙結(jié)構(gòu)。在該方法中通常采用規(guī)則的球形顆粒[13]，其示意圖如圖1所示，也有學者基于真實粗粒土顆粒的形態(tài)來進行堆積[14]。在堆積完成后，則可以獲得結(jié)構(gòu)中任意顆粒的位置及所占據(jù)的空間，即獲得多孔介質(zhì)中固相骨架及孔隙的空間結(jié)構(gòu)。離散元重構(gòu)法對于砂土等顆粒本身接近球體并且具有較好連通性的粗粒土多孔介質(zhì)材料具有較好的適用性，對于黏土、巖石等低連通性多孔介質(zhì)材料則適用性較差[15]。

圖1 離散元重構(gòu)法構(gòu)建的多孔介質(zhì)空間結(jié)構(gòu)Fig.1 Discrete element method to simulate the spatial structure of porous media

1.1.2序列圖像重構(gòu)法

序列圖像重構(gòu)法是基于成像技術(shù)獲取巖土材料樣本的一系列二維掃描圖像，再通過計算機圖形學方法將二維圖像重構(gòu)為三維數(shù)字圖像[16]，該技術(shù)不僅適用于粗粒土同時也適用于頁巖等巖石樣本，巖石樣品重構(gòu)后的三維數(shù)字圖像也稱為數(shù)字巖芯[17-18]。目前常見的成像方法主要包括計算機斷層掃描(CT)成像技術(shù)和聚焦離子束掃描電鏡(FIB-SEM)成像技術(shù)[19]。CT成像技術(shù)是通過發(fā)射端發(fā)射X射線，穿過巖土介質(zhì)時不同類型介質(zhì)會對射線產(chǎn)生不同程度的衰減，由接收端接收后還原得到各位置的衰減系數(shù)，生成該掃描面的CT掃描圖像[20]，典型的粗粒土CT掃描圖像如圖2(a)所示；FIB-SEM成像技術(shù)則是通過鎵離子束對巖土介質(zhì)樣本進行連續(xù)切割，同時在電子束下成像的技術(shù)，該方法相比于CT技術(shù)具有更高的分辨率，更適用于頁巖等孔隙尺寸在納米級的巖土介質(zhì)的研究[21]，采用FIB-SEM法獲得的頁巖二維數(shù)字圖像如圖2(b)所示[22]。

通過上述方法獲得的二維圖像通常含有噪聲信息不能直接應用，因此需要經(jīng)過濾波、閉運算、二值化等過程進行預處理提高圖像精度。經(jīng)過預處理的二維數(shù)字圖像則可以經(jīng)過規(guī)格化、像素值插值以及三維圖像分割等步驟獲得重構(gòu)的三維多孔介質(zhì)空間結(jié)構(gòu)[23]。圖2(c)展示了一組粗粒土的CT圖像重構(gòu)結(jié)果，圖2(d)展示了一組頁巖的FIB-SEM圖像重構(gòu)結(jié)果[22]。

1.2 基于孔隙重構(gòu)模型的滲流分析

基于上述方法建立多孔介質(zhì)孔隙重構(gòu)模型后，可以對其中孔隙中的滲流問題進行分析和計算，該問題實為復雜邊界條件下滲流方程的求解問題。目前常用的方法包括計算流體力學方法(CFD)和格子玻爾茲曼方法(LBM).

CFD是流體力學數(shù)值模擬的常用方法，該方法基于流體的連續(xù)介質(zhì)假設(shè)，通常以歐拉方程和納維-斯托克斯(N-S)方程作為控制方程，基于包括有限差分、有限體積、有限單元、有限分析和邊界元等離散方法對控制方程進行求解[16]。ZHANG et al[25]基于CT圖像對透水混凝土內(nèi)部的孔隙結(jié)構(gòu)進行了重構(gòu)，并對滲流過程進行了模擬和分析；SANEMATSU et al[26]基于CT圖像獲得了砂巖樣本內(nèi)部的孔隙結(jié)構(gòu)，并采用有限元法求解了其中的流場，基于流場計算結(jié)果采用拉格朗日粒子追蹤法對其中納米顆粒的運移和滯留行為進行分析。由于CFD模型的計算精度很大程度上依賴于網(wǎng)格劃分精細程度，因此對于復雜孔隙結(jié)構(gòu)的計算代價通常較大。

圖2 基于序列圖像重構(gòu)法獲取的多孔介質(zhì)孔隙結(jié)構(gòu)Fig.2 Spatial structure of porous media obtained by sequential image reconstruction method

LBM方法是一種介于流體微觀分子動力學與宏觀流體力學分析之間的計算方法，該方法將流體離散成介尺度的流體粒子，并通過粒子的平衡態(tài)分布函數(shù)特性計算其宏觀滲流力學特征，對于復雜邊界的計算求解問題有較大的優(yōu)勢[27]，目前被廣泛應用于多孔介質(zhì)滲流問題細觀尺度研究。針對球堆積孔隙重構(gòu)模型，劉一飛等[13]采用LBM研究了多孔介質(zhì)粒徑級配對土體滲透系數(shù)的影響；LONG和HILPERT[28]采用LBM對其中膠體的收集過程進行了模擬。劉振宇和吳慧英[29]基于CT圖像對砂巖樣本中的孔隙結(jié)構(gòu)進行了重構(gòu)，并基于LBM對其中氣液兩相滲流過程進行了研究。研究表明LBM方法能夠?qū)碗s邊界條件下的多種多孔介質(zhì)滲流問題進行有效求解。

2 孔隙網(wǎng)絡(luò)模型

基于多孔介質(zhì)孔隙重構(gòu)模型進行滲流分析時，需要針對復雜的三維邊界條件求解滲流控制方程，計算代價通常較大。為了提高計算效率，研究者們提出了孔隙網(wǎng)絡(luò)模型的概念[30]?？紫毒W(wǎng)絡(luò)模型的核心思想是將孔隙結(jié)構(gòu)簡化為孔隙和孔喉的組合[31]，通常將孔隙結(jié)構(gòu)中膨大的空腔定義為孔隙，將其中狹長的通道定義為孔喉，通?？梢圆捎们蝮w代表孔隙、用圓柱體代表孔喉[32]，其示意圖如圖3所示，也有研究采用立方體和棱柱體分別代表孔隙和孔喉，以描述在邊角處的毛細作用[17]。經(jīng)過上述簡化，孔隙結(jié)構(gòu)分解為大量規(guī)則的孔隙和孔喉，從而形成由孔隙和孔喉組成的網(wǎng)絡(luò)，可以在此基礎(chǔ)上進行滲流過程分析。

圖3 孔隙網(wǎng)絡(luò)模型示意圖Fig.3 Schematic diagram of pore network model

2.1 孔隙網(wǎng)絡(luò)模型的建立

在多孔介質(zhì)孔隙重構(gòu)模型的基礎(chǔ)上，通過一定的幾何規(guī)則對孔隙和孔喉進行區(qū)分，該過程稱為孔隙結(jié)構(gòu)提取。目前常用的孔隙結(jié)構(gòu)提取方法包括四面體剖分法、居中軸線法和最大球法等。

2.1.1四面體剖分法

四面體剖分法適用于粗粒土，該方法將孔隙結(jié)構(gòu)模型中所有固相顆?？醋髌淝蛐幕蝾w粒質(zhì)心的代表點，采用Delaunay四面體法對整個模型中的全部代表點進行四面體剖分，該剖分方法可以保證任意四面體的外接球不包含其余的點，以此保證所得的四面體均為最接近頂點所組成的四面體[33]，圖1所示的顆粒堆積體所對應的四面體剖分結(jié)果如圖4(a)所示。

在完成四面體剖分后，針對每一個四面體，可以確定位于該四面體頂點處的四個顆粒的最大內(nèi)切球作為孔隙，進而獲得孔隙的半徑和位置信息，如圖4(b)所示。對于規(guī)則球體堆積得到的孔隙結(jié)構(gòu)模型，可以直接利用固相顆粒的半徑和位置關(guān)系確定其最大內(nèi)切球，而對于不規(guī)則固相顆粒堆積體或由序列圖像重構(gòu)得到的孔隙結(jié)構(gòu)模型，由于其中固相顆粒形狀不規(guī)則，則對其最大內(nèi)切球的確定需要利用極值搜索方法進行，該過程應遵循如下原則：①內(nèi)切球的球心應位于剖分得到的四面體內(nèi)；②當內(nèi)切球與任一位置顆粒接觸時，即得到該處的內(nèi)切球大??；③內(nèi)切球的位置與半徑應取為實數(shù)[23]。

圖4 基于四面體剖分法的孔隙結(jié)構(gòu)提取Fig.4 Pore structure extraction based on tetrahedral subdivision method

在孔隙網(wǎng)絡(luò)模型中，相鄰孔隙間利用3個固相顆粒圍成的通道進行滲流，因此由這3個固相顆粒形成的通道即為孔隙間的孔喉。BRYANT et al[33]提出了孔喉面積的等效計算方法，如圖4(c)所示，其中rc表示3個球形成的橫截面內(nèi)的內(nèi)切圓半徑，re表示與截面中空隙部分滲流面積相等的圓半徑，將二者的均值定義為等效的孔喉半徑，基于上述方法可以確定孔喉的半徑。

2.1.2居中軸線法

居中軸線將多孔介質(zhì)孔隙結(jié)構(gòu)看成形狀不規(guī)則的連通管道，并將管道中軸線相連所構(gòu)成的網(wǎng)絡(luò)做為該孔隙結(jié)構(gòu)的居中軸線，通?？梢圆捎每s減算法[34]或烈火模擬算法[35]得到，圖5展示了一組基于真實巖心樣本所提取得到的居中軸線圖[36]。基于居中軸線進行合理的分割和簡化可以得到孔隙網(wǎng)絡(luò)模型，通常將中軸線之間的節(jié)點定義為孔隙，將中軸線上孔隙空間的局部最小區(qū)域定義為孔喉[37]。AL-RAOUSH和WILLSON[38]細化了基于居中軸線法的孔隙識別方法，將居中軸線連線的交點作為孔隙中心，并采用膨脹法在該中心放置逐步增大的圓球，直到圓球接觸巖石骨架，基于該方法確定了孔隙的半徑。趙秀才[39]對居中軸線法進行了進一步的優(yōu)化，對提取得到的枝節(jié)路徑和冗余邊界節(jié)點進行了處理。

圖5 LINDQUIST et al[36]采用居中軸線法提出的孔隙網(wǎng)絡(luò)模型Fig.5 Pore network extracted via medial axis method by LINDQUIST et al[36]

2.1.3最大球法

最大球法最早由SINGH和MOHANTY[40]提出，后續(xù)學者對該方法進行了進一步的發(fā)展和完善[41-43]。對于孔隙結(jié)構(gòu)模型中孔隙空間中的任意一點，可以尋找以該點為球心并且與固相骨架相切的最大內(nèi)切球，在去除包含在企圖內(nèi)切球中的冗余球體后，則可以獲得能夠描述孔隙空間的全部無冗余內(nèi)切球的集合，并將其中局部半徑最大的內(nèi)切球半徑定義為孔隙半徑，將兩個孔隙中局部半徑最小的內(nèi)切球半徑定義為孔喉半徑[44]，其示意圖如圖6所示。

2.2 基于孔隙網(wǎng)絡(luò)模型的滲流分析

在采用上述方法對孔隙結(jié)構(gòu)模型中的孔隙和孔喉進行區(qū)分后，則可以將邊界形狀復雜的孔隙結(jié)構(gòu)模型轉(zhuǎn)化幾何形狀規(guī)則的孔隙網(wǎng)絡(luò)模型。根據(jù)圖1所示球體堆積模型的孔隙結(jié)構(gòu)，可以得到對應的孔隙網(wǎng)絡(luò)模型如圖3所示。

圖6 采用最大球法提取得到的孔隙和孔喉Fig.6 Pore body and pore throat extracted by maximum sphere method

利用孔隙網(wǎng)絡(luò)模型，可以簡化多孔介質(zhì)孔隙結(jié)構(gòu)滲流問題的計算分析。對于單相滲流問題，由于多孔介質(zhì)中的滲流過程主要由模型中狹窄的通道控制，因此通常假設(shè)孔隙網(wǎng)絡(luò)模型中的單個孔隙中壓強處處相同，流入每個孔隙的流量應符合質(zhì)量平衡；孔喉兩端連接的孔隙壓強可以不同，其中的滲流過程可以用Hagen-Poseuille管流公式計算[45]，從而使研究者免于對復雜邊界條件下的滲流方式進行求解，極大提高了計算效率。

基于上述研究思路，研究者們基于孔隙網(wǎng)絡(luò)模型開展了多孔介質(zhì)中單相滲流、兩相流及溶質(zhì)運移等問題的研究。GAO et al[45]構(gòu)建了等直徑球體簡單立方堆積條件下所對應的孔隙網(wǎng)絡(luò)模型，并基于該方法研究了不同粒徑條件下多孔介質(zhì)孔隙結(jié)構(gòu)的滲透系數(shù)，計算結(jié)果與文獻中的試驗結(jié)果較為接近；GONG和PIRI[46]構(gòu)建了基于CT掃描圖片重構(gòu)的巖心結(jié)構(gòu)所對應的孔隙網(wǎng)絡(luò)模型，并在飽和及非飽和條件下對該巖心中的溶質(zhì)運移問題進行了研究；楊永飛等[47]基于數(shù)字巖心構(gòu)建孔隙網(wǎng)絡(luò)模型，研究了頁巖試樣中的孔隙結(jié)構(gòu)特征及油相流動能力。

3 等效孔隙網(wǎng)絡(luò)模型

孔隙網(wǎng)絡(luò)模型孔隙及孔喉基于實際孔隙結(jié)構(gòu)建立，其拓撲結(jié)構(gòu)較為復雜。針對多孔介質(zhì)孔隙結(jié)構(gòu)中的孔隙尺寸、孔喉尺寸、孔隙與孔喉的連接情況、孔隙率等數(shù)據(jù)進行統(tǒng)計分析，可以獲得不同類型多孔介質(zhì)材料孔隙結(jié)構(gòu)的特征參數(shù)取值范圍。一般認為，當多孔介質(zhì)孔隙結(jié)構(gòu)的特征參數(shù)相同時，其滲流特征具有相似性。基于上述思路提出等效孔隙網(wǎng)絡(luò)模型，該模型與孔隙網(wǎng)絡(luò)模型的區(qū)別在于該模型的建立不依賴于具體的多孔介質(zhì)樣本，而是基于描述多孔介質(zhì)孔隙結(jié)構(gòu)特征的特征參數(shù)，將孔隙排列在一定規(guī)則的網(wǎng)格頂點上，并生成孔隙之間連接的孔喉，從而實現(xiàn)對孔隙結(jié)構(gòu)模型拓撲結(jié)構(gòu)的簡化。

3.1 等效孔隙網(wǎng)絡(luò)模型的構(gòu)建

構(gòu)建等效孔隙網(wǎng)絡(luò)模型，首先需要對描述多孔介質(zhì)孔隙結(jié)構(gòu)特征的參數(shù)進行定義和分類統(tǒng)計，之后將根據(jù)特征參數(shù)的取值，生成不同類型的等效孔隙網(wǎng)絡(luò)模型。本文基于球形孔隙及圓柱形孔喉的等效孔隙網(wǎng)絡(luò)模型的相關(guān)特征進行介紹，也有學者采用立方體及棱柱體作為孔隙和孔喉進行模擬[48]。

3.1.1描述多孔介質(zhì)孔隙結(jié)構(gòu)的特征參數(shù)

基于球形孔隙和圓柱形孔喉的基本假設(shè)，可以用如下參數(shù)來描述多孔介質(zhì)孔隙。

1) 孔隙半徑：多孔介質(zhì)中膨大的空腔為孔隙，當用球體對孔隙進行描述時，對應的球體半徑則記為孔隙半徑。

2) 孔喉半徑：多孔介質(zhì)中狹長的通道為孔喉，當用圓柱體對孔喉進行描述時，對應的圓柱體半徑記為孔喉半徑。

3) 孔喉長度：描述孔喉的圓柱體的長度記為孔喉長度。

4) 孔隙率：多孔介質(zhì)中孔隙和孔喉的總體積為孔隙體積，孔隙體積占多孔介質(zhì)總體積(孔隙體積+固相體積)的比例為孔隙率。

5) 配位數(shù)：配位數(shù)是指一個孔隙連接的孔喉數(shù)量。

在多孔介質(zhì)孔隙結(jié)構(gòu)中，每一個孔隙、孔喉的尺寸以及配位數(shù)均可能不同，因此除孔隙率這一宏觀參數(shù)外，孔隙半徑、孔喉半徑、孔喉長度及配位數(shù)在多孔介質(zhì)中均有一定的分布特征。對于頁巖等具有較強各向異性的多孔介質(zhì)，其配位數(shù)的分布在不同方向上也存在顯著差異。

張鵬偉[15]通過文獻調(diào)研對不同類型巖土介質(zhì)孔隙特征參數(shù)進行了總結(jié)，指出對于砂土、砂巖及頁巖3種巖土介質(zhì)，其孔隙尺寸由毫米級別逐漸過渡到納米級別，其孔隙連通性也逐漸由高向低轉(zhuǎn)變，文獻中所述的上述3種巖土介質(zhì)的平均配位數(shù)范圍如表1所示。

表1 不同巖土介質(zhì)的平均配位數(shù)取值范圍Table 1 Value range of average coordination number of different porous media

3.1.2等效孔隙網(wǎng)絡(luò)模型的構(gòu)建

基于上述孔隙特征參數(shù)的統(tǒng)計結(jié)果，研究者們建立了不同類型的等效孔隙網(wǎng)絡(luò)模型。FATT[49]最早基于毛細管模型提出了采用二維蜂窩狀網(wǎng)格、規(guī)則方形網(wǎng)格、雙重六邊形網(wǎng)格以及三重六邊形網(wǎng)格的二維等效孔隙網(wǎng)絡(luò)模型，并利用上述模型對基質(zhì)吸力曲線進行了模擬。后來的研究者將上述思路拓展到了三維情況[50-51]，例如REEVES和CELIA[52]將孔隙排列在立方體網(wǎng)格的格點上建立了三維等效孔隙網(wǎng)絡(luò)模型，上述文獻中所報道的孔隙網(wǎng)絡(luò)模型的配位數(shù)是恒定的，取決于所采用的網(wǎng)格類型。也有學者將等效孔隙網(wǎng)絡(luò)模型稱為非重建模型[16]或規(guī)則拓撲孔隙網(wǎng)絡(luò)模型[18]。

RAOOF和HASSANIZADEH[53]提出了最大配位數(shù)為26的等效孔隙網(wǎng)絡(luò)模型，將孔隙排列在立方體網(wǎng)格的格點上，每一個孔隙最多與周圍的26個孔隙相連，并且基于配位數(shù)的正態(tài)分布隨機生成每個孔隙與周圍孔隙的連接關(guān)系，更充分地考慮了孔隙連通的多向性和無序性。張鵬偉[15]基于最大配位數(shù)為26的等效孔隙網(wǎng)絡(luò)模型，通過考慮配位數(shù)在不同方向的連接比例，建立了能夠反映各向異性的等效孔隙網(wǎng)絡(luò)模型，并對不同巖土介質(zhì)的滲透率進行了計算，驗證了該模型的有效性，并將模型應用到了頁巖氣開采和二氧化碳注入頁巖過程中的水氣兩相流模擬中。上述等效孔隙網(wǎng)絡(luò)的結(jié)構(gòu)示意圖如圖7所示。

圖7 等效孔隙網(wǎng)絡(luò)模型示意圖Fig.7 Schematic diagram of equivalent pore network model

3.2 基于孔隙網(wǎng)絡(luò)模型的多相滲流和物質(zhì)輸移問題分析

基于等效孔隙網(wǎng)絡(luò)模型，研究者們對多孔介質(zhì)中的氣液滲流、溶質(zhì)運移、膠體運移等問題進行了研究。

3.2.1滲流分析

基于頁巖等巖體的孔隙結(jié)構(gòu)統(tǒng)計參數(shù)的取值范圍，ZHANG et al[54]構(gòu)建了等效孔隙網(wǎng)絡(luò)模型描述頁巖孔隙連通情況及低滲透性，同時綜合考慮了氣體在納米尺度孔隙中的Knudsen擴散及滑移流作用，研究了不同壓力條件下頁巖氣滲流過程，結(jié)果表明頁巖氣的表觀滲透率與孔隙壓力顯著相關(guān)，隨著頁巖儲層壓力的降低滑移流和克努森擴散逐漸占主導地位，頁巖氣流動存在明顯的指進效應，如圖8所示；進一步考慮頁巖的各向異性進行計算分析，通過與試驗結(jié)果的比較，進一步驗證了等效孔隙網(wǎng)絡(luò)模型在滲流模擬上的有效性[55]；此外，ZHANG et al[56]還基于等效孔隙網(wǎng)絡(luò)模型研究了二氧化碳地質(zhì)封存和頁巖氣開發(fā)耦合過程中二氧化碳與甲烷的競爭吸附過程。

3.2.2溶質(zhì)運移分析

RAOOF et al[57]發(fā)展等效孔隙網(wǎng)絡(luò)模型描述多孔介質(zhì)中溶質(zhì)運移過程，考慮溶質(zhì)隨著孔隙水的對流作用，以及由于多孔介質(zhì)復雜孔隙流場引起的機械彌散作用，并研究了飽和度對溶質(zhì)運移過程的影響[58]；張興昊等[59]在此基礎(chǔ)上考慮了溶質(zhì)分子擴散作用，計算結(jié)果表明孔隙結(jié)構(gòu)參數(shù)對水動力彌散過程有顯著影響，在高流速條件下水動力彌散系數(shù)與機械彌散系數(shù)接近，在低流速條件下水動力彌散系數(shù)與分子擴散系數(shù)接近，而在中等流速的過渡階段，水動力彌散系數(shù)小于機械彌散系數(shù)與分子擴散系數(shù)的和，如圖9所示，其原因是分子擴散作用導致溶質(zhì)在低流速區(qū)域的運移。

3.2.3膠體運移分析

近年來，等效孔隙網(wǎng)絡(luò)模型也逐漸被應用到膠體運移和滯留的研究領(lǐng)域。YANG和BALHOFF[60]建立了帶有收縮管孔喉的等效孔隙網(wǎng)絡(luò)模型，并利用粒子追蹤的方法研究了其中膠體的運動過程；SEETHA et al[61]建立了帶有圓柱形孔喉的等效孔隙網(wǎng)絡(luò)模型，并提出了孔喉尺度沉積速率的表達式；LIN et al[62]發(fā)展了能描述膠體運移和滯留過程的等效孔隙網(wǎng)絡(luò)模型，模型中考慮了多孔介質(zhì)固相表面粗糙度和化學異質(zhì)性對表面沉積過程的影響作用和對膠體出流和滯留曲線的影響，結(jié)果表明膠體的收集效率系數(shù)、粘附效率系數(shù)、最大沉積面積比等描述膠體運動的孔隙參數(shù)在多孔介質(zhì)存在一定的空間分布特征[62]；通過訓練神經(jīng)網(wǎng)絡(luò)模型實現(xiàn)對膠體的收集效率系數(shù)的快速計算，并基于等效孔隙網(wǎng)絡(luò)模型建立了從孔隙參數(shù)到宏觀運移參數(shù)的升尺度計算方法[63]，此外，在等效孔隙網(wǎng)絡(luò)模型計算分析中考慮了表面沉積、篩濾作用及粒橋3種滯留作用，研究在不同膠體粒徑、注入濃度及流速條件下膠體的宏觀出流及滯留曲線的變化趨勢，結(jié)果表明僅考慮表面沉積作用產(chǎn)生均勻的滯留曲線，而同時考慮粒橋作用則可以重現(xiàn)試驗中膠體在入口處大量滯留的超指數(shù)型滯留曲線，如圖10所示[64]。

圖8 頁巖多流態(tài)計算結(jié)果及指進效應[54]Fig.8 Calculation results of shale multi-flow states and fingering effect[54]

D為水動力彌散系數(shù)，D*為機械彌散系數(shù)，D′為分子擴散系數(shù)圖9 不同孔隙流速條件下對應的彌散系數(shù)Fig.9 Corresponding dispersion coefficient under different pore velocity

4 結(jié)論與展望

本文將多孔介質(zhì)孔隙結(jié)構(gòu)模型分為孔隙重構(gòu)模型、孔隙網(wǎng)絡(luò)模型和等效孔隙網(wǎng)絡(luò)模型，分別介紹了3種模型的概念、建立方法以及在滲流分析領(lǐng)域的應用。從孔隙結(jié)構(gòu)模型到孔隙網(wǎng)絡(luò)模型再到等效孔隙網(wǎng)絡(luò)模型，模型對多孔介質(zhì)孔隙結(jié)構(gòu)的描述逐漸由真實性向等效性發(fā)展，同時計算效率則顯著提高?？紫督Y(jié)構(gòu)模型對于研究滲流過程的微觀機理具有重要作用，目前研究集中于多孔介質(zhì)中滲流分析，包括氣體和液體滲流、溶質(zhì)運移和膠體運移等問題的數(shù)值模擬。

圖10 膠體粒徑對表面沉積和粒橋的影響[64]Fig.10 Effect of colloid size on surface deposition and particle bridging[64]

目前研究主要針對具體的巖土材料樣本構(gòu)建對應的孔隙結(jié)構(gòu)模型。由于巖土材料孔隙結(jié)構(gòu)復雜多樣，今后的研究中應重視對不同類別巖土材料的多孔介質(zhì)孔隙結(jié)構(gòu)特征的統(tǒng)計分析，基于統(tǒng)計參數(shù)的等效性建立更具有普適意義的模型，提高計算分析效率，深入認識多孔介質(zhì)滲流和物質(zhì)輸移過程孔隙尺度的物理機理。同時，隨著孔隙結(jié)構(gòu)模型發(fā)展和計算效率提升，建立微觀機理與宏觀特征之間的內(nèi)在聯(lián)系和定量關(guān)系，進而基于多孔介質(zhì)材料特征實現(xiàn)宏觀滲流和物質(zhì)輸移行為的預測分析是未來的重要研究方向。