邢蘭昌， 齊淑英， 李 馨， 魏 偉， 韓維峰

（1.中國石油大學(xué)（華東）控制科學(xué)與工程學(xué)院，山東青島266580；2.中國石油勘探開發(fā)研究院新能源研究所，河北廊坊065007）

0 引 言

天然氣水合物（簡稱水合物）是水與天然氣（主要成分為甲烷）在高壓低溫條件下形成的類冰狀結(jié)晶物質(zhì)［1］。自然界中水合物廣泛分布于內(nèi)陸湖的深水環(huán)境與凍土帶和海底0.1 ～1.0 km以下的大陸斜坡和深海盆地的沉積物［2-4］，其在沉積物中以多種分布形式存在，宏觀上有分散型、塊狀、脈狀等；微觀上有懸浮、接觸、膠結(jié)等，多種空間分布形式導(dǎo)致含水合物沉積物等物理性質(zhì)呈現(xiàn)各向異性［5-6］，從而為刻畫含水合物沉積物的性質(zhì)以及評價含水合物飽和度帶來了困難。

為了深入研究含水合物沉積物的電學(xué)特性以及電學(xué)參數(shù)與含水合物飽和度之間的關(guān)系，研究者對野外鉆探獲取的巖心或者人工合成的樣品進行測試。但是，由于水合物通常賦存于固結(jié)程度較弱的沉積物中，而且水合物穩(wěn)定性對溫度壓力條件的依賴性高，所以野外鉆探取樣技術(shù)非常復(fù)雜、取樣難度大且成本高。實驗室內(nèi)制備人工樣品并進行測試是一種有效的研究手段。在實驗過程中已逐步采用多種監(jiān)測技術(shù)以盡可能直觀地獲取水合物的空間分布特征［7-10］，如CT、ECT等。但是CT技術(shù)要求樣品小，ERT 技術(shù)分辨率受到限制。此外，開展水合物模擬實驗測試耗時較長、可模擬的實驗條件有限、探測信息較少（如無法獲知物理場參數(shù)的空間分布信息），巖石物理數(shù)值模擬為克服上述困難提供了一條可行的途徑。

本文基于流場-電場耦合的有限元數(shù)值計算方法研究含水合物模擬巖心的電學(xué)響應(yīng)特性。通過流場計算獲得多孔介質(zhì)中油／水空間分布來模擬水合物／水在巖心中的微觀分布，通過電場計算獲得含水合物／水模擬巖心的電學(xué)特性，最終形成一套針對含水合物模擬巖心開展巖電數(shù)值模擬實驗的方法。

1 巖電數(shù)值模擬實驗

首先需要獲得具有不同含水飽和度的模擬巖心樣品，然后對巖心樣品進行測試獲得電學(xué)特性參數(shù)。在實際巖心中建立不同含水飽和度的常用方法有油／氣驅(qū)法、離心法、烘／風(fēng)干法、半滲透隔板法和自吸增水法等［11-14］。

本研究數(shù)值模擬實驗中采用驅(qū)替法來獲得具有不同含水飽和度的模擬巖心。考慮到油與水合物電學(xué)性質(zhì)和微觀分布狀態(tài)的相似性，利用水驅(qū)油過程所形成的油／水分布來模擬水合物／水在模擬巖心孔隙中的微觀分布狀態(tài)?；贑OMSOL Multiphysics 平臺建立含油水模擬巖心的流場-電場耦合數(shù)值仿真模型。

圖1 所示為帶有平板電極和環(huán)形電極的模擬巖心。首先將模擬巖心進行油飽和，然后固定出口壓力，以恒定的流速從入口處注入水，對模擬巖心中的油進行驅(qū)替。通過調(diào)節(jié)模擬巖心和流體的性質(zhì)以及驅(qū)替時間來控制油水兩相在孔隙空間中的分布狀態(tài)和含油／水飽和度；在驅(qū)替過程中采用上述兩種類型的電極對處于不同狀態(tài)的模擬巖心進行電學(xué)參數(shù)測量，進而分析電學(xué)參數(shù)與含油飽和度之間的關(guān)系。

圖1 兩端帶有平板和環(huán)形電極的模擬巖心

2 耦合數(shù)值模型的建立

2.1 幾何結(jié)構(gòu)

建立如圖2 所示的長18 mm、寬9 mm的二維多孔介質(zhì)模型來模擬柱塞狀巖心，其中骨架顆粒在多孔介質(zhì)內(nèi)以等徑圓形交叉排列，圓的直徑設(shè)定為1 mm。模型入口位于x ＝0 mm處，出口位于x ＝18 mm處。

