全文快速搜索:   高级搜索

  中国石油大学学报(自然科学版)  2019, Vol. 43 Issue (6): 77-87  DOI:10.3969/j.issn.1673-5005.2019.06.009
0

引用本文 [复制中英文]

董长银, 闫切海, 李彦龙, 等. 天然气水合物储层颗粒级尺度微观出砂数值模拟[J]. 中国石油大学学报(自然科学版), 2019, 43(6): 77-87. DOI: 10.3969/j.issn.1673-5005.2019.06.009.
[复制中文]
DONG Changyin, YAN Qiehai, LI Yanlong, et al. Numerical simulation of sand production based on a grain scale microcosmic model for natural gas hydrate reservoir[J]. Journal of China University of Petroleum (Edition of Natural Science), 2019, 43(6): 77-87. DOI: 10.3969/j.issn.1673-5005.2019.06.009.
[复制英文]

基金项目

国家自然科学基金项目(51774307);国家重点研发计划(2017YFC0307304);国家海洋地质调查二级项目(DD20190231);中国石油重大专项(ZD2019-184)

作者简介

董长银(1976-), 男, 教授, 博士, 博士生导师, 研究方向为油气开采工程、出砂与防砂完井、复杂多相流动。E-mail:dongcy@upc.edu.cn

文章历史

收稿日期:2019-05-15
天然气水合物储层颗粒级尺度微观出砂数值模拟
董长银1, 闫切海1, 李彦龙2, 徐鸿志3, 周玉刚1, 尚校森1, 陈强2, 宋洋1     
1. 非常规油气开发教育部重点实验室(中国石油大学(华东)), 山东青岛 266580;
2. 中国地质调查局青岛海洋地质研究所, 山东青岛 266071;
3. 中国石油集团工程技术研究有限公司, 天津 345000
摘要: 构建颗粒级尺度的天然气水合物储层地层砂颗粒结构模型, 在此基础上进行天然气水合物储层颗粒级尺度微观出砂过程和形态模拟, 研究出砂机制和出砂形态及其影响因素。颗粒结构模型根据地层砂粒度分布曲线随机确定颗粒直径、圆球度系数及位置倾角, 根据测井曲线利用随机方法表征颗粒胶结强度和水合物饱和度的随机非均质性分布和粒间应力。使用水合物强度比例系数和影响系数表征水合物分解对胶结强度的影响。根据构建的颗粒剥落判别模型, 结合水合物分解模型, 建立颗粒剥落孔道延伸扩展模拟及出砂形态模拟方法。结果表明:对于胶结强度较低的高水合物饱和度储层, 出砂首先从孔壁边界上开始, 沿胶结弱面不均匀扩展形成类蚯蚓洞形态, 最终总体呈现前端类蚯蚓洞与后端连续垮塌的复合出砂形态; 水合物分解会大幅度降低储层强度, 生产过程中出砂前沿与分解前沿非常接近; 随着初期生产延续, 最大剥落出砂粒径由16 μm降低为6 μm, 剥落出砂速度呈上升趋势; 水合物饱和度及其分解对出砂形态和出砂速度有明显影响; 水合物饱和度越低, 其分解对储层总体内聚强度的影响越小, 越不容易出砂; 饱和度为0.4、0.3和0.18时的砂粒剥落速度分别约降低到平均水合物饱和度0.48对应的砂粒剥落速度的80%、69%和58%。
关键词: 天然气水合物    出砂模拟    出砂预测    颗粒级尺度    出砂形态    
Numerical simulation of sand production based on a grain scale microcosmic model for natural gas hydrate reservoir
DONG Changyin1 , YAN Qiehai1 , LI Yanlong2 , XU Hongzhi3 , ZHOU Yugang1 , SHANG Xiaosen1 , CHEN Qiang2 , SONG Yang1     
1. Key Laboratory of Unconventional Oil & Gas Development(China University of Petroleum (East China)), Ministry of Education, Qingdao 266580, China;
2. Qingdao Institute of Marine Geology, China Geological Survey, Qingdao 266071, China;
3. CNPC Engineering Technology Research Company Limited, Tianjin 345000, China
Abstract: A sand particle structure model was built in terms of grain size scales for natural gas hydrate reservoir, and the microcosmic sand production process and morphology of gas hydrate reservoir were simulated, in which the sand production mechanism, sand production morphology and the influencing factors were studied. In the particle structure model, the particle diameter, sphericity coefficient and position inclination angle were randomly determined according to the particle size distribution of the hydrate formation sands. Particle cementation strength, hydrate saturation, inter-grain stress and their random heterogeneous distribution were characterized according to the well logging curve using a stochastic method. A hydrate strength ratio coefficient and its influence coefficient were used to characterize the effect of hydrate decomposition on the bonding strength of sands. Based on the constructed particle exfoliation discriminant model, and combined with a hydrate decomposition model. Numerical simulations of particle exfoliation channel expansion and sand production morphology were conducted. The results show that, for reservoirs with high hydrate saturation and low cementation strength, sand production can start from the boundary of the pore wall, then it can expand unevenly along the weak cementation surface to form a wormhole like path, and finally it can form a composite sand production morphology, in which the front end looks like a wormhole and the back end will continuously collapses. The decomposition of hydrate can greatly reduce the strength of the reservoir, and the sand production front and the decomposition front become very close during the production process. With the continuation of initial production, the maximum exfoliated sand particle size can be reduced from 16 μm to 6 μm, and the speed of exfoliated sand production showed an upward trend. Hydrate saturation and its decomposition have a significant impact on the size distribution of the sands produced and sand production rate. The lower the hydrate saturation is, the less impact of its decomposition has on the overall cohesion of the reservoir, and the less possible for sand out. Compared with an average hydrate saturation of 0.48, for hydrate saturation of 0.4, 0.3 and 0.18, the exfoliation rate of sand particles can be decreased to about 80%, 69% and 58% of average value, respectively.
Keywords: natural gas hydrate    sand production simulation    sand production prediction    grain size scale    sand production morphology    

