徐 源，張歡歡，邢蘭昌?，魏 偉，韓維峰，魏周拓，鄧少貴

基于電-力-聲多物理場(chǎng)耦合數(shù)值模型的含水合物多孔介質(zhì)聲速和衰減特性研究*

徐 源1，張歡歡1，邢蘭昌1?，魏 偉2，韓維峰2，魏周拓3，鄧少貴3

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

在物理模擬實(shí)驗(yàn)中對(duì)水合物微觀(guān)賦存模式和飽和度進(jìn)行準(zhǔn)確控制和評(píng)價(jià)尚存在技術(shù)困難，僅依賴(lài)實(shí)驗(yàn)技術(shù)研究含水合物沉積物聲學(xué)特性、建立儲(chǔ)層參數(shù)解釋模型存在局限性。采用基于有限元的數(shù)字巖石物理技術(shù)，針對(duì)懸浮、接觸和膠結(jié)三種典型的水合物微觀(guān)賦存模式分別建立多孔介質(zhì)的三維電?力?聲多物理場(chǎng)耦合模型，考察了微觀(guān)賦存模式和水合物飽和度對(duì)多孔介質(zhì)聲速和衰減的影響規(guī)律，對(duì)比了聲速數(shù)值模擬與理論模型計(jì)算結(jié)果，建立了聲波衰減參數(shù)與水合物飽和度之間的關(guān)系式。研究結(jié)果表明：（1）對(duì)于三種水合物賦存模式，由于水合物相比孔隙水具有更高的彈性模量，多孔介質(zhì)的聲速隨著水合物飽和度的增大而增大；水合物的存在導(dǎo)致聲波在傳播過(guò)程中遇到更多不連續(xù)的聲阻抗界面，聲衰減隨著水合物飽和度的增大而近似線(xiàn)性增大；（2）懸浮和接觸賦存模式條件下，水合物飽和度對(duì)多孔介質(zhì)的聲速和衰減影響規(guī)律基本一致；對(duì)于相同的水合物飽和度，膠結(jié)模式條件下含水合物多孔介質(zhì)具有更高的聲速和更小的聲衰減；（3）通過(guò)合理選擇參數(shù)值，利用權(quán)重方程與Lee改進(jìn)的Biot-Gassmann Theory（BGTL）模型估算的含懸浮和接觸模式水合物多孔介質(zhì)的聲速較為準(zhǔn)確；通過(guò)等效介質(zhì)理論模型C計(jì)算的含膠結(jié)模式水合物多孔介質(zhì)的聲速更為準(zhǔn)確。研究結(jié)果可為獲取復(fù)雜條件下含水合物沉積物的聲學(xué)特性提供數(shù)值建模方法，為基于聲波測(cè)井?dāng)?shù)據(jù)的水合物儲(chǔ)層精細(xì)評(píng)價(jià)提供理論支撐。

天然氣水合物；數(shù)字巖石物理；微觀(guān)賦存模式；水合物飽和度；聲速；聲衰減

0 引 言

天然氣水合物是由水和天然氣分子在低溫高壓環(huán)境下結(jié)合形成的與冰相似的化合物，其廣泛分布于海洋沉積物和凍土區(qū)域中[1-4]，是一種極具開(kāi)發(fā)前景且儲(chǔ)量巨大的新型清潔能源[5-7]。地球物理測(cè)井是準(zhǔn)確識(shí)別和評(píng)價(jià)水合物儲(chǔ)層的重要手段，其中聲波測(cè)井是評(píng)價(jià)水合物飽和度（h）的重要方法之一[8-9]。研究含水合物沉積物的聲學(xué)響應(yīng)特性，深入探究水合物飽和度、微觀(guān)賦存模式等因素對(duì)沉積物聲波傳播速度及其能量衰減的影響規(guī)律，將為聲波測(cè)井?dāng)?shù)據(jù)解釋、水合物儲(chǔ)層識(shí)別與評(píng)價(jià)提供理論基礎(chǔ)。

天然氣水合物在儲(chǔ)層中具有多種微觀(guān)賦存模式，不同的微觀(guān)賦存模式會(huì)對(duì)沉積物的聲學(xué)、電學(xué)、力學(xué)、熱學(xué)等物性參數(shù)產(chǎn)生不同程度的影響[10-13]。為了準(zhǔn)確評(píng)價(jià)儲(chǔ)層中水合物的飽和度，在研究水合物飽和度與沉積物中聲速和衰減之間的定量關(guān)系時(shí)，必須考慮水合物微觀(guān)賦存模式對(duì)沉積物聲速及衰減的影響規(guī)律。國(guó)內(nèi)外有許多學(xué)者采用室內(nèi)模擬實(shí)驗(yàn)方法研究了含水合物（模擬）沉積物的聲學(xué)特性[14-20]。文獻(xiàn)[14-17]研究了水合物的存在對(duì)不同類(lèi)型沉積物（固結(jié)、松散）聲速的影響。PRIEST等[18-19]研究了不同體系（過(guò)量氣、定量氣、溶解氣）中水合物的生成對(duì)沉積物聲速的影響規(guī)律。SAHOO等[21]和LIU等[22]分別研究了聲波頻率和水合物飽和度對(duì)沉積物聲波衰減特性的影響規(guī)律。雖然研究者提出了制備具有不同水合物飽和度、不同微觀(guān)分布模式的含水合物沉積物樣品的實(shí)驗(yàn)方法，但是對(duì)于所測(cè)試的整個(gè)樣品而言，實(shí)驗(yàn)中所得到的飽和度、微觀(guān)分布模式是否準(zhǔn)確、是否具有重復(fù)性、是否方便控制，仍然是值得進(jìn)一步探討的問(wèn)題。而數(shù)值模擬方法能夠準(zhǔn)確、方便地控制水合物飽和度、微觀(guān)賦存模式等影響因素，同時(shí)數(shù)值模擬方法具有重復(fù)性好、時(shí)間成本低的顯著優(yōu)勢(shì)。

