含水合物沉積物孔隙流動特性數(shù)值模擬

張吉東，劉曉惠，邊培明，吳強，殷振元

(1.清華大學清華大學深圳國際研究生院，廣東深圳，518055；2.山東科技大學能源與礦業(yè)工程學院，山東青島，266590；3.深圳市百勤石油技術(shù)有限公司，廣東深圳，518055；4.黑龍江科技大學安全工程學院，黑龍江哈爾濱，150022)

天然氣水合物廣泛分布于海底的沉積物層及永久凍土層中，儲量豐富且具有高能量儲存能力，燃燒時只釋放CO2和H2O，是未來清潔能源的重要組成部分[1]。目前，中國[2]、日本[3]等國家均對天然氣水合物進行了試采工作，但由于儲層力學性質(zhì)及其流動的復雜性，天然氣水合物的商業(yè)化開采仍是世界性難題。在含水合物沉積物中，滲透率是流體流動性能的關(guān)鍵指標，決定了水合物的生長、分布及甲烷氣體的開采效率。為此，了解流體在孔隙中的滲透特性對于準確預(yù)測產(chǎn)液行為至關(guān)重要?？紤]到通過室內(nèi)物理實驗確定滲透率特性影響因素需要復雜的測試設(shè)備和精確控制的環(huán)境[4]，本文基于數(shù)值模擬手段研究含水合物沉積物儲層的流動特性。

沉積物孔隙結(jié)構(gòu)是控制滲透率演化規(guī)律的首要因素，而孔隙結(jié)構(gòu)變化則受多種因素的影響，孔隙內(nèi)固相水合物生成和分解是影響孔隙結(jié)構(gòu)的重要因素。目前，學者們普遍將沉積物中的水合物的微觀結(jié)構(gòu)分為表層覆蓋、顆粒膠結(jié)與孔隙填充等形態(tài)[5]；同時，伴隨著水合物飽和度增加，水合物形態(tài)通常從表層覆蓋與孔隙填充型逐漸向顆粒膠結(jié)型轉(zhuǎn)化[6-7]。沉積物粒徑是影響滲透率演化的另一因素，宋永臣等[8-9]開展的沉積物滲透率實驗結(jié)果表明，沉積物的絕對滲透率隨粒徑增大而升高，但由于實驗沉積物砂粒種類較少，并未全面分析黏土至粗砂多粒徑范圍內(nèi)滲透率變化的差異。LIU等[10]發(fā)現(xiàn)南海神狐海域含甲烷水合物沉積物屬于泥質(zhì)細粉砂型儲層。當前對真實儲層粒徑分布的含水合物非均質(zhì)沉積物儲層的滲透率演化規(guī)律研究較少，因此，有必要利用數(shù)值方法深入分析粒徑導致的滲透率變化。張永超等[11]基于核磁共振(NMR)結(jié)果，發(fā)現(xiàn)水合物主要存在于孔徑較大的孔隙中，水合物的生成具有明顯的非均質(zhì)性；蔣觀利等[12]基于CT 實驗，發(fā)現(xiàn)甲烷水合物生成早期結(jié)構(gòu)非常疏松，呈絮狀，后期逐漸密實；李晨安等[13]通過以Al2O3球為多孔介質(zhì)的CT 掃描實驗，發(fā)現(xiàn)水合物生成早期主要位于氣-水界面，孔隙內(nèi)水合物密度較低，隨著水合物飽和度增加，多孔介質(zhì)孔隙結(jié)構(gòu)完全被水合物占據(jù)，整個生成過程中存在明顯的非均質(zhì)性分布。結(jié)合上述分析，當前人工合成的沉積物儲層砂質(zhì)顆粒粒徑通常較單一，無法有效呈現(xiàn)真實含水合物沉積物粒徑范圍廣、差異大的特點；同時，當前砂質(zhì)顆?？臻g非均質(zhì)性對滲透率演化規(guī)律影響的研究較少，非均質(zhì)性對滲透率變化的機理仍不明確。為此，本文在孔隙尺度探究水合物飽和度、砂質(zhì)顆粒粒徑和砂質(zhì)顆?？臻g非均質(zhì)性3種因素對滲透率演化規(guī)律的影響。