初步研究与试采实践表明, 出砂是制约天然气水合物储层高效安全开采的关键问题之一[1-3]。对于未来科学合理的水合物储层防控砂决策, 准确的出砂模拟与出砂规律预测是关键基础。区别于传统的石油与天然气储层出砂特点, 天然气水合物储层出砂受水合物相变分解影响较大, 使其出砂模拟更加复杂, 出砂规律预测也更加困难[4-7]。传统石油与天然气储层开采过程中的出砂预测主要基于储层岩石宏观破坏和出砂临界条件, 或基于连续介质和流体力学理论[8-11], 流体对储层物性的影响几乎可以忽略。关于复杂相变条件下的天然气水合物开采中的出砂模拟和预测研究目前处于初步阶段, 其研究方法主要有两类:一类是基于传统岩石力学的宏观破坏机制, 本构力学模型涉及线-弹性模型、弹-塑性模型、连续介质模型等[4-7],初步揭示了水合物储层的出砂规律, 但难以从微观上阐明不同胶结类型和程度水合物储层岩石颗粒剥离和运移动态过程的差异[12-16];第二类方法主要使用商业软件进行数值模拟计算, 其基本原理是建立温度-水力-应力-化学耦合出砂模拟模型, 使用TOUGH+HYDRATE、FLAC3D、PFC及ABAQUS等大型商业软件, 建立软件之间的接口, 进行求解或二次开发, 实现天然气水合物储层开采过程中的生产动态模拟和出砂规律预测[17-20], 此类方法考虑的影响因素全面, 但其模拟计算过程复杂并且需要借助于庞大的第三方模拟软件系统, 使模拟方法和手段不便于工程灵活应用。笔者考虑储层岩石物性、水合物饱和度的宏观分布规律及地层砂粒度特征, 初步构建颗粒级尺度的天然气水合物储层地层砂颗粒结构模型, 进行天然气水合物储层颗粒级尺度微观出砂过程和形态模拟, 研究出砂形态及其影响因素, 为水合物储层的出砂规律预测探索一条新的途径。

1 天然气水合物储层颗粒级尺度建模 1.1 基本颗粒网格模型

天然气水合物储层颗粒结构建模的目的是构建储层颗粒级尺度的微观结构, 作为降压开采条件下颗粒出砂模拟的基础和依据。为了使建立的模型能够更加接近真实储层的微观结构, 微观结构模型的构建及物性表征需要考虑地层砂粒径分布、砂粒形状的不规则性、储层胶结类型、储层岩石强度分布、储层地应力分布、泥质含量、水合物赋存状态、水合物的饱和度分布。

构建以颗粒为对象单元的模型(particle-objective model, POM模型)和以网格为对象单元的模型(grid-objective model, GOM模型)进行二维水合物储层颗粒骨架结构建模。POM模型如图 1(a)所示, 使用颗粒位置(x, y)和半径r描述颗粒的基本位置和尺寸属性。POM模型可以模拟颗粒的自然沉积过程, 颗粒不规则排列形成非均匀网格, 其生成可以采用类似文献[22]中的方法。该模型计算复杂, 但与实际情况比较相符。GOM模型则以均匀网格为对象单元(图 1(b)), 网格划分为均匀正方形, 每个网格中填充一个地层砂砂粒; 网格尺寸由地层砂粒径确定。该方法相对简单且计算速度快。

图 1 水合物储层二维颗粒结构基本模型示意图 Fig.1 Sketch map of basic 2-D particle distribution model of gas-hydrate reservoir

水合物储层的微观出砂过程受颗粒尺寸、形状和位置的影响, 在颗粒级尺度颗粒结构建模中需要确定每一个网格中颗粒的粒径和形状等特征参数。表达储层地层砂粒径分布特征的主要资料是地层砂粒度分布曲线。利用随机选取方法确定每个网格的颗粒粒径, 但总体粒径分布规律严格符合地层砂粒度分布筛析曲线。由于地层砂筛析曲线表达的是累积质量百分比与粒径的关系, 使用时须将某一粒径范围内的颗粒质量转化为颗粒数量:

$ {N_i} = \frac{{{M_{\rm{g}}}{W_i}}}{{{\rho _{\rm{g}}}\left( {\frac{4}{3}{\rm{ \mathsf{ π} }}d_i^3} \right)}}. $ (1)

式中, Mg为砂样总质量, kg; ρg为地层砂密度, kg/m3; di为第i组分地层砂粒径, m; Wi为第i组分砂粒径所占质量分数; Ni为第i组分地层砂的颗粒数目。

实际地层砂颗粒的形状系数各异, 为了简化模型并降低计算量, 预设8种颗粒圆球度系数a在0.65~1.0均匀分布。建模时以随机分布的方式刻画每一个网格对象的颗粒圆球度, 将每一个网格颗粒的圆球度在8种圆球度中随机选取确定。

1.2 储层颗粒间强度及水合物饱和度非均质性表征 1.2.1 水合物储层颗粒间微观胶结强度模型

图 2(a)~(d)为天然气水合物储层颗粒微观胶结的4种物理模式, GR指岩石颗粒,GH指水合物颗粒。当多孔介质中的水合物饱和度超过25%时, 胶结方式趋向于骨架颗粒支撑模式和颗粒包裹模式[24]。天然气水合物储层中, 分解前的水合物本身在颗粒间起胶结作用。典型研究表明, 降压法开采时水合物储层岩石强度是其分解后强度的3~9倍, 水合物分解后使水合物储层胶结强度大幅度降低[25]。在进行微观出砂模拟时, 须考虑水合物的分解对储层颗粒间微观胶结强度的影响作用。