數(shù)值模擬方法在巖石物理領(lǐng)域具有廣泛的應(yīng)用，逐漸形成了數(shù)字巖石物理技術(shù)。利用數(shù)字巖石物理技術(shù)不僅能夠?qū)崿F(xiàn)巖石孔隙尺度建模，也能夠研究微觀(guān)因素對(duì)巖石聲速、電阻率、滲透率等物理參數(shù)的影響[23]。針對(duì)巖石聲學(xué)特性的數(shù)值模擬方法主要包括旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分方法、聲格子方法和有限元方法。SAENGER等[24]采用旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分方法模擬出聲波在各向異性、黏彈性介質(zhì)中的傳播過(guò)程。DEL VALLE-GARCíA等[25]利用聲格子方法研究了聲波在多孔介質(zhì)中的傳播過(guò)程和規(guī)律。MAKARYNSKA等[26]利用有限元方法研究了不同含水飽和度巖石的聲學(xué)特性。ZHU等[27]利用有限差分方法實(shí)現(xiàn)了油氣儲(chǔ)層介質(zhì)彈性波場(chǎng)的數(shù)值模擬。孫建孟等[28]、劉向君等[29]采用有限元方法計(jì)算了數(shù)字巖心的彈性參數(shù)。YANG等[30]采用旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分方法模擬了地震波在海底天然氣水合物沉積地層中的傳播。卜慶濤等[31]基于實(shí)驗(yàn)條件對(duì)水合物模擬儲(chǔ)層開(kāi)展了聲學(xué)特性數(shù)值模擬，發(fā)現(xiàn)水合物儲(chǔ)層聲速的實(shí)驗(yàn)與數(shù)值模擬結(jié)果對(duì)應(yīng)關(guān)系良好。目前，基于數(shù)字巖石物理技術(shù)研究微觀(guān)賦存模式、水合物飽和度等因素對(duì)含水合物沉積物聲學(xué)特性的影響規(guī)律鮮有報(bào)道。

本文采用基于有限元的數(shù)字巖石物理技術(shù)，針對(duì)懸浮、接觸和膠結(jié)三種水合物微觀(guān)賦存模式分別建立多孔介質(zhì)的三維電?力?聲多物理場(chǎng)耦合模型，研究了微觀(guān)賦存模式和水合物飽和度對(duì)多孔介質(zhì)聲速和衰減的影響規(guī)律。研究結(jié)果可為含水合物多孔介質(zhì)聲學(xué)特性數(shù)值模擬提供建模方法，為聲波測(cè)井在水合物儲(chǔ)層精細(xì)評(píng)價(jià)中的應(yīng)用提供理論依據(jù)和模型支撐。

1 含水合物多孔介質(zhì)多物理場(chǎng)有限元建模

1.1 幾何結(jié)構(gòu)與材料

采用多物理場(chǎng)耦合計(jì)算平臺(tái)（COMSOL Multiphysics）來(lái)建立含水合物多孔介質(zhì)的有限元數(shù)值模型。數(shù)值模型的幾何結(jié)構(gòu)包含壓電晶片和多孔介質(zhì)兩部分。如圖1所示，兩塊壓電晶片位于多孔介質(zhì)兩側(cè)，分別用于模擬發(fā)射和接收聲波。兩塊壓電晶片的幾何尺寸和材料相同，晶片厚度為1 mm，邊長(zhǎng)為1.624 mm。水飽和多孔介質(zhì)為正方體，邊長(zhǎng)為1.624 mm。幾何結(jié)構(gòu)包含64個(gè)大小相同的球形石英砂顆粒，半徑為0.21 mm，石英砂顆粒之間規(guī)則相切排布。多孔介質(zhì)的孔隙度為相互連通的孔隙體積與總體積的比值，本文所構(gòu)建的多孔介質(zhì)幾何結(jié)構(gòu)的孔隙體積約為1.8 mm3，孔隙度為0.42。

依據(jù)ECKER[32]提出的懸浮、接觸和膠結(jié)三種典型的水合物微觀(guān)賦存模式（如圖2所示），構(gòu)建含水合物多孔介質(zhì)模型的三維幾何結(jié)構(gòu)。圖3所示為三種水合物微觀(guān)賦存模式條件下的含水合物多孔介質(zhì)的三維立體結(jié)構(gòu)。含懸浮模式水合物多孔介質(zhì)截面中，代表水合物的9個(gè)球體均勻分布于孔隙空間中，水合物與石英砂顆粒不接觸。含接觸模式水合物多孔介質(zhì)截面中，代表水合物的9個(gè)球體與石英砂顆粒相接觸（相交）。含膠結(jié)模式水合物多孔介質(zhì)截面中，代表水合物的薄層包裹在石英砂顆粒外表面，相鄰石英砂顆粒表面的水合物層相互接觸（相交），將石英砂顆粒膠結(jié)在一起。

圖1 水飽和多孔介質(zhì)數(shù)值模型幾何結(jié)構(gòu)

圖2 沉積物中水合物的三種微觀(guān)賦存模式（改自文獻(xiàn)[32]）

圖3 含不同微觀(guān)賦存模式水合物的多孔介質(zhì)模型三維立體結(jié)構(gòu)：（a）懸浮模式；（b）接觸模式；（c）膠結(jié)模式

為得到不同的水合物飽和度，對(duì)于懸浮和接觸模式，保持代表水合物的球體個(gè)數(shù)不變，改變球體的直徑；對(duì)于膠結(jié)模式，改變石英砂顆粒表面水合物層的厚度。受到三維模型幾何結(jié)構(gòu)的限制（保持水合物為單一球形顆粒），懸浮模式條件下模型所模擬的水合物飽和度存在上限，當(dāng)水合物與石英砂相切時(shí)水合物飽和度為23% ～ 24%?？紤]到地層中水合物飽和度過(guò)高時(shí)，水合物的賦存模式傾向于其他模式（如接觸模式），因此選取含2%、4%、6%、8%、10%、12%、14%、16%、18%、20%、22%等11個(gè)水合物飽和度的多孔介質(zhì)進(jìn)行數(shù)值求解。懸浮與接觸模式條件下，以上11種水合物飽和度對(duì)應(yīng)的水合物顆粒半徑的范圍為68.47 ～ 152.14 μm。膠結(jié)模式條件下對(duì)應(yīng)的水合物層厚度范圍為1.16 ～ 13.18 μm。