當前格子玻爾茲曼方法(LBM)、孔隙網(wǎng)絡(luò)模型(PNM)和有限單元法(FEM)是研究水合物生成過程中的滲透率演化規(guī)律中常用的方法。HOU 等[14]基于LBM 方法從水合物飽和度、礦物顆粒排列、水合物生長及分布形態(tài)等方面分析了水合物對多孔介質(zhì)滲透率變化的影響；ZHANG等[15]基于LBM方法建立了考慮分解動力學、氣體流動和傳熱機理的孔隙尺度框架模型，研究了甲烷水合物分解過程中非等溫多重物理化學過程。但LBM 方法無法表征砂粒粒徑，并不能有效分析粒徑因素對滲透率演化規(guī)律的影響。MAHABADI 等[16]建立了孔隙網(wǎng)絡(luò)模型，模擬在初始飽和度為10.0%～60.0%時的水合物分解，計算相對水滲透率和相對氣滲透率，并修正了STONE方程的擬合參數(shù)。雖然PNM方法可以很好地描述多孔介質(zhì)宏觀尺度屬性，但由于簡化了多孔介質(zhì)，因此丟失了一些水合物分布形態(tài)的信息。由于上述2種方法并不能有效研究砂質(zhì)顆粒粒徑及水合物形貌對含水合物沉積物滲透率變化的影響，而基于有限單元法的COMSOL Multiphysics 在含水合物沉積物砂質(zhì)顆粒和水合物形貌構(gòu)建方面更為有效，本文采用基于有限單元法的COMSOL Multiphysics模擬器，求解含水合物沉積物滲透率演化。

本文重點討論砂質(zhì)顆粒粒徑、水合物飽和度和砂質(zhì)顆?？臻g非均質(zhì)性對含水合物沉積物滲透率的影響。首先，采用基于有限單元法的COMSOL Multiphysics分析上述因素對滲透率演化規(guī)律的影響，其中，非均質(zhì)幾何模型顆粒及位置遵循正態(tài)分布并通過MATLAB 程序生成；其次，為驗證模型的正確性，將數(shù)值模擬結(jié)果與實驗數(shù)據(jù)進行對比分析，并通過修正Tokyo模型對本文數(shù)值模擬結(jié)果進行擬合，分析均質(zhì)模型與非均質(zhì)模型下修正Tokyo模型參數(shù)的差異。

1 幾何模型及計算方法

1.1 砂質(zhì)顆粒粒徑劃分

為研究砂質(zhì)顆粒粒徑對水合物沉積物滲透率的影響，參考文獻[17]，將幾何模型顆粒劃分為4種不同范圍的粒徑尺寸，并選擇范圍上限作為模擬粒徑，如表1所示。

表1 砂質(zhì)顆粒粒徑劃分[17]Table 1 Sand particle size classification[17]

1.2 含水合物沉積物幾何模型

在含水合物沉積物儲層中，天然氣水合物呈現(xiàn)3種不同的形態(tài)：表層覆蓋、孔隙填充和顆粒膠結(jié)。圖1所示為含水合物沉積物CT 圖像。由圖1(c)和(d)可見，表層覆蓋水合物通常在實驗前期即在孔隙中生成；由圖1(b)可知孔隙填充和顆粒膠結(jié)水合物主要在實驗中后期被觀察到，且為水合物存儲的主要形態(tài)。本文主要研究孔隙填充型水合物的生成對滲透率演化規(guī)律的影響，同時，數(shù)值模擬均以球形表征沉積物固體顆粒和孔隙填充型水合物形狀。在模型中，將砂質(zhì)顆粒包裹，且將未被水合物填充的空間定義為孔隙，砂質(zhì)顆粒間相互連通相對狹窄的空間定義為孔喉。

圖1 含水合物沉積物CT圖像[18]Fig.1 CT images of hydrate-bearing sediments[18]

為探究砂質(zhì)顆?？臻g非均質(zhì)性對滲透率演化規(guī)律的影響，分別建立均質(zhì)幾何模型與非均質(zhì)幾何模型，如圖2所示。由圖2(a)可見：砂質(zhì)顆粒(深灰色)與水合物(紅色)均勻分布，砂質(zhì)顆粒間距ds=a×D50/2，其中a為砂質(zhì)顆粒粒徑與顆粒間距的比。由于模型砂質(zhì)顆粒及水合物幾何關(guān)系的限制，當球形水合物直徑與砂質(zhì)顆粒間距相等時，水合物飽和度達到最大值，但由于孔隙未被完全填充及孔喉部分無水合物的原因，本文水合物飽和度SH最大值為50.0%。由圖2(b)可見：沉積物中D50=125.0 μm顆粒(深灰色)體積占全部顆粒的體積分數(shù)為70.0%，D50=62.5 μm顆粒(淺灰色)體積占全部顆粒體積的30.0%，D50=5.0 μm 水合物(紅色)體積占孔隙體積定義為飽和度SH，2種砂質(zhì)顆粒粒徑及水合物尺寸均遵循正態(tài)分布，位置則是在保證互不重疊基礎(chǔ)上的隨機分布；非均質(zhì)幾何模型初始孔隙度φ及邊界長度L與D50=125.0 μm 的均質(zhì)幾何模型的一致：φ=35.0%，L=0.89 mm，W=0.38 mm。