图 2 天然气水合物颗粒胶结模式及微观粒间强度模型 Fig.2 Sand cementation mode of gas hydrate reservoir and micro strength model

为了定量表征水合物分解对颗粒间胶结强度的影响, 构建如图 2(e)所示的天然气水合物储层颗粒间微观胶结强度模型。天然气水合物储层颗粒间微观胶结强度由两部分组成:一是原始水合物胶结强度, 其对总胶结强度的贡献比例系数用α表示; 二是泥质胶结物胶结强度(类似于传统石油天然气储层的胶结强度), 其对总胶结强度的贡献比例系数为(1-α)。

图 2(e)所示的模型中, 当天然气水合物分解后, 水合物胶结强度消失, 其强度比例系数降低为0。同时, 为了表征水合物分解成水后对原始胶结强度的影响, 定义系数β为水合物分解对原始胶结物胶结强度的降低系数(简称水合物分解强度降低系数)。根据上述定义和假设, 有如下关系:

$ {S_0} = {S_{\rm{h}}}\alpha + {S_{\rm{c}}}\left( {1 - \alpha } \right), $ (2)
$ {S_1} = {S_{\rm{c}}}\left( {1 - \alpha } \right)\beta . $ (3)

式中, S0为颗粒间原始总胶结强度, MPa; Sh为颗粒间原始水合物胶结强度, MPa; Sc为颗粒间原始胶结物胶结强度, MPa; S1为水合物分解后颗粒间剩余强度, MPa。

原始水合物胶结强度占总胶结强度的比例系数(简称水合物强度系数)α与水合物储层孔隙度、水合物饱和度和泥质含量有关, 也与水合物的微观赋存状态有关。水合物储层的孔隙中, 泥质含量越高, 水合物饱和度越低, 则水合物强度系数α越小。简化起见, 暂忽略系数α与水合物赋存状态的关系, 其基本关系表达式为

$ \alpha = \frac{{{S_{{\rm{CH}}}}\varphi }}{{{S_{{\rm{CH}}}}\varphi + \left( {1 - \varphi } \right){S_{{\rm{VCL}}}}}}. $ (4)

式中, SGH为天然气水合物饱和度; SVCL为泥质含量; φ为孔隙度。

上述模型描述了天然气水合物的存在(由水合物饱和度表征)对储层颗粒间胶结强度变化的影响。系数αβ均为0~1, 具体数值由水合物储层物性确定; 当水合物饱和度为0时, 水合物强度系数α为0, 则储层胶结强度仅由胶结物贡献, 即与传统石油与天然气储层相同。当水合物分解强度降低系数β为0时, 表示水合物分解将使颗粒间胶结强度完全丧失; 当β为1时, 则表示水合物分解影响水合物本身的胶结强度, 而对胶结物胶结强度无影响。实际应用时, 系数α由式(4)确定, 系数β由经验确定。上述水合物储层颗粒间微观胶结强度模型可以直观描述水合物储层的颗粒间胶结强度的原始状态及其随水合物分解的变化规律。

1.2.2 非均质性表征方法

图 1所示的水合物储层颗粒级微观结构模型中, 需要确定每个颗粒的胶结强度和水合物饱和度。对于天然气水合物储层的垂直井, 沿井筒垂直方向储层的强度和水合物分布规律可由测井资料直接或解释获得[18], 其结果可以直接表征储层物性在纵向上的非均质性; 横向(平面)上的强度和饱和度分布规律及其非均质性则由正态分布随机函数确定, 但其均值取同一深度测井资料上对应的数值。

对于给定的深度h, 假设其对应图 1(b)所示的网格模型中的第(i, j)个网格的强度或水合物饱和度, 确定方法如下:

$ {K_{i{\rm{min}}}} = {K_{i0}}\left( {1 - {X_{\rm{s}}}} \right), $ (5)
$ {K_{i\max }} = {K_{i0}}\left( {1 + {X_{\rm{s}}}} \right), $ (6)
$ f\left( {i,j} \right) = \frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}} \sigma }}\exp \left[ { - \frac{{{{\left( {\frac{j}{N} - \mu } \right)}^2}}}{{2{\sigma ^2}}}} \right], $ (7)
$ {K_{ij}} = {K_{i\min }} + \frac{{{K_{i\max }} - {K_{i\min }}}}{{{f_{\max }} - {f_{\min }}}}\left[ {f\left( {i,j} \right) - {f_{\min }}} \right]. $ (8)

式中, ij为网格纵向和横向序号; N为横向网格数; Xs为储层强度或水合物饱和度非均质系数; Ki0为根据测井资料插值获得的纵向第i个网格序列的强度或水合物饱和度; KiminKimax分别为纵向第i个网格位置处横向网格强度或水合物饱和度非均质分布的最小值和最大值; f(i, j)为第(i, j)网格上的强度或水合物饱和度正态分布系数; Kij为第(i, j)网格上的强度或水合物饱和度; μ为正态分布均值参数, 建议取μ=0.5;σ为正态随机分布系数, 建议取σ=0.2。

储层强度或水合物饱和度非均质系数Xs用来表征储层强度或水合物饱和度在横向(平面)上的非均质性, 该系数决定了非均质波动范围的最低值和最高值。系数Xs的确定可以参考同层段由测井资料获得的储层强度或水合物饱和度的纵向非均质性,

$ {X_{\rm{s}}} = \frac{{{K_{\min }} + {K_{\max }}}}{{2\bar s}}. $ (9)