在分別定義石英砂、孔隙水、水合物和壓電晶片等材料的基礎(chǔ)上，根據(jù)表1設(shè)定多孔介質(zhì)中各組分的聲速和密度等參數(shù)。壓電晶片材料為鋯鈦酸鉛（PZT-5H），密度為7 500 kg/m3。

表1 多孔介質(zhì)各組分聲學(xué)特性參數(shù)[33]

1.2 物理場(chǎng)控制方程與邊界條件

采用靜電、固體力學(xué)、壓力聲學(xué)多物理場(chǎng)耦合方法實(shí)現(xiàn)數(shù)值建模方案。采用聲?結(jié)構(gòu)邊界耦合壓力聲學(xué)和固體力學(xué)兩個(gè)物理場(chǎng)，基于壓電效應(yīng)實(shí)現(xiàn)靜電和固體力學(xué)兩個(gè)物理場(chǎng)的耦合。

靜電物理場(chǎng)中控制方程為高斯方程，如式（1）所示。

式中：為電位移，C/m2；為電場(chǎng)強(qiáng)度，V/m；為電勢(shì)，V；v為自由電荷密度，C/m3。

固體力學(xué)物理場(chǎng)控制方程為運(yùn)動(dòng)方程和壓電材料本構(gòu)方程，計(jì)算得到位移、應(yīng)力和應(yīng)變值，方程如式（2）所示。

式中：為壓電晶片材料的密度，kg/m3；為位移場(chǎng)，m；V為體積力，N/m3；tt為壓電晶片振動(dòng)加速度，m/s2；為應(yīng)力，Pa；為應(yīng)變，無(wú)量綱；E與分別為壓電晶片材料的彈性矩陣（Pa）和耦合矩陣（C/m2）；0為真空介電常數(shù)，8.854 × 10?12F/m；S為壓電晶片材料的相對(duì)介電常數(shù)，無(wú)量綱。

壓力聲學(xué)物理場(chǎng)控制方程為時(shí)域波動(dòng)方程，如式（3）所示。

式中：0為聲波傳播介質(zhì)的密度，kg/m3；為介質(zhì)的聲速，m/s；t為總聲壓，Pa。

以上三個(gè)物理場(chǎng)的邊界條件設(shè)定如下：①對(duì)于靜電物理場(chǎng)，將用作發(fā)射聲波的壓電晶片內(nèi)端面定義為電勢(shì)邊界，電勢(shì)為矩形波調(diào)制諧波后的電壓脈沖，脈沖幅值0為100 V，頻率為；將用作接收聲波的壓電晶片內(nèi)端面定義為懸浮電位邊界（floating potential）；將兩塊壓電晶片的外端面定義為接地邊界。②對(duì)于固體力學(xué)物理場(chǎng)，將位移場(chǎng)初始值設(shè)為0，將兩塊壓電晶片外端面定義為輥支承邊界，防止其沿厚度方向移動(dòng)。③對(duì)于壓力聲學(xué)物理場(chǎng)，為保證聲波在抵達(dá)多孔介質(zhì)幾何外邊界時(shí)沿法向產(chǎn)生的反射最少，將多孔介質(zhì)的外邊界定義為平面輻射波邊界。

聲?結(jié)構(gòu)（本模型中的多孔介質(zhì)?壓電晶片）邊界條件如式（4）所示。

式中：為邊界上的單位法向量，無(wú)量綱；A為邊界上單位面積所受的力，Pa。

1.3 網(wǎng)格劃分與求解器

使用自由四面體網(wǎng)格對(duì)壓電晶片幾何結(jié)構(gòu)進(jìn)行網(wǎng)格劃分。采用三種數(shù)量級(jí)別的網(wǎng)格、網(wǎng)格尺寸及數(shù)量如表2所示。利用穩(wěn)態(tài)（靜態(tài)）求解器對(duì)以上描述靜電、固體力學(xué)物理場(chǎng)的微分方程進(jìn)行求解。

表2 壓電晶片部分網(wǎng)格尺寸及數(shù)量

使用自由四面體網(wǎng)格對(duì)多孔介質(zhì)幾何結(jié)構(gòu)進(jìn)行網(wǎng)格劃分。采用多種數(shù)量級(jí)別的網(wǎng)格、網(wǎng)格尺寸及數(shù)量如表3所示。使用與瞬態(tài)求解器（廣義方法）求解上述靜電、固體力學(xué)、壓力聲學(xué)物理場(chǎng)中的微分方程，選擇并行稀疏直接求解方法（multi-frontal massively parallel sparse direct solver, MUMPS）。在獲取多孔介質(zhì)聲速的計(jì)算模型中，將求解時(shí)長(zhǎng)設(shè)定為2 μs，時(shí)間步長(zhǎng)為1 ns；在獲取多孔介質(zhì)聲衰減的計(jì)算模型中，將求解時(shí)長(zhǎng)設(shè)定為20 μs，時(shí)間步長(zhǎng)為10 ns。使用組合參數(shù)掃描方法對(duì)表3中不同數(shù)量級(jí)別的網(wǎng)格依次進(jìn)行計(jì)算。

表3 多孔介質(zhì)部分網(wǎng)格尺寸及數(shù)量

1.4 網(wǎng)格選取和首波到達(dá)時(shí)刻提取方法

