基于流場-電場耦合的含水合物模擬巖心電學特性數值實驗

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

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

0 引 言

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

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

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

1 巖電數值模擬實驗

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

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

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

圖1 兩端帶有平板和環形電極的模擬巖心

2 耦合數值模型的建立

2.1 幾何結構

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

圖2 二維多孔介質模型

2.2 邊界條件和材料物性參數

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

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

2.3 物理模型和求解方法

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

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

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

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

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

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

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

2.4 網格剖分與選擇

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

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

圖3 各截面處電勢均值隨網格數變化曲線

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

3.1 壁面接觸角

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

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

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

3.2 界面張力

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

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

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

3.3 油水黏度比

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

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

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

4 飽和度對電場響應特性的影響

4.1 電極形狀的選擇

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

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

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

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

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

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

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

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

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

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

4.3 含流體飽和度的影響

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

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

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

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

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

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

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

5 結 語

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

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

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