0 引言
天然气水合物在储层中具有多种微观赋存模式,不同的微观赋存模式会对沉积物的声学、电学、力学、热学等物性参数产生不同程度的影响[10,11,12,13]。为了准确评价储层中水合物的饱和度,在研究水合物饱和度与沉积物中声速和衰减之间的定量关系时,必须考虑水合物微观赋存模式对沉积物声速及衰减的影响规律。国内外有许多学者采用室内模拟实验方法研究了含水合物（模拟）沉积物的声学特性[14,15,16,17,18,19,20]。文献[14,15,16,17] 研究了水合物的存在对不同类型沉积物（固结、松散）声速的影响。PRIEST等[18,19]研究了不同体系（过量气、定量气、溶解气）中水合物的生成对沉积物声速的影响规律。SAHOO等[21]和LIU等[22]分别研究了声波频率和水合物饱和度对沉积物声波衰减特性的影响规律。虽然研究者提出了制备具有不同水合物饱和度、不同微观分布模式的含水合物沉积物样品的实验方法,但是对于所测试的整个样品而言,实验中所得到的饱和度、微观分布模式是否准确、是否具有重复性、是否方便控制,仍然是值得进一步探讨的问题。而数值模拟方法能够准确、方便地控制水合物饱和度、微观赋存模式等影响因素,同时数值模拟方法具有重复性好、时间成本低的显著优势。
数值模拟方法在岩石物理领域具有广泛的应用,逐渐形成了数字岩石物理技术。利用数字岩石物理技术不仅能够实现岩石孔隙尺度建模,也能够研究微观因素对岩石声速、电阻率、渗透率等物理参数的影响[23]。针对岩石声学特性的数值模拟方法主要包括旋转交错网格有限差分方法、声格子方法和有限元方法。SAENGER等[24]采用旋转交错网格有限差分方法模拟出声波在各向异性、黏弹性介质中的传播过程。DEL VALLE-GARCÍA等[25]利用声格子方法研究了声波在多孔介质中的传播过程和规律。MAKARYNSKA等[26]利用有限元方法研究了不同含水饱和度岩石的声学特性。ZHU等[27]利用有限差分方法实现了油气储层介质弹性波场的数值模拟。孙建孟等[28]、刘向君等[29]采用有限元方法计算了数字岩心的弹性参数。YANG等[30]采用旋转交错网格有限差分方法模拟了地震波在海底天然气水合物沉积地层中的传播。卜庆涛等[31]基于实验条件对水合物模拟储层开展了声学特性数值模拟,发现水合物储层声速的实验与数值模拟结果对应关系良好。目前,基于数字岩石物理技术研究微观赋存模式、水合物饱和度等因素对含水合物沉积物声学特性的影响规律鲜有报道。
本文采用基于有限元的数字岩石物理技术,针对悬浮、接触和胶结三种水合物微观赋存模式分别建立多孔介质的三维电-力-声多物理场耦合模型,研究了微观赋存模式和水合物饱和度对多孔介质声速和衰减的影响规律。研究结果可为含水合物多孔介质声学特性数值模拟提供建模方法,为声波测井在水合物储层精细评价中的应用提供理论依据和模型支撑。
1 含水合物多孔介质多物理场有限元建模
1.1 几何结构与材料
采用多物理场耦合计算平台（COMSOL Multiphysics）来建立含水合物多孔介质的有限元数值模型。数值模型的几何结构包含压电晶片和多孔介质两部分。如图1所示,两块压电晶片位于多孔介质两侧,分别用于模拟发射和接收声波。两块压电晶片的几何尺寸和材料相同,晶片厚度h为1 mm,边长L为1.624 mm。水饱和多孔介质为正方体,边长为1.624 mm。几何结构包含64个大小相同的球形石英砂颗粒,半径r为0.21 mm,石英砂颗粒之间规则相切排布。多孔介质的孔隙度为相互连通的孔隙体积与总体积的比值,本文所构建的多孔介质几何结构的孔隙体积约为1.8 mm3,孔隙度φ为0.42。
Fig. 1 Geometry of numerical model of water-saturated porous media图1 水饱和多孔介质数值模型几何结构 |
Fig. 2 Three micro-distribution modes of hydrate in sediments （modified from reference [32]）图2 沉积物中水合物的三种微观赋存模式（改自文献[32]） |
Fig. 3 Three-dimensional cubic geometries of hydrate-bearing porous media with different micro-distribution modes: (a) suspending mode; (b) contact mode; (c) cementation mode图3 含不同微观赋存模式水合物的多孔介质模型三维立体结构：（a）悬浮模式;（b）接触模式;（c）胶结模式 |
为得到不同的水合物饱和度,对于悬浮和接触模式,保持代表水合物的球体个数不变,改变球体的直径;对于胶结模式,改变石英砂颗粒表面水合物层的厚度。受到三维模型几何结构的限制（保持水合物为单一球形颗粒）,悬浮模式条件下模型所模拟的水合物饱和度存在上限,当水合物与石英砂相切时水合物饱和度为23% ~ 24%。考虑到地层中水合物饱和度过高时,水合物的赋存模式倾向于其他模式（如接触模式）,因此选取含2%、4%、6%、8%、10%、12%、14%、16%、18%、20%、22%等11个水合物饱和度的多孔介质进行数值求解。悬浮与接触模式条件下,以上11种水合物饱和度对应的水合物颗粒半径r的范围为68.47 ~ 152.14 μm。胶结模式条件下对应的水合物层厚度范围为1.16 ~ 13.18 μm。
在分别定义石英砂、孔隙水、水合物和压电晶片等材料的基础上,根据表1设定多孔介质中各组分的声速和密度等参数。压电晶片材料为锆钛酸铅（PZT-5H）,密度ρ为7 500 kg/m3。
1.2 物理场控制方程与边界条件
采用静电、固体力学、压力声学多物理场耦合方法实现数值建模方案。采用声-结构边界耦合压力声学和固体力学两个物理场,基于压电效应实现静电和固体力学两个物理场的耦合。
静电物理场中控制方程为高斯方程,如式（1）所示。
$\left\{ \begin{align} & \nabla \cdot D={{\rho }_{\text{v}}} \\ & E=-\nabla U \\ \end{align} \right.$ (1)
式中：D为电位移,C/m2;E为电场强度,V/m;U为电势,V;ρv为自由电荷密度,C/m3。
固体力学物理场控制方程为运动方程和压电材料本构方程,计算得到位移、应力和应变值,方程如式（2）所示。
$\left\{ \begin{align} & \rho {{\mathbf{u}}_{\text{tt}}}=\nabla \cdot \mathbf{T}+{{\mathbf{F}}_{\text{V}}} \\ & {{\mathbf{u}}_{\text{tt}}}=\frac{{{\partial }^{2}}\mathbf{u}}{\partial {{t}^{2}}} \\ & \mathbf{T}={{\mathbf{C}}^{\text{E}}}\cdot \mathbf{S}-\mathbf{e}\cdot \mathbf{E} \\ & \mathbf{D}=\mathbf{e}\cdot \mathbf{S}+{{\varepsilon }_{0}}{{\mathbf{K}}^{\text{S}}}\cdot \mathbf{E} \\ \end{align} \right.$ (2)
式中：ρ为压电晶片材料的密度,kg/m3;u为位移场,m;FV为体积力,N/m3;utt为压电晶片振动加速度,m/s2;T为应力,Pa;S为应变,无量纲;CE与e分别为压电晶片材料的弹性矩阵（Pa）和耦合矩阵（C/m2）;ɛ0为真空介电常数,8.854 × 10-12 F/m;KS为压电晶片材料的相对介电常数,无量纲。
压力声学物理场控制方程为时域波动方程,如式（3）所示。
$\nabla \cdot \left( \frac{-\nabla {{p}_{\text{t}}}}{{{\rho }_{0}}} \right)\text{+}\frac{1}{{{\rho }_{0}}{{c}^{2}}}\frac{{{\partial }^{2}}{{p}_{\text{t}}}}{\partial {{t}^{2}}}=0$ (3)
式中：ρ0为声波传播介质的密度,kg/m3;c为介质的声速,m/s;pt为总声压,Pa。
以上三个物理场的边界条件设定如下：①对于静电物理场,将用作发射声波的压电晶片内端面定义为电势边界,电势U为矩形波调制谐波后的电压脉冲,脉冲幅值U0为100 V,频率为f;将用作接收声波的压电晶片内端面定义为悬浮电位边界（floating potential）;将两块压电晶片的外端面定义为接地边界。②对于固体力学物理场,将位移场初始值设为0,将两块压电晶片外端面定义为辊支承边界,防止其沿厚度方向移动。③对于压力声学物理场,为保证声波在抵达多孔介质几何外边界时沿法向产生的反射最少,将多孔介质的外边界定义为平面辐射波边界。
声-结构（本模型中的多孔介质-压电晶片）边界条件如式（4）所示。
$\left\{ \begin{align} & -\mathbf{n}\cdot \left( -\frac{\nabla {{p}_{\text{t}}}}{\rho } \right)=-\mathbf{n}\cdot {{\mathbf{u}}_{\text{tt}}} \\ & {{\mathbf{F}}_{\text{A}}}={{p}_{\text{t}}}\mathbf{n} \\ \end{align} \right.$ (4)
式中：n为边界上的单位法向量,无量纲;FA为边界上单位面积所受的力,Pa。
1.3 网格划分与求解器
使用自由四面体网格对压电晶片几何结构进行网格划分。采用三种数量级别的网格、网格尺寸及数量如表2所示。利用稳态（静态）求解器对以上描述静电、固体力学物理场的微分方程进行求解。
Table 2 Sizes and counts of meshes for piezoelectric wafer表2 压电晶片部分网格尺寸及数量 |
| 网格编号 | 最大尺寸 / μm | 最小尺寸 / μm | 网格数量 |
|---|---|---|---|
| 1 | 200 | 40 | 5 724 |
| 2 | 80 | 30 | 37 971 |
| 3 | 60 | 20 | 91 274 |
Table 3 Sizes and counts of meshes for porous media表3 多孔介质部分网格尺寸及数量 |
| 网格编号 | 最大尺寸 / μm | 最小尺寸 / μm | 网格数量 |
|---|---|---|---|
| 1 | 200 | 50 | 96 672 |
| 2 | 100 | 40 | 150 653 |
| 3 | 80 | 30 | 237 535 |
| 4 | 60 | 20 | 588 993 |
| 5 | 50 | 15 | 818 336 |
| 6 | 30 | 10 | 3 061 727 |
1.4 网格选取和首波到达时刻提取方法
使用表2所示三种数量级别的网格测试了压电晶片中心频率的模型计算值对网格的依赖性,根据测试结果选取合适的网格尺寸进行后续研究。根据计算结果得到三种网格下压电晶片振动位移随激励电信号频率的变化关系,最大振动位移对应的频率即为压电晶片中心频率。结果表明,压电晶片的中心频率在三种数量级别网格下均为642.5 kHz,由于网格数量成倍增加,可认为基于这三种网格所得到的结果均为正确值。为减少计算时间,针对压电晶片选择单元数量为5 724的网格级别,即最大网格尺寸为200 μm,最小网格尺寸为40 μm。
Fig. 4 Received acoustic waves for different meshes图4 不同网格条件下接收的声波 |
Fig. 5 Variation of the arrival time for the first trough of acoustic wave with mesh counts图5 首个声波波谷到达时刻随网格数量的变化 |
Fig. 6 Received acoustic waves of different samples图6 不同样品接收的声波 |
Fig. 7 Variation of Δt with the arrival time for the first trough of acoustic waves图7 Δt随首个波谷到达时刻变化 |
可以将图7所示Δt随t的变化趋势划分为三个区间,两个转折点对应的首个波谷到达时刻分别为1.256 μs和1.457 μs。针对三个区间内Δt随t的变化规律,分别采用线性函数进行拟合,拟合结果如式（5）所示。
$\left\{ \begin{align} & \Delta t=-0.2192t+0.8759\quad \ \ 0<t\le 1.256 \text{ }\!\!\mu\!\!\text{ s} \\ & \Delta t=-0.02377t+0.6366\,\ \,\,1.256<t\le 1.457 \text{ }\!\!\mu\!\!\text{ s} \\ & \Delta t=~0.03876t+0.5456\ \quad \,\,t>1.457 \text{ }\!\!\mu\!\!\text{ s} \\ \end{align} \right.$ (5)
采用水饱和多孔介质时,接收声波的首个波谷到达时刻t的模型计算值为1.515 μs,结合上述线性拟合式[式（5）],可以得到水饱和多孔介质的Δt为0.604 μs,因此其首波到达时刻的真实值t0为0.911 μs。已知多孔介质样品沿声波传播方向的长度L为1.624 mm,因此能够计算出水饱和多孔介质的声速Vp为1 782.7 m/s。在后续的模型求解数据分析中,采用相同的方法对所获取的首波到达时刻t进行校正,从而得到真实的首波时刻t0和多孔介质声速。
2 含水合物多孔介质声速特性
Fig. 8 Received acoustic waves (time period of 2 μs): (a) suspending mode; (b) contact mode; (c) cementation mode图8 接收的声波信号（时间长度为2 μs）：（a）悬浮模式;（b）接触模式;（c）胶结模式 |
Fig. 9 Sound velocities of porous media variations with hydrate saturation图9 多孔介质声速随水合物饱和度变化 |
（1）三种赋存模式对应的含水合物多孔介质的声速均随水合物饱和度增大而增大。随着水合物饱和度的增大,更多的孔隙水转化为水合物,由于水合物比孔隙水具有更高的体积模量,因此多孔介质整体的体积模量和剪切模量增大,从而使多孔介质具有更高的声速。
（3）含胶结模式水合物多孔介质中水合物饱和度由0增大到2%时,声速出现大幅增长现象,由1 782.7 m/s增长为1 956.6 m/s;水合物饱和度由2%增大到22%时,声速随水合物饱和度的增长速率趋于稳定,饱和度为22%时声速达到2 343.4 m/s。在相同水合物饱和度下,含胶结模式水合物多孔介质的声速显著高于含悬浮模式和接触模式水合物多孔介质[16]。胶结模式条件下的水合物使得多孔介质的体积模量和剪切模量大幅增长,因此多孔介质声速大幅提升。随着水合物饱和度进一步增大,水合物饱和度变化成为声速变化的主导因素,因此声速增长速率趋于稳定（与悬浮和接触模式类似）。
Fig. 10 Results comparison of numerical simulation and theoretical models for sound velocity: (a) suspending mode; (b) contact mode; (c) cementation mode图10 声速数值模拟与理论模型计算结果对比：（a）悬浮模式;（b）接触模式;（c）胶结模式 |
分析图10可知：
（1）对于悬浮模式（图10a）,声速数值模拟结果低于时间平均方程计算结果,高于伍德方程计算结果;权重方程（W取1.9、N取0.8）计算出的声速与数值模拟结果具有较高的一致性;BGTL模型（G取0.9、n取1.4）声速计算结果与数值模拟结果基本一致;等效介质理论模型A计算出的声速随水合物饱和度增长速率略低于数值模拟结果。
（2）对于接触模式（图10b）,声速数值模拟结果低于时间平均方程计算结果,高于伍德方程计算结果;权重方程（W取1.9、N取1）计算出的声速与数值模拟结果具有较高的一致性;BGTL模型（G取1、N取1.5）计算结果与数值模拟结果一致性较好;等效介质理论模型B计算出的声速随水合物饱和度增长速率略低于数值模拟结果。
（3）对于胶结模式（图10c）,声速数值模拟结果低于时间平均方程计算结果,高于伍德方程计算结果;权重方程（W取1.7、N取2.9）计算出的声速与数值模拟结果差异最小;BGTL模型（G取1、n取0.6）计算结果在水合物饱和度大于10%时与数值模拟结果几乎一致,但在小于10%饱和度区间内,BGTL模型计算结果显著高于数值模拟结果;等效介质理论模型C计算结果与数值模拟结果高度一致。
（4）对于三种赋存模式对应的含水合物多孔介质,在所研究的水合物饱和度区间内,基于时间平均方程的多孔介质声速预测值偏高,而伍德方程预测值偏低[39];权重方程和BGTL模型对预测含悬浮和接触模式水合物多孔介质声速具有优势;等效介质理论模型C对预测含胶结模式多孔介质声速具有独特优势。
3 含水合物多孔介质声衰减特性
声波在多孔介质中传播时,在声阻抗不连续的界面处会产生折射与反射现象,从而导致声波能量衰减。本节讨论多孔介质中水合物饱和度和赋存模式对声衰减的影响。经过数值求解,得到针对含悬浮模式、接触模式和胶结模式水合物多孔介质所接收到的声波信号,如图11所示。研究声波信号的衰减时,将信号的振幅（以最大电压表示）作为考察指标。基于信号振幅的降低定义了用于量化声波能量衰减的参数Ma,Ma由公式（6）计算得到。
${{M}_{\text{a}}}=\frac{{{U}_{0}}-U}{{{U}_{0}}}$ (6)
Fig. 11 Received acoustic waves (20 μs): (a) suspending mode; (b) contact mode; (c) cementation mode图11 接收的声波信号（20 μs）：（a）悬浮模式;（b）接触模式;（c）胶结模式 |
Fig. 12 Variations of wave attenuation Ma with hydrate saturation图12 声波衰减参数Ma随水合物饱和度变化 |
分析图12可知：
（1）三种水合物赋存模式对应的含水合物多孔介质声衰减参数Ma随着水合物饱和度的增大都呈现出近似线性增长的规律[22]。与水饱和多孔介质相比较,水合物的存在导致声波在传播过程中遇到更多不连续的声阻抗界面,从而加剧声波能量的衰减。采用线性函数对两者关系进行拟合,得到式（7）。
$\left\{ \begin{align} & {{M}_{\text{a}}}=0.540{{S}_{\text{h}}}+0.675 悬浮模式 \quad \\ & {{M}_{\text{a}}}=0.539{{S}_{\text{h}}}+0.675 接触模式 \quad \\ & {{M}_{\text{a}}}=0.175{{S}_{\text{h}}}+0.675 胶结模式 \quad \\ \end{align} \right.$. (7)
（2）含悬浮和接触模式水合物多孔介质的声衰减参数Ma随水合物饱和度的变化规律基本一致,这与以上声速随水合物饱和度变化规律类似,分析其原因为悬浮和接触模式条件下的水合物对多孔介质内部声阻抗结构影响的差异性不显著。
（3）在相同的水合物饱和度条件下,含胶结模式水合物多孔介质的声衰减参数Ma显著低于悬浮和接触模式。与悬浮和接触模式相比较,胶结模式条件下多孔介质内部组分界面出现差异。主要差异在于：前两种模式中,水直接与石英砂接触,而在胶结模式中水合物与石英砂接触,同时水与石英砂不接触。考虑到水合物与石英砂之间的声阻抗差异小于水与石英砂,因此胶结模式条件下多孔介质声衰减参数Ma较低,即衰减较小。
4 结论
基于电-力-声多物理场耦合建模方法,针对悬浮、接触和胶结三种水合物微观赋存模式,分别建立了含水合物多孔介质的有限元数值模型,分析了微观赋存模式和水合物饱和度对多孔介质声速及衰减的影响规律,对比了声速数值模拟与理论模型计算结果,建立了声波衰减参数与水合物饱和度之间的关系式。得到如下结论：
（1）对于三种水合物赋存模式,由于水合物相比孔隙水具有更高的弹性模量,多孔介质的声速随着水合物饱和度的增大而增大;水合物的存在导致声波在传播过程中遇到更多不连续的声阻抗界面,声衰减随着水合物饱和度的增大而近似线性增大。
（2）悬浮和接触赋存模式条件下,水合物饱和度对多孔介质的声速和衰减影响规律基本一致;对于相同的水合物饱和度,胶结模式条件下含水合物多孔介质具有更高的声速和更小的声衰减。
（3）通过合理选择参数值,利用权重方程与BGTL模型估算的含悬浮和接触模式水合物多孔介质的声速较为准确;通过等效介质理论模型C计算的含胶结模式水合物多孔介质的声速更为准确。