使用表2所示三種數(shù)量級(jí)別的網(wǎng)格測(cè)試了壓電晶片中心頻率的模型計(jì)算值對(duì)網(wǎng)格的依賴(lài)性，根據(jù)測(cè)試結(jié)果選取合適的網(wǎng)格尺寸進(jìn)行后續(xù)研究。根據(jù)計(jì)算結(jié)果得到三種網(wǎng)格下壓電晶片振動(dòng)位移隨激勵(lì)電信號(hào)頻率的變化關(guān)系，最大振動(dòng)位移對(duì)應(yīng)的頻率即為壓電晶片中心頻率。結(jié)果表明，壓電晶片的中心頻率在三種數(shù)量級(jí)別網(wǎng)格下均為642.5 kHz，由于網(wǎng)格數(shù)量成倍增加，可認(rèn)為基于這三種網(wǎng)格所得到的結(jié)果均為正確值。為減少計(jì)算時(shí)間，針對(duì)壓電晶片選擇單元數(shù)量為5724的網(wǎng)格級(jí)別，即最大網(wǎng)格尺寸為200 μm，最小網(wǎng)格尺寸為40 μm。

圖4 不同網(wǎng)格條件下接收的聲波

將激勵(lì)信號(hào)頻率設(shè)定為642.5 kHz，使用表3中不同級(jí)別網(wǎng)格，計(jì)算了水飽和多孔介質(zhì)接收聲波的首個(gè)波谷到達(dá)時(shí)刻，測(cè)試了模型計(jì)算值對(duì)網(wǎng)格的依賴(lài)性。圖4所示為接收聲波信號(hào)，信號(hào)首個(gè)波谷到達(dá)時(shí)刻隨網(wǎng)格數(shù)量變化情況如圖5所示。結(jié)果表明，當(dāng)網(wǎng)格數(shù)量達(dá)到588 993后，首個(gè)波谷到達(dá)時(shí)刻的模型計(jì)算值已經(jīng)不隨網(wǎng)格數(shù)量的增加而變化。為了平衡計(jì)算精度和計(jì)算時(shí)間，選擇單元數(shù)量為588 993的網(wǎng)格級(jí)別，對(duì)應(yīng)的最大和最小網(wǎng)格尺寸分別為60 μm和20 μm。

圖5 首個(gè)聲波波谷到達(dá)時(shí)刻隨網(wǎng)格數(shù)量的變化

從接收聲波信號(hào)中準(zhǔn)確提取首波到達(dá)時(shí)刻是獲取多孔介質(zhì)聲速的前提。將整個(gè)多孔介質(zhì)幾何結(jié)構(gòu)設(shè)定為同種材料（即賦予相同的聲速值），結(jié)合多孔介質(zhì)幾何尺寸能夠計(jì)算出首波到達(dá)時(shí)刻理論值。為了建立首波到達(dá)時(shí)刻提取方法，首先賦予多孔介質(zhì)中各組分相同的聲速和密度值，然后保持密度不變，在1 500 ～ 4 500 m/s范圍內(nèi)設(shè)定一系列材料聲速值進(jìn)行數(shù)值求解。經(jīng)過(guò)數(shù)值求解得到樣品的接收聲波信號(hào)如圖6所示，首先從接收信號(hào)中提取波形的首個(gè)波谷到達(dá)時(shí)刻，然后計(jì)算其與首波到達(dá)時(shí)刻理論值0的差值Δ=–0。圖7所示為Δ隨的變化情況。

可以將圖7所示Δ隨的變化趨勢(shì)劃分為三個(gè)區(qū)間，兩個(gè)轉(zhuǎn)折點(diǎn)對(duì)應(yīng)的首個(gè)波谷到達(dá)時(shí)刻分別為1.256 μs和1.457 μs。針對(duì)三個(gè)區(qū)間內(nèi)Δ隨的變化規(guī)律，分別采用線(xiàn)性函數(shù)進(jìn)行擬合，擬合結(jié)果如式（5）所示。

采用水飽和多孔介質(zhì)時(shí)，接收聲波的首個(gè)波谷到達(dá)時(shí)刻t的模型計(jì)算值為1.515 μs，結(jié)合上述線(xiàn)性擬合式［式（5）］，可以得到水飽和多孔介質(zhì)的Δt為0.604 μs，因此其首波到達(dá)時(shí)刻的真實(shí)值t0為0.911 μs。已知多孔介質(zhì)樣品沿聲波傳播方向的長(zhǎng)度L為1.624 mm，因此能夠計(jì)算出水飽和多孔介質(zhì)的聲速Vp為1 782.7 m/s。在后續(xù)的模型求解數(shù)據(jù)分析中，采用相同的方法對(duì)所獲取的首波到達(dá)時(shí)刻t進(jìn)行校正，從而得到真實(shí)的首波時(shí)刻t0和多孔介質(zhì)聲速。

圖7 Δt隨首個(gè)波谷到達(dá)時(shí)刻變化

2 含水合物多孔介質(zhì)聲速特性

圖8所示為懸浮模式、接觸模式和膠結(jié)模式條件下含水合物多孔介質(zhì)接收的聲波信號(hào)。利用所接收的聲波信號(hào)計(jì)算出首個(gè)波谷到達(dá)時(shí)刻，從而得到三種微觀(guān)賦存模式下不同水合物飽和度的多孔介質(zhì)樣品的聲速。含水合物多孔介質(zhì)樣品的聲速隨水合物飽和度的變化情況如圖9所示。

分析圖8和圖9可知：

（1）三種賦存模式對(duì)應(yīng)的含水合物多孔介質(zhì)的聲速均隨水合物飽和度增大而增大。隨著水合物飽和度的增大，更多的孔隙水轉(zhuǎn)化為水合物，由于水合物比孔隙水具有更高的體積模量，因此多孔介質(zhì)整體的體積模量和剪切模量增大，從而使多孔介質(zhì)具有更高的聲速。