圖2 基于CT圖像重構(gòu)的空間均質(zhì)與非均質(zhì)幾何模型(孔隙填充型)Fig.2 Constructed homogeneous and heterogeneous geometric models(pore-filling type)

1.3 計算方法

本文模型采用二維模型，其中幾何模型面積S1表示為

式中：L和W分別為模型的長度和寬度，μm；nx和ny分別為幾何模型x和y方向的顆粒數(shù)目；D50為砂質(zhì)顆粒直徑，μm；ds為砂質(zhì)顆粒間距，μm。

沉積物顆粒面積S2及水合物面積S3均基于COMSOL Multiphysics中的域積分求得，因此，含水合物幾何模型初始孔隙度φ0及飽和度SH分別表示為：

數(shù)值模擬中將含水合物沉積物砂質(zhì)顆粒、水合物顆粒表面及模型上下邊界均設(shè)置為無滑移邊界，即|u|= 0。同時在模型的左邊界設(shè)置入口壓力Pu；右邊界設(shè)置出口壓力Pd。模型中的流動區(qū)域基于網(wǎng)格劃分中的邊界層來劃分，模型網(wǎng)格采用超細化單元，模型邊界條件及網(wǎng)格示意如圖2(a)所示。在上述模型中，進行如下假設(shè)：

1)流體為不可壓縮流體；

2)流體流動是由壓差驅(qū)動的單相流動；

3)忽略重力的影響。

根據(jù)達西定律，COMSOL Multiphysics計算流體滲透率的表達式為：

式中：k為滲透率，m2；μ為黏度，Pa·s；Qg為氣體流量，m3/s；A為樣品橫截面積，m2；Pu和Pd分別為上游壓力和下游壓力，Pa。為確保幾何模型流體流態(tài)為層流，數(shù)值模擬采用的上下游壓差為500.0 Pa，雷諾數(shù)Re為1.27。同時，流速u=Qg/A，表示為幾何模型出口端的平均速度。數(shù)值模型輸入?yún)?shù)及數(shù)值模擬的計算流程分別如表2和圖3所示。

圖3 數(shù)值模擬計算流程Fig.3 Flow chart of numerical modelling process

表2 數(shù)值模擬參數(shù)Table 2 Parameters of numerical simulation

2 含水合物沉積物滲透率的影響因素

基于建立的幾何模型與數(shù)值模擬計算流程，探究砂質(zhì)顆粒粒徑、水合物飽和度及砂質(zhì)顆?？臻g非均質(zhì)性3種因素對含水合物沉積物滲透率演化的影響。

2.1 砂質(zhì)顆粒粒徑

不同粒徑下含水合物沉積物絕對滲透率與飽和度的關(guān)系如圖4(a)所示。由圖4(a)可見：在無水合物情況下，黏土顆粒至中砂顆粒4種沉積物絕對滲透率分別為1.06×10-16，1.04×10-13，1.65×10-12和6.41×10-12m2，顆粒粒徑對滲透率變化影響顯著，中砂型(D50=500.0 μm)水合物沉積物絕對滲透率是黏土型(D50=2.0 μm)沉積物的6×104倍；為量化粒徑對滲透率的影響，將飽和度SH=41.1%時的各粒徑含水合物沉積物絕對滲透率與相對滲透率進行對比分析，如圖4(b)所示。由圖4(b)可見：在一定飽和度下，隨著沉積物顆粒粒徑增長，沉積物絕對滲透率呈非線性上升趨勢，不同粒徑間的相對滲透率差異很小；其中粉砂型(D50=62.5 μm)沉積物滲透率至中砂型(D50=500.0 μm)沉積物絕對滲透率以3.0～4.0 倍的增長率逐漸上升，而黏土型(D50=2.0 μm)與粉砂型(D50=62.5 μm)沉積物絕對滲透率則相差962 倍，D50=62.5 μm 時相鄰粒徑滲透率變化加劇，因此，粉砂(D50=62.5 μm)為影響滲透率變化速率的拐點粒徑。

圖4 顆粒粒徑對水合物沉積物滲透率演化規(guī)律的影響Fig.4 Effect of particle size on permeability of hydratebearing sediments

