天然气水合物非成岩出砂储层砂流数值模拟

万庭辉，张可倪，李占钊，王静丽，于彦江

(1.中国地质调查局广州海洋地质调查局，广州 510760；2.东北石油大学石油工程学院，大庆 163318；3.暨南大学地下水与地球科学研究院，广州 510632)

天然气水合物分布广、埋藏浅、清洁无污染、储量巨大，被视为油气领域最具潜力的清洁能源。苏联、加拿大、美国、日本、中国等国家已先后在冻土层或海域成功进行了天然气水合物试验性开采。2017年5月18日，由中国地质调查局组织实施的南海北部神狐海域水合物试采成功，实现了持续产气60 d和产气总量超过30万m3的世界纪录。2020年3月，由中国地质调查局组织实施的中国海域天然气水合物第二轮试采取得成功，并超额完成目标任务。本轮试采1个月产气总量86.14万m3、日均产气量2.87万m3，是第一轮60 d产气总量的2.8倍。攻克了深海浅软地层水平井钻采核心技术，实现了从探索性试采向试验性试采的重大跨越[1-12]。

天然水合物开采过程中，水合物相变引起的储层应力行为的改变，可能会导致一系列地质风险问题，如井壁失稳、井眼出砂、地层垮塌以及海床滑坡等。苏联麦索亚哈陆域水合物田(1967年)[13-14]、加拿大麦肯齐陆域水合物试采(2007年和2008年)[15-17]、阿拉斯加北坡陆域水合物试采(2008年和2012年)[17-20]和日本南海海槽海域水合物试采(2013年和2017年)[18-20]，这些试采工程均以降压开采为主，大部分地区的开采过程中均有不同程度的出砂情况[20-22]；尤其是2013年和2017年日本南海海槽进行的两次海域天然气水合物试采均因井筒大量出砂而被迫中止，这表明出砂问题严重影响天然气水合物开采过程中的安全和效率，全面了解“砂”的行为，以期实现有效控砂，是水合物安全高效、长期可控开采的关键。

数值模拟是了解天然气水合物开采储层出砂的有效方法。在模型研究方面，Ranjith等[23]综述了2014年以前开展的与出砂问题相关的模型研究。认为出砂的基本机理包括剪切和拉伸破坏、临界压力梯度、压力下降和侵蚀。目前为止，连续体和离散元法都已应用在出砂模型中。其中连续体法包括基于塑性理论、混合理论及微观力学建立起来的方法。大部分方法都考虑了流体力学或剪切破坏机制，有些将它们结合起来。Uchida等[24-25]提出了一个水合物储层出砂模型，综合描述了水合物的出砂特征，包括颗粒的分离、运移、沉积变形和水合物的离解。这些过程是传热-流动-力学的耦合过程。颗粒脱离引起应力降低，沉积物剪切变形引起颗粒脱离，颗粒流动改变多相流体压力和温度分布。研究发现，在不同的操作方法中，减小降压率是降低出砂量的最有效方法。Yan等[26]利用多场耦合理论系统分析了影响井周围分解分布、塑性变形区及出砂规律的影响因素，结合TOUGH+HYDRATE和ABAQUS来分析出砂的多场耦合。Ning等[27]构建了基于连续-离散介质理论的水合物开采储层响应与出砂数值模型，将连续介质模型获取的应力、流速和水合物饱和度等数据作为边界和初始条件传递给离散元模型，进而分析出砂机理、规律并预测出砂量。

目前中外水合物开采出砂数值模拟研究大多是建立温度-水力-应力-化学耦合出砂模拟模型，使用TOUGH+HYDRATE水合物模拟器获取开采时流场如压力、饱和度、产气产水速率等参数变化以及温度场如温度和热量变化规律，然后导入应力场模拟软件FLAC3D，获取地层中的应力和位移演化情况，结合颗粒流软件PFC3D、工程模拟有限元软件ABAQUS等[28-39]，对水合物储层出砂机制和出砂规律进行研究。此类研究方法考虑影响因素全面，虽然有通过力学耦合描述“砂”的行为，但都没有涉及“砂”的流动过程。针对泥质粉砂型天然气水合物储层，水合物分解后储层多孔介质中“砂”流动过程的数值模拟，现通过通用储层渗流模拟器GPSFLOW，实现水、气(甲烷)、砂三相三组分模拟，把砂当作一种非牛顿流体，建立宾汉流体型砂流模型，把“砂”当成宾汉流体进行模拟，并用理想模型进行验证。

1 软件介绍及研究方法

1.1 软件简介

GPSFLOW(general-purpose subsurface flow simulator)软件是一套考虑等温复杂介质系统中油、气、水三相的牛顿和非牛顿流体多相流动的商业数值模拟程序，允许用户设定液相组分为牛顿流体或非牛顿流体。GPSFLOW假设非牛顿流体流动符合幂律、宾汉流体等流变模型，黏度随着剪切速率和孔隙流速的变化而变化。程序可应用于一维、二维或三维的、各向异性的、非均质孔隙裂缝系统。