（2）含懸浮和接觸模式水合物多孔介質(zhì)的聲速隨水合物飽和度的增大而增大的規(guī)律基本一致，聲速增長(zhǎng)速率穩(wěn)定，與DUCHKOV等[17]采用模擬實(shí)驗(yàn)方法得到的聲速與飽和度之間的關(guān)系一致。當(dāng)水合物飽和度為22%時(shí)，兩種賦存模式下多孔介質(zhì)的聲速分別達(dá)到2 090.1 m/s和2 117.3 m/s。參見(jiàn)圖2和圖3，接觸模式可通過(guò)將懸浮模式下的水合物顆粒運(yùn)移到與顆粒接觸而得到，因此水合物飽和度變化所引起的多孔介質(zhì)彈性模量的變化差異不大。

（3）含膠結(jié)模式水合物多孔介質(zhì)中水合物飽和度由0增大到2%時(shí)，聲速出現(xiàn)大幅增長(zhǎng)現(xiàn)象，由1 782.7 m/s增長(zhǎng)為1 956.6 m/s；水合物飽和度由2%增大到22%時(shí)，聲速隨水合物飽和度的增長(zhǎng)速率趨于穩(wěn)定，飽和度為22%時(shí)聲速達(dá)到2 343.4 m/s。在相同水合物飽和度下，含膠結(jié)模式水合物多孔介質(zhì)的聲速顯著高于含懸浮模式和接觸模式水合物多孔介質(zhì)[16]。膠結(jié)模式條件下的水合物使得多孔介質(zhì)的體積模量和剪切模量大幅增長(zhǎng)，因此多孔介質(zhì)聲速大幅提升。隨著水合物飽和度進(jìn)一步增大，水合物飽和度變化成為聲速變化的主導(dǎo)因素，因此聲速增長(zhǎng)速率趨于穩(wěn)定（與懸浮和接觸模式類(lèi)似）。

圖8 接收的聲波信號(hào)（時(shí)間長(zhǎng)度為2 μs）：（a）懸浮模式；（b）接觸模式；（c）膠結(jié)模式

圖9 多孔介質(zhì)聲速隨水合物飽和度變化

采用時(shí)間平均方程[34]、伍德方程[35]、權(quán)重方程[36]、Lee改進(jìn)的Biot-Gassmann（Biot-Gassmann Theory modified by Lee, BGTL）模型[37]和等效介質(zhì)理論模型[32,38]等對(duì)多孔介質(zhì)的聲速進(jìn)行計(jì)算，計(jì)算時(shí)采用的介質(zhì)聲速、密度、剪切模量和體積模量等參數(shù)參考表1。將三種微觀(guān)賦存模式對(duì)應(yīng)的含水合物多孔介質(zhì)聲速數(shù)值模擬結(jié)果與理論模型計(jì)算結(jié)果進(jìn)行了對(duì)比，如圖10所示。其中基于等效介質(zhì)理論模型A、B和C的聲速計(jì)算結(jié)果分別與含懸浮、接觸、膠結(jié)模式水合物多孔介質(zhì)數(shù)值模擬結(jié)果進(jìn)行了比較。

圖10 聲速數(shù)值模擬與理論模型計(jì)算結(jié)果對(duì)比：（a）懸浮模式；（b）接觸模式；（c）膠結(jié)模式

分析圖10可知：

（1）對(duì)于懸浮模式（圖10a），聲速數(shù)值模擬結(jié)果低于時(shí)間平均方程計(jì)算結(jié)果，高于伍德方程計(jì)算結(jié)果；權(quán)重方程（取1.9、取0.8）計(jì)算出的聲速與數(shù)值模擬結(jié)果具有較高的一致性；BGTL模型（取0.9、取1.4）聲速計(jì)算結(jié)果與數(shù)值模擬結(jié)果基本一致；等效介質(zhì)理論模型A計(jì)算出的聲速隨水合物飽和度增長(zhǎng)速率略低于數(shù)值模擬結(jié)果。

（2）對(duì)于接觸模式（圖10b），聲速數(shù)值模擬結(jié)果低于時(shí)間平均方程計(jì)算結(jié)果，高于伍德方程計(jì)算結(jié)果；權(quán)重方程（取1.9、取1）計(jì)算出的聲速與數(shù)值模擬結(jié)果具有較高的一致性；BGTL模型（取1、取1.5）計(jì)算結(jié)果與數(shù)值模擬結(jié)果一致性較好；等效介質(zhì)理論模型B計(jì)算出的聲速隨水合物飽和度增長(zhǎng)速率略低于數(shù)值模擬結(jié)果。

（3）對(duì)于膠結(jié)模式（圖10c），聲速數(shù)值模擬結(jié)果低于時(shí)間平均方程計(jì)算結(jié)果，高于伍德方程計(jì)算結(jié)果；權(quán)重方程（取1.7、取2.9）計(jì)算出的聲速與數(shù)值模擬結(jié)果差異最??；BGTL模型（取1、取0.6）計(jì)算結(jié)果在水合物飽和度大于10%時(shí)與數(shù)值模擬結(jié)果幾乎一致，但在小于10%飽和度區(qū)間內(nèi)，BGTL模型計(jì)算結(jié)果顯著高于數(shù)值模擬結(jié)果；等效介質(zhì)理論模型C計(jì)算結(jié)果與數(shù)值模擬結(jié)果高度一致。

（4）對(duì)于三種賦存模式對(duì)應(yīng)的含水合物多孔介質(zhì)，在所研究的水合物飽和度區(qū)間內(nèi)，基于時(shí)間平均方程的多孔介質(zhì)聲速預(yù)測(cè)值偏高，而伍德方程預(yù)測(cè)值偏低[39]；權(quán)重方程和BGTL模型對(duì)預(yù)測(cè)含懸浮和接觸模式水合物多孔介質(zhì)聲速具有優(yōu)勢(shì)；等效介質(zhì)理論模型C對(duì)預(yù)測(cè)含膠結(jié)模式多孔介質(zhì)聲速具有獨(dú)特優(yōu)勢(shì)。

3 含水合物多孔介質(zhì)聲衰減特性

