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

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

0 引 言

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

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

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

1 巖電數(shù)值模擬實(shí)驗(yàn)

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

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

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

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

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

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

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

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

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

采用COMSOL中CFD 模塊模擬多孔介質(zhì)孔隙中水驅(qū)替油的兩相流動(dòng)過(guò)程。數(shù)值模型采用速度入口和壓力出口，在入口處設(shè)置注入水流速，在出口處設(shè)置壓力值。初始狀態(tài)時(shí)，多孔介質(zhì)模型處于油飽和狀態(tài)，驅(qū)替時(shí)在入口處注入不同流速（如7.81、31.24 mm／s等）的水來(lái)驅(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模塊計(jì)算多孔介質(zhì)區(qū)域的電場(chǎng)參數(shù)。模型中的電極材料為銅，將油視為不導(dǎo)電的物質(zhì)，其電導(dǎo)率設(shè)置為0 S／m，相對(duì)介電常數(shù)設(shè)置為2.2，水的電導(dǎo)率設(shè)置為0.05 S／m，相對(duì)介電常數(shù)為80。在實(shí)際巖電實(shí)驗(yàn)中所采用的激勵(lì)電壓通常為幅值1 ～5 V的正弦波信號(hào)，因此數(shù)值模型中在電極間施加幅值為1 V的正弦波交流電壓。

2.3 物理模型和求解方法

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

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

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

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

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

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

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

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

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

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

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

3 微觀因素對(duì)相分布特性的影響

3.1 壁面接觸角

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

圖4 壁面接觸角不同時(shí)的相分布圖

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

3.2 界面張力

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

圖5 不同界面張力時(shí)相分布圖

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

3.3 油水黏度比

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

圖6 不同黏度比時(shí)相分布圖

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

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

4.1 電極形狀的選擇

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

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

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

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

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

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

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

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

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

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

4.3 含流體飽和度的影響

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

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

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

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

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

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

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

5 結(jié) 語(yǔ)

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

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

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