式中s, 为同层位测井资料获得的强度或水合物饱和度平均值; KminKmax分别为同层位测井资料获得的强度或水合物饱和度的最小值和最大值。

1.3 颗粒间微应力表征

颗粒间应力是影响颗粒剥落过程的重要因素。垂直井三维柱坐标下储层单元的宏观主应力由垂向主应力、径向应力和切向应力表示, 3个主应力可根据密度测井资料计算得到[11, 17]。对于二维模型, 垂向主应力即为纵向主应力, 径向应力表现为横向主应力。对于某特定地层砂颗粒而言, 其通过胶结物与周围颗粒连接在一起, 固定在固体骨架上, 这样的胶结连接称为连接键, 每个连接键有一定的连接强度, 可以同时传递压力和剪切力。由接触力学, 颗粒间的相互作用力可分解为通过相邻两个颗粒几何中心的法向接触力Fn和在切平面内的切向作用力Fs。颗粒间地应力的表征考虑颗粒间的相互作用及应力传递。

对于二维岩体单元, 在宏观垂向应力σv和径向应力σr作用下, 边界颗粒均匀承担外界应力, 其受力平衡方程为

$ {\sigma _{{{\rm{n}}_1}}} + {\tau _{{{\rm{s}}_1}}} + {\sigma _{{{\rm{n}}_2}}} + {\tau _{{{\rm{s}}_2}}} + \cdots + {\sigma _{{\rm{n}}n}} + {\tau _{{\rm{s}}n}} = \frac{{{\sigma _{\rm{v}}}}}{m}. $ (10)

式中, m为所研究部分边界颗粒的数目; σni为第i个颗粒的法向应力; τsi为第i个边界颗粒的切向应力。

以连接键上的法向力和剪切力为未知量X=[σn1+τs1, …, σnn+τsn]T, S=[σg1, σg2, …, σgn]T, 其中n为连接键个数, 对所有颗粒列力学平衡方程为