1.2 研究方法

非牛顿流体的表观黏度通常取决于孔隙流速或压力梯度。不同非牛顿流体在孔隙介质中流动的流变函数不同。图1给出了在孔隙介质流动中经常遇到的非牛顿流体的剪切应力和剪切速率间的典型关系。

τ为剪切应力；γ为剪切速率；μ为流体动力黏性系数(即黏度)；μb为宾汉流体的塑性黏度系数；μa为假塑性流体的塑性黏度系数；A、B、C、D分别为流体的拟合曲线

在天然气水合物储层中，水合物作为胶结物能有效增加储层强度，水合物分解后，储层胶结程度降低，地层应力重新分布，当有效压力梯度达到临界值时，孔径较小的砂体颗粒被水合物分解产生的水气携带并开始流动，表现为储层出砂；当有效压力梯度不满足条件时，微小砂粒度又会静止沉淀，反过来影响储层孔隙度、渗透率等，这一现象和非牛顿流体中宾汉流体的流动特性非常相似。通过模拟器GPSFLOW，在软件中将“非水”组分变为“砂”，并且把“砂”当成宾汉流体模拟。在多相流达西定律的基础上，软件结合砂流的非牛顿流动过程，砂、水与气的混合流动过程，描述水合物分解后储层中“砂”的运移过程。

1.2.1 控制方程

对于一个包含气、水、砂三种物质成分的三相流体等温系统，任意流动区域的三相物质平衡方程如下。

ρgVg)+qg

(1)

(2)

(3)

β相(β=g为气，β=w为水，β=s为砂)的流动服从多相流达西定律，即

(4)

三种流体相在任意时刻的孔隙中均满足如下约束条件：

Sw+Ss+Sg=1

(5)

砂作为特殊流体，假设其相对渗透率总是1.0，与水相和气相之间不存在虹吸现象，即毛细压力为0。水气之间按常规的两相流处理其毛细压和相对渗透率的计算。

1.2.2 方程离散

GPSFLOW程序利用有限体积法对连续性方程进行空间离散，用一系列积分有限差分方程表达油气水的物质平衡方程；利用全隐式的稳定性和时间步长较长的特点，或利用自适应隐式方法(AIM)可以加快模拟和降低空间存储要求的特点去求解这些离散非线性方程，对于简单的问题，也可以采用隐式压力显式饱和度(IMPES)方法，以提高计算速度；一般通过确定有限子域或网格块的属性来描述流体和岩石的热力学性质；通过有限差分逼近法去计算通过连通网格块表面部分的质量流量；应用牛顿-拉甫森迭代程序求解这些离散的非线性的有限差分物质平衡方程。

对上述质量守恒方程用积分有限差分方法把连续性方程空间离散化，用向后差分方法实现该方程的时间离散化。写成残差形式为

(6)

(7)

(8)

式中：i=1,2,…,N；m可以是时间点n或n+1，如果m=n，即单元i是应用隐式压力显式饱和度法(IMPES)，如果m=n+1，即单元i是应用全隐式方法；N为网格的总数量；n为前时间点；n+1为有待解决的现有时间点；Vi为成分i的体积;R为组分在当前时间点质量守恒计算的残余量；λ为β相流度，定义为流体的相对渗透率与黏度之比；γ为导流系数；ψ为网格的总相压力；Q为源汇项；Δt为时间步长;ηi包含临近块(j)组或单元i的结点，这样直接包含了单元i；下脚标ij+1/2为单元i和j分界面的一个恰当平均值，与β相流度(相对渗透率与黏度之比)有关。

流动边界的传导系数定义为

(9)

(10)

式中：Aij为连通单元i和j的共同界面面积；di为从单元i的中心到单元i和j交界面的距离；dj为从单元j的中心到单元i和j交界面的距离；kij+1/2为沿着单元i和j连通处的平均绝对渗透率；Di为单元i中心的深度；单元i的汇点/源点项定义为

(11)

应用牛顿-拉甫森迭代法去求解一个流动系统的方程，每个单元都有油气水三成分的物质平衡方程，表现为3×N个耦合非线性方程。每个单元选择3个主要变量(x1、x2、x3)，分别是水压力、含砂饱和度和含气饱和度，如表1所示。

表1 主要变量的选择

根据这3个主要变量，利用牛顿-拉甫森迭代方法得

[xk,p+1-xk,p]=0

(12)

式(12)中：指数k=1，2，3分别代表主要变量x1、x2、x3的下角标1，2，3；p为迭代点。式(12)可表示为

k=1,2,3，i=1,2,…,N

(13)

迭代中产生的主要变量的增量表达式为

δxk,p+1=xk,p+1-xk,p

(14)

式(13)代表3×N个未知δxk,p+1的一系列3×N个线性方程。

本程序采用这样一个数值算法去建立式(13)的雅可比矩阵：对于全隐式项来说，用数值微分法去求解雅可比矩阵。对于隐式压力显式饱和法来说，用更简单的半解析方法去求解雅可比矩阵。形成的线性方程可通过迭代法去求解。