圖2 二維多孔介質(zhì)模型

2.2 邊界條件和材料物性參數(shù)

采用COMSOL中CFD 模塊模擬多孔介質(zhì)孔隙中水驅(qū)替油的兩相流動過程。數(shù)值模型采用速度入口和壓力出口，在入口處設(shè)置注入水流速，在出口處設(shè)置壓力值。初始狀態(tài)時，多孔介質(zhì)模型處于油飽和狀態(tài)，驅(qū)替時在入口處注入不同流速（如7.81、31.24 mm／s等）的水來驅(qū)替多孔介質(zhì)內(nèi)的油。模型中設(shè)定水和油的典型物性參數(shù)為：水密度ρwater＝998.2 kg／m3、水黏度μwater＝1 mPa·s，油密度ρoil＝880 kg／m3，油黏度μoil＝0.02 Pa·s。

采用COMSOL中的AC／DC模塊計算多孔介質(zhì)區(qū)域的電場參數(shù)。模型中的電極材料為銅，將油視為不導(dǎo)電的物質(zhì)，其電導(dǎo)率設(shè)置為0 S／m，相對介電常數(shù)設(shè)置為2.2，水的電導(dǎo)率設(shè)置為0.05 S／m，相對介電常數(shù)為80。在實際巖電實驗中所采用的激勵電壓通常為幅值1 ～5 V的正弦波信號，因此數(shù)值模型中在電極間施加幅值為1 V的正弦波交流電壓。

2.3 物理模型和求解方法

水與油在模擬巖心孔隙中流動時，兩相之間存在復(fù)雜的動態(tài)相界面。數(shù)值模型中采用的界面追蹤方法主要有相場法和水平集法［15-16］。其中相場法采用連續(xù)的變量來模擬流體流動的不連續(xù)現(xiàn)象，能夠捕捉兩相流動過程中的相界面，對于尺寸較小的幾何模型（微流體仿真），相場法比水平集法對界面識別的精確度更高，因此本研究的數(shù)值模型采用相場法。通過求解Cahn-Hillard方程實現(xiàn)對兩相界面的跟蹤［17］，即

式中：G 為化學(xué)勢能；γ 為控制Cahn-Hillard 方程擴散尺度的遷移系數(shù)；φ為相場參數(shù)，表示油時相場參數(shù)的值為1，表示水時相場參數(shù)的值為－1，表示兩相界面時相場參數(shù)的值為0。

模擬巖心中流體流動通道尺寸小、流速低，通過計算雷諾數(shù)得知模擬巖心中的流態(tài)為層流，因此選用層流模型對流場參數(shù)進行計算。在微小流道中通常不考慮重力的影響，此時黏性力和慣性力支配著流體在模擬巖心中的流動狀態(tài)。假設(shè)水和油均為不可壓縮流體，水與油互不相溶，由流體力學(xué)理論得到描述流動過程的動量守恒方程和質(zhì)量守恒方程［18］：

式中：v 為流體的速度矢量，m／s；ρ 為流體的密度，kg／m3；I為單位矩陣，τ 為黏性切應(yīng)力張量，Pa；Fst為單位體積的界面張力矢量，N／m。

在分析流場與電場耦合問題時，考慮到模型中所采用的電極激勵電壓幅值較低（如1 V），不足以使油相和水相電離而影響相分布，因此不考慮電場對油水相分布的影響，數(shù)值模型中的物理場耦合為流場對電場影響的單向耦合分析。忽略電場中磁效應(yīng)的影響，由電荷守恒方程得［19］：

式中：σ為電導(dǎo)率，S／m；qv為體電荷密度，C／m3；u 為流速，m／s；E為電場強度，V／m。

采用MUMPS（Multi-frontal Massively Parallel Sparse direct Solver）算法、瞬態(tài)求解器對流場和電場控制方程進行求解，將多孔介質(zhì)內(nèi)的流場和電場進行耦合時選用層流物理場接口、相場物理場接口和電流物理場接口。

2.4 網(wǎng)格剖分與選擇

由于模擬巖心中流體通道不規(guī)則，所以采用自由剖分的三角形網(wǎng)格。網(wǎng)格劃分的精細程度影響計算結(jié)果的準確度和所需的計算時間，因此需要討論網(wǎng)格數(shù)的影響。針對圖2 所示的多孔介質(zhì)模型，共考察了6種具有不同單元數(shù)量的網(wǎng)格，即6 452、12 066、14 408、20 496、28 232 和38 858。以多孔介質(zhì)模型的入口作為起點（見圖2），截取多孔介質(zhì)x ＝1.5、4.5、7.5、10.5、13.5 mm 5 個截面，由于多孔介質(zhì)內(nèi)相分布的變化將導(dǎo)致電場中電勢分布的變化，因此選取不同截面處電勢的均值來討論數(shù)值模擬結(jié)果對網(wǎng)格數(shù)量的依賴性。