上述結(jié)果表明，沉積物中的黏土體積分數(shù)對水合物儲層的滲透率有較大影響，儲層黏土體積分數(shù)將會直接影響到水合物開采的效率。根據(jù)文獻[19]，我國南海神狐海域某區(qū)塊沉積物顆粒屬于泥質(zhì)粉砂型[20]。

2.2 水合物飽和度

為量化研究飽和度對滲透率的影響，以各粒徑無水合物時的滲透率作為初始滲透率，對不同粒徑下的滲透率進行歸一化處理，如圖5所示，由圖5可見：相對滲透率均隨飽和度呈非線性下降的趨勢；隨著粒徑增加，不同飽和度間的滲透率差異不大，這也許是均質(zhì)模型本身的局限性(如隨飽和度增加，迂曲度變化較小等)導致的，同時，本文模擬結(jié)果與XU 等[21]的模擬結(jié)果一致。從圖5可以看出：隨著飽和度增加，滲透率間的差異逐漸增大；當飽和度介于0～33.9%時，含水合物沉積物滲透率隨飽和度呈線性降低近1個數(shù)量級，而當飽和度高于33.9%時，含水合物沉積物滲透率則呈指數(shù)快速下降近3個數(shù)量級。因此，SH=33.9%為飽和度對滲透率影響的拐點，當飽和度高于此拐點時，含水合物沉積物滲透率將快速下降。

圖5 飽和度對水合物沉積物滲透率演化規(guī)律的影響Fig.5 Effect of hydrate saturation on permeability of hydrate-bearing sediments

為探究飽和度對滲透率的影響，繪制不同飽和度下幾何模型的流速場，如圖6所示。由圖6可見：隨著水合物飽和度逐漸增加，沉積物顆粒與水合物間通道逐漸變窄；當飽和度從12.5%增至41.1%時，沉積物顆粒與水合物間孔喉半徑R1與沉積物顆粒間孔喉半徑R2的比(Rr=R1/R2)從1.64 降至0.35，含水合物沉積物流動通道堵塞并生成新的孔喉通道，致使沉積物孔隙通道流速降低。隨著飽和度逐漸增加，兩孔喉流速的比值(vr=v1/v2)從0.64(SH=12.5%)升至3.12(SH=41.1%)，孔喉半徑R1逐漸成為流速控制的孔隙通道，這種流速與孔喉呈非線性下降趨勢，導致飽和度SH高于33.9%時，沉積物滲透率急劇下降。

圖6 不同水合物飽和度下孔喉及孔喉流速場演化規(guī)律Fig.6 Evolution of pore throat and associated flow field against hydrate saturation

基于上述現(xiàn)象，我們推斷在水合物沉積物生成過程中，若前期表層覆蓋型水合物率先生成，孔喉通道被迅速堵塞。儲層滲透率將快速下降，導致生成水合物飽和度較低。而當孔隙填充型水合物占主導地位時，儲層滲透率下降較慢，儲層水合物飽和度較高；最后，當孔隙填充型水合物逐漸生長而與顆粒間生成新的孔喉時，滲透率將進一步降低。因此，含水合物沉積物儲層滲透率的演化規(guī)律與水合物生長模式密切相關(guān)。

2.3 砂質(zhì)顆?？臻g非均質(zhì)性

基于圖1(a)的CT 圖像可以看出，儲層顆粒及水合物粒徑及位置均為隨機分布。為更加準確地描述水合物生成過程中儲層滲透率的演化規(guī)律，以D50=125.0 μm 的顆粒體積分數(shù)為70.0%，D50=62.5 μm 的顆粒體積分數(shù)為30.0%組成非均質(zhì)幾何模型，并將其與D50=125.0 μm的均質(zhì)幾何模型進行對比分析，探究儲層顆粒、水合物分布非均質(zhì)性對滲透率演化的影響，其中，均質(zhì)幾何模型初始滲透率k0=3.96×10-12m2，非均質(zhì)幾何模型初始滲透率k0=1.35×10-12m2，2 種模型的相對滲透率結(jié)果如圖7所示。

圖7 水合物沉積物(空間均質(zhì)與非均質(zhì)分布)滲透率隨飽和度演化規(guī)律及模型驗證Fig.7 Relationship between relative permeability of hydrate-bearing sediments with hydrate saturation and model validation

