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

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

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

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

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

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

1 幾何模型及計算方法

1.1 砂質顆粒粒徑劃分

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

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

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

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

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

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

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

1.3 計算方法

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

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

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

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

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

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

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ù)值模擬計算流程，探究砂質顆粒粒徑、水合物飽和度及砂質顆?？臻g非均質性3種因素對含水合物沉積物滲透率演化的影響。

2.1 砂質顆粒粒徑

不同粒徑下含水合物沉積物絕對滲透率與飽和度的關系如圖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

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

2.2 水合物飽和度

為量化研究飽和度對滲透率的影響，以各粒徑無水合物時的滲透率作為初始滲透率，對不同粒徑下的滲透率進行歸一化處理，如圖5所示，由圖5可見：相對滲透率均隨飽和度呈非線性下降的趨勢；隨著粒徑增加，不同飽和度間的滲透率差異不大，這也許是均質模型本身的局限性(如隨飽和度增加，迂曲度變化較小等)導致的，同時，本文模擬結果與XU 等[21]的模擬結果一致。從圖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ī)律與水合物生長模式密切相關。

2.3 砂質顆?？臻g非均質性

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

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

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

本文中，由于均質幾何模型最大飽和度有一定適用范圍(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ù)值模擬計算值。

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

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

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

3 結論

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

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

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

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