圖3 顯示了不同截面處電勢均值隨網(wǎng)格數(shù)變化的曲線。由圖可知：當(dāng)網(wǎng)格數(shù)小于20 496 時，電勢均值隨著網(wǎng)格數(shù)的增加而變化，但是當(dāng)網(wǎng)格數(shù)從20 496 增至38 858 時，各截面處電勢均值趨于穩(wěn)定。綜合考慮計算結(jié)果與計算時間，最終選擇網(wǎng)格數(shù)為28 232。

圖3 各截面處電勢均值隨網(wǎng)格數(shù)變化曲線

3 微觀因素對相分布特性的影響

3.1 壁面接觸角

潤濕性影響模擬巖心孔隙中油／水分布和流動特性，進而影響模擬巖心的電學(xué)響應(yīng)特性。潤濕性可通過壁面接觸角θ來表示。當(dāng)θ ＜π／2 時，壁面親水；當(dāng)θ ＞π／2 時，壁面親油。圖4 為壁面接觸角不同時模擬巖心內(nèi)的油水相分布圖（藍色代表水，紅色代表油）。

圖4 壁面接觸角不同時的相分布圖

分析圖4 可知：當(dāng)模擬巖心顆粒具有親水壁面時（θ ＜π／2），水以薄膜的狀態(tài)分布于顆粒表面，注入水易沿著顆粒表面爬行，由于毛管力強，殘余油被水圈閉在孔隙中不能流動，其多以簇團狀的形態(tài)存在；當(dāng)顆粒具有親油壁面時（θ ＞π／2），油受到顆粒表面附著力的牽制，油在顆粒表面的流動受到阻礙，弱化了水對顆粒表面油膜的驅(qū)替作用，使水在孔隙中形成竄流，在孔隙中多處形成以柱狀形態(tài)存在的殘余油。

3.2 界面張力

界面張力為兩不相混溶的液體接觸界面產(chǎn)生的力，其大小與毛細管滯留阻力成正比例關(guān)系，影響油和水在模擬巖心孔隙中的微觀分布形態(tài)。將油和水之間界面張力γ12的變化范圍設(shè)置為1 ～3 mN／m，圖5 為界面張力不同時多孔介質(zhì)內(nèi)油水相分布圖（藍色代表水，紅色代表油）。

圖5 不同界面張力時相分布圖

分析圖5 可知：油水界面張力越大，水的驅(qū)替阻力越大；界面張力的增大使得將大油滴分散成小油滴所需的分散功升高，導(dǎo)致油的不連續(xù)性增強，水包裹油的現(xiàn)象越明顯。

3.3 油水黏度比

油水黏度比的不同會引起驅(qū)替過程中主導(dǎo)力的變化，進而影響模擬巖心中油水的流動特性和相空間分布。將油水黏度比M的變化范圍設(shè)定為10 ～50，圖6所示為油水黏度比不同時模擬巖心內(nèi)油水相分布圖（藍色代表水，紅色代表油）。

圖6 不同黏度比時相分布圖

分析圖6 可知：隨著油相黏度的增大，水驅(qū)替的阻力逐漸增大，模擬巖心孔隙中的油不連續(xù)性逐漸增強，相同驅(qū)替時間時殘余油飽和度隨之增加；由于模擬巖心上下不滲透邊界的影響，驅(qū)替過程中邊界處的流量大于中部的流體流量，導(dǎo)致油大部分存在于模擬巖心的中部；隨著油水黏度比的增大，油更傾向于以簇團狀的形式存在。

4 飽和度對電場響應(yīng)特性的影響

4.1 電極形狀的選擇

除了模擬巖心內(nèi)的流體性質(zhì)和空間相分布以外，電極的形狀和位置也影響電場參數(shù)的變化特性。

采用指標參數(shù)F 來衡量模擬巖心內(nèi)飽和單相水時的電場分布的均勻性［20］，F(xiàn)表示為

式中：Emax和Emean分別為模擬巖心內(nèi)部電場強度的最大值和平均值，V／m。

采用靈敏度S來衡量模擬巖心內(nèi)測量區(qū)域的靈敏場分布［21，22］，S表示為：

式中：Ewater是模擬巖心內(nèi)充滿水時兩電極之間的平均電場強度，V／m；E是驅(qū)替過程中某一狀態(tài)下兩電極之間的電場強度，V／m。