$ \left( {\begin{array}{*{20}{c}} {{a_{11}}}& \cdots &{{a_{1n}}}\\ \vdots & \ddots & \vdots \\ {{a_{{\rm{n1}}}}}& \cdots &{{a_{nn}}} \end{array}} \right) \times \left( {\begin{array}{*{20}{c}} {{\sigma _{{{\rm{n}}_1}}}}& + &{{\tau _{{{\rm{s}}_1}}}}\\ {}& \vdots &{}\\ {{\sigma _{{{\rm{n}}_n}}}}& + &{{\tau _{{{\rm{s}}_n}}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{\sigma _{{{\rm{g}}_1}}}}\\ \vdots \\ {{\sigma _{{{\rm{g}}_n}}}} \end{array}} \right). $ (11)

式中, aij为系数, 当两个颗粒之间相互连接时为1, 不连接时取0;σgi为边界颗粒应力, 当颗粒不为边界颗粒时σgi取0。式(11)可写成以下两种的形式:

$ \mathit{\boldsymbol{AX}} = \mathit{\boldsymbol{S}}, $ (12)
$ \mathit{\boldsymbol{X}} = {\mathit{\boldsymbol{A}}^{ - 1}}\mathit{\boldsymbol{S}}. $ (13)

A为非奇异矩阵时, A-1为A的逆矩阵; 当A为奇异矩阵时, A-1A的广义逆矩阵。使用上述模型可以得到岩体在胶结状态下颗粒间的法向应力与切向应力分布。

2 颗粒级尺度微观出砂形态模拟方法 2.1 颗粒剥落条件判别模型

以水合物储层流通通道壁面上的胶结骨架砂为例进行受力分析, 如图 3所示。假设骨架砂颗粒等效直径为dg, 在壁面上被临近的颗粒所挤压。在应力作用下骨架砂颗粒剥落时受到的阻力包括骨架砂颗粒之间胶结产生的内聚力、骨架砂颗粒因剥落运动与周围颗粒之间的摩擦力、骨架砂颗粒与剥落运动方向的砂粒之间的抗拉强度; 骨架砂颗粒剥落时的动力为各种应力在壁面上造成的应力合力, 包括地应力、孔隙压力和井底压力等。

图 3 骨架砂颗粒在流通通道壁面受力示意图 Fig.3 Force analysis of particle on rock body surface

骨架砂颗粒剥落时受到的阻力为

$ {F_{\rm{r}}} = {\rm{ \mathsf{ π} }}{\left( {\frac{{{d_{\rm{g}}}}}{2}} \right)^2}\left[ {4{S_0} + \mu \left( {2\left( {{\sigma _{\rm{ \mathsf{ θ} }}} - {p_{\rm{p}}}} \right) + 2\left( {{\sigma _{\rm{z}}} - {p_{\rm{p}}}} \right)} \right) + {T_0}} \right]. $ (14)

式中, Fr为骨架砂颗粒剥落时受到的阻力, N; dg为骨架砂颗粒的粒径, μm; S0为骨架砂颗粒间的内聚强度, MPa; μ为骨架砂颗粒间的摩擦系数; σθσz分别为孔道壁上的切向应力和轴向应力, MPa; pp为骨架砂颗粒间的孔隙压力, MPa; T0为骨架砂颗粒间的抗拉强度, MPa。

骨架砂颗粒剥落时受到的动力为

$ {F_{\rm{m}}} = {\rm{ \mathsf{ π} }}{\left( {\frac{{{d_{\rm{g}}}}}{2}} \right)^2}\left( {{\sigma _{\rm{r}}} + {p_{\rm{p}}} - {p_{\rm{w}}}} \right). $ (15)

式中, Fm为骨架砂颗粒剥落时受到的动力, N; σr为孔道壁上的径向应力, MPa; pp为骨架砂颗粒间的孔隙压力, MPa; pw为孔道中的井底压力, MPa。

在应力作用下骨架砂颗粒剥落的条件为Fm>Fr

2.2 水合物相变分解及其饱和度变化模型

水合物储层降压开采过程中, 储层可分为完全分解区、分解区和未分解区3个区域。在完全分解区内, 水合物完全分解, 只存在水和气体; 在分解区内, 除了水和气体外, 还存在正处于分解状态的水合物; 而在未分解区, 水合物还未开始分解, 处于原始饱和度状态。参考Kim-Bishnoi、Naval Goel并结合李淑霞等[24]的研究, 水合物分解速率为

$ {u_{\rm{h}}} = \frac{{{\rm{d}}{c_{\rm{h}}}}}{{{\rm{d}}t}} = {K_0}{A_{{\rm{hs}}}}{S_{\rm{w}}}{S_{\rm{h}}}{\varphi ^2}{p_{\rm{e}}}\exp \left( { - \frac{{\Delta E}}{{RT}}} \right)\left( {1 - \frac{1}{K}} \right). $ (16)

式中, μh为水合物分解速率, mol/(m3·d); ch为水合物浓度, mol/m3; K0为水合物分解表观速率因子, 1.071×1013, mol/(d·kPa·m2); Ahs为水合物球形颗粒表面积, m2; pe为水合物三相平衡压力, kPa; ΔE为分解活化能, 7330, J/mol; R为理想气体常数, 8.3144, J/(mol·K); T为温度, K; K为相平衡值。

三相平衡压力pe和相平衡值K分别为

$ {p_{\rm{e}}} = 1.15\exp \left( {49.3185 - \frac{{9459}}{T} \times {{10}^{ - 3}}} \right), $ (17)
$ K = \frac{{{k_1}}}{p}\exp \left( {\frac{{{k_2}}}{{\left( {T - 273.15} \right) - {k_3}}}} \right). $ (18)

式中, p为压力, kPa; k1k2k3由试验数据拟合得到。

在分解初始时刻, 压力、温度及各相饱和度分别为p0jT0jSw0jSCHjSg0j, 则i时刻j网格的水合物饱和度为

$ {S_{{\rm{CH}}\left( {i,j} \right)}} = {S_{{\rm{CH}}\left( {0,j} \right)}} - \frac{{\sum\limits_{i = 0}^{i - 1} {{u_{{\rm{h}}\left( {i,j} \right)}}\Delta t} }}{{{v_{{\rm{p}}j}}}}, $ (19)
$ \begin{array}{l} {u_{{\rm{h}}\left( {i,j} \right)}} = {K_0}{A_{{\rm{hs}}}}{S_{{\rm{w}}\left( {i - 1,j} \right)}}{S_{{\rm{CH}}\left( {i - 1,j} \right)}}\varphi _j^2{p_{\rm{e}}}\left( {{T_{\left( {i - 1,j} \right)}}} \right) \times \\ \exp \left( { - \frac{{\Delta E}}{{R{T_{\left( {i - 1,j} \right)}}}}} \right)\left[ {1 - \frac{1}{{K\left( {{p_{\left( {i - 1,j} \right)}},{T_{\left( {i - 1,j} \right)}}} \right)}}} \right]. \end{array} $ (20)
2.3 出砂过程与形态模拟方法

裸眼井壁或射孔孔眼附近的应力集中和流速条件都促使砂粒剥落产出首先从壁面上开始[25]图 4展示了带初始孔眼的天然气水合物储层单元颗粒级尺度微观出砂模拟的基本过程。图 4(a)中, 原始孔眼周围形成初始边界; 由于边界上各个颗粒的粒径、强度、水合物饱和度、微应力及周围流速不同, 其是否达到剥落的条件亦有差异。首先判断所有边界颗粒是否达到剥落条件, 达到条件的颗粒, 剥落形成出砂, 并形成如图 4(b)所示的新边界。在新边界上, 砂粒的微流场、应力等发生改变, 继续判断最容易剥落产出的砂粒模拟产出过程, 如图 4(c)所示。每一次的边壁颗粒脱落模拟, 相当于出砂空洞完成一次扩展; 以此反复, 可以完成整个颗粒剥落出砂过程以及出砂形态的模拟。

图 4 微观出砂模拟过程及形态扩展示意图 Fig.4 Process of sand fall-down and cavity extension

基于微观出砂过程模拟, 对于出砂临界条件预测, 只要模拟获得初始出砂的位置和生产条件(生产压差、流速/产量), 即可得到出砂临界条件; 对于出砂粒径, 则通过实时对全部已产出的砂粒直径进行数学统计即可获得; 对于出砂范围动态预测, 只需实时监控出砂孔洞的前沿位置即可; 对于累积出砂量, 则对应微观过程模拟中的出砂孔洞容积; 产出流体的含砂体积分数即为实时出砂速度与流体产出流量之比。

3 出砂形态模拟结果分析 3.1 颗粒骨架建模结果

利用本文中提出的颗粒级尺度微观出砂模拟方法, 研制开发自主产权的水合物储层出砂模拟器软件, 对典型的天然气水合物储层进行出砂模拟。模拟使用的基础数据:储层深度为1 381~1 399 m, 层厚为18 m, 井眼直径为219.07 mm, 裸眼完井, 储层压力为14.55 MPa, 储层温度为15.1 ℃, 储层泥质含量为20.660 3%, 岩石泊松比为0.238 6, 岩石弹性模量为148.6 MPa, 岩石内聚强度为0.31 MPa, 储层平均孔隙度为36.81%, 原始垂直主应力15.234 MPa, 平均水合物饱和度为0.48, 地层砂粒度中值为0.015 mm, 天然气相对密度为0.557 2, 生产压差为3.849 MPa, 原始最小水平主应力为13.695 MPa, 原始最大水平主应力为15.953 MPa。根据测井资料计算得到的该储层岩石内聚强度和水合物饱和度分布曲线如图 5所示。

图 5 测井得到的储层岩石内聚力及水合物饱和度纵向分布 Fig.5 Rock coheresive strength and gas hydrate saturation distribution from well logging

设定储层强度和水合物饱和度非均质系数均为0.15, 并认为其遵循正态分布, 在纵向上服从测井资料计算得到的强度和水合物饱和度分布规律。根据储层微观结构建模方法, 构建得到颗粒总胶结强度分布和水合物饱和度分布如图 6所示。

图 6 储层微观结构建模案例结果图像 Fig.6 Micro 2-D model of gas hydrate resorvior for case study

图 6(a)为按照总内聚强度分布输出的二维储层结构建模结果图像。红色深色区域表示较高的胶结强度, 绿色浅色区域表示较低的胶结强度。例如在井深1 384和1 386~1 387 m附近, 为高强度分布区域; 而在1 387和1 385 m附近为低强度区域。在井深纵向方向上, 总胶结强度的分布具有随机特性, 但其总体特征受测井资料获得内聚强度分布规律控制。图 6(b)为按照原始水合物饱和度分布输出的二维储层结构建模结果图像。相类似, 在井深纵向方向上水合物饱和度的随机分布总体特征受测井资料获得水合物饱和度分布规律控制。

3.2 颗粒级尺度出砂过程和形态模拟结果

为了探究微观出砂过程, 利用基础数据首先构建厚度1.5 mm的微单元进行颗粒级尺度微观出砂模拟。模拟得到的出砂过程及形态如图 7所示。图 7中红色区域为砂粒剥落区域, 每个细小方格即为一个地层砂颗粒位置。根据模拟结果, 天然气水合物微单元的水合物分解首先从孔壁开始, 出砂也从孔壁边界开始。在孔壁边界上, 颗粒的粒径和位置、胶结强度、应力和流体流速不同, 根据颗粒剥落判别条件, 最容易达到剥落条件的颗粒首先剥落产出, 如图 7(a)所示。颗粒剥落形成新的壁面边界并改变原有的流场和应力场; 在新的壁面边界上, 容易脱落的砂粒继续剥落产出, 出砂孔道向前扩展延伸, 如图 7(b)~(d)所示。

图 7 微单元颗粒级微观出砂形态模拟结果 Fig.7 Particle-scale sanding simulation image of micro rock unit

由于单元颗粒模型的非均质性, 出砂孔道为非均匀不规则扩展, 总体而言是沿胶结弱面向前延伸, 最终形成如图 7(d)所示的类蚯蚓洞形态的出砂孔道。水合物分解后储层出砂前沿在微观层面产生众多细小的出砂蚯蚓洞, 降低了储层整体强度和各部分之间链接, 最终造成宏观层面上储层呈现连续垮塌式出砂。孔道的具体位置、形态、深度等除了与生产条件有关外, 更与原始胶结强度、水合物饱和度、地层砂粒径等及其随机分布规律相关。

3.3 储层尺度出砂形态模拟结果

利用基础数据和图 6的颗粒微观建模结果, 设定水合物分解强度影响系数为0.85, 井底流压为3.85 MPa, 利用自主研发的模拟软件进行储层尺度的出砂过程和形态模拟, 得到如图 8所示的模拟结果。

图 8 储层尺度出砂过程和形态模拟结果 Fig.8 Resersoir-scale simulation image of sanding process and pattern

图 8清晰地表达了非均质天然气水合物储层的地层砂粒的剥落和出砂形态的演变过程。随着水合物分解前沿推进, 以裸眼井壁为边界的出砂剖面从边界上不均匀扩散(图 8(a)), 然后出砂亏空孔道以类蚯蚓洞的形态扩展(图 8(b)); 随着水合物分解及出砂前沿推进, 蚯蚓洞逐步延伸扩展, 相互接触, 造成连续垮塌式剥落和出砂。图 8的案例中最终呈现前端类蚯蚓洞与后端连续垮塌的复合出砂形态。整个储层从井筒到地层分为出砂区、水合物完全分解区、水合物分解过渡区、未分解区等几个区域。由于本例中储层胶结强度总体较低, 水合物饱和度较高, 水合物分解造成储层强度大幅度降低, 出砂前沿基本接近水合物分解前沿。

图 9为出砂前沿达到4.0 m时的储层胶结强度分布。图 9(a)出现连续垮塌式剥落出砂, 在剥落出砂区域的地层强度消失; 在已分解区域, 储层强度降低为原始胶结物胶结强度的85%(由水合物分解强度影响系数确定); 在分解过渡区, 总强度与水合物饱和度有关, 从胶结物胶结强度逐步过渡到未分解区的原始储层内聚强度。图 9(b)(c)中由于水合物分解对胶结物的影响, 分解区的胶结物胶结强度低于未分解区的原始胶结物胶结强度。

图 9 模拟储层胶结强度分布 Fig.9 Simulated rock strength distribution during sanding process

图 10为模拟得到的生产动态模拟及出砂速度结果。随着生产的进行, 水合物分解前沿向深部推进, 在有限的生产时间和空间范围内(本例中为18 d), 水合物分解面积增加, 使产水量和产气量增加; 生产18 d时, 分解前沿推进到约5 m, 出砂前沿与分解前沿十分接近(图 10(b))。本例中储层胶结物胶结强度平均约为0.15 MPa, 水合物胶结强度平均约为0.2 MPa, 在总内聚强度中占主要作用。水合物分解除了使水合物胶结强度消失, 还部分影响胶结物胶结强度。因此本例中水合物分解即意味着砂粒剥落出砂。从图 10(c)可以看出, 随着生产继续, 出砂粒径呈逐步减小趋势, 由初始的16 μm降低为6 μm; 这是由于随着分解和出砂前沿推进, 地应力集中现象减弱, 流体流动剖面面积增加, 流速降低, 促使砂粒剥落的条件减弱, 只有相对细小的颗粒达到剥落条件。与此同时, 由于出砂半径增加, 砂粒产出剖面增加, 出砂速度呈上升趋势。需要指出的是, 此处的出砂速度指的是砂粒的剥落速度。

图 10 生产动态模拟及出砂预测结果 Fig.10 Simulation reaults of production performance and sand production rate

为了探究水合物饱和度对储层强度和出砂规律的影响, 保持其他条件不变, 仅将平均水合物饱和度降低到0.4、0.3和0.18, 模拟得到水合物饱和度0.3条件下的出砂形态图像以及不同水合物饱和度下的剥落出砂速度, 如图 11所示。

图 11 不同水合物饱和度下的模拟结果 Fig.11 Simulation reaults of different gas hydrate saturation

根据图 11(a), 在储层平均水合物饱和度约为0.3条件下, 与图 8(d)相比, 总体出砂趋势明显减弱, 总体出砂形态呈现更加明显的类蚯蚓洞形态。随着储层平均水合物饱和度的降低, 砂粒剥落速度也明显降低(图 11(b)); 当储层初始平均水合物饱和度分别降低为0.4、0.3和0.18时, 砂粒剥落速度分别约降低到平均水合物饱和度0.48对应的砂粒剥落速度的80%、69%和58%。水合物饱和度越低, 水合物分解速率越低, 产气量和产水量越低, 流体对地层砂粒的冲刷携带作用越弱; 另外, 储层水合物饱和度越低, 水合物分解对储层总体内聚强度的影响越小, 即水合物分解区的储层强度越高, 相同生产条件下越不容易出砂。

4 结论

(1) 对于胶结强度较低的高水合物饱和度储层, 其微观出砂前沿首先从孔壁边界上不均匀扩散, 出砂孔洞沿胶结弱面延伸形成类蚯蚓洞形态; 随着逐步延伸扩展, 造成连续垮塌式剥落出砂。总体出砂形态呈现前端类蚯蚓洞与后端连续垮塌的复合出砂形态。

(2) 对于总体内聚强度较低且水合物胶结强度占主要地位的水合物储层, 水合物分解会大幅度降低储层强度, 水合物分解即意味着砂粒剥落出砂, 生产过程中出砂前沿与分解前沿非常接近。在生产初期随着生产延续, 最大剥落出砂粒径由16 μm降低为6 μm, 剥落出砂速度呈上升趋势。

(3) 水合物的饱和度及其分解对出砂形态和出砂速度有明显影响。水合物饱和度越低, 气水产量越低导致对砂粒的冲刷携带作用减弱; 低饱和度储层的水合物分解对储层总体内聚强度的影响较小, 相同生产条件下不容易出砂。

参考文献
[1]
付亚荣. 可燃冰研究现状及商业化开采瓶颈[J]. 石油钻采工艺, 2018, 40(1): 68-80.
FU Yarong. Research status of combustible ice and the bottle neck of its commercial exploitation[J]. Oil Drilling & Production Technology, 2018, 40(1): 68-80.
[2]
卢静生, 李栋梁, 何勇, 等. 天然气水合物开采过程中出砂研究现状[J]. 新能源进展, 2017, 5(5): 394-402.
LU Jingsheng, LI Dongliang, HE Yong, et al. Research status of sand production during the gas hydrate exploitation process[J]. Advances in New and Renewable Energy, 2017, 5(5): 394-402. DOI:10.3969/j.issn.2095-560X.2017.05.011
[3]
董长银, 钟奕昕, 武延鑫, 等. 水合物储层高泥质细粉砂筛管挡砂机制及控砂可行性评价试验[J]. 中国石油大学学报(自然科学版), 2018, 42(6): 79-87.
DONG Changyin, ZHONG Yixin, WU Yanxin, et al. Experimental study on sand retention mechanisms and feasibility evaluation of sand control for gas hydrate reservoirs with highly clayey fine sands[J]. Journal of China University of Petroleum (Edition of Natural Science), 2018, 42(6): 79-87. DOI:10.3969/j.issn.1673-5005.2018.06.009
[4]
赵景芳, 宋林松, 吉飞, 等. 天然气水合物降压开采储层出砂数值模拟[J]. 中国海上油气, 2019, 31(2): 116-124.
ZHAO Jingfang, SONG Linsong, JI Fei, et al. Numerical simulation of reservoir sand production during gas hydrate depressurizing production[J]. China Offshore Oil and Gas, 2019, 31(2): 116-124.
[5]
程远方, 沈海超, 赵益忠, 等. 天然气水合物藏开采物性变化的流固耦合研究[J]. 石油学报, 2010, 31(4): 607-611.
CHENG Yuanfang, SHEN Haichao, ZHAO Yizhong, et al. Study on fluid-solid coupling of physical variation of gas hydrate reservoirs during natural gas development[J]. Acta Petrolei Sinica, 2010, 31(4): 607-611.
[6]
UCHIDA S, KLAR A, YAMAMOTO K. Sand production model in gas hydrate-bearing sediments[J]. International Journal of Rock Mechanics and Mining Sciences, 2016, 86: 303-316. DOI:10.1016/j.ijrmms.2016.04.009
[7]
JUNG J W, JANG J, SANTAMARINA J C, et al. Gas production from hydrate-bearing sediments: the role of fine particles[J]. Energy & Fuels, 2012, 26(1): 480-487.
[8]
MOHAMAD-HUSSEIN A, NI Q. Numerical modeling of onset and rate of sand production in perforated wells[J]. Journal of Petroleum Exploration and Production Technology, 2018(8): 1255-1271.
[9]
MORITA N, WHITFILL D L, MASSIE I, et al. Realistic sand-production prediction: numerical approach[J]. SPE Production Engineering, 1989, 4(1): 15-24. DOI:10.2118/16989-PA
[10]
董长银, 周崇, 钟奕昕, 等. 中等强度砂岩饱水力学参数变化试验及动态出砂规律[J]. 中国石油大学学报(自然科学版), 2017, 41(6): 108-116.
DONG Changyin, ZHOU Chong, ZHONG Yixin, et al. Experimental study on mechanical properties of medium strength sandstone saturated with water and its dynamic sanding behaviors[J]. Journal of China University of Petroleum(Edition of Natural Science), 2017, 41(6): 108-116. DOI:10.3969/j.issn.1673-5005.2017.06.013
[11]
董长银, 张清华, 崔明月, 等. 复杂条件下疏松砂岩油藏动态出砂预测研究[J]. 石油钻探技术, 2015, 43(6): 80-86.
DONG Changyin, ZHANG Qinghua, CUI Mingyue, et al. A dynamic sanding prediction model for unconsolidated sandstone reservoirs with complicated production conditions[J]. Petroleum Drilling Techniques, 2015, 43(6): 80-86.
[12]
PINKERT S, GROZIC J L H. Failure mechanisms in cemented hydrate-bearing sands[J]. Journal of Chemical & Engineering Data, 2015, 60(2): 376-382.
[13]
KAJIYAMA S, HYODO M, NAKATA Y, et al. Shear behaviour of methane hydrate bearing sand with various particle characteristics and fines[J]. Soils & Foundations, 2017, 57(2): 176-193.
[14]
UCHIDA S, KLAR A, YAMAMOTO K. Sand production modeling of the 2013 Nankai offshore gas production test: 1st International Conference on Energy Geotechnics, 29-31 August, 2016, Kiel, Germany[C].London: Taylor & Francis Group, 2016: 451-458.
[15]
KLAR A, UCHIDA S, SOGA K, et al. Explicitly coupled thermal flow mechanical formulation for gas-hydrate sediments[J]. SPE Journal, 2013, 18(2): 196-206. DOI:10.2118/162859-PA
[16]
LI Y, SONG Y, LIU W, et al. A new strength criterion and constitutive model of gas hydrate-bearing sediments under high confining pressures[J]. Journal of Petroleum Science & Engineering, 2013, 109(5): 45-50.
[17]
宁伏龙, 吴能友, 李实, 等. 基于常规测井方法估算原位水合物储集层力学参数[J]. 石油勘探与开发, 2013, 40(4): 507-512.
NING Fulong, WU Nengyou, LI Shi, et al. Estimation of in-situ mechanical properties of gas hydrate-bearing sediments by well logging[J]. Petroleum Exploration and Development, 2013, 40(4): 507-512.
[18]
COLLETT T, BAHK J J, BAKER R, et al. Methane hydrates in nature-current knowledge and challenges[J]. Journal of Chemical & Engineering Data, 2015, 60(2): 319-329.
[19]
RUTQVIST J, MORIDIS G J, GROVER T, et al. Coupled multiphase fluid flow and wellbore stability analysis associated with gas production from oceanic hydrate-bearing sediments[J]. Journal of Petroleum Science & Engineering, 2012, 92/93(4): 65-81.
[20]
刘乐乐, 鲁晓兵, 张旭辉. 天然气水合物分解区演化数值分析[J]. 石油学报, 2014, 35(5): 942-951.
LIU Lele, LU Xiaobing, ZHANG Xuhui. Numerical analysis on evolution of natural gas hydrate decomposition region in hydrate-bearing sediment[J]. Acta Petrolei Sinica, 2014, 35(5): 942-951.
[21]
薄启炜, 董长银, 张琪, 等. 砾石充填层孔喉结构可视化模拟[J]. 石油勘探与开发, 2003, 30(4): 108-110.
BO Qiwei, DONG Changyin, ZHANG Qi, et al. Visual simulation of porous structure in packed gravels[J]. Petroleum Exploration and Development, 2003, 30(4): 108-110. DOI:10.3321/j.issn:1000-0747.2003.04.035
[22]
胡高伟, 李承峰, 业渝光, 等. 沉积物孔隙空间天然气水合物微观分布观测[J]. 地球物理学报, 2014, 57(5): 1675-1682.
HU Gaowei, LI Chengfeng, YE Yuguang, et al. Observation of gas hydrate distribution in sediment pore space.ChineseJ[J]. Earth and Planetary Physics, 2014, 57(5): 1675-1682.
[23]
颜荣涛, 韦昌富, 魏厚振, 等. 水合物形成对含水合物砂土强度影响[J]. 岩土工程学报, 2012, 34(7): 1234-1239.
YAN Rongtao, WEI Changfu, WEI Houzhen, et al. Effect of hydrate formation on mechanical strength of hydrate-bearing sand[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(7): 1234-1239.
[24]
陈月明, 李淑霞. 天然气水合物开采理论与技术[M]. 青岛: 中国石油大学出版社, 2011: 20-30.
[25]
董长银, 李怀文. 油气井防砂技术[M]. 2版. 北京: 中国石化出版社, 2017: 15-40.