聲波在多孔介質(zhì)中傳播時(shí)，在聲阻抗不連續(xù)的界面處會(huì)產(chǎn)生折射與反射現(xiàn)象，從而導(dǎo)致聲波能量衰減。本節(jié)討論多孔介質(zhì)中水合物飽和度和賦存模式對(duì)聲衰減的影響。經(jīng)過(guò)數(shù)值求解，得到針對(duì)含懸浮模式、接觸模式和膠結(jié)模式水合物多孔介質(zhì)所接收到的聲波信號(hào)，如圖11所示。研究聲波信號(hào)的衰減時(shí)，將信號(hào)的振幅（以最大電壓表示）作為考察指標(biāo)?；谛盘?hào)振幅的降低定義了用于量化聲波能量衰減的參數(shù)a，a由公式（6）計(jì)算得到。

從圖11所示接收聲波信號(hào)中，提取三種賦存模式對(duì)應(yīng)的含水合物多孔介質(zhì)在不同水合物飽和度條件下的信號(hào)振幅，計(jì)算得到多孔介質(zhì)聲波衰減參數(shù)Ma。圖12展示了含水合物多孔介質(zhì)聲波衰減參數(shù)Ma隨水合物飽和度變化的情況。

圖12 聲波衰減參數(shù)Ma隨水合物飽和度變化

分析圖12可知：

（1）三種水合物賦存模式對(duì)應(yīng)的含水合物多孔介質(zhì)聲衰減參數(shù)a隨著水合物飽和度的增大都呈現(xiàn)出近似線(xiàn)性增長(zhǎng)的規(guī)律[22]。與水飽和多孔介質(zhì)相比較，水合物的存在導(dǎo)致聲波在傳播過(guò)程中遇到更多不連續(xù)的聲阻抗界面，從而加劇聲波能量的衰減。采用線(xiàn)性函數(shù)對(duì)兩者關(guān)系進(jìn)行擬合，得到式（7）。

（2）含懸浮和接觸模式水合物多孔介質(zhì)的聲衰減參數(shù)a隨水合物飽和度的變化規(guī)律基本一致，這與以上聲速隨水合物飽和度變化規(guī)律類(lèi)似，分析其原因?yàn)閼腋『徒佑|模式條件下的水合物對(duì)多孔介質(zhì)內(nèi)部聲阻抗結(jié)構(gòu)影響的差異性不顯著。

（3）在相同的水合物飽和度條件下，含膠結(jié)模式水合物多孔介質(zhì)的聲衰減參數(shù)a顯著低于懸浮和接觸模式。與懸浮和接觸模式相比較，膠結(jié)模式條件下多孔介質(zhì)內(nèi)部組分界面出現(xiàn)差異。主要差異在于：前兩種模式中，水直接與石英砂接觸，而在膠結(jié)模式中水合物與石英砂接觸，同時(shí)水與石英砂不接觸?？紤]到水合物與石英砂之間的聲阻抗差異小于水與石英砂，因此膠結(jié)模式條件下多孔介質(zhì)聲衰減參數(shù)a較低，即衰減較小。

4 結(jié) 論

基于電?力?聲多物理場(chǎng)耦合建模方法，針對(duì)懸浮、接觸和膠結(jié)三種水合物微觀(guān)賦存模式，分別建立了含水合物多孔介質(zhì)的有限元數(shù)值模型，分析了微觀(guān)賦存模式和水合物飽和度對(duì)多孔介質(zhì)聲速及衰減的影響規(guī)律，對(duì)比了聲速數(shù)值模擬與理論模型計(jì)算結(jié)果，建立了聲波衰減參數(shù)與水合物飽和度之間的關(guān)系式。得到如下結(jié)論：

（1）對(duì)于三種水合物賦存模式，由于水合物相比孔隙水具有更高的彈性模量，多孔介質(zhì)的聲速隨著水合物飽和度的增大而增大；水合物的存在導(dǎo)致聲波在傳播過(guò)程中遇到更多不連續(xù)的聲阻抗界面，聲衰減隨著水合物飽和度的增大而近似線(xiàn)性增大。

（2）懸浮和接觸賦存模式條件下，水合物飽和度對(duì)多孔介質(zhì)的聲速和衰減影響規(guī)律基本一致；對(duì)于相同的水合物飽和度，膠結(jié)模式條件下含水合物多孔介質(zhì)具有更高的聲速和更小的聲衰減。

（3）通過(guò)合理選擇參數(shù)值，利用權(quán)重方程與BGTL模型估算的含懸浮和接觸模式水合物多孔介質(zhì)的聲速較為準(zhǔn)確；通過(guò)等效介質(zhì)理論模型C計(jì)算的含膠結(jié)模式水合物多孔介質(zhì)的聲速更為準(zhǔn)確。

[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] YE J L, QIN X W, XIE W W, et al. The second natural gas hydrate production test in the South China Sea[J]. China geology, 2020, 3(2): 197-209. DOI: 10.31035/ cg2020043.

[3] LIU J W, LI X S. Recent advances on natural gas hydrate exploration and development in the South China Sea[J]. Energy & fuels, 2021, 35(9): 7528-7552. DOI: 10.1021/ acs.energyfuels.1c00494.

[4] 蘇正, 劉麗華. 南海北部陸坡天然氣水合物成藏特征研究進(jìn)展分析[J]. 新能源進(jìn)展, 2020, 8(1): 35-41. DOI: 10.3969/j.issn.2095-560X.2020.01.006.

[5] 周守為, 李清平, 呂鑫, 等. 天然氣水合物開(kāi)發(fā)研究方向的思考與建議[J]. 中國(guó)海上油氣, 2019, 31(4): 1-8. DOI: 10.11935/j.issn.1673-1506.2019.04.001.

[6] LI J F, YE J L, QIN X W, et al. The first offshore natural gas hydrate production test in South China Sea[J]. China geology, 2018, 1(1): 5-16. DOI: 10.31035/cg2018003.