綜合考慮平板電極和環(huán)形電極電場均勻性和靈敏度，以下僅對電場均勻度和靈敏度均較高的平板電極進行討論。

4.2 基于圖像處理的飽和度獲取

通過求解數(shù)值模型獲得不同驅(qū)替時刻時模擬巖心中油水的相空間分布圖，然后采用圖像法對相分布圖像進行處理，獲得不同時刻時模擬巖心中的含水／油飽和度。通常采用閾值法對相分布二維圖像進行處理。

選定RGB（紅色、綠色、藍色）顏色空間，然后根據(jù)各顏色分量軸的變化范圍設(shè)定合理的閾值，當(dāng)?shù)? 個顏色分量值≥150 時設(shè)定為藍色，其余情況為紅色。通過統(tǒng)計各種顏色的像素數(shù)量即可得到模擬巖心中的含水／油飽和度。圖7 所示為一典型水驅(qū)油過程中含油飽和度So隨時間的變化圖。

圖7 模擬巖心中含油飽和度隨時間的變化

4.3 含流體飽和度的影響

模擬巖心初始狀態(tài)為油飽和，此時平板電極與模擬巖心形成一個等效電容器。初始狀態(tài)時油為連續(xù)相，隨著實驗進程的推進水逐漸進入模擬巖心并將部分油從出口驅(qū)出，在某一時刻模擬巖心內(nèi)部的水成為連續(xù)相，此時兩平板電極與模擬巖心形成一個等效電阻器。以下將實驗過程分為兩個階段進行討論，即油為連續(xù)相和水為連續(xù)相。

在油為模擬巖心中連續(xù)相的階段，等效電容器的容抗值受到兩個因素的影響，即頻率和電容。在電極與模擬巖心尺寸一定的前提下，電容取決于等效介電常數(shù)的大小，即受到含油／水飽和度的直接影響。模擬巖心的等效容抗值XC表示如下：

式中：f 為電場頻率，Hz；C 為電容，F(xiàn)。在驅(qū)替實驗的某一時刻，由于油水分布以及含油／水飽和度為一確定值，因此電容C已經(jīng)確定，則容抗僅受到電場頻率f的影響。

圖8 為電場頻率1 Hz ～1 MHz范圍內(nèi)，模擬巖心中含油飽和度分別為1.000、0.920、0.847、0.755 時的容抗值變化圖。分析圖8 并結(jié)合式（7）可知：隨著頻率的升高，模擬巖心的容抗降低；在電場頻率一定的條件下，含油飽和度越高則模擬巖心等效介電常數(shù)越小，導(dǎo)致電容降低和容抗升高，通過測量容抗值即可計算得到模擬巖心中的含油飽和度。

圖8 不同含油飽和度時容抗值隨頻率變化圖

在水為模擬巖心中連續(xù)相的階段，等效電阻器的電阻值主要取決于等效電導(dǎo)率，即取決于含油／水飽和度的大小。圖9 為電場頻率1 Hz ～1 MHz范圍內(nèi)，模擬巖心中含油飽和度分別為0.589、0.553、0.420、0.357、0.255、0 時的電阻值變化圖。分析圖9 可知：電場頻率對模擬巖心電阻值影響很??；當(dāng)頻率一定時，含油飽和度越高則電阻值越大，通過測量電阻值即可計算得到模擬巖心中的含油飽和度。

圖9 不同含油飽和度時電阻值隨頻率變化圖

5 結(jié) 語

本研究提出了一種基于流場-電場耦合的含水合物模擬巖心巖電數(shù)值實驗方法，包括模擬巖心建立和電學(xué)參數(shù)計算兩部分。

在構(gòu)建多孔介質(zhì)幾何結(jié)構(gòu)的基礎(chǔ)上進行流場和電場耦合計算。利用水驅(qū)油過程所形成的油／水分布來模擬水合物／水的微觀分布，通過改變壁面接觸角、界面張力、油水黏度比、驅(qū)替時間等因素實現(xiàn)對相空間分布和含水合物飽和度的調(diào)節(jié)，進而建立模擬巖心。利用平板電極在模擬巖心兩側(cè)施加1 Hz ～1 MHz頻率的交流電場，通過電場計算得到模擬巖心的等效容抗和等效電阻，從而獲得寬頻率范圍的電學(xué)響應(yīng)特性。

所提出的模擬巖心巖電數(shù)值仿真方法為含天然氣水合物和含油水巖石物理性質(zhì)的研究提供了新的途徑，為物理模擬實驗中電學(xué)傳感器的參數(shù)優(yōu)化、電學(xué)測試數(shù)據(jù)的分析提供了有力工具。