• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于激電法評價含水合物沉積物滲透率數(shù)值模擬研究*

    2023-08-31 08:52:50邢蘭昌張歡歡韓維峰楊金秀葛新民
    新能源進展 2023年4期
    關(guān)鍵詞:虛部激電水合物

    邢蘭昌,王 碩,張歡歡,魏 偉,韓維峰,楊金秀,葛新民

    基于激電法評價含水合物沉積物滲透率數(shù)值模擬研究*

    邢蘭昌1,?,王 碩1,張歡歡1,魏 偉2,韓維峰2,楊金秀3,葛新民3

    (1. 中國石油大學(xué)(華東) 控制科學(xué)與工程學(xué)院,山東 青島 266580;2. 中國石油勘探開發(fā)研究院 新能源研究所,河北 廊坊 065007;3. 中國石油大學(xué)(華東) 地球科學(xué)與技術(shù)學(xué)院,山東 青島 266580)

    含水合物沉積物樣品滲透率的實驗測量過程中存在成本高、周期長、水合物相態(tài)變化等問題,水合物開采過程中缺乏可對儲層滲透率動態(tài)變化進行實時監(jiān)測的技術(shù)。提出基于激電響應(yīng)評價含水合物沉積物滲透率的新方法,構(gòu)建含水合物多孔介質(zhì)的流場與電場有限元數(shù)值模型,探討水合物飽和度、骨架顆粒尺寸等因素對復(fù)電導(dǎo)率和滲透率的影響規(guī)律及機理,建立基于復(fù)電導(dǎo)率參數(shù)的含水合物多孔介質(zhì)滲透率評價模型。研究結(jié)果表明:(1)隨著多孔介質(zhì)骨架顆粒粒徑及等效孔徑的增大,復(fù)電導(dǎo)率譜的弛豫時間增大,復(fù)電導(dǎo)率虛部極大值所對應(yīng)頻率降低,雙電層極化主導(dǎo)頻率范圍內(nèi)的復(fù)電導(dǎo)率虛部增大,多孔介質(zhì)有效滲透率增大;(2)隨著孔隙填充型水合物的飽和度增加,雙電層極化主導(dǎo)的頻率范圍內(nèi)復(fù)電導(dǎo)率虛部減小,界面極化主導(dǎo)的頻率范圍內(nèi)復(fù)電導(dǎo)率虛部增大,含水合物多孔介質(zhì)有效滲透率及其變化率持續(xù)減小;(3)雙電層極化主導(dǎo)的頻率范圍內(nèi)復(fù)電導(dǎo)率虛部與滲透率之間存在確定性關(guān)系,在飽和水條件下選取的固定頻率處復(fù)電導(dǎo)率虛部與滲透率呈冪函數(shù)關(guān)系,結(jié)合飽和與非飽和條件下復(fù)電導(dǎo)率之間的關(guān)系以及歸一化滲透率與飽和度之間的關(guān)系可建立基于復(fù)電導(dǎo)率虛部的含水合物多孔介質(zhì)有效滲透率評價模型。研究結(jié)果可為開發(fā)基于激電法的含水合物沉積物樣品滲透率實驗測量技術(shù)以及天然氣水合物儲層滲透率動態(tài)監(jiān)測技術(shù)提供理論與模型基礎(chǔ)。

    天然氣水合物;滲透率;激發(fā)極化;復(fù)電導(dǎo)率;水合物飽和度;數(shù)值模擬

    0 引 言

    天然氣水合物儲量巨大,具有能量密度高、清潔無污染等優(yōu)點,被認(rèn)為是繼頁巖氣、煤層氣、致密氣之后最具開采潛力的戰(zhàn)略性接續(xù)能源[1-2]。我國已完成多次海域天然氣水合物試采,為將來實施“生產(chǎn)性試采”奠定了堅實的基礎(chǔ)。為了最終實現(xiàn)海域天然氣水合物的商業(yè)化開發(fā),仍需解決一系列的難題,如具有開采潛力水合物資源的評估、水合物甜點的精準(zhǔn)探測、水合物經(jīng)濟高效開采等[3]。我國南海天然氣水合物大多賦存于未固結(jié)成巖的粉砂、黏土質(zhì)粉砂、粉砂質(zhì)黏土等沉積物中,水合物儲層呈現(xiàn)高泥質(zhì)含量、低滲透性和高非均質(zhì)性等特點[4]。降壓試采結(jié)果表明,水合物儲層滲流能力差、天然氣產(chǎn)量低、儲層改造困難,尚未達(dá)到商業(yè)化開采要求。天然氣水合物降壓開采涉及熱量傳遞、水合物相變分解、氣?水兩相滲流、儲層變形等多個相互耦合的物理/化學(xué)過程,儲層的溫度、孔隙壓力、水合物/水/天然氣含量等持續(xù)發(fā)生變化。水合物分解吸熱和氣流焦耳?湯姆遜效應(yīng)等綜合作用會引起儲層局部溫度降低,可能導(dǎo)致地層水結(jié)冰以及水合物二次生成,此外水合物分解產(chǎn)出的水促使黏土發(fā)生體積膨脹,繼而堵塞流體流動通道,最終導(dǎo)致地層滲透率顯著降低。可見地層中水合物分解所引起的上述系列變化嚴(yán)重制約了天然氣產(chǎn)出效率的提升[5-7]。掌握天然氣水合物開采過程中儲層滲透率的實時變化規(guī)律對制定科學(xué)、高效的水合物開發(fā)方案,評價儲層改造效果,以及現(xiàn)場實時調(diào)整和動態(tài)優(yōu)化開發(fā)策略具有重要的意義。

    針對天然氣水合物儲層的滲透性評價問題,研究人員已經(jīng)開展了理論分析、數(shù)值模擬、實驗室測試以及現(xiàn)場原位測試等工作[5-14]。測試含水合物沉積物滲透率的方法可分為流動測試法、非流動測試法以及兩者相結(jié)合的方法。流動測試法又可分為穩(wěn)態(tài)法和非穩(wěn)態(tài)法,其特點在于測試過程涉及單相或多相流體的滲流過程,屬于侵入式、直接式測量方式。典型的非流動測試法包括聲波法、電阻率法、基于X射線計算機斷層成像(X-ray computed tomography,X-CT)技術(shù)以及基于核磁共振(nuclear magnetic resonance, NMR)技術(shù)的測試法,其特點在于測試過程中不涉及流體滲流過程,屬于非侵入式、間接式測量方式。非流動測試法對滲透率模型的可靠性要求較高。由于水合物對環(huán)境條件非常敏感,采用侵入式測量方式的流動測試法在實施過程中容易受到水合物分解、二次形成的影響,從而引入不可預(yù)測性測量誤差;而采用非流動測試法可以有效地降低上述影響,其中部分方法(如NMR測井法)在儲層原位滲透率評價方面更具優(yōu)勢。但是,由于NMR測井存在作業(yè)成本高、探測深度較淺等問題,在天然氣水合物開采過程中難以對儲層滲透率進行實時、長期的遠(yuǎn)探測。

    本文提出基于激電法定量評價天然氣水合物儲層滲透率的新思路。激電法,又稱為激發(fā)極化(induced polarization, IP)法,是一種以巖礦石的激電效應(yīng)、激電特性差異為物理基礎(chǔ)的電性勘探方法[15]。激電法可分為時間域激電法和頻率域激電法,前者在斷電后測量瞬態(tài)響應(yīng)隨時間的衰減特性,后者在供電期間記錄穩(wěn)態(tài)響應(yīng)隨頻率的變化特性。在一定頻率范圍(如1 mHz ~ 100 Hz)內(nèi)開展掃頻測量可以得到復(fù)電導(dǎo)率譜,通過分析復(fù)電導(dǎo)率譜來研究地層的激電特性,稱為譜激電法(spectral induced polarization, SIP)或復(fù)電導(dǎo)率(complex conductivity, CC)法。SIP或CC法是一類低頻地電方法,通過將SIP響應(yīng)參數(shù)與多孔介質(zhì)的結(jié)構(gòu)特征參數(shù)相關(guān)聯(lián),可實現(xiàn)對多孔介質(zhì)滲透率的定量評價,其基本假設(shè)為激電效應(yīng)主導(dǎo)的電極化(如低頻電化學(xué)極化)過程的長度尺度與控制滲透率的幾何長度尺度密切相關(guān)[16-18]。目前,低頻段電學(xué)阻抗譜方法在天然氣水合物相關(guān)實驗研究中得到了初步應(yīng)用。TZIRITA等[19]測試了四氫呋喃水合物及其與砂和高嶺石混合物的阻抗(50 ~ 1 000 Hz),實驗數(shù)據(jù)展示了水合物體系電學(xué)參數(shù)的頻率依賴特性。DU FRANE等[20]測試了甲烷水合物的阻抗譜(20 Hz ~ 2 MHz)并確定了水合物的電導(dǎo)率,隨后對水合物與石英砂和玻璃珠混合物(不含孔隙水)進行了測試[21],但是在數(shù)據(jù)分析時把低頻(如低于100 kHz)部分作為噪聲過大的無效數(shù)據(jù)進行了舍棄。邢蘭昌等[22]測試了四氫呋喃水合物生成和分解過程中的阻抗譜(10 mHz ~ 8 MHz),通過分析等效電路元件參數(shù)的變化規(guī)律探討了水合物生成和分解過程的阻抗譜參數(shù)特征。李棟梁等[23]利用以甲烷為主的混合氣為模擬氣在砂巖樣品中合成了水合物并測試了其電阻抗(0.05 kHz ~ 200 kHz),結(jié)果顯示含水合物巖石的介電常數(shù)在1 kHz ~ 100 kHz頻段呈現(xiàn)頻散特性。為了實現(xiàn)基于低頻(SIP頻率范圍)電阻抗譜對松散沉積物中水合物飽和度進行定量評價的目標(biāo),邢蘭昌等[24-29]開展了系列研究工作,設(shè)計了復(fù)電阻率/復(fù)電導(dǎo)率與聲學(xué)參數(shù)的聯(lián)合探測(ultrasound combined with electrical impedance, UCEI)系統(tǒng),在開展水合物模擬實驗測試的基礎(chǔ)上提取了用于構(gòu)建電聲聯(lián)合巖石物理模型的有效特征參數(shù),以寬頻復(fù)電導(dǎo)率譜及其導(dǎo)出參數(shù)為基礎(chǔ),結(jié)合泥質(zhì)修正阿爾奇公式、Simandoux公式、Cole-Cole模型、頻散度模型等建立了一系列水合物飽和度計算模型,形成了基于SIP的多孔介質(zhì)中水合物飽和度評價方法。

    在提出利用激電法定量評價含水合物沉積物滲透率新思路的基礎(chǔ)上,討論激電效應(yīng)參數(shù)與滲透率的關(guān)系模型,建立含水合物多孔介質(zhì)的流場與電場有限元數(shù)值模型,研究水合物飽和度、骨架顆粒尺寸等因素對復(fù)電導(dǎo)率和滲透率的影響規(guī)律及機理,基于此建立基于復(fù)電導(dǎo)率參數(shù)的含水合物多孔介質(zhì)滲透率評價模型,提出有待進一步解決的問題。本研究可為將來開發(fā)基于激電法的含水合物沉積物樣品滲透率實驗測量新技術(shù)以及水合物儲層滲透率動態(tài)監(jiān)測技術(shù)提供理論與模型基礎(chǔ)。

    1 基于激電法評價滲透率基本原理

    1.1 多孔介質(zhì)復(fù)電導(dǎo)率與滲透率

    激電響應(yīng)主要受到巖石的巖性所控制,直接感知礦物與流體之間界面的極化效應(yīng)[30]。通過對多孔介質(zhì)進行SIP測量可以獲得多孔介質(zhì)的復(fù)電導(dǎo)率,復(fù)電導(dǎo)率綜合了孔隙水導(dǎo)電、礦物顆粒表面導(dǎo)電、礦物顆粒?流體界面處電化學(xué)極化的貢獻。顆粒表面電導(dǎo)率可表示為一個復(fù)數(shù),其實部描述了電荷沿著顆粒?流體界面的電遷移現(xiàn)象,虛部刻畫了電荷在顆粒?流體界面處的極化過程。多孔介質(zhì)骨架顆粒表面電導(dǎo)率、內(nèi)比表面積等參數(shù)與孔隙尺寸密切相關(guān),而孔隙尺寸是用于預(yù)測滲透率的關(guān)鍵參數(shù)[31]。

    如式(1)所示,多孔介質(zhì)復(fù)電導(dǎo)率可以采用模值與相位角或者實部與虛部的形式來表示。

    式中:*為復(fù)電導(dǎo)率,S/m;||為復(fù)電導(dǎo)率模值,S/m;為相角,rad;′為復(fù)電導(dǎo)率實部,S/m;″為復(fù)電導(dǎo)率虛部,S/m。

    對于飽和水溶液的多孔介質(zhì),其復(fù)電導(dǎo)率實部可表示為

    式中:w為孔隙水的電導(dǎo)率,S/m;為地層因子,結(jié)合阿爾奇公式可將表示為

    在多孔介質(zhì)材料中,流體在壓力差的作用下發(fā)生運移,通??刹捎脻B透率來定量表示流體流動的難易程度。以達(dá)西定律為基礎(chǔ)可將多孔介質(zhì)中流體流動的量表示為

    當(dāng)多孔介質(zhì)的孔隙被單相流體所飽和時(如上述飽和水溶液的多孔介質(zhì)),式(6)所定義的滲透率為多孔介質(zhì)的絕對滲透率;當(dāng)多孔介質(zhì)孔隙中含有多相流體時,每相流體的滲透率稱為該相流體的有效(絕對)滲透率;有效滲透率與絕對滲透率的比值稱為相對滲透率或歸一化滲透率[33]。

    1.2 多孔介質(zhì)滲透率計算模型

    1.2.1 采用幾何尺度參數(shù)的模型

    將多孔介質(zhì)的不規(guī)則孔隙空間簡化為等徑平行毛細(xì)管束,則對于半徑為、長度為的根毛細(xì)管束,利用Hagen-Poiseuille方程可以得到總流量為

    對比式(6)和式(7)可得滲透率的計算式

    實際巖石多孔介質(zhì)的孔隙空間具有不規(guī)則性,流體實際流動路徑a比定義壓力梯度的直線距離更長。于是可引入迂曲度參數(shù)a= (a/)2,用實際路徑a替代式(7)中的直線距離可得

    (10)

    上述推導(dǎo)過程均立足于毛細(xì)管截面為圓形的假設(shè),對于非圓形截面的情況,將式(11)拓展為[34-37]

    式中:eff為等效水力半徑;s為形狀因子。

    對于具有單一粒徑尺寸的球形顆粒多孔介質(zhì)材料,其滲透率計算式可表示為

    式中:為顆粒直徑,m。CARMAN[38]指出對于均勻球形顆粒常數(shù)p= 4.8 ± 0.3,通常近似為5。

    式中:por為內(nèi)比表面積。考慮到巖石內(nèi)比表面積的分形特征,PAPE等[40]提出如下滲透率計算模型:

    上述滲透率模型中僅考慮代表幾何長度尺度的參數(shù)(包含電/水力迂曲度相等的情況),為了實現(xiàn)利用激電參數(shù)評價滲透率的目的,需要采用等效的地球物理長度尺度參數(shù)[44]。以下將基于SIP參數(shù)的滲透率評價模型分為兩類,第一類模型采用激電效應(yīng)的作用強度(以復(fù)電導(dǎo)率虛部或極化率來表示)參數(shù)來評價滲透率,第二類模型采用激電效應(yīng)的演化速度(以復(fù)電導(dǎo)率譜弛豫時間表示)參數(shù)來評價滲透率。

    1.2.2 采用激電效應(yīng)強度參數(shù)的模型

    基于對激電效應(yīng)強度與比表面積存在關(guān)聯(lián)的認(rèn)識,B?RNER等[30]對泥質(zhì)砂、粉砂和黏土進行了復(fù)電導(dǎo)率測量,證明了復(fù)電導(dǎo)率虛部與比表面積存在線性相關(guān)性,進而得到以下滲透系數(shù)()計算模型:

    式中:為擬合常數(shù)(下文相同);′′(0)為角頻率為0時復(fù)電導(dǎo)率的虛部,S/m;為滲透系數(shù),m/s。滲透系數(shù)與滲透率的關(guān)系為=/,其中:為流體的密度,kg/m3;為重力加速度,m/s2。

    SLATER等[45]研究了滲透系數(shù)模型對松散多孔介質(zhì)材料的適用性,結(jié)果表明地層因子變化范圍較小,因此將其忽略,僅利用復(fù)電導(dǎo)率虛部對滲透系數(shù)進行計算。

    針對B?rner模型適用巖性條件范圍較窄的問題,SLATER等[46]提出了如下滲透系數(shù)模型:

    根據(jù)水力半徑與電荷密度的關(guān)系、復(fù)電導(dǎo)率虛部與電荷密度的關(guān)系以及地層因子與孔隙度的關(guān)系,REVIL[47]建立了滲透率、地層因子與復(fù)電導(dǎo)率虛部之間的關(guān)系模型。

    WELLER等[31]測量了從不同地區(qū)獲取的純砂巖、泥質(zhì)砂巖樣品的復(fù)電導(dǎo)率,利用復(fù)電導(dǎo)率虛部和地層因子來計算滲透率。

    考慮到比表面積、復(fù)電導(dǎo)率虛部以及表面電導(dǎo)率等參數(shù)之間的比例關(guān)系,可將這些參數(shù)統(tǒng)一表示為,進而得到以下滲透率計算通式。

    1.2.3 采用激電效應(yīng)速度參數(shù)的模型

    SCOTT等[48]的研究結(jié)果表明,特征孔隙半徑與弛豫時間之間存在較強關(guān)聯(lián)。通過對未固結(jié)砂進行測試,BINLEY等[49]建立了弛豫時間與滲透系數(shù)之間的冪律關(guān)系式;此外,KEMNA等[50]的研究結(jié)果也表明弛豫時間與滲透率之間存在冪律關(guān)系:

    ZISSER等[51]建立了弛豫時間、地層因子與滲透率之間的關(guān)系模型。

    2 含水合物沉積物電場和流場數(shù)值模型建立

    2.1 數(shù)值建模方案

    1.2節(jié)中所引入的多孔介質(zhì)滲透率評價模型僅適用于多孔介質(zhì)被水溶液完全飽和的條件,若要對非飽和條件下多孔介質(zhì)的有效絕對滲透率進行評價,則需進一步考慮含水飽和度的影響。對于含水合物多孔介質(zhì)而言,可將其看作孔隙中部分水被水合物所替代的情況,即非飽和條件下的多孔介質(zhì)包括骨架顆粒、水合物和孔隙水三部分。

    基于多物理場耦合計算平臺COMSOL Multiphysics建立含水合物沉積物電場和流場的有限元數(shù)值模型。首先,建立含水合物多孔介質(zhì)的幾何結(jié)構(gòu),采用二維幾何結(jié)構(gòu)以節(jié)省數(shù)值計算量,采用圓面代表多孔介質(zhì)中的骨架顆粒和水合物顆粒;然后,計算并設(shè)定含水合物多孔介質(zhì)中各相的電參數(shù),如帶有雙電層的骨架顆粒、水合物顆粒、孔隙水;最后,在剖分網(wǎng)格的基礎(chǔ)上對電場和流場控制方程進行數(shù)值求解,最終計算不同條件下的復(fù)電導(dǎo)率和滲透率。

    2.2 幾何結(jié)構(gòu)與材料參數(shù)

    圖1所示為數(shù)值建模中所采用的代表性幾何結(jié)構(gòu),其中水合物以懸浮模式賦存于孔隙空間中。石英砂顆粒直徑為48 μm時,多孔介質(zhì)模型的長度為0.864 mm、寬度為0.432 mm。通過調(diào)節(jié)水合物顆粒的尺寸以實現(xiàn)對不同水合物飽和度條件的模擬。

    圖1 二維多孔介質(zhì)幾何結(jié)構(gòu)

    為了對流場和電場模型進行求解,需要對材料的特性參數(shù)進行設(shè)定。孔隙水的密度和黏度分別為1 000 kg/m3和0.001 Pa?s。采用式(28)對孔隙水的電導(dǎo)率進行計算,孔隙水相對介電常數(shù)設(shè)定為80。水合物的電導(dǎo)率和相對介電常數(shù)分別設(shè)定為1 × 10?5S/m和60[20,54-55]。采用文獻[27]中的方法對帶有雙電層的石英砂顆粒的等效電導(dǎo)率和等效介電常數(shù)進行計算。

    式中:f為孔隙水電導(dǎo)率,S/m;f為孔隙水鹽度,mol/L。

    2.3 流場/電場控制方程及邊界條件

    在流場模型中,不可壓縮流體的連續(xù)性方程為

    式中:為速度矢量,m/s。不可壓縮流體的動量守恒方程為

    式中:為壓力,Pa;為單位矩陣;為體積力矢量,N/m3;上標(biāo)T表示轉(zhuǎn)置運算。通過對以上流動控制方程進行數(shù)值求解,可以得到多孔介質(zhì)出口邊界處的質(zhì)量流量,進而可得到滲流速度。通過聯(lián)立式(31)和式(32)即可計算得到多孔介質(zhì)的絕對滲透率。

    式中:為滲流速度,m/s;為絕對滲透率,m2;?為壓力梯度,Pa/m;為多孔介質(zhì)的長度,m。

    在電場模型中,交變電場作用下總電流密度為傳導(dǎo)電流密度c與位移電流密度d之和。數(shù)值模型中求解的電場控制方程如下

    式中:為電導(dǎo)率,S/m;為電場強度矢量,V/m;為電位移矢量,C/m2;0為真空中的絕對介電常數(shù),取值為8.854 × 10?12F/m;為相對介電常數(shù)。

    在流場模型中,將流體入口以及出口均設(shè)置為壓力邊界,入口與出口之間的壓差設(shè)定為1 000 Pa。在電場模型中,對模型施加交流電場,即在模型左端設(shè)置正弦交流電壓= sin(),頻率范圍為1 × 103~ 1 × 106Hz,將模型右端接地。

    2.4 網(wǎng)格剖分與數(shù)值求解

    有限元網(wǎng)格單元的數(shù)量直接影響數(shù)值計算的穩(wěn)定性、模型求解所需的時間和算力、模型解的準(zhǔn)確度。以下利用模型解的網(wǎng)格依賴性檢驗來獲取合理的有限元網(wǎng)格單元數(shù)量。采用自由三角形網(wǎng)格將幾何區(qū)域離散化。流場求解采用廣義最小殘差法(generalized minimum residual method, GMRES);電場計算采用并行稀疏直接法(multi-frontal massivelyparallel sparse direct solver, MUMPS)。以石英砂顆粒直徑為48 μm為例,飽和水多孔介質(zhì)的復(fù)電導(dǎo)率虛部和絕對滲透率隨網(wǎng)格單元數(shù)變化的曲線如圖2所示。復(fù)電導(dǎo)率虛部和絕對滲透率的數(shù)值隨著單元數(shù)的增加而趨于穩(wěn)定,當(dāng)網(wǎng)格數(shù)量達(dá)到2.41 × 106時,網(wǎng)格單元最大和最小尺寸分別為0.70 μm和8.64 × 10?3μm。

    3 含水合物多孔介質(zhì)滲透率評價模型建立

    3.1 石英砂粒徑對復(fù)電導(dǎo)率/滲透率影響

    圖3和圖4分別展示了水合物飽和度為30%條件下石英砂顆粒粒徑不同時多孔介質(zhì)的復(fù)電導(dǎo)率虛部頻散曲線和有效絕對滲透率變化曲線。

    分析圖3可知:①復(fù)電導(dǎo)率虛部呈現(xiàn)隨頻率升高先增大后減小、繼而再增大的總體變化趨勢,雙電層極化和界面極化(Maxwell-Wagner效應(yīng))分別在低頻段(低于虛部取極大值的頻率)和高頻段(高于虛部取極小值的頻率)占據(jù)主導(dǎo)地位,在兩者之間的中頻段雙電層極化強度隨著頻率升高而減弱、界面極化隨著頻率升高而增強;②復(fù)電導(dǎo)率虛部極大值所對應(yīng)的頻率隨著粒徑的增大而降低,即弛豫時間隨著粒徑的增大而增大,已有的理論與實驗結(jié)果也表明雙電層極化弛豫時間與顆粒粒徑的平方成正比[56-57];③在低于虛部極大值所對應(yīng)頻率的頻率段范圍內(nèi)(雙電層極化主導(dǎo)頻率范圍),復(fù)電導(dǎo)率虛部隨著粒徑增大而增大,其原因在于弛豫時間增大引起顆粒表面電導(dǎo)率虛部增大,從而使得多孔介質(zhì)復(fù)電導(dǎo)率虛部增大[58]。

    圖4 有效絕對滲透率隨石英砂顆粒粒徑的變化趨勢

    分析圖4可知,含水合物多孔介質(zhì)的有效絕對滲透率隨粒徑增大而增大,在雙對數(shù)坐標(biāo)平面中有效絕對滲透率與石英砂顆粒粒徑的平方值呈現(xiàn)近似線性關(guān)系(2= 0.996)。當(dāng)孔隙度保持不變時,顆粒粒徑增大的同時顆粒間孔隙尺寸增大,多孔介質(zhì)樣品尺寸不變的條件下則孔隙數(shù)量減小。孔隙直徑增大引起過流斷面面積的增大、孔隙數(shù)量減少導(dǎo)致繞流路徑縮短,兩者均使得滲透率增大[59-62]。

    3.2 水合物飽和度對復(fù)電導(dǎo)率/滲透率影響

    圖5和圖6分別展示了粒徑為48 μm時不同飽和度條件下多孔介質(zhì)的復(fù)電導(dǎo)率虛部頻散曲線和有效絕對滲透率變化曲線。

    圖5 水合物飽和度不同時復(fù)電導(dǎo)率虛部的頻散曲線

    圖6 水合物飽和度不同時有效絕對滲透率曲線

    分析圖5可知,隨著測試頻率的升高,復(fù)電導(dǎo)率虛部隨著水合物飽和度的增加呈現(xiàn)先減小后增大的趨勢。在雙電層極化作用主導(dǎo)的低頻范圍,顆粒表面電導(dǎo)率虛部隨水合物飽和度增加而減小,而顆粒表面電導(dǎo)率近似等于多孔介質(zhì)復(fù)電導(dǎo)率虛部[參見式(5)];在界面極化作用主導(dǎo)的高頻范圍,水合物飽和度增加致使界面極化作用增強,導(dǎo)致多孔介質(zhì)復(fù)電導(dǎo)率虛部增大。

    分析圖6可知,隨著水合物飽和度的升高,多孔介質(zhì)的有效絕對滲透率呈現(xiàn)下降趨勢。水合物飽和度的增加使得孔隙空間減小,導(dǎo)致流體可流通孔隙的等效半徑減小[63-64],從而引起滲透率的下降。當(dāng)水合物飽和度較低時(如圖中20%以下),多孔介質(zhì)有效絕對滲透率下降速率相對較高,這與公開發(fā)表的研究結(jié)果相一致[65-67]。對于水合物以孔隙填充微觀模式賦存于多孔介質(zhì)中的情況,水合物的形成導(dǎo)致流動通道迂曲度的增加。由于水合物更容易在大孔隙中形成且大孔隙對滲透率的貢獻度相對更高,因此在水合物生成初期(水合物飽和度較低時),水合物飽和度的增加所引起的有效絕對滲透率下降速率相對更高。隨著水合物飽和度的逐步升高,水合物所占據(jù)小孔隙的數(shù)量逐漸增多,但是小孔隙對滲透率影響相對較小,因此有效絕對滲透率的下降速率有所降低。

    3.3 基于復(fù)電導(dǎo)率的滲透率評價模型

    以上探討了水合物飽和度、骨架顆粒粒徑對多孔介質(zhì)復(fù)電導(dǎo)率和滲透率的影響規(guī)律。在此基礎(chǔ)上,首先建立水飽和條件下基于復(fù)電導(dǎo)率參數(shù)的絕對滲透率評價模型,然后引入水合物飽和度的影響,從而獲得含水合物多孔介質(zhì)的有效絕對滲透率評價模型。

    圖7展示了水飽和(即無水合物)條件下多孔介質(zhì)的絕對滲透率與復(fù)電導(dǎo)率虛部之間的關(guān)系。

    圖7 多孔介質(zhì)絕對滲透率與復(fù)電導(dǎo)率虛部的關(guān)系

    分析圖7可知,采用冪函數(shù)對絕對滲透率與復(fù)電導(dǎo)率虛部之間的關(guān)系進行擬合可取得良好的效果,擬合模型如式(34)所示(2= 0.993)。

    圖8 基于復(fù)電導(dǎo)率虛部的絕對滲透率模型計算結(jié)果

    Fig. 8 Calculation results from the absolute permeability model based on the imaginary part of complex conductivity

    DAI等[63]提出了以水合物飽和度為變量的歸一化滲透率計算模型,如式(35)所示。同時,非水飽和條件下的復(fù)電導(dǎo)率虛部可表示為水飽和條件下復(fù)電導(dǎo)率虛部和水合物飽和度的函數(shù)[27,68],如式(36)所示。結(jié)合式(34)~ 式(36)可得式(37)所示的有效滲透率計算模型。為了提高式(37)中有效滲透率的計算準(zhǔn)確度,對常數(shù)1進行重新擬合,在式中表示為2。

    圖9所示為多孔介質(zhì)骨架顆粒粒徑為24 μm時,通過式(37)計算得到的有效絕對滲透率與其參考值之間的對比圖。其中,常數(shù)和2的取值分別為3.0和305。圖中可見,在水合物飽和度低于40%范圍內(nèi),有效絕對滲透率的計算誤差均在±15%以內(nèi)。

    圖9 基于復(fù)電導(dǎo)率虛部的有效絕對滲透率模型計算結(jié)果

    4 討 論

    4.1 黏土礦物對電場/流場的影響

    我國海域天然氣水合物主要賦存于未固結(jié)成巖的粉砂、黏土質(zhì)粉砂、粉砂質(zhì)黏土等沉積物中。其中南海神狐海域水合物儲層中多為富含有孔蟲等古生物化石的沉積物,沉積物主要成分的粒徑低于63 μm[4,14]。黏土具有比表面積大、陽離子交換能力強的特點,從而對復(fù)電導(dǎo)率的實部和虛部影響特別顯著,同時部分黏土礦物(如蒙脫石)遇水后發(fā)生體積膨脹,也會顯著降低沉積物的滲透率。

    水合物分解產(chǎn)生純水,從而為黏土礦物的吸水和體積膨脹提供了物質(zhì)條件,推測可知地層中水合物分解過程中存在因水合物飽和度降低引起的滲透率升高、黏土礦物體積膨脹引起的滲透率降低等相互對立的復(fù)雜情況。同時,黏土礦物的類型(如蒙脫石、伊利石、高嶺石等)、黏土含量、黏土的微觀分布特征(如分散狀、層狀、包裹狀、層狀等)均對復(fù)電導(dǎo)率的實部和虛部產(chǎn)生顯著影響[69-70]。綜上所述,黏土礦物對電場和流場參數(shù)的影響效應(yīng)顯著、影響規(guī)律復(fù)雜,需要深入探討存在黏土礦物條件下的含水合物沉積物的滲透率與復(fù)電導(dǎo)率之間的定量關(guān)系。

    4.2 水合物賦存模式對電場/流場的影響

    儲層中的天然氣水合物以固體相態(tài)存在,呈現(xiàn)出多種賦存模式。水合物賦存模式可概括為顆粒排擠型、孔隙浸潤型兩大類[4]。其中,顆粒排擠型主要包含裂隙填充、脈狀、層狀等,孔隙浸潤型可分為孔隙填充、接觸膠結(jié)、顆粒包裹、骨架支撐、局部聚集等形式[71]。由于天然氣水合物對溫度和壓力等環(huán)境條件非常敏感,水合物分解過程中易出現(xiàn)二次生成的現(xiàn)象,在水合物分解/二次生成等過程中水合物的賦存模式可能發(fā)生動態(tài)轉(zhuǎn)變。

    有研究表明,在水合物飽和度相同的條件下,水合物賦存模式的差異也會導(dǎo)致沉積物的孔隙特征尺寸及其統(tǒng)計分布發(fā)生變化,進而引起沉積物復(fù)電導(dǎo)率、滲透率的變化。據(jù)此可以推測,水合物賦存模式必然影響復(fù)電導(dǎo)率與滲透率之間的關(guān)系,因此需要深入研究水合物賦存模式對復(fù)電導(dǎo)率虛部、滲透率的影響規(guī)律與機理。

    4.3 滲透率相關(guān)激電響應(yīng)特征參數(shù)的提取

    多孔介質(zhì)的孔隙度和孔隙表面積(與孔隙尺寸負(fù)相關(guān))是控制多孔介質(zhì)滲透率的關(guān)鍵參數(shù)。有研究人員將直流電導(dǎo)率與多孔介質(zhì)滲透率相關(guān)聯(lián)[18]。沉積物的直流電導(dǎo)率來源于兩類導(dǎo)電機制,即孔隙流體導(dǎo)電、骨架顆粒與流體界面處的雙電層導(dǎo)電(即表面電導(dǎo))??紫读黧w導(dǎo)電機制受到孔隙流體化學(xué)性質(zhì)以及連通孔隙度的影響,而表面電導(dǎo)與沉積物內(nèi)比表面積有關(guān)。對于表面電導(dǎo)足夠小并可以忽略的情況(如孔隙水電導(dǎo)率較高),直流電導(dǎo)率與滲透率之間存在正相關(guān)關(guān)系;然而在表面電導(dǎo)不可忽略的條件下(如較低的孔隙水電導(dǎo)率、含黏土的沉積物等),則表面電導(dǎo)主導(dǎo)多孔介質(zhì)直流電導(dǎo)率的變化,此時滲透率隨著直流電導(dǎo)率的上升而呈現(xiàn)降低的趨勢。由以上分析可見,單獨利用多孔介質(zhì)的直流電導(dǎo)率來預(yù)測滲透率存在顯著的局限性。

    從激電效應(yīng)中可以獲得復(fù)電導(dǎo)率虛部、弛豫時間等參數(shù),如1.2節(jié)中所述,研究人員提出了多個基于激電響應(yīng)參數(shù)的滲透率預(yù)測模型,但是這些模型的推廣應(yīng)用遇到困難。在實際應(yīng)用中,利用激電法測試得到的數(shù)據(jù)往往受到電磁耦合效應(yīng)的干擾,如何從激電/電磁耦合響應(yīng)中分離出真實的激電響應(yīng)數(shù)據(jù)是一個具有挑戰(zhàn)性的課題。此外,多孔介質(zhì)的激電響應(yīng)受到若干物理性、化學(xué)性、結(jié)構(gòu)性因素以及測試所采用電場頻率的影響,因此激電響應(yīng)本身具有極高的復(fù)雜度。由上可見,即使能夠獲取真實的激電響應(yīng)數(shù)據(jù),針對含有多種賦存模式水合物的含黏土沉積物,采用何種激電響應(yīng)參數(shù)來建立其與滲透率之間的關(guān)系模型以及如何提取出這些參數(shù)仍然是值得進一步深入研究的問題。

    5 結(jié) 論

    提出基于激電法定量評價含水合物沉積物滲透率的新思路,構(gòu)建了含水合物多孔介質(zhì)的流場與電場有限元數(shù)值模型,探討了水合物飽和度、骨架顆粒尺寸等對復(fù)電導(dǎo)率和滲透率的影響規(guī)律及機理,建立了基于復(fù)電導(dǎo)率參數(shù)的含水合物多孔介質(zhì)滲透率評價模型,提出了值得深入研究的問題。得到以下結(jié)論:

    (1)建立含水合物沉積物滲透率評價模型的激電響應(yīng)參數(shù)可分為兩類,即描述激電效應(yīng)強度的復(fù)電導(dǎo)率虛部或極化率、描述激電效應(yīng)演化速度的復(fù)電導(dǎo)率譜弛豫時間;

    (2)多孔介質(zhì)骨架顆粒粒徑增大、等效孔徑增大,引起復(fù)電導(dǎo)率譜弛豫時間的延長、復(fù)電導(dǎo)率虛部極大值所對應(yīng)頻率的降低、雙電層極化主導(dǎo)頻率范圍內(nèi)復(fù)電導(dǎo)率虛部的增大以及多孔介質(zhì)滲透率的增大;

    (3)隨著孔隙填充型水合物的飽和度增加,雙電層極化主導(dǎo)的頻率范圍內(nèi)復(fù)電導(dǎo)率虛部減小、界面極化主導(dǎo)的頻率范圍內(nèi)復(fù)電導(dǎo)率虛部增大,多孔介質(zhì)滲透率及其變化率持續(xù)減小;

    (4)雙電層極化主導(dǎo)的頻率范圍內(nèi)復(fù)電導(dǎo)率虛部與滲透率之間存在確定性關(guān)系;在飽和水條件下選取的固定頻率處的復(fù)電導(dǎo)率虛部與滲透率呈冪函數(shù)關(guān)系,結(jié)合飽和與非飽和條件下復(fù)電導(dǎo)率之間的關(guān)系以及歸一化滲透率與飽和度之間的關(guān)系可建立基于復(fù)電導(dǎo)率虛部的含水合物多孔介質(zhì)有效滲透率評價模型。

    為了開發(fā)基于激電法的天然氣水合物儲層滲透率動態(tài)監(jiān)測新技術(shù),需要深入研究黏土礦物類型及含量、水合物微觀賦存模式、骨架顆粒/孔隙尺寸分布、游離氣飽和度及微觀分布、孔隙水及其與骨架顆粒界面化學(xué)性質(zhì)等因素對沉積物低頻交變電場響應(yīng)和滲流特性的影響規(guī)律及機理,設(shè)計含水合物沉積物復(fù)電導(dǎo)率譜與滲透率聯(lián)合測試實驗裝置,在實驗的基礎(chǔ)上進一步優(yōu)化激電響應(yīng)數(shù)據(jù)處理與特征參數(shù)提取方法以及滲透率評價模型。

    [1] LI X S, XU C G, ZHANG Y, et al. Investigation into gas production from natural gas hydrate: a review[J]. Applied energy, 2016, 172: 286-322. DOI: 10.1016/j.apenergy. 2016.03.101.

    [2] 李清平, 周守為, 趙佳飛, 等. 天然氣水合物開采技術(shù)研究現(xiàn)狀與展望[J]. 中國工程科學(xué), 2022, 24(3): 214-224. DOI: 10.15302/J-SSCAE-2022.03.022.

    [3] 孫金聲, 程遠(yuǎn)方, 秦緒文, 等. 南海天然氣水合物鉆采機理與調(diào)控研究進展[J]. 中國科學(xué)基金, 2021, 35(6): 940-951. DOI: 10.16262/j.cnki.1000-8217.2021.06.014.

    [4] 寧伏龍, 梁金強, 吳能友, 等. 中國天然氣水合物賦存特征[J]. 天然氣工業(yè), 2020, 40(8): 1-24. DOI: 10.3787/ j.issn.1000-0976.2020.08.001.

    [5] REN X W, GUO Z Y, NING F L, et al. Permeability of hydrate-bearing sediments[J]. Earth-science reviews, 2020, 202: 103100. DOI: 10.1016/j.earscirev.2020.103100.

    [6] 劉樂樂, 萬義釗, 李承峰, 等. 天然氣水合物儲層有效絕對滲透率現(xiàn)場測試進展[J]. 海洋地質(zhì)前沿, 2022, 38(11): 40-55. DOI: 10.16028/j.1009-2722.2022.232.

    [7] 曾家明, 李棟梁, 梁德青, 等. 天然氣水合物儲層滲透率研究進展[J]. 新能源進展, 2021, 9(1): 25-34. DOI: 10.3969/j.issn.2095-560X.2021.01.004.

    [8] LI G, XU Z L, LI X S, et al. Permeability investigation and hydrate migration of hydrate–bearing silty sands and silt[J]. Journal of natural gas science and engineering, 2021, 89: 103891. DOI: 10.1016/j.jngse.2021.103891.

    [9] CHOI J H, MYSHAKIN E M, LEI L, et al. An experimental system and procedure of unsteady-state relative permeability test for gas hydrate-bearing sediments[J]. Journal of natural gas science and engineering, 2020, 83: 103545. DOI: 10.1016/j.jngse.2020.103545.

    [10] MAHABADI N, DAI S, SEOL Y, et al. Impact of hydrate saturation on water permeability in hydrate-bearing sediments[J]. Journal of petroleum science and engineering,2019, 174: 696-703. DOI: 10.1016/J.PETROL.2018.11.084.

    [11] PAN L L, LEI L, SEOL Y. Pore-scale influence of methane hydrate on permeability of porous media[J]. Journal of natural gas science and engineering, 2021, 87: 103758. DOI: 10.1016/j.jngse.2020.103758.

    [12] WU Z R, LIU W G, ZHENG J N, et al. Effect of methane hydrate dissociation and reformation on the permeability of clayey sediments[J]. Applied energy, 2020, 261: 114479. DOI: 10.1016/j.apenergy.2019.114479.

    [13] ZHANG Y C, LI C F, MA J S, et al. Investigating the effective permeability evolution as a function of hydrate saturation in the hydrate-bearing sands using a kinetic-theory-based pore network model[J]. Computers and geotechnics, 2022, 150: 104930. DOI: 10.1016/j.compgeo. 2022.104930.

    [14] 秦緒文, 陸程, 王平康, 等. 中國南海天然氣水合物開采儲層水合物相變與滲流機理: 綜述與展望[J]. 中國地質(zhì), 2022, 49(3): 749-769. DOI: 10.12029/gc20220306.

    [15] 何繼善. 頻率域電法的新進展[J]. 地球物理學(xué)進展, 2007, 22(4): 1250-1254. DOI: 10.3969/j.issn.1004-2903. 2007.04.035.

    [16] STINGACIU L R, WEIHERMüLLER L, HABER-POHLMEIER S, et al. Determination of pore size distribution and hydraulic properties using nuclear magnetic resonance relaxometry: a comparative study of laboratory methods[J]. Water resources research, 2010, 46(11): W11510. DOI: 10.1029/2009WR008686.

    [17] HUBBARD S S, LINDE N. Hydrogeophysics[EB/OL]. 2011. [2023-03-27]. http://escholarship.org/uc/item/11c8s8d4.

    [18] SLATER L. Near surface electrical characterization of hydraulic conductivity: from petrophysical properties to aquifer geometries—a review[J]. Surveys in geophysics, 2007, 28(2/3): 169-197. DOI: 10.1007/s10712-007-9022-y.

    [19] TZIRITA A. A study of electrical and thermal properties and their use to detect natural gas hydrates in ocean sediments[D]. College Station: Texas A&M University. 1992.

    [20] DU FRANE W L, STERN L A, WEITEMEYER K A, et al.Electrical properties of polycrystalline methane hydrate[J]. Geophysical research letters, 2011, 38(9): L09313. DOI: 10.1029/2011GL047243.

    [21] DU FRANE W L, STERN L A, CONSTABLE S, et al. Electrical properties of methane hydrate + sediment mixtures[J]. Journal of geophysical research: solid earth, 2015, 120(7): 4773-4783. DOI: 10.1002/2015jb011940.

    [22] 邢蘭昌, 陳強, 劉昌嶺. 基于電化學(xué)阻抗譜測試方法研究四氫呋喃水合物的生成和分解過程[J]. 巖礦測試, 2015, 34(6): 704-711. DOI: 10.15898/j.cnki.11-2131/td. 2015.06.016.

    [23] 李棟梁, 盧靜生, 梁德青. 祁連山凍土區(qū)天然氣水合物形成對巖芯電阻率及介電常數(shù)的影響[J]. 新能源進展, 2016, 4(3): 179-183. DOI: 10.3969/j.issn.2095-560X. 2016.03.003.

    [24] 邢蘭昌, 祁雨, 朱泰, 等. 含甲烷水合物沉積物電–聲響應(yīng)特性聯(lián)合探測: 裝置開發(fā)與實驗研究[J]. 新能源進展, 2018, 6(2): 119-129. DOI: 10.3969/j.issn.2095-560X.2018.02.006.

    [25] XING L C, ZHU T, NIU J L, et al. Development and validation of an acoustic-electrical joint testing system for hydrate-bearing porous media[J]. Advances in mechanical engineering, 2020, 12(3): 1-11. DOI: 10.1177/ 1687814020908981.

    [26] 邢蘭昌, 牛佳樂, 魏偉, 等. 基于寬頻復(fù)電阻率的含黏土沉積物中水合物飽和度計算方法[J]. 新能源進展, 2020, 8(4): 251-257. DOI: 10.3969/j.issn.2095-560X. 2020.04.001.

    [27] XING L C, QI S Y, XU Y, et al. Numerical study on complex conductivity characteristics of hydrate-bearing porous media[J]. Journal of natural gas science and engineering, 2021, 95: 104145. DOI: 10.1016/j.jngse. 2021.104145.

    [28] XING L C, NIU J L, ZHANG S L, et al. Experimental study on hydrate saturation evaluation based on complex electrical conductivity of porous media[J]. Journal of petroleum science and engineering, 2022, 208: 109539. DOI: 10.1016/j.petrol.2021.109539.

    [29] 徐源, 張歡歡, 邢蘭昌, 等. 基于電-力-聲多物理場耦合數(shù)值模型的含水合物多孔介質(zhì)聲速和衰減特性研究[J].新能源進展, 2022, 10(5): 400-409. DOI: 10.3969/j.issn.2095-560X.2022.05.002.

    [30] B?RNER F D, SCHOPPER J R, WELLER A. Evaluation of transport and storage properties in the soil and groundwater zone from induced polarization measurements[J]. Geophysical prospecting, 1996, 44(4): 583-601. DOI: 10.1111/j.1365-2478.1996.tb00167.x.

    [31] WELLER A, SLATER L, BINLEY A, et al. Permeability prediction based on induced polarization: insights from measurements on sandstone and unconsolidated samples spanning a wide permeability range[J]. Geophysics, 2015, 80(2): D161-D173. DOI: 10.1190/geo2014-0368.1

    [32] VINEGAR H J, WAXMAN M H. Induced polarization of shaly sands[J]. Geophysics, 1984, 49(8): 1267-1287. DOI: 10.1190/1.1441755.

    [33] 蔡建超, 夏宇軒, 徐賽, 等. 含水合物沉積物多相滲流特性研究進展[J]. 力學(xué)學(xué)報, 2020, 52(1): 208-223. DOI: 10.6052/0459-1879-19-362.

    [34] ZHANG Z Y, WELLER A. Fractal dimension of pore-space geometry of an Eocene sandstone formation[J]. Geophysics, 2014, 79(6): D377-D387. DOI: 10.1190/ geo2014-0143.1.

    [35] KOSTEK S, SCHWARTZ L M, JOHNSON D L. Fluid permeability in porous media: comparison of electrical estimates with hydrodynamical calculations[J]. Physical review B, 1992, 45(1): 186-195. DOI: 10.1103/ PhysRevB.45.186.

    [36] THOMPSON A H, KATZ A J, KROHN C E. The microgeometry and transport properties of sedimentary rock[J]. Advances in physics, 1987, 36(5): 625-694. DOI: 10.1080/00018738700101062.

    [37] SCHEIDEGGER A E. The physics of flow through porous media[M]. 3rd ed. Toronto: University of Toronto Press, 1974. DOI: 10.3138/j.ctvfrxmtw.

    [38] CARMAN P C. Fluid flow through granular beds[J]. Chemical engineering research and design, 1997, 75(S1): S32-S48. DOI: 10.1016/S0263-8762(97)80003-2.

    [39] GUéGUEN Y, PALCIAUSKAS V. Introduction to the Physics of Rocks[M]. Princeton: Princeton University Press, 1994.

    [40] PAPE H, RIEPE L, SCHOPPER J R. A pigeon-hole model for relating permeability to specific surface[J]. The log analyst, 1982, 23(1): 5-13.

    [41] JOHNSON D L, KOPLIK J, SCHWARTZ L M. New pore-size parameter characterizing transport in porous media[J]. Physical review letters, 1986, 57(20): 2564-2567. DOI: 10.1103/PhysRevLett.57.2564.

    [42] AVELLANEDA M, TORQUATO S. Rigorous link between fluid permeability, electrical conductivity, and relaxation times for transport in porous media[J]. Physics of fluids A: fluid dynamics, 1991, 3(11): 2529-2540. DOI: 10.1063/1.858194.

    [43] REVIL A, CATHLES III L M. Permeability of shaly sands[J]. Water resources research, 1999, 35(3): 651-662. DOI: 10.1029/98WR02700.

    [44] FIANDACA G, MAURYA P K, BALBARINI N, et al. Permeability estimation directly from logging-while-drilling induced polarization data[J]. Water resources research, 2018, 54(4): 2851-2870. DOI: 10.1002/2017WR022411.

    [45] SLATER L, LESMES D P. Electrical-hydraulic relationships observed for unconsolidated sediments[J]. Water resources research, 2002, 38(10): 1213. DOI: 10.1029/2001WR001075.

    [46] SLATER L D, GLASER D R. Controls on induced polarization in sandy unconsolidated sediments and application to aquifer characterization[J]. Geophysics, 2003, 68(5): 1547-1558. DOI: 10.1190/1.1620628.

    [47] REVIL A. Spectral induced polarization of shaly sands: influence of the electrical double layer[J]. Water resources research, 2012, 48(2): W02517. DOI: 10.1029/ 2011WR011260.

    [48] SCOTT J B T, BARKER R D. Determining pore-throat size in Permo-Triassic sandstones from low-frequency electrical spectroscopy[J]. Geophysical research letters, 2003, 30(9): 1450. DOI: 10.1029/2003GL016951.

    [49] BINLEY A, SLATER L D, FUKES M, et al. Relationship between spectral induced polarization and hydraulic properties of saturated and unsaturated sandstone[J]. Water resources research, 2005, 41(12): W12417. DOI: 10.1029/2005WR004202.

    [50] KEMNA A, MüNCH H M, TITOV K, et al. Relation of SIP relaxation time of sands to salinity, grain size and hydraulic conductivity[C]//Near Surface 2005-11th European Meeting of Environmental and Engineering Geophysics. Houten: European Association of Geoscientists & Engineers, 2005: cp-13-00047. DOI: 10.3997/2214-4609-pdb.13.P054.

    [51] ZISSER N, KEMNA A, NOVER G. Relationship between low-frequency electrical properties and hydraulic permeability of low-permeability sandstones[J]. Geophysics, 2010, 75(3): E131-E141. DOI: 10.1190/1.3413260.

    [52] REVIL A, FLORSCH N. Determination of permeability from spectral induced polarization in granular media[J]. Geophysical journal international, 2010, 181(3): 1480-1498. DOI: 10.1111/j.1365-246X.2010.04573.x.

    [53] REVIL A, KOCH K, HOLLIGER K. Is it the grain size or the characteristic pore size that controls the induced polarization relaxation time of clean sands and sandstones?[J]. Water resources research, 2012, 48(5): W05602. DOI: 10.1029/2011WR011561.

    [54] WAITE W F, SANTAMARINA J C, CORTES D D, et al. Physical properties of hydrate-bearing sediments[J]. Reviews of geophysics, 2009, 47(4): RG4003. DOI: 10.1029/2008RG000279.

    [55] HAUKALID K, FOLGER? K, BARTH T, et al. Hydrate formation in water-in-crude oil emulsions studied by broad-band permittivity measurements[J]. Energy & fuels, 2017, 31(4): 3793-3803. DOI: 10.1021/acs.energyfuels. 6b03416.

    [56] TITOV K, KOMAROV V, TARASOV V, et al. Theoretical and experimental study of time domain-induced polarization in water-saturated sands[J]. Journal of applied geophysics, 2002, 50(4): 417-433. DOI: 10.1016/S0926-9851(02)00168-4.

    [57] LEROY P, REVIL A, KEMNA A, et al. Complex conductivity of water-saturated packs of glass beads[J]. Journal of colloid and interface science, 2008, 321(1): 103-117. DOI: 10.1016/j.jcis.2007.12.031.

    [58] REVIL A. Effective conductivity and permittivity of unsaturated porous materials in the frequency range 1 mHz–1GHz[J]. Water resources research, 2013, 49(1): 306-327. DOI: 10.1029/2012WR012700.

    [59] SHEPHERD R G. Correlations of permeability and grain size[J]. Groundwater, 1989, 27(5): 633-638. DOI: 10.1111/j.1745-6584.1989.tb00476.x.

    [60] GHASSEMI A, PAK A. Pore scale study of permeability and tortuosity for flow through particulate media using Lattice Boltzmann method[J]. International journal for numerical and analytical methods in geomechanics, 2011, 35(8): 886-901. DOI: 10.1002/nag.932.

    [61] 喬太斌, 楊玉雙, 李如如, 等. 多孔介質(zhì)中固體體積分?jǐn)?shù)與顆粒尺度對流體絕對滲透率的影響[J]. 山西大學(xué)學(xué)報(自然科學(xué)版), 2017, 40(1): 92-99. DOI: 10.13451/ j.cnki.shanxi.univ(nat.sci.).2017.01.013.

    [62] 楊斌, 徐曾和, 楊天鴻, 等. 高水力梯度條件下顆粒堆積型多孔介質(zhì)滲流規(guī)律試驗研究[J]. 巖土力學(xué), 2018, 39(11): 4017-4024. DOI: 10.16285/j.rsm.2017.0643.

    [63] DAI S, SEOL Y. Water permeability in hydrate-bearing sediments: a pore-scale study[J]. Geophysical research letters, 2014, 41(12): 4176-4184. DOI: 10.1002/2014GL060535.

    [64] 李世龍, 李剛, 魏納, 等. 含甲烷水合物的石英砂滲透率實驗和分形模型對比研究[J]. 新能源進展, 2020, 8(4): 264-271. DOI: 10.3969/j.issn.2095-560X.2020.04.003.

    [65] CHEN Y, SUN B J, CHEN L T, et al. Simulation and observation of hydrate phase transition in porous medium via microfluidic application[J]. Industrial & engineering chemistry research, 2019, 58(12): 5071-5079. DOI: 10.1021/acs.iecr.9b00168.

    [66] LUO Y, SUN Y, LI L, et al. Image-based pore-network modeling of two-phase flow in hydrate-bearing porous media[J]. Energy, 2022, 252: 124044. DOI: 10.1016/j. energy.2022.124044.

    [67] XU J C, BU Z W, LI H Y, et al. Pore-scale flow simulation on the permeability in hydrate-bearing sediments[J]. Fuel, 2022, 312: 122681. DOI: 10.1016/j.fuel.2021.122681.

    [68] DENG Y P, SHI X Q, REVIL A, et al. Complex conductivity of oil-contaminated clayey soils[J]. Journal of hydrology, 2018, 561: 930-942. DOI: 10.1016/j. jhydrol.2018.04.055.

    [69] XING L C, ZHANG H H, WANG S, et al. Pore-scale modelling on complex-conductivity responses of hydrate-bearing clayey sediments: implications for evaluating hydrate saturation and clay content[J]. Geoenergy science and engineering, 2023, 221: 211356. DOI: 10.1016/j. geoen.2022.211356.

    [70] XING L C, ZHANG S L, ZHANG H H, et al. Saturation estimation with complex electrical conductivity for hydrate-bearing clayey sediments: an experimental study[J]. Journal of ocean university of China, 2023. DOI: 10.1007/s11802-023-5492-x.

    [71] DAI S, SANTAMARINA J C, WAITE W F, et al. Hydrate morphology: physical properties of sands with patchy hydrate saturation[J]. Journal of geophysical research: solid earth, 2012, 117(B11): B11205. DOI: 10.1029/ 2012jb009667.

    Numerical Study on Permeability Evaluation for Hydrate-Bearing Sediments Based on Induced Polarization Principle

    XING Lanchang1,?, WANG Shuo1, ZHANG Huanhuan1, WEI Wei2,HAN Weifeng2, YANG Jinxiu3, GE Xinmin3

    (1. College of Control Science and Engineering, China University of Petroleum (East China), Qingdao 266580, Shandong, China;2. Department of Alternative Energy, PetroChina Research Institute of Petroleum Exploration & Development, Langfang 065007, Hebei, China; 3. School of Geosciences, China University of Petroleum (East China), Qingdao 266580, Shandong, China)

    There are significant problems such as high cost, long period and phase-state change of hydrate in the experimental process for measuring the permeability of hydrate-bearing samples. Moreover, real-time monitoring technologies for the dynamic evolution of reservoir permeability in hydrate exploitation processes are still lacking. A new method for evaluating the permeability of hydrate-bearing sediments based on induced polarization (IP) responses was proposed, and finite-element numerical models of flow field and electrical field for hydrate-bearing porous media were constructed. The effects and corresponding mechanisms of hydrate saturation and skeleton particle size on the complex conductivity and permeability were analyzed. Finally, a permeability evaluation model of hydrate-bearing porous media based on complex conductivity was established. It was demonstrated that: (1) with the increase of the skeleton particle size and the equivalent pore diameter of porous medium, the relaxation time of the complex-conductivity spectrum increases, the frequency corresponding to the maximum imaginary complex conductivity decreases, the imaginary complex conductivity within the frequency range of electrical-double-layer (EDL) polarization increases, and the effective permeability of the porous medium increases; (2) with the increase of the pore-filling-hydrate saturation, the imaginary complex conductivity within the frequency range of EDL polarization decreases, the imaginary complex conductivity within the frequency range of interface polarization increases, and the effective permeability and its changing rate of the hydrate-bearing porous medium decrease continuously; (3) there is a definite relationship between the imaginary complex conductivity and permeability within the frequency range of EDL polarization; the imaginary complex conductivity at a fixed frequency under water-saturated conditions exhibits a power-law relationship with permeability; combining the relationship between the complex conductivity under saturated/unsaturated conditions and that between the normalized permeability and saturation, an evaluation model for the effective permeability of hydrate-bearing porous media can be established based on the imaginary complex conductivity. This study can provide a theoretical and model basis for the development of an experimental measurement technology for the permeability of hydrate-bearing samples based on IP principle and a dynamic monitoring technology for the permeability of natural gas hydrate reservoirs.

    natural gas hydrate; permeability; induced polarization; complex conductivity; hydrate saturation; numerical simulation

    2095-560X(2023)04-0320-13

    TK01;P631

    A

    10.3969/j.issn.2095-560X.2023.04.004

    2022-12-21

    2023-03-27

    中石油“十四五”前瞻性基礎(chǔ)性重大科技項目(2021DJ4901);國家留學(xué)基金項目(202106455003);中央高?;究蒲袠I(yè)務(wù)費專項資金項目(20CX05005A);中國石油科技創(chuàng)新基金項目(2018D-5007-0214);山東省自然科學(xué)基金項目(ZR2019MEE095)

    邢蘭昌,E-mail:xinglc@upc.edu.cn

    邢蘭昌, 王碩, 張歡歡, 等. 基于激電法評價含水合物沉積物滲透率數(shù)值模擬研究[J]. 新能源進展, 2023, 11(4): 320-332.

    : XING Lanchang, WANG Shuo, ZHANG Huanhuan, et al. Numerical study on permeability evaluation for hydrate-bearing sediments based on induced polarization principle[J]. Advances in new and renewable energy, 2023, 11(4): 320-332.

    邢蘭昌(1983-),男,博士,副教授,碩士生導(dǎo)師,主要從事天然氣水合物與多相流相關(guān)檢測理論與方法、多物理場耦合數(shù)值模擬方法、智能感知與檢測技術(shù)、計算機測控系統(tǒng)研究。

    王 碩(1997-),女,碩士研究生,主要從事天然氣水合物與多相流相關(guān)檢測及數(shù)值模擬方法的研究。

    張歡歡(1997-),女,碩士研究生,主要從事天然氣水合物與多相流相關(guān)檢測及數(shù)值模擬方法的研究。

    猜你喜歡
    虛部激電水合物
    格點量子色動力學(xué)數(shù)據(jù)的虛部分布與信號改進*
    兩類特殊多項式的復(fù)根虛部估計
    大功率激電測深方法在豫西董家埝銀礦床勘查中的應(yīng)用
    氣井用水合物自生熱解堵劑解堵效果數(shù)值模擬
    高頻大地電磁測深與激電中梯在金礦勘查中的應(yīng)用研究
    大功率激電測量在冀北溫家營—馬家溝銀多金屬礦勘查中的應(yīng)用
    激電聯(lián)合剖面在判斷矽卡巖型礦床礦體產(chǎn)狀中的應(yīng)用
    例談復(fù)數(shù)應(yīng)用中的計算兩次方法
    熱水吞吐開采水合物藏數(shù)值模擬研究
    天然氣水合物保壓轉(zhuǎn)移的壓力特性
    av女优亚洲男人天堂| 国产片特级美女逼逼视频| 极品教师在线视频| 久久久久久久久久成人| 日韩一本色道免费dvd| 精品国产一区二区久久| av一本久久久久| 丰满少妇做爰视频| 好男人视频免费观看在线| 日韩一区二区三区影片| 国产成人午夜福利电影在线观看| 久久久午夜欧美精品| 亚州av有码| 精品卡一卡二卡四卡免费| 亚洲精品456在线播放app| 国产精品人妻久久久久久| 三级国产精品欧美在线观看| 国产亚洲91精品色在线| 男人和女人高潮做爰伦理| 久久久久久久久久久丰满| 精品少妇黑人巨大在线播放| 观看av在线不卡| 亚洲情色 制服丝袜| 国产黄色视频一区二区在线观看| 少妇的逼水好多| 久久99一区二区三区| 色吧在线观看| 日韩免费高清中文字幕av| 亚洲不卡免费看| 国产在线视频一区二区| 久久国产乱子免费精品| 青青草视频在线视频观看| 久久精品夜色国产| 中文资源天堂在线| 国产一区二区三区av在线| 日韩强制内射视频| 热re99久久精品国产66热6| 午夜免费男女啪啪视频观看| 精品国产一区二区三区久久久樱花| 少妇丰满av| 日本av手机在线免费观看| 国产精品国产三级国产专区5o| 午夜免费观看性视频| 三级国产精品欧美在线观看| 日韩强制内射视频| 欧美 日韩 精品 国产| 高清视频免费观看一区二区| 精品久久久久久久久亚洲| 人人澡人人妻人| 国产视频内射| 国产成人freesex在线| 欧美精品一区二区免费开放| 91久久精品国产一区二区成人| 国产精品久久久久久久久免| 国产黄片美女视频| av免费在线看不卡| 能在线免费看毛片的网站| 极品少妇高潮喷水抽搐| 欧美三级亚洲精品| 一边亲一边摸免费视频| 少妇熟女欧美另类| 一个人看视频在线观看www免费| 日产精品乱码卡一卡2卡三| 国产精品三级大全| 建设人人有责人人尽责人人享有的| 久久久精品免费免费高清| 亚洲av欧美aⅴ国产| av又黄又爽大尺度在线免费看| 一本久久精品| 精品人妻一区二区三区麻豆| 国产精品一区二区性色av| 777米奇影视久久| 99久国产av精品国产电影| 国产高清三级在线| a级毛片在线看网站| 国产亚洲午夜精品一区二区久久| 久久99热6这里只有精品| 欧美+日韩+精品| av播播在线观看一区| 亚洲情色 制服丝袜| 男人狂女人下面高潮的视频| 美女脱内裤让男人舔精品视频| 国产成人91sexporn| 一本久久精品| 黄色欧美视频在线观看| 亚洲精品视频女| 久久ye,这里只有精品| 国产日韩欧美视频二区| 五月开心婷婷网| 日韩不卡一区二区三区视频在线| 久久国内精品自在自线图片| 日日摸夜夜添夜夜添av毛片| 欧美精品一区二区免费开放| av在线播放精品| 三上悠亚av全集在线观看 | 777米奇影视久久| av女优亚洲男人天堂| 欧美变态另类bdsm刘玥| 免费观看a级毛片全部| 最近中文字幕高清免费大全6| 亚洲国产欧美日韩在线播放 | 日韩欧美精品免费久久| 黑人巨大精品欧美一区二区蜜桃 | 亚洲国产日韩一区二区| 国内揄拍国产精品人妻在线| 亚洲四区av| 美女中出高潮动态图| 22中文网久久字幕| 亚洲情色 制服丝袜| 在线亚洲精品国产二区图片欧美 | 日韩制服骚丝袜av| 免费观看的影片在线观看| 日韩视频在线欧美| 日本色播在线视频| 一本大道久久a久久精品| 成人美女网站在线观看视频| 边亲边吃奶的免费视频| 天堂8中文在线网| 久久久久久久大尺度免费视频| 成人免费观看视频高清| 亚州av有码| 乱系列少妇在线播放| 国产精品国产av在线观看| 欧美亚洲 丝袜 人妻 在线| 在线观看美女被高潮喷水网站| 极品教师在线视频| 在线看a的网站| 一级片'在线观看视频| 午夜av观看不卡| 亚洲中文av在线| 久久久亚洲精品成人影院| 日本wwww免费看| a级毛片免费高清观看在线播放| 美女福利国产在线| 三级国产精品片| 午夜影院在线不卡| 国产欧美亚洲国产| 日产精品乱码卡一卡2卡三| 高清av免费在线| 插阴视频在线观看视频| 日韩电影二区| 少妇人妻一区二区三区视频| 日韩免费高清中文字幕av| 99久久精品一区二区三区| 日产精品乱码卡一卡2卡三| 精品少妇久久久久久888优播| 国产爽快片一区二区三区| 免费观看av网站的网址| 免费大片18禁| 免费黄频网站在线观看国产| 国产极品天堂在线| 丰满少妇做爰视频| 高清视频免费观看一区二区| 在线观看三级黄色| 美女视频免费永久观看网站| 美女脱内裤让男人舔精品视频| 亚洲综合色惰| 好男人视频免费观看在线| 国产免费又黄又爽又色| 交换朋友夫妻互换小说| 久久影院123| 午夜精品国产一区二区电影| 交换朋友夫妻互换小说| 亚洲婷婷狠狠爱综合网| 久久毛片免费看一区二区三区| 久久国产乱子免费精品| 国产精品一区www在线观看| 另类亚洲欧美激情| 草草在线视频免费看| 老司机影院成人| 午夜免费男女啪啪视频观看| 色吧在线观看| 亚洲国产欧美日韩在线播放 | 亚洲精品久久午夜乱码| 国产色婷婷99| 欧美日韩精品成人综合77777| 一级毛片我不卡| av专区在线播放| 国产真实伦视频高清在线观看| 最后的刺客免费高清国语| 免费人成在线观看视频色| 中文欧美无线码| 人妻人人澡人人爽人人| 高清视频免费观看一区二区| 九色成人免费人妻av| 丝袜在线中文字幕| 美女大奶头黄色视频| 嫩草影院新地址| 久久久a久久爽久久v久久| 汤姆久久久久久久影院中文字幕| av不卡在线播放| 亚洲精品一二三| 亚洲欧洲国产日韩| 人妻系列 视频| 三级经典国产精品| 大片电影免费在线观看免费| 欧美成人午夜免费资源| 男的添女的下面高潮视频| 黄色怎么调成土黄色| 乱码一卡2卡4卡精品| 亚洲av国产av综合av卡| 99久久精品国产国产毛片| 日韩在线高清观看一区二区三区| 日韩大片免费观看网站| 亚洲国产欧美日韩在线播放 | 一级av片app| 国产熟女欧美一区二区| 亚洲情色 制服丝袜| 国产亚洲av片在线观看秒播厂| 91午夜精品亚洲一区二区三区| 国产成人a∨麻豆精品| 久久精品熟女亚洲av麻豆精品| 日韩中字成人| 免费黄色在线免费观看| av.在线天堂| 免费大片18禁| 男人添女人高潮全过程视频| 亚洲欧美成人综合另类久久久| 亚洲人与动物交配视频| 成人亚洲欧美一区二区av| 精品一区在线观看国产| 亚洲av在线观看美女高潮| 午夜日本视频在线| 国产免费一区二区三区四区乱码| 国产av一区二区精品久久| 日韩精品有码人妻一区| 色视频www国产| 欧美xxⅹ黑人| 中文字幕av电影在线播放| av专区在线播放| 亚洲怡红院男人天堂| 久久精品夜色国产| 三级国产精品欧美在线观看| 午夜老司机福利剧场| 国产成人一区二区在线| 精品久久国产蜜桃| 国产精品麻豆人妻色哟哟久久| av播播在线观看一区| 男人和女人高潮做爰伦理| 看免费成人av毛片| 国产亚洲午夜精品一区二区久久| 老女人水多毛片| 国产一区二区三区av在线| 成年人午夜在线观看视频| 免费大片黄手机在线观看| 国产精品欧美亚洲77777| 久久精品国产鲁丝片午夜精品| 久久99精品国语久久久| 国产精品嫩草影院av在线观看| 精品人妻熟女毛片av久久网站| xxx大片免费视频| a级毛色黄片| 99re6热这里在线精品视频| 国产精品秋霞免费鲁丝片| 国产日韩一区二区三区精品不卡 | 久久久亚洲精品成人影院| 亚洲熟女精品中文字幕| 亚洲欧美日韩卡通动漫| 观看免费一级毛片| 日本黄大片高清| 少妇人妻久久综合中文| 精品少妇内射三级| 精品99又大又爽又粗少妇毛片| 国产av精品麻豆| 男女国产视频网站| 51国产日韩欧美| 在线播放无遮挡| 国产精品无大码| 久久韩国三级中文字幕| 日本wwww免费看| 久久久精品94久久精品| 狂野欧美激情性bbbbbb| 精品视频人人做人人爽| 天天躁夜夜躁狠狠久久av| 一级毛片电影观看| 国产无遮挡羞羞视频在线观看| 内射极品少妇av片p| 一级毛片电影观看| 超碰97精品在线观看| 大香蕉久久网| 亚洲av.av天堂| 色婷婷久久久亚洲欧美| 国产伦精品一区二区三区四那| 国产极品天堂在线| 亚洲欧美日韩卡通动漫| 国产日韩欧美亚洲二区| 中文天堂在线官网| av在线播放精品| 精品久久久噜噜| 2022亚洲国产成人精品| 久久97久久精品| a级片在线免费高清观看视频| 黑人高潮一二区| 中文字幕免费在线视频6| 婷婷色麻豆天堂久久| 成人无遮挡网站| 在线精品无人区一区二区三| 久久午夜综合久久蜜桃| 国产av一区二区精品久久| 国产精品人妻久久久影院| 成年av动漫网址| 免费高清在线观看视频在线观看| av国产精品久久久久影院| 在线观看人妻少妇| 国产永久视频网站| 偷拍熟女少妇极品色| 午夜日本视频在线| 国产精品久久久久久精品电影小说| 九九在线视频观看精品| 王馨瑶露胸无遮挡在线观看| 性高湖久久久久久久久免费观看| 亚洲精品久久午夜乱码| 久久鲁丝午夜福利片| 久久综合国产亚洲精品| 日韩成人伦理影院| 国产精品女同一区二区软件| 午夜免费观看性视频| 国产一区亚洲一区在线观看| 免费黄频网站在线观看国产| 爱豆传媒免费全集在线观看| 亚洲av成人精品一区久久| 99热网站在线观看| 成年av动漫网址| 69精品国产乱码久久久| 丝袜在线中文字幕| 3wmmmm亚洲av在线观看| 如何舔出高潮| 国产午夜精品久久久久久一区二区三区| 丝袜脚勾引网站| 三级国产精品欧美在线观看| av女优亚洲男人天堂| 制服丝袜香蕉在线| 99久久精品国产国产毛片| 人体艺术视频欧美日本| 精品久久久精品久久久| 欧美日韩综合久久久久久| 啦啦啦视频在线资源免费观看| 亚洲无线观看免费| 成人18禁高潮啪啪吃奶动态图 | 亚洲av在线观看美女高潮| 欧美bdsm另类| 国产深夜福利视频在线观看| www.色视频.com| 乱人伦中国视频| 国产精品秋霞免费鲁丝片| 狠狠精品人妻久久久久久综合| 日本色播在线视频| 国产91av在线免费观看| 国内揄拍国产精品人妻在线| 久久精品国产a三级三级三级| 最黄视频免费看| 国产无遮挡羞羞视频在线观看| 欧美亚洲 丝袜 人妻 在线| 亚洲欧美日韩另类电影网站| 这个男人来自地球电影免费观看 | 人妻一区二区av| 亚洲国产精品一区三区| 十八禁网站网址无遮挡 | 黑丝袜美女国产一区| 国产免费视频播放在线视频| 欧美最新免费一区二区三区| 99九九在线精品视频 | 欧美日韩av久久| 精品熟女少妇av免费看| 日本av手机在线免费观看| 18禁动态无遮挡网站| 久久久国产精品麻豆| 男女国产视频网站| 丁香六月天网| 欧美亚洲 丝袜 人妻 在线| 亚洲va在线va天堂va国产| 国产av精品麻豆| 能在线免费看毛片的网站| 最近中文字幕2019免费版| 大片免费播放器 马上看| 看十八女毛片水多多多| 久久午夜福利片| 欧美最新免费一区二区三区| 在线观看免费高清a一片| 国产成人精品福利久久| 黄色毛片三级朝国网站 | 简卡轻食公司| 国产亚洲5aaaaa淫片| 日韩一区二区三区影片| av免费观看日本| 纯流量卡能插随身wifi吗| 又大又黄又爽视频免费| 久久久久久久久久久丰满| 久久久久久久久大av| 国产乱人偷精品视频| 99久久中文字幕三级久久日本| 一个人免费看片子| 午夜激情久久久久久久| 国产av一区二区精品久久| 国产高清国产精品国产三级| 麻豆成人av视频| 啦啦啦中文免费视频观看日本| 免费看光身美女| 啦啦啦中文免费视频观看日本| 欧美激情极品国产一区二区三区 | 成人18禁高潮啪啪吃奶动态图 | 夜夜看夜夜爽夜夜摸| 国产精品熟女久久久久浪| 边亲边吃奶的免费视频| 国产视频内射| 少妇丰满av| 十八禁高潮呻吟视频 | 97超碰精品成人国产| 国产日韩一区二区三区精品不卡 | 国产精品久久久久久久电影| 丝袜脚勾引网站| 亚洲一区二区三区欧美精品| 91精品国产九色| 欧美bdsm另类| 国产色爽女视频免费观看| 日韩三级伦理在线观看| 在线观看免费视频网站a站| 91久久精品国产一区二区三区| 人妻一区二区av| 亚洲,一卡二卡三卡| 看免费成人av毛片| 男男h啪啪无遮挡| av天堂中文字幕网| 麻豆成人午夜福利视频| 国产精品偷伦视频观看了| 欧美日韩亚洲高清精品| 极品人妻少妇av视频| 乱码一卡2卡4卡精品| 国产成人精品婷婷| 中文资源天堂在线| 91精品伊人久久大香线蕉| 国产精品久久久久久久电影| 简卡轻食公司| 日产精品乱码卡一卡2卡三| 久热久热在线精品观看| www.色视频.com| 日本欧美视频一区| 一级毛片电影观看| 伊人久久精品亚洲午夜| 插阴视频在线观看视频| 久久免费观看电影| 婷婷色av中文字幕| 日本wwww免费看| 免费少妇av软件| 黑人猛操日本美女一级片| 亚洲,欧美,日韩| 黄色毛片三级朝国网站 | 国产日韩欧美视频二区| 久久久久精品久久久久真实原创| 日本91视频免费播放| 免费久久久久久久精品成人欧美视频 | 中文精品一卡2卡3卡4更新| 日韩av不卡免费在线播放| 国产成人aa在线观看| 国产午夜精品久久久久久一区二区三区| 精品99又大又爽又粗少妇毛片| 亚洲综合色惰| 日韩精品免费视频一区二区三区 | 日本黄大片高清| 国产中年淑女户外野战色| 欧美日韩在线观看h| 狠狠精品人妻久久久久久综合| 最近最新中文字幕免费大全7| 精品一区二区三卡| 亚洲情色 制服丝袜| 亚洲经典国产精华液单| 国模一区二区三区四区视频| 亚洲av日韩在线播放| 人妻系列 视频| 国内揄拍国产精品人妻在线| 少妇精品久久久久久久| 男女国产视频网站| 女性生殖器流出的白浆| 男人爽女人下面视频在线观看| 亚洲精品自拍成人| 蜜桃久久精品国产亚洲av| 97超视频在线观看视频| av视频免费观看在线观看| 啦啦啦中文免费视频观看日本| 插逼视频在线观看| 韩国av在线不卡| 少妇人妻 视频| 亚洲国产毛片av蜜桃av| 我要看日韩黄色一级片| 女性生殖器流出的白浆| av女优亚洲男人天堂| 亚洲美女黄色视频免费看| 亚洲精品乱码久久久v下载方式| 亚洲一级一片aⅴ在线观看| 亚洲国产毛片av蜜桃av| 欧美性感艳星| 黑丝袜美女国产一区| 国产成人aa在线观看| 观看av在线不卡| 欧美精品高潮呻吟av久久| 激情五月婷婷亚洲| 亚洲经典国产精华液单| 国产美女午夜福利| 精华霜和精华液先用哪个| 国产高清三级在线| 欧美性感艳星| 一级毛片黄色毛片免费观看视频| 一级黄片播放器| 亚洲av成人精品一二三区| 日韩精品免费视频一区二区三区 | 国产精品成人在线| 蜜桃久久精品国产亚洲av| 日韩av免费高清视频| 欧美+日韩+精品| 最近手机中文字幕大全| 国产精品福利在线免费观看| 97在线人人人人妻| 欧美另类一区| 免费av不卡在线播放| 下体分泌物呈黄色| av免费观看日本| 精品久久久噜噜| 天美传媒精品一区二区| 寂寞人妻少妇视频99o| 色视频www国产| 成人国产av品久久久| 亚洲欧美一区二区三区黑人 | 久久久久视频综合| 欧美激情国产日韩精品一区| 国模一区二区三区四区视频| 男人舔奶头视频| 99久久人妻综合| 成人特级av手机在线观看| 久久久国产精品麻豆| 国产精品99久久99久久久不卡 | 99久久综合免费| 欧美精品一区二区大全| 99热这里只有精品一区| 夫妻午夜视频| 国产精品一区二区三区四区免费观看| 国产精品一区www在线观看| av免费观看日本| 亚洲三级黄色毛片| 中文字幕av电影在线播放| 亚洲av欧美aⅴ国产| 久久久久久人妻| 97在线视频观看| 女性生殖器流出的白浆| 久久久久国产精品人妻一区二区| 免费播放大片免费观看视频在线观看| 久热这里只有精品99| 久久 成人 亚洲| 三级国产精品欧美在线观看| 欧美日韩亚洲高清精品| 一级a做视频免费观看| 免费大片黄手机在线观看| 下体分泌物呈黄色| 亚洲怡红院男人天堂| kizo精华| 国产精品久久久久成人av| 亚洲欧美精品专区久久| 国产一区有黄有色的免费视频| 黄色怎么调成土黄色| 高清黄色对白视频在线免费看 | 午夜激情久久久久久久| 五月开心婷婷网| 久久久久国产精品人妻一区二区| 少妇人妻久久综合中文| 午夜福利,免费看| 国产精品熟女久久久久浪| 精品久久久久久久久亚洲| 日韩欧美一区视频在线观看 | 一级,二级,三级黄色视频| 高清在线视频一区二区三区| 偷拍熟女少妇极品色| 大话2 男鬼变身卡| 国产爽快片一区二区三区| 亚洲精品国产色婷婷电影| 国产有黄有色有爽视频| 国产精品伦人一区二区| 一级a做视频免费观看| 欧美最新免费一区二区三区| 韩国高清视频一区二区三区| 桃花免费在线播放| 99re6热这里在线精品视频| 亚洲国产色片| 菩萨蛮人人尽说江南好唐韦庄| 亚洲国产精品999| 亚洲精品日本国产第一区| av福利片在线| 18+在线观看网站| 欧美 亚洲 国产 日韩一| 成人黄色视频免费在线看| 99久久精品一区二区三区| 欧美日韩国产mv在线观看视频| 免费高清在线观看视频在线观看| 99热这里只有是精品50| 91久久精品国产一区二区成人| 亚洲av国产av综合av卡| 99热这里只有是精品50| 少妇猛男粗大的猛烈进出视频| 亚洲激情五月婷婷啪啪| tube8黄色片| 最近中文字幕高清免费大全6| 男女啪啪激烈高潮av片| 国产日韩欧美亚洲二区| 久久精品久久精品一区二区三区| 亚洲激情五月婷婷啪啪| 亚洲国产毛片av蜜桃av| 午夜福利视频精品| 国产精品成人在线| 久久婷婷青草| 一本久久精品| 99热国产这里只有精品6| 国产午夜精品一二区理论片| 一区二区三区乱码不卡18| 少妇人妻一区二区三区视频| 成人午夜精彩视频在线观看| 丰满少妇做爰视频| 五月伊人婷婷丁香|