[7] KOH D Y, KANG H, LEE J W, et al. Energy-efficient natural gas hydrate production using gas exchange[J]. Applied energy, 2016, 162: 114-130. DOI: 10.1016/j. apenergy.2015.10.082.

[8] 鐘廣法, 張迪, 趙巒嘯. 大洋鉆探天然氣水合物儲(chǔ)層測(cè)井評(píng)價(jià)研究進(jìn)展[J]. 天然氣工業(yè), 2020, 40(8): 25-44. DOI: 10.3787/j.issn.1000-0976.2020.08.002.

[9] SHANKAR U, RIEDEL M. Gas hydrate saturation in the Krishna-Godavari basin from P-wave velocity and electrical resistivity logs[J]. Marine and petroleum geology, 2011, 28(10): 1768-1778. DOI: 10.1016/j. marpetgeo.2010.09.008.

[10] KOU X, LI X S, WANG Y, et al. Effects of gas occurrence pattern on distribution and morphology characteristics of gas hydrates in porous media[J]. Energy, 2021, 226: 120401. DOI: 10.1016/j.energy.2021.120401.

[11] 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.

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

[13] 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.

[14] WINTERS W J, WAITE W F, MASON D H, et al. Methane gas hydrate effect on sediment acoustic and strength properties[J]. Journal of petroleum science and engineering, 2007, 56(1/3): 127-135. DOI: 10.1016/j.petrol. 2006.02.003.

[15] YUN T S, FRANCISCA F M, SANTAMARINA J C, et al. Compressional and shear wave velocities in uncemented sediment containing gas hydrate[J]. Geophysical research letters, 2005, 32(10): L10609. DOI: 10.1029/2005GL022607.

[16] HU G W, YE Y G, ZHANG J, et al. Acoustic properties of hydrate-bearing unconsolidated sediments measured by the bender element technique[J]. Chinese journal of geophysics, 2012, 55(6): 635-647. DOI: 10.1002/cjg2. 1758.

[17] DUCHKOV A D, DUCHKOV A A, PERMYAKOV M E, et al. Acoustic properties of hydrate-bearing sand samples: laboratory measurements (setup, methods, and results)[J]. Russian geology geophysics, 2017, 58(6): 727-737. DOI: 10.1016/j.rgg.2016.09.029.

[18] PRIEST J A, BEST A I, CLAYTON C R I. A laboratory investigation into the seismic velocities of methane gas hydrate-bearing sand[J]. Journal of geophysical research: solid earth, 2005, 110(B4): B04102. DOI: 10.1029/ 2004JB003259.

[19] PRIEST J A, REES E V L, CLAYTON C R I. Influence of gas hydrate morphology on the seismic velocities of sands[J]. Journal of geophysical research: solid earth, 2009, 114(B11): B11205. DOI: 10.1029/2009JB006284.

[20] REN S R, LIU Y J, LIU Y X, et al. Acoustic velocity and electrical resistance of hydrate bearing sediments[J]. Journal of petroleum science and engineering, 2010, 70(1/2): 52-56. DOI: 10.1016/j.petrol.2009.09.001.

[21] SAHOO S K, NORTH L J, MARíN-MORENO H, et al. Laboratory observations of frequency-dependent ultrasonic P-wave velocity and attenuation during methane hydrate formation in Berea sandstone[J]. Geophysical journal international, 2019, 219(1): 713-723. DOI: 10.1093/gji/ggz311.

[22] LIU Z C, NING F L, HU G W, et al. Characterization of seismic wave velocity and attenuation and interpretation of tetrahydrofuran hydrate-bearing sand using resonant column testing[J]. Marine and petroleum geology, 2020, 122: 104620. DOI: 10.1016/j.marpetgeo.2020.104620.

[23] 趙建鵬, 陳惠, 李寧, 等. 三維數(shù)字巖心技術(shù)巖石物理應(yīng)用研究進(jìn)展[J]. 地球物理學(xué)進(jìn)展, 2020, 35(3): 1099-1108. DOI: 10.6038/pg2020DD0486.

[24] SAENGER E H, BOHLEN T. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid[J]. Geophysics, 2004, 69(2): 583-591. DOI: 10.1190/1.1707078.

[25] DEL VALLE-GARCíA R, SáNCHEZ-SESMA F J. Rayleigh waves modeling using an elastic lattice model[J]. Geophysical research letters, 2003, 30(16): 1866. DOI: 10.1029/2003GL017600.

[26] MAKARYNSKA D, GUREVICH B, CIZ R, et al. Finite element modelling of the effective elastic properties of partially saturated rocks[J]. Computers & geosciences, 2008, 34(6): 647-657. DOI: 10.1016/j.cageo.2007.06.009.

[27] ZHU X, MCMECHAN G A. Numerical simulation of seismic responses of poroelastic reservoirs using Biot theory[J]. Geophysics, 1991, 56(3): 328-339. DOI: 10.1190/1.1443047.

[28] 孫建孟, 閆國(guó)亮, 姜黎明, 等. 基于數(shù)字巖心研究流體性質(zhì)對(duì)裂縫性低滲透儲(chǔ)層彈性參數(shù)的影響規(guī)律[J]. 中國(guó)石油大學(xué)學(xué)報(bào)(自然科學(xué)版), 2014, 38(3): 39-44. DOI: 10.3969/j.issn.1673-5005.2014.03.006.

[29] 劉向君, 朱洪林, 梁利喜. 基于微CT技術(shù)的砂巖數(shù)字巖石物理實(shí)驗(yàn)[J]. 地球物理學(xué)報(bào), 2014, 57(4): 1133- 1140. DOI: 10.6038/cjg20140411.

[30] YANG J J, HE B S, ZHANG J Z. Multicomponent seismic forward modeling of gas hydrates beneath the seafloor[J]. Applied geophysics, 2014, 11(4): 418-428. DOI: 10.1007/s11770-014-0465-x.

