页岩气藏两相流固耦合无网格数值模拟.pdf

收稿日期:改回日期:基金项目:国家自然科学基金“考虑微观力作用的纳微米孔隙流动理论及模拟研究”()中国石油西南油气田公司科研项目“页岩气多尺度非线性渗流规律研究及应用”()作者简介:吴建发()男教授级高级工程师 年毕业于西南石油学院石油工程专业 年毕业于西南石油大学油气田开发工程专业获博士学位现主要从事页岩气开发技术研究工作:/页岩气藏两相流固耦合无网格数值模拟吴建发朱维耀张德良陈 震吴天鹏(中国石油西南油气田分公司四川 成都 北京科技大学北京)摘要:针对页岩气藏受两相复杂流动与流固耦合作用影响导致产能预测难度大的问题将渗流场分为主改造区、次改造区和未改造区结合多流态统一输运模型及等效连续非均匀介质物理模型建立多尺度气水两相流动流固耦合数学模型并利用无网格广义有限差分法对数学模型进行编程求解 研究结果表明:无网格广义有限差分法能适应不同计算域情况避免传统差分网格与非结构网格耦合所导致的运算不稳定压降前缘在未改造区的传播仅为 左右固体位移主要发生在次改造区与未改造区的交界区域忽略应力场的影响气井投产初期产气量被显著高估到生产后期这种影响逐渐下降低初始含水饱和度的页岩气藏更有益于开发储层优选时应重点考虑含水饱和度 研究成果对流固耦合数值模拟及非常规油气藏产能预测具有指导意义关键词:页岩气藏无网格法数值模拟多尺度渗流流固耦合中图分类号:文献标识码:文章编号:()(.):.:引 言页岩气是中国天然气开发的重要组成部分近年来占比达到 左右 页岩气藏储层地质条件复杂而两相流动与流固耦合是其复杂性的主要根源可导致页岩气井产量预测存在较大不确定 第 期吴建发等:页岩气藏两相流固耦合无网格数值模拟 性 目前对页岩气产量的预测通常仅考虑含水对气体流动的单方面影响而不是多尺度条件下的两相流关于流固耦合问题大多采用经验公式或本构方程但为进一步提高预测的准确性应建立完整的流固耦合数学模型 因此页岩气井产量预测需同时建立考虑两相流以及流固耦合的页岩气生产数学模型但这将导致模型具有非常复杂的非线性 传统基于网格的算法一般采用有限差分或有限体积法求解渗流场、有限元法求解应力场每个时间步中算法和网格的频繁更换使计算效率和稳定性严重下降且非常依赖于网格划分质量 近年来无网格法由于其相对的简单性和稳定性引起广泛关注并已经被用于页岩气藏模拟 无网格法不但被证实具有足够解决非线性问题的能力而且可达到与传统方法相同的精度 因此该文建立了一套同时考虑两相流、流固耦合以及基质缝网多区耦合的页岩气生产数学模型并利用无网格广义有限差分法对其进行求解以期有效模拟页岩气藏的流动过程 页岩气藏物理模型 考虑纳微米效应的页岩气多尺度流动页岩气的流动贯穿纳微米级孔隙、微裂缝、人工裂缝呈现多尺度运移的特征其中基质内的运移为受纳微米效应所影响的非线性流动而裂缝内的流动则为线性达西流 依据 等所推导的气体传输多尺度多流态统一输运模型页岩气在多尺度介质内的传输需要考虑如式()所示的表观渗透率:()式中:为表观渗透率为多孔介质的有效渗透率 为气相流体黏度为气体的努森扩散系数/为流体压力 为滑移系数通常被指定为考虑到水的存在页岩气的运动方程需考虑相对渗透率:()式中:为气相流体流速/为气相的相对渗透率 为 个渗流截面间的距离页岩气在基质内存在吸附解吸行为采用 等温吸附模型进行描述:()式中:为单位质量页岩的等温吸附量/为 体积/为 压力 “人工气藏”物理模型水力压裂后的页岩储层是具有多尺度复杂裂隙缝网结构的多孔岩石介质其缝网形态的数学描述是页岩气开采储层介质模型的关键 根据已有研究成果将页岩气的流动分为 大区域:主改造区、次改造区、未改造区(图 图中蓝色曲线表示各区渗透率变化)页岩气由未改造区流入改造区再由改造区流入水平井筒形成页岩气储层完整的流动体系 物理模型中改造区内的流动为线性达西流未改造区中的流动为考虑纳微米效应影响的非线性流 水平井筒内的压降相对很小通常对产量的影响不超过 因此考虑为无限导流即忽略井筒内的压降各射孔点处的井底流压均相等图 储层压裂分区示意图 为更好地模拟改造区与未改造区的耦合使各种物理量等效地分布在模型上采用 等提出的非均匀等效连续介质模型渗透率的平面分布如图 所示 该模型本质上是复合分区模型如何将显式的裂缝转换为等效的连续介质其关键是要得到主改造区的等效渗透率()该文采用 种基本缝网形态的等效形式将不同裂缝分布形态划分为羽状缝、簇状缝、网状缝和树状缝不同形态下的等效计算数学模型见文献 特 种 油 气 藏第 卷50040302010?/mD图 连续等效渗透率分布 页岩气藏流固耦合渗流数学模型 气水两相渗流控制方程气水两相的渗流数学模型以连续性方程为基础描述渗流控制数学模型中某微元的两相连续性方程为:()()()式中:为散度算子、分别为气、水的渗流速度矢量/、分别为气、水的密度/、分别为单位体积岩体所含气、水的质量/、分别为气、水质量源汇强度/()为时间气体的运动方程见式()水的运动方程除遵循达西定律外还需要引入毛管力即:()()式中:为水相流体的相对渗透率为水相流体黏度为气水之间的毛管压力是含水饱和度 的函数考虑气相与水相的饱和度之和恒为 且气体存在吸附质量、通过单位体积多孔介质内的流体质量来表征即:()()()式中:为页岩孔隙度为含水饱和度为标态条件下的甲烷密度/为页岩密度/将式()、()代入式()中并考虑气相存在吸附解吸的过程由此可以得到气水两相的渗流控制方程:()()()()()()()式中:为梯度算子()表示计算标量函数 的梯度 应力场控制方程孔隙压力的变化导致有效应力发生变化从而引起气藏岩石变形因此流固耦合关系是依靠有效应力原理来体现的 考虑到页岩储层具有甲烷吸附所引起的系统应变的特征得到页岩储层应力场控制方程为:()()()()()式中:、为、方向上的位移量 为页岩的杨氏弹性模量 为页岩的泊松比 为 系数为岩石的体积模量即所受压力与体积应变之比为吸附应变系数/、分别为二维情况下的空间坐标由于物理模型采用等效连续介质因而弹性模量和泊松比也随区域位置的不同而改变 引入非均匀分布的权重因子 弹性模量 和泊松比 的分布为:()()()式中:、分别为模型主改造区和未改造区的弹性模量、分别模型为主改造区和未改造区的泊松比 第 期吴建发等:页岩气藏两相流固耦合无网格数值模拟 流固耦合问题有效应力的增大将导致页岩储层的孔隙裂隙空间减小引起有效孔隙度和渗透率下降 当孔隙度从初始状态经过应变达到另一状态时其关系式可表示为:()式中:为储层初始状态条件下的孔隙度为体积应变且有:()式中:为网格体积为网格体积的变化量、分别为 和 方向上的线应变渗透率遵循与孔隙度相关的三次方关系:(/)()式中:为储层初始状态的渗透率除孔隙度、渗透率变化外气水的物理性质包括密度和黏度也会随压力的变化而改变 目前已有不少研究采用经验公式来拟合气水两相的密度、黏度数值为确保计算准确搜集整理了压力为 、温度为 的气水密度和黏度数据在数值模拟的全程加以调用 数学模型的求解 无网格广义有限差分法数值模拟中无网格法以“点簇”的概念代替网格 设某待求解偏微分方程的形式如下:()()式中:为待求量为方程系数()为有关 的任意初等函数在计算域中随机选择一个节点(设编号为)假设在围绕点 的一定距离内(距离范围越大计算越精确但效率降低)共有 个邻点那么对未知量 在点 处做泰勒展开可在点 处形成如下残差函数:()()()()式中:()为残差函数 为点 的邻点的编号取值为、和 为点 与点 分别在 和 方向上的距离为点 关于点 的广义有限差分的权重函数由文献获得当 取得准确值时残差函数()就等于零根据偏导数的意义此时()的偏导数应取极小值 记未知量在点 处的偏导数向量为:()式中:为矩阵的转置符号为未知量在点 处的偏导数向量将残差函数()转化为:()()()()()()()()()()式中:为对角矩阵 为求解特性矩阵与 和 的值有关联立式()()经系列运算得到:()式中:、分别为整理得到的左端系数矩阵与右端向量具体按照下式计算:、()因而对于原偏微分方程式()可将其左端直接转化为待求量的线性组合 若记原方程左端的系数矩阵为 则 ()求解该方程即获各节点处的数值结果 渗流场和应力场数学模型的求解 渗流场的求解由于渗流场方程包含对时间的微分同时系数具有非线性因而对时间项做差分 将所有偏微分置于方程左端其余项置于方程右端式()、()化为:()特 种 油 气 藏第 卷()式中:、均为与 时刻有关的系数显然根据广义有限差分法将式()、()中关于压力 偏导数的系数矩阵分别与 相乘即可将压力 偏导数转化为线性代数和的形式 将待求量的线性和全部置于左端化为:()式中:、为方程的压力系数、为方程的饱和度系数、分别为 个方程的常数项需要注意的是中包含关于毛管力 的偏导数和 若按通常数值微分的方法求取即按照导数的数学定义进行计算那么可能会由于数值逼近的截断和摄入误差导致结果并不可靠 然而通过广义有限差分方法可更准确地求取这 个偏微分具体是通过将 个偏微分项的系数矩阵与 相乘再与所有节点的毛管力 矩阵相乘即得到偏导数和的值 应力场的求解式()中应力场方程组的待求量不仅包括 和 也包括压力 个待求量的偏导数存在 套系数矩阵同样与 相乘即可将式()转换为线性模式 将待求量全部置于左端其余项置于右端式()化为:()式中:、均为方程系数、为方程的常数项此外常数项、中所包含的偏导数、同样通过广义有限差分法进行计算 至此该流固耦合问题可以实现求解 模型验证及示例模拟分析 模型验证为验证文中模型通过两相流及流固耦合的 个案例进行验证 页岩气水两相流的案例来自文献的研究 文献采用了有限元与有限体积结合的方法对页岩气水两相流进行了模拟计算文中考虑模型中的气、水两相流并将流固耦合模型蜕化为与文献相同的应力敏感经验公式对计算结果进行验证(图)计算所采用的参数均来自文献 由图 可知文中计算结果与文献的研究结果基本一致最大误差小于 图 与文献的计算结果对比 流固耦合的案例来自文献的研究 文献采用了现有的有限元计算方法库模拟了流固 第 期吴建发等:页岩气藏两相流固耦合无网格数值模拟 耦合条件下的页岩气单相流动 文中考虑模型蜕化为气相单向流即取消水相流动的控制方程其余特性与现有模型相同与文献的计算结果进行对比验证(图)计算所采用的参数均来自文献 由图 可知文中计算结果与文献的研究结果基本一致最大误差小于 通过以上 种模型的对比验证表明文中提出的模型具有较高的可靠性图 与文献的计算结果对比 示例计算根据提出的数值模型及求解方法进行示例计算 模型规模为 初始条件如下:压力初始值为 人工改造区的含水饱和度最高为 基质区的含水饱和度为 储层温度为 边界条件如下:射孔点的井底流压为 外边界无流动 根据中国南方海相龙马溪组页岩气储层物性模拟的基本参数参考文献分别模拟不同流固耦合参数影响下压力及产气量变化曲线模拟计算 得到压力、饱和度和位移的分布结果如图 所示 由图 可知:时压力扰动波及改造区所有区域前缘已突破到未改造区但在未改造区的传播非常有限仅约 左右 由图 可知:主改造区的高含水随着生产进行而下降说明开井后的压裂液返排次改造区内的含水饱和度上升表明气相具有更好流动性随着压力下降开始进入强改造区而水则滞留在渗透率较低的区域未改造区含水饱和度上升的区域范围较小证实了压力波及前缘的传播距离在未改造地带十分有限 由图 可知:位移的变化量随着生产进行逐渐增加主改造区的压降最大但处于应力平衡状态没有发生较大的位移位移主要发生在次改造区与未改造区的交界区域因为页岩基质具有更显著的应力敏感性同时压降前缘仅能波及到这一区域35525201510?/MPa30a?10db?1000d图 压力场的模拟结果 0 8.0 4.0 7.0 6.0 5.?a?10db?1000d图 含水饱和度的模拟结果 文中求解方法具有较高的稳定性要求条件仅为求解域是光滑函数计算过程中没有遇到不收敛的情况 因此该方法对不同节点排布的适应性很强 特 种 油 气 藏第 卷14064281012?/cma?10db?1000d图 储层岩石位移的模拟结果 影响因素分析 流固耦合的影响模拟计算了考虑和不考虑流固耦合时的 产气量、产水量(图)由图 可知:流固耦合作用明显降低了生产井的产气量 累计产气量相差最高可达 但随着生产进行这种差异逐渐降低不考虑应力场的耦合作用产水期大幅缩短产水量显著增大说明流固耦合作用导致储层弹性变形对产水的影响更加显著图 流固耦合对页岩储层累计产量的影响 初始含水饱和度的影响模拟计算了不同初始含水饱和度(记为)条件下 的产气量和产水量(图)由图 可知:产气量随着含水饱和度的升高而下降且下降幅图 初始含水饱和度 对页岩储层累计产量的影响 第 期吴建发等:页岩气藏两相流固耦合无网格数值模拟 度基本一致而产水量则随之上升且上升幅度逐渐增大 主要是因为毛管力对气、水的作用机制不同:气体的渗流控制方程并未包含毛管力但水的渗流控制方程中毛管力却是流动阻力 当含水饱和度较低时毛管力较高对水产生了巨大阻碍当含水饱和度逐渐升高毛管力随之减小水的流动阻力显著减少因此累计产水量发生大幅增加说明初始含水饱和度较低的页岩气藏更有益于气体生产储层优选时应优先考虑 结 论()采用无网格广义有限差分法能够较好地适应不同计算域情况收敛性好、计算稳定可对非线性方程中的非线性系数(如毛管力)进行准确的线性转换代替传统的数值微分方法()压降扰动能够波及改造区所有区域压降前缘能突破未改造区但在未改造区的传播非常有限 页岩基质具有更显著的应力敏感性固体位移主要发生在次改造区与未改造区的交界区域压降前缘仅能波及到该区域()流固耦合作用非常显著忽略应力场的耦合作用时在投产初期产气量会被显著高估最高可达 生产后期影响程度逐渐下降同时产水期大幅缩短产水量显著增大 水在基质区的流动性很差在主改造区的流动性较强()当含水饱和度较低时毛管力较高水相流动阻力较大当含水饱和度逐渐升高毛管力随之减小水相流动阻力显著减少产水量大幅增加初始含水饱和度是页岩气藏开发过程中重点考虑内容参考文献:邹才能赵群丛连铸等 中国页岩气开发进展、潜力及前景 天然气工业():():.:():亓倩 页岩气储层多级压裂水平井多场耦合非线性渗流理论研宄 北京:北京科技大学 :():():.:.:编辑 姜广义