從圖7可以看出：2 種模型的滲透率變化速率有很大差異，非均質(zhì)模型在隨飽和度增加過程中變化速率逐漸減小，當SH=20.0%時，相對滲透率已有一個數(shù)量級的變化；均質(zhì)模型滲透率變化速率則隨飽和度增加逐漸增大，當SH=38.1%時，相對滲透率降低1個數(shù)量級。上述結(jié)果表明：非均質(zhì)性會顯著降低初始滲透率且滲透率對飽和度變化更加敏感。同時，將MINGAWA 等[22]、LIANG等[23]、DELLI等[24]、KUMAR等[25]、LEI等[26]的實驗數(shù)據(jù)與數(shù)值模擬結(jié)果進行對比，發(fā)現(xiàn)相對滲透率實驗數(shù)據(jù)位于2條曲線之間，其中，均質(zhì)模型曲線為上邊界，非均質(zhì)模型曲線為下邊界。

本文中，由于均質(zhì)幾何模型最大飽和度有一定適用范圍(SHmax=50.0%)，根據(jù)圖6所示流速場可以看出，當SH=41.1%時，流動通道幾乎被完全堵塞，而孔隙內(nèi)仍有部分空間未填充水合物；同時，其他學者實驗數(shù)據(jù)中水合物最大飽和度約為60.0%[25]。因此，對Tokyo 模型中SH項添加系數(shù)β來量化實際飽和度與理想化SH=100.0%間的差異，修正后的模型表示為

采用絕對平均相差偏差kAARD評價修正的Tokyo模型與數(shù)值模擬擬合相對滲透率(k/k0)的效果，kAARD表示為

式中：為修正Tokyo 模型計算值；yi為數(shù)值模擬計算值。

在均質(zhì)模型中，隨水合物飽和度增加，當飽和度高于33.9%時含水合物沉積物絕對滲透率有突變的現(xiàn)象，如圖5所示。為此，以SH=33.9%為分界點，利用式(5)進行分段擬合。由于水合物分布形態(tài)的差異，非均質(zhì)模型中并未出現(xiàn)滲透率突變的現(xiàn)象，為此，非均質(zhì)模型曲線未進行分段擬合，2 個模型的擬合結(jié)果如圖7所示，均質(zhì)與非均質(zhì)Tokyo模型參數(shù)如表3所示。

表3 修正的Tokyo模型擬合參數(shù)Table 3 Fitting parameters of corrected Tokyo models

由于非均質(zhì)模型中孔隙填充水合物粒徑僅5 μm，且位置隨機分布，無限定飽和度適用范圍，Tokyo模型擬合效果良好，而均質(zhì)模型由于幾何模型最大飽和度有一定適用范圍，設(shè)定β=2的修正模型擬合效果更好。當考慮水合物在孔喉處生長時應(yīng)適當增加N，但原因仍不明確。在本文中，非均質(zhì)模型中N高于均質(zhì)模型，因此，除水合物在孔喉處生長外，N升高與模型的非均質(zhì)性有關(guān)。

3 結(jié)論

1)含水合物沉積物顆粒粒徑是影響滲透率演化的重要因素，中砂型(D50=500.0 μm)與黏土型(D50=2.0 μm)沉積物絕對滲透率相差5 個數(shù)量級；在一定飽和度下，沉積物絕對滲透率隨砂質(zhì)顆粒粒徑增大呈非線性上升趨勢，且不同粒徑間的相對滲透率差異很??；砂質(zhì)顆粒粒徑下降至拐點(D50=62.5 μm)后，絕對滲透率急劇下降，粉砂水合物沉積物絕對滲透率是黏土水合物沉積物的962倍(SH=41.1%)，儲層黏土體積分數(shù)將會直接影響到水合物開采的效率。

2)含水合物沉積物相對滲透率均隨飽和度呈非線性下降的趨勢，儲層流動通道逐漸變窄，兩孔喉流速比vr與孔喉半徑比Rr呈非線性下降的趨勢，SH=33.9%為飽和度對滲透率影響的拐點，當飽和度高于此拐點時，含水合物沉積物滲透率快速下降。

3)模型非均質(zhì)性會顯著降低滲透率且滲透率變化對飽和度更加敏感，非均質(zhì)模型中修正Tokyo模型的指數(shù)N比均質(zhì)模型的大，N升高與非均質(zhì)性有關(guān)；現(xiàn)有的相對滲透率實驗數(shù)據(jù)通常位于2條曲線之間，其中均質(zhì)模型曲線為上邊界，非均質(zhì)模型曲線為下邊界。

4)由于實驗測試結(jié)果間相差較大，后續(xù)將進行實驗探究，進一步分析均質(zhì)幾何模型與非均質(zhì)幾何模型的差異。