• <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)移的壓力特性
    亚洲最大成人中文| 淫妇啪啪啪对白视频| 女人精品久久久久毛片| 性欧美人与动物交配| 成年女人毛片免费观看观看9| 精品久久久久久成人av| 18禁黄网站禁片午夜丰满| 99久久国产精品久久久| 国产欧美日韩一区二区精品| 国产麻豆成人av免费视频| 中文字幕人妻熟女乱码| 日本三级黄在线观看| 国产亚洲精品一区二区www| 国产精品99久久99久久久不卡| 亚洲国产精品999在线| 久久精品人人爽人人爽视色| 热99re8久久精品国产| 亚洲在线自拍视频| 丝袜在线中文字幕| 亚洲精品国产精品久久久不卡| 午夜福利视频1000在线观看 | 精品一品国产午夜福利视频| av网站免费在线观看视频| 性少妇av在线| 无限看片的www在线观看| 波多野结衣高清无吗| 99在线人妻在线中文字幕| 亚洲人成电影观看| 黑人欧美特级aaaaaa片| 国产一区二区三区视频了| 午夜免费观看网址| 午夜福利视频1000在线观看 | 天堂√8在线中文| 精品乱码久久久久久99久播| 久久国产精品人妻蜜桃| 精品日产1卡2卡| 老汉色∧v一级毛片| 日韩中文字幕欧美一区二区| 天天躁夜夜躁狠狠躁躁| 久久人人97超碰香蕉20202| 久久精品国产亚洲av香蕉五月| 18禁美女被吸乳视频| 在线av久久热| 99国产精品一区二区三区| 午夜免费成人在线视频| 欧美日韩乱码在线| 中出人妻视频一区二区| 色婷婷久久久亚洲欧美| 久久亚洲精品不卡| av视频在线观看入口| 国产97色在线日韩免费| 757午夜福利合集在线观看| 国产野战对白在线观看| 久久久精品欧美日韩精品| 美女 人体艺术 gogo| 看免费av毛片| 成人国产一区最新在线观看| 日韩欧美免费精品| 性欧美人与动物交配| 亚洲专区字幕在线| 亚洲国产精品999在线| 免费女性裸体啪啪无遮挡网站| 精品一区二区三区四区五区乱码| 女人爽到高潮嗷嗷叫在线视频| 亚洲精品在线观看二区| а√天堂www在线а√下载| 老司机午夜福利在线观看视频| 91九色精品人成在线观看| 日本免费一区二区三区高清不卡 | 亚洲情色 制服丝袜| 精品少妇一区二区三区视频日本电影| 亚洲成av片中文字幕在线观看| 午夜福利视频1000在线观看 | 国产精品自产拍在线观看55亚洲| 精品久久蜜臀av无| 又紧又爽又黄一区二区| 一二三四社区在线视频社区8| 两性夫妻黄色片| 久久国产精品人妻蜜桃| 久久人人97超碰香蕉20202| 国产精品日韩av在线免费观看 | 99re在线观看精品视频| 久久精品91无色码中文字幕| 亚洲欧美日韩高清在线视频| 日韩欧美国产一区二区入口| 亚洲av第一区精品v没综合| 久久久久国产一级毛片高清牌| 久久精品亚洲精品国产色婷小说| 国产一区在线观看成人免费| 久久中文字幕一级| 成人18禁在线播放| 国产亚洲精品久久久久久毛片| 亚洲激情在线av| 亚洲五月婷婷丁香| 色综合婷婷激情| 乱人伦中国视频| 女性生殖器流出的白浆| 中文字幕色久视频| 国产精品1区2区在线观看.| 久久精品aⅴ一区二区三区四区| 最近最新免费中文字幕在线| 50天的宝宝边吃奶边哭怎么回事| 韩国精品一区二区三区| 亚洲av熟女| 国产成人啪精品午夜网站| 黑人巨大精品欧美一区二区mp4| 人人妻人人爽人人添夜夜欢视频| 三级毛片av免费| 女人爽到高潮嗷嗷叫在线视频| 啦啦啦免费观看视频1| 亚洲五月天丁香| 亚洲天堂国产精品一区在线| 成人国产一区最新在线观看| 免费人成视频x8x8入口观看| 黄色a级毛片大全视频| 悠悠久久av| 亚洲中文av在线| 免费不卡黄色视频| 精品福利观看| 国产黄a三级三级三级人| bbb黄色大片| 国产亚洲精品一区二区www| 日本欧美视频一区| 香蕉国产在线看| 欧美最黄视频在线播放免费| 亚洲国产毛片av蜜桃av| 丝袜人妻中文字幕| 热re99久久国产66热| 涩涩av久久男人的天堂| 日本 欧美在线| xxx96com| 国产亚洲精品第一综合不卡| 久久久久国内视频| 看片在线看免费视频| 嫩草影视91久久| 欧美日韩亚洲综合一区二区三区_| 老司机午夜福利在线观看视频| 在线观看66精品国产| 精品乱码久久久久久99久播| 性少妇av在线| 欧美精品啪啪一区二区三区| 午夜免费成人在线视频| 亚洲va日本ⅴa欧美va伊人久久| 日韩高清综合在线| 给我免费播放毛片高清在线观看| 欧美成狂野欧美在线观看| 最新在线观看一区二区三区| 国产99久久九九免费精品| 欧美精品亚洲一区二区| 午夜福利高清视频| 此物有八面人人有两片| 国产成年人精品一区二区| 搡老妇女老女人老熟妇| 日韩高清综合在线| 午夜精品久久久久久毛片777| 日本vs欧美在线观看视频| 国产精品 欧美亚洲| 91九色精品人成在线观看| 老司机靠b影院| 国产主播在线观看一区二区| 国产成人精品在线电影| 香蕉国产在线看| 高清在线国产一区| 免费无遮挡裸体视频| 精品久久久久久,| 欧美日韩中文字幕国产精品一区二区三区 | 操美女的视频在线观看| 一级毛片女人18水好多| 国产精品亚洲一级av第二区| 黄色 视频免费看| 欧美黄色片欧美黄色片| 性色av乱码一区二区三区2| 欧美色欧美亚洲另类二区 | 中文字幕人成人乱码亚洲影| 精品人妻在线不人妻| 一夜夜www| 免费在线观看黄色视频的| 中文字幕最新亚洲高清| 看片在线看免费视频| 国产精品免费视频内射| 国产成人系列免费观看| 人人妻人人爽人人添夜夜欢视频| 亚洲一码二码三码区别大吗| 大香蕉久久成人网| 人人妻人人澡人人看| 久久午夜综合久久蜜桃| 国产精品1区2区在线观看.| 在线观看一区二区三区| 国产精品爽爽va在线观看网站 | 好男人电影高清在线观看| 亚洲第一av免费看| 国产精品av久久久久免费| 日韩有码中文字幕| 国产精品一区二区免费欧美| 久久香蕉精品热| 国产精品久久久av美女十八| a在线观看视频网站| 波多野结衣巨乳人妻| 久久狼人影院| aaaaa片日本免费| 一级黄色大片毛片| 亚洲熟妇中文字幕五十中出| 亚洲人成电影观看| 别揉我奶头~嗯~啊~动态视频| 精品国产美女av久久久久小说| 亚洲精品国产精品久久久不卡| 亚洲精品久久国产高清桃花| 国产97色在线日韩免费| 真人做人爱边吃奶动态| 欧美一区二区精品小视频在线| 亚洲熟妇熟女久久| 99香蕉大伊视频| 看免费av毛片| 久久狼人影院| 免费不卡黄色视频| 少妇粗大呻吟视频| 18禁美女被吸乳视频| 黄色丝袜av网址大全| 又紧又爽又黄一区二区| 两个人视频免费观看高清| 欧美久久黑人一区二区| 国产精品久久久av美女十八| 日韩有码中文字幕| 在线观看免费视频日本深夜| 大型av网站在线播放| 可以在线观看毛片的网站| 午夜免费成人在线视频| 男人的好看免费观看在线视频 | 一个人观看的视频www高清免费观看 | 精品电影一区二区在线| 18美女黄网站色大片免费观看| 久久 成人 亚洲| 在线观看一区二区三区| 久久久久久亚洲精品国产蜜桃av| 一区二区三区高清视频在线| 午夜福利在线观看吧| 丝袜美足系列| 男女床上黄色一级片免费看| 亚洲第一电影网av| 美女免费视频网站| 日韩一卡2卡3卡4卡2021年| 一二三四在线观看免费中文在| 久久久久久久午夜电影| 91老司机精品| 免费少妇av软件| av福利片在线| 身体一侧抽搐| 久久久久国内视频| 亚洲av五月六月丁香网| 香蕉丝袜av| 国内久久婷婷六月综合欲色啪| 此物有八面人人有两片| 99香蕉大伊视频| 亚洲av五月六月丁香网| а√天堂www在线а√下载| 一级毛片精品| 亚洲熟妇中文字幕五十中出| 最新美女视频免费是黄的| 村上凉子中文字幕在线| 人人妻人人澡人人看| 在线观看免费视频日本深夜| 日本 欧美在线| 一本综合久久免费| 亚洲av日韩精品久久久久久密| 中文字幕最新亚洲高清| 两个人免费观看高清视频| 亚洲欧美激情综合另类| 身体一侧抽搐| 久久九九热精品免费| 天天一区二区日本电影三级 | 成人永久免费在线观看视频| 亚洲熟妇熟女久久| 精品国产乱子伦一区二区三区| av欧美777| 两人在一起打扑克的视频| 亚洲精品国产区一区二| 男人的好看免费观看在线视频 | 一本综合久久免费| 久久久精品欧美日韩精品| 人妻丰满熟妇av一区二区三区| 一a级毛片在线观看| 一进一出抽搐动态| 午夜福利高清视频| 别揉我奶头~嗯~啊~动态视频| 在线观看66精品国产| 一夜夜www| 亚洲国产欧美网| www.999成人在线观看| 免费久久久久久久精品成人欧美视频| 日本欧美视频一区| 久久精品国产99精品国产亚洲性色 | 欧美黑人欧美精品刺激| 欧美另类亚洲清纯唯美| 午夜福利欧美成人| 精品欧美一区二区三区在线| av天堂在线播放| 欧美人与性动交α欧美精品济南到| 欧美日韩精品网址| 99国产精品一区二区蜜桃av| 欧美精品啪啪一区二区三区| 一夜夜www| 久久精品国产综合久久久| 国产免费av片在线观看野外av| 国内精品久久久久久久电影| 此物有八面人人有两片| 国产麻豆69| 国产精品一区二区精品视频观看| 日韩视频一区二区在线观看| 色播在线永久视频| 国产精品日韩av在线免费观看 | 国产av精品麻豆| 久久青草综合色| 两个人免费观看高清视频| 成年女人毛片免费观看观看9| 亚洲伊人色综图| 免费在线观看黄色视频的| 妹子高潮喷水视频| 两性夫妻黄色片| 欧美成人午夜精品| 亚洲国产欧美一区二区综合| 日韩精品中文字幕看吧| 亚洲一码二码三码区别大吗| 搞女人的毛片| 黄色视频,在线免费观看| 国产成人免费无遮挡视频| 欧美激情 高清一区二区三区| 欧美日韩亚洲国产一区二区在线观看| 一级a爱片免费观看的视频| 97超级碰碰碰精品色视频在线观看| 丝袜在线中文字幕| 69精品国产乱码久久久| 老司机靠b影院| 91国产中文字幕| 麻豆成人av在线观看| 欧美成人性av电影在线观看| 欧美 亚洲 国产 日韩一| 久久精品人人爽人人爽视色| 亚洲精华国产精华液的使用体验 | 日本精品一区二区三区蜜桃| 成人美女网站在线观看视频| АⅤ资源中文在线天堂| 免费在线观看日本一区| 久久久久久国产a免费观看| 亚洲在线观看片| 亚洲av免费在线观看| 亚洲人成网站在线播| 国产精品无大码| aaaaa片日本免费| 欧美+亚洲+日韩+国产| 欧美xxxx性猛交bbbb| 日本免费一区二区三区高清不卡| 国产麻豆成人av免费视频| 欧美三级亚洲精品| 日本欧美国产在线视频| 永久网站在线| 好男人在线观看高清免费视频| 亚洲中文字幕日韩| 亚洲av免费高清在线观看| 亚洲电影在线观看av| 国产精品久久视频播放| 欧美日韩综合久久久久久 | 人妻制服诱惑在线中文字幕| 日本a在线网址| 国产免费男女视频| 日本免费a在线| av在线亚洲专区| 99九九线精品视频在线观看视频| 亚洲成av人片在线播放无| 色综合站精品国产| 在线观看舔阴道视频| 国产精品日韩av在线免费观看| 99久久精品一区二区三区| 天堂√8在线中文| 免费在线观看日本一区| 成年女人永久免费观看视频| 天堂av国产一区二区熟女人妻| 久久亚洲精品不卡| 俺也久久电影网| 男女之事视频高清在线观看| 尤物成人国产欧美一区二区三区| 麻豆久久精品国产亚洲av| 最近中文字幕高清免费大全6 | 男女视频在线观看网站免费| 欧美色视频一区免费| 大又大粗又爽又黄少妇毛片口| 一级a爱片免费观看的视频| 亚洲人成伊人成综合网2020| 国产毛片a区久久久久| 看十八女毛片水多多多| 搡老妇女老女人老熟妇| 欧美潮喷喷水| 成人国产麻豆网| 韩国av在线不卡| 在线a可以看的网站| 欧美性猛交黑人性爽| 亚洲人成网站高清观看| 成人国产麻豆网| aaaaa片日本免费| 亚洲中文字幕一区二区三区有码在线看| 亚洲七黄色美女视频| 日日夜夜操网爽| 日韩精品有码人妻一区| xxxwww97欧美| 12—13女人毛片做爰片一| 校园人妻丝袜中文字幕| 国产精品野战在线观看| 日韩欧美一区二区三区在线观看| 无人区码免费观看不卡| 亚洲国产高清在线一区二区三| 深爱激情五月婷婷| 亚洲av成人精品一区久久| 五月伊人婷婷丁香| 无人区码免费观看不卡| 91久久精品电影网| 亚洲最大成人手机在线| 看十八女毛片水多多多| 国产高清有码在线观看视频| 伊人久久精品亚洲午夜| 国产精品日韩av在线免费观看| 尤物成人国产欧美一区二区三区| 精品久久久久久久久亚洲 | 亚洲电影在线观看av| 最后的刺客免费高清国语| 国产精品98久久久久久宅男小说| 久久久久久久午夜电影| 国产 一区 欧美 日韩| 精品久久久久久久久亚洲 | 国产男靠女视频免费网站| 国产真实乱freesex| 国产熟女欧美一区二区| 男女之事视频高清在线观看| 亚洲中文日韩欧美视频| 十八禁网站免费在线| 亚洲国产精品久久男人天堂| 亚洲精品日韩av片在线观看| 狂野欧美白嫩少妇大欣赏| 国产极品精品免费视频能看的| 欧美成人a在线观看| 欧美绝顶高潮抽搐喷水| 亚洲无线在线观看| 亚洲人成网站高清观看| 99久久精品热视频| 国产精品伦人一区二区| 在线观看av片永久免费下载| a级毛片a级免费在线| 国产精品免费一区二区三区在线| 午夜视频国产福利| 久9热在线精品视频| 黄色丝袜av网址大全| 成年版毛片免费区| 色综合亚洲欧美另类图片| 精品久久久久久久末码| 老女人水多毛片| 精品一区二区三区人妻视频| a级毛片a级免费在线| 亚洲欧美激情综合另类| 老司机深夜福利视频在线观看| 搡老岳熟女国产| 久久天躁狠狠躁夜夜2o2o| 成人综合一区亚洲| 一级a爱片免费观看的视频| 午夜免费激情av| 91在线观看av| 搡老岳熟女国产| 好男人在线观看高清免费视频| 日日摸夜夜添夜夜添小说| 亚洲自拍偷在线| 久久久久国内视频| 国产av麻豆久久久久久久| 热99在线观看视频| 熟女人妻精品中文字幕| 尾随美女入室| 国产精品久久久久久精品电影| 能在线免费观看的黄片| 免费av不卡在线播放| 成人特级av手机在线观看| 久久久久久久久久久丰满 | 在线观看免费视频日本深夜| 日韩欧美精品免费久久| 亚洲一区二区三区色噜噜| 别揉我奶头~嗯~啊~动态视频| 嫩草影视91久久| 国产69精品久久久久777片| 精品国产三级普通话版| 久久久久久久久久久丰满 | 老司机午夜福利在线观看视频| 国产高清激情床上av| 国产大屁股一区二区在线视频| 日本精品一区二区三区蜜桃| 少妇的逼水好多| 日韩强制内射视频| 国产伦在线观看视频一区| 校园人妻丝袜中文字幕| 99视频精品全部免费 在线| 亚洲四区av| 亚洲专区中文字幕在线| 最新在线观看一区二区三区| 最新中文字幕久久久久| 天堂网av新在线| 亚洲最大成人av| 内射极品少妇av片p| av视频在线观看入口| 亚洲成人免费电影在线观看| 国产人妻一区二区三区在| 精品人妻1区二区| 欧美成人免费av一区二区三区| 一区福利在线观看| 中文字幕久久专区| 国产精品日韩av在线免费观看| av在线亚洲专区| 夜夜爽天天搞| 亚洲一区二区三区色噜噜| 国产亚洲欧美98| 国内揄拍国产精品人妻在线| 亚洲在线自拍视频| 亚洲精品国产成人久久av| 国产精品亚洲美女久久久| 午夜老司机福利剧场| 人妻夜夜爽99麻豆av| 欧美色欧美亚洲另类二区| 日韩中字成人| 亚洲欧美日韩高清在线视频| 嫩草影院精品99| h日本视频在线播放| 国产精品嫩草影院av在线观看 | 一区二区三区免费毛片| 12—13女人毛片做爰片一| 性欧美人与动物交配| 久久精品国产亚洲av香蕉五月| 精品一区二区三区av网在线观看| 最近视频中文字幕2019在线8| .国产精品久久| 国产淫片久久久久久久久| 国内精品美女久久久久久| 舔av片在线| 中文在线观看免费www的网站| 男人舔奶头视频| 校园人妻丝袜中文字幕| 久久久久免费精品人妻一区二区| 日本精品一区二区三区蜜桃| 免费看av在线观看网站| 日韩人妻高清精品专区| 亚洲黑人精品在线| 尾随美女入室| 女生性感内裤真人,穿戴方法视频| 老司机福利观看| 男女啪啪激烈高潮av片| 一夜夜www| 欧美最黄视频在线播放免费| 国产av在哪里看| 69av精品久久久久久| 日日摸夜夜添夜夜添小说| 九九久久精品国产亚洲av麻豆| 最近在线观看免费完整版| 婷婷亚洲欧美| 亚洲av日韩精品久久久久久密| 亚洲在线自拍视频| av在线观看视频网站免费| 国产精品久久视频播放| 午夜免费男女啪啪视频观看 | 欧美最新免费一区二区三区| or卡值多少钱| 国产色婷婷99| 欧美bdsm另类| 亚洲人成伊人成综合网2020| 成年女人看的毛片在线观看| 两个人视频免费观看高清| 欧美黑人巨大hd| 老师上课跳d突然被开到最大视频| 村上凉子中文字幕在线| 人妻少妇偷人精品九色| 国产精品亚洲一级av第二区| 免费搜索国产男女视频| 亚洲无线观看免费| 国产一区二区亚洲精品在线观看| 最近最新中文字幕大全电影3| 亚洲精品一区av在线观看| www.www免费av| 国产精品亚洲美女久久久| 十八禁网站免费在线| 久久久久精品国产欧美久久久| 日韩一本色道免费dvd| 欧美xxxx黑人xx丫x性爽| 两性午夜刺激爽爽歪歪视频在线观看| 国产免费av片在线观看野外av| 91精品国产九色| 此物有八面人人有两片| 好男人在线观看高清免费视频| 久久精品91蜜桃| 在线观看午夜福利视频| 亚洲经典国产精华液单| 91在线精品国自产拍蜜月| 尾随美女入室| 91久久精品国产一区二区三区| 久久久久久久亚洲中文字幕| 少妇被粗大猛烈的视频| bbb黄色大片| 麻豆久久精品国产亚洲av| 久久久久国内视频| 欧美区成人在线视频| 日本 欧美在线| 国产精品98久久久久久宅男小说| 校园人妻丝袜中文字幕| 国产伦一二天堂av在线观看| 听说在线观看完整版免费高清| 亚洲av中文字字幕乱码综合| 亚洲熟妇熟女久久| 精品一区二区三区av网在线观看| 午夜a级毛片| 蜜桃久久精品国产亚洲av| 久久这里只有精品中国| 欧美日韩亚洲国产一区二区在线观看| 国产午夜精品久久久久久一区二区三区 | 国产v大片淫在线免费观看| 日韩欧美精品v在线|