1.2.3 宾汉流体流变模型

宾汉流体是一种黏塑性材料，在低应力下表现为刚性，在高应力下像黏性流体一样流动。砂流将按非牛顿流体类型中的宾汉流体进行处理。由于其切应力与剪切速率呈线性关系(图1)，因此在高应力下宾汉流体具有牛顿流体性质[40-41]。

(15)

式(15)中：τ为剪应力，Pa；dv/dy为剪切速率，s-1；η为运动黏性系数，Pa·s；τ0为屈服应力或屈服值，Pa。在数值模拟过程中，对宾汉流体的有效压力梯度进行修正比描述其表观黏度有更高的运算效率[42]，宾汉流体服从的达西定律为

(16)

(17)

式中：μb为宾汉流体的塑性黏度系数，与砂质量分数相关，Pa·s；∇Φe为有效压力梯度，Pa/m；k为绝对渗透率，D；kr为相对渗透率；v为达西流速，m/s；G为使得宾汉流体流动的最小压力梯度，Pa/m[43]，满足：

(18)

式(18)中：α为经验系数，或者拟合参数。

1.2.4 流体性质计算

模型中储层流动系统由水、气和砂3个组分组成。流体性质计算主要包括流体的密度、黏度、相对渗透率和毛细压等，这些参数一般是温度与压力的函数，气体(甲烷)参数可通过PVT表输入进行插值，模型暂时只考虑等温情况。砂的黏度和密度被认为是常数，不受温压变化影响，相对渗透率总是1，毛细压为0。

2 软件应用

2.1 模型建立

建立均值理想模型，模拟开采区域长宽均为1 000 m，厚10 m；模型四周为恒压边界，开采井井筒半径为0.15 m，位于模型中心处。均值储层原始渗透率k=2 mD，孔隙度Φ=0.41，储层初始压力15 MPa，储层初始温度14.1 ℃；假定水合物分解后储层多孔介质中砂饱和度为Ss=0.1，气饱和度为Sg=0.2，其余为水。砂流启动压力梯度为1 MPa/m，其中残余砂饱和度为0.067，理想模型储层特征参数如表2所示。

表2 理想模型储层特征参数

理想模型xy平面网格数量为3 293个(图2所示)，z方向上分为10个网格。整个模拟过程所要求解的方程数为32 930个，近井筒周围网格进行了加密处理。

图2 网格划分示意图

2.2 模拟结果

上述理想模型以井底流压Pwf=0.5 MPa开采10 d后，随即以井底流压Pwf=2 MPa开采50 d，该模型开采10、30、50、60 d的压力和砂饱和度场图分别如图3和图4所示。

图3 压力场图

图4 砂相饱和度场图

模拟开采期间，砂的流动呈现宾汉流体的流动特征，砂的产出速率和累积产砂量如图5所示。开采初期(0～10 d)，压力迅速传递，生产压差小于启动压力梯度，砂相不能流动，产砂速率为0 m3/d，累积产砂0 m3；开采后期(10～60 d)，生产压差大于启动压力梯度，砂相持续流动，由于井周砂相饱和度不断增加，井周渗透率不断减小，产砂速率由最初的0.6 m3/d，逐渐减小至0.000 01 m3/d，累积产砂1.62 m3。

图5 产砂速率及累积产砂

2.3 效果分析

为探究泥质粉砂型天然气水合物储层，水合物分解后储层多孔介质的宏观砂流运移规律，结合通用储层渗流模拟器GPSFLOW，创新性地提出把“砂”当作非牛顿流体，通过建立宾汉流体型砂流模型，实现水、气(甲烷)、砂三相三组分模拟，并用理想模型进行了验证。

模拟结果表明该砂流模拟方法简单高效，并可提供试采场地规模、长时间尺度的砂流运移规律可视化显示；后续通过耦合水合物开采渗流模型，结合试采场地天然气水合物储层层间非均质地质模型、各小层砂流启动压力梯度、降压条件等实际生产数据，未来有望实现对水合物开采过程中场地规模砂流运移规律的动态监测，有力支撑天然气水合物试采工程现场决策。

3 结论

探讨了“砂”在非成岩储层的流动过程，使用通用储层渗流模拟器GPSFLOW实现宾汉型砂流的数值模拟。模拟器充分考虑了“砂”的非牛顿流动过程，可实现三维非结构网格、非均质孔隙/裂隙介质中的“砂”运移情况的模拟，通过隐式积分有限差分求解模型，具有良好的稳定性和收敛性。提出的砂流模拟方法简单高效，后续通过耦合水合物开采渗流模型，有望实现对水合物开采过程中砂流的准确模拟。通过设置不同岩性、砂的性质、降压条件等，可以分析水合物开采潜力、砂的运移规律、地层压力分布，预测地层出砂情况和地层稳定性，为实际试采工程提供技术支撑，预测出砂情况和评估潜在风险。