[31] 卜慶濤, 劉圣彪, 胡高偉, 等. 含水合物沉積物聲學(xué)特性——實(shí)驗(yàn)?zāi)M與數(shù)值模擬的對(duì)比分析[J]. 海洋地質(zhì)前沿, 2020, 36(9): 56-67. DOI: 10.16028/j.1009-2722. 2020.079.

[32] ECKER C. Seismic characterization of methane hydrate structures[M]. Stanford University?ProQuest Dissertations Publishing, 1998, 9837195.

[33] LEE M W, COLLETT T S. Pore- and fracture-filling gas hydrate reservoirs in the gulf of Mexico gas hydrate joint industry project leg II green canyon 955 H well[J]. Marine and petroleum geology, 2012, 34(1): 62-71. DOI: 10.1016/j.marpetgeo.2011.08.002.

[34] WYLLIE M R J, GREGORY A R, GARDNER G H F. An experimental investigation of factors affecting elastic wave velocities in porous media[J]. Geophysics, 1958, 23(3): 459-493. DOI: 10.1190/1.1438493.

[35] PEARSON C F, HALLECK P M, MCGUIRE P L, et al. Natural gas hydrate deposits: a review of in situ properties[J]. The journal of physical chemistry, 1983, 87(21): 4180-4185. DOI: 10.1021/j100244a041.

[36] WOOD W T, STOFFA P L, SHIPLEY T H. Quantitative detection of methane hydrate through high-resolution seismic velocity analysis[J]. Journal of geophysical research: solid earth, 1994, 99(B5): 9681-9695. DOI: 10.1029/94jb00238.

[37] LEE M W. Biot–Gassmann theory for velocities of gas hydrate-bearing sediments[J]. Geophysics, 2002, 67(6): 1711-1719. DOI: 10.1190/1.1527072.

[38] HELGERUD M B, DVORKIN J, NUR A, et al. Elastic-wave velocity in marine sediments with gas hydrates: Effective medium modeling[J]. Geophysical research letters, 1999, 26(13): 2021-2024. DOI: 10.1029/1999GL900421.

[39] 雷亞妮, 吳時(shí)國(guó), 孫金, 等. 海洋天然氣水合物綜合地球物理識(shí)別方法評(píng)述[J]. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2022, 53(3): 864-878. DOI: 10.11817/j.issn.1672-7207. 2022.03.009.

Sound Velocity and Attenuation Characteristics of Hydrate-Bearing Porous Media Based on Numerical Model of Electrical-Mechanical-Acoustic Multi-Physics-Field Coupling

XU Yuan1, ZHANG Huan-huan1, XING Lan-chang1, WEI Wei2, HAN Wei-feng2, WEI Zhou-tuo3, DENG Shao-gui3

(1. College of Information and Control 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 still technical difficulties in accurately controlling and evaluating the micro-distribution mode and saturation of hydrate in physical simulation experiments. Limitations exist in the experimental technologies for investigating the acoustic characteristics of hydrate-bearing sediments and establishing interpretation models of reservoir parameters. Three-dimensional numerical models were established for hydrate-bearing porous media based on digital rock physics technology. For the three kinds of hydrate micro-distribution modes (suspension, contact, and cementation), finite-element models were established individually based on the method of electrical-mechanical-acoustic multi-physics-field coupling. The effects of micro-distribution mode and hydrate saturation on sound velocity and attenuation of porous media were examined. The results of sound velocity from the numerical and theoretical models were compared. The relationship between acoustic attenuation parameters and hydrate saturation was established. It was demonstrated that: (1) for the three micro-distribution modes, the sound velocity of porous media increased with the increase of hydrate saturation due to the higher elastic modulus of hydrate compared with pore water; the existence of hydrate caused the sound wave to encounter more discontinuous acoustic impedance interfaces in the process of propagation, thus the sound attenuation increased approximately linearly with the increase of hydrate saturation; (2) under the conditions of suspension and contact modes, the influence of hydrate saturation on sound velocity and attenuation in porous media was quite similar; under the condition of cementation mode, the sound velocity was higher and the sound attenuation was lower for the porous media with the same hydrate saturation; (3) based on a reasonable selection of parameters, the sound velocity of hydrate porous media with suspension and contact modes of hydrate estimated by the weighted equation and Biot-Gassmann Theory modified by Lee (BGTL) model was more accurate; the sound velocity of porous media containing cemented hydrate calculated by the effective medium theoretical model C was more accurate. This study may provide a numerical modeling method for obtaining the acoustic characteristics of hydrate-bearing sediments under complex conditions, and provide a theoretical support for the fine evaluation of hydrate reservoirs based on acoustic logging data.

natural gas hydrate; digital rock physics; micro-distribution mode; hydrate saturation; sound velocity; sound attenuation

2095-560X（2022）05-0400-10

TK01；P631；P315

A

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

2022-04-29

2022-07-25

中石油重大科技項(xiàng)目（ZD2019-184-001）；中央高?；究蒲袠I(yè)務(wù)費(fèi)專(zhuān)項(xiàng)資金項(xiàng)目（20CX05005A）；中國(guó)石油科技創(chuàng)新基金項(xiàng)目（2018D-5007-0214）；山東省自然科學(xué)基金項(xiàng)目（ZR2019MEE095）

邢蘭昌，E-mail：xinglc@upc.edu.cn

徐 源（1998-），男，碩士研究生，主要從事天然氣水合物儲(chǔ)層物性參數(shù)檢測(cè)技術(shù)研究。

張歡歡（1997-），女，碩士研究生，主要從事天然氣水合物儲(chǔ)層物性參數(shù)檢測(cè)技術(shù)研究。

邢蘭昌（1983-），男，博士，副教授，主要從事天然氣水合物、CO2等多相流相關(guān)檢測(cè)理論與方法、多物理場(chǎng)耦合數(shù)值模擬方法、智能感知與檢測(cè)技術(shù)、計(jì)算機(jī)測(cè)控系統(tǒng)研究工作。