CN108133080B - 一种考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法 - Google Patents

Abstract

本发明涉及一种考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法。一种考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法，包括如下步骤：（1）根据吸附气、自由气质量守恒方程，分别获得吸附气浓度和自由气浓度控制方程；根据裂缝中气相和水相质量守恒方程，获得裂缝中气浓度和水饱和度控制方程；（2）通过有限元方法，建立吸附气浓度、自由气浓度、裂缝气浓度和裂缝水饱和度有限元数值方程；（3）对上述有限元数值方程求解后得到任意位置任意时刻的吸附气浓度、自由气浓度、裂缝气浓度和裂缝水饱和度。本发明在常规离散裂缝模型基础上有四方面的改进，使其适应于页岩气的数值模拟。

Description

一种考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法

技术领域

本发明涉及一种考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法。

背景技术

近年来，页岩气产业在美国乃至全球都得到大力发展，数值模拟方法是油气产量预测的关键技术。目前，双重介质模型在页岩气藏数值模拟中应用较多，但是实际页岩气藏并不是任何区域都存在裂缝，这无疑夸大了裂缝的作用。双重介质模型中裂缝的方位和网格的走向相关，限制了实际裂缝方位的有效模拟。吸附气在页岩气藏中既可以传输还可以脱附转化为自由气，吸附气的传输一般通过表面扩散来研究，目前的研究通过将表面扩散的影响转化为对页岩气渗透率的影响，同时大多数学者通过兰格缪尔吸附模型，将吸附气加入到自由气控制方程中的源相，以表征吸附气解吸后对自由气的补充。但这两种方法无法表征任意位置和时间的吸附气的表面扩散情况和浓度的变化。

鄂尔多斯盆地中生界页岩气藏压力较低，在大型水力压裂后，半数以上气井的返排率在50％以下，裂缝中存在大量的压裂液未能返排出，对页岩气井的生产造成一定的影响，这类井产气量明显降低，并伴有压裂液的产出，但现有页岩气数值模拟方法中很少研究产出液的影响。

数值模拟技术是指导页岩气藏开发的关键技术。由于页岩中具有较多的天然裂缝，并且在压裂过程中产生压裂裂缝，且裂缝是页岩中流体传输的重要通道，裂缝的长度、开度、方位等会对页岩气的产量造成很大的影响。因此裂缝的模型是页岩气数值模拟中非常关键的步骤。目前，评价裂缝性油气藏的数值模拟方法主要有多重介质模型和离散裂缝模型。

虽然离散裂缝模型可以很好的解决复杂裂缝网络的数值模拟问题，但页岩气藏除了具有复杂裂缝外，页岩气藏存在大量的吸附气，吸附气对页岩气的产量贡献比较大，因此吸附气在孔隙壁面的传输，以及吸附气和自由气之间的耦合都制约页岩气数值模拟结果的准确性。目前，离散裂缝模型大多应用在常规油气的数值模拟预测中，在页岩气藏中应用很少，即使有应用也是在特殊情况下：如糜利栋提出的页岩气离散裂缝模型，只是针对单相气问题研究，但现有页岩气的产水问题比较严重，不容忽视；其次，吸附气的以源相的方式加入到自由气的控制方程中，这样不利于研究任意位置任意时间下吸附气的量化表征；再者，常规气藏中，自由气直接和孔壁接触，默认自由气在孔隙壁面的传输速度为0，页岩气藏中，自由气和吸附气接触，吸附气在吸附气浓度差下沿着孔隙壁面传输，因此自由气在与吸附气接触的边界处的传输速度不为0，这意味这吸附气的传输速度对自由气的传输造成一定的影响。离散裂缝模型不能解决以上三点问题。

发明内容

本发明旨在针对上述问题，提出一种改进型的离散裂缝模型，裂缝中存在气水两相传输，考虑吸附气对自由气传输的影响，将吸附气采用单独的控制方程表征并和自由气控制方程耦合，与页岩气的实际传输特征更接近的考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法。

本发明的技术方案在于：

一种考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法，模型假设如下：(1)考虑基质孔隙表面吸附气的传输，基质孔隙中存在单相自由气；(2)人工裂缝中存在气水两相流体，不考虑天然裂缝的影响；(3)气相为单组份，基质中气相微可压缩，裂缝中气水均微可压缩，基质和裂缝不可压缩，粘度恒定；(4)吸附气和自由气不能瞬时达到平衡；(5)气藏在生产过程中温度保持不变；(6)不考虑构造对气藏的影响，基质均质、等厚、各向同性；包括如下步骤：

(1)根据吸附气、自由气质量守恒方程，分别获得吸附气浓度和自由气浓度控制方程；根据裂缝中气相和水相质量守恒方程，获得裂缝中气浓度和水饱和度控制方程；

①吸附气以源相的方式加入到自由气的控制方程中，这样不利于研究任意位置任意时间下吸附气的量化表征；气体的吸附过程和解吸过程是同时进行的，自由气分子吸附在孔隙表面的同时，一部分吸附在孔隙表面的吸附气分子在脱离孔隙表面；单位时间内吸附在单位固体体积上的吸附量通过兰格缪尔吸附模型表示

β_ad＝k_a(1-θ)ρ_m (1)

其中

单位时间内在单位固体体积上解吸量可用Henry模型来表示

β_de＝k_dq_m (3)

单位时间内的净解吸附速率为：

β_net＝k_a(1-θ)ρ_m-k_dq_m (4)

β_net表示单位时间内的净解吸附速率，得到吸附气质量守恒方程：

其中

将式(4)和式(6)带入式(5)中，得到吸附气浓度控制方程

②自由气直接和孔壁接触，默认自由气在孔隙壁面的传输速度为0，页岩气藏中，自由气和吸附气接触，吸附气在吸附气浓度差下沿着孔隙壁面传输，因此自由气在与吸附气接触的边界处的传输速度不为0。基于此，根据毛细管层流物理模型得到考虑吸附气影响的自由气质量通量J_m，如下式，括号中第二项为吸附气的影响

其中，J_m为考虑表面扩散传输对自由气影响的气体质量通量

自由气质量守恒方程

将式(9)和式(10)带入(8)得到自由气浓度控制方程

③离散裂缝模型大多应用在常规油气的数值模拟预测中，在页岩气藏中应用很少，即使有应用也是在特殊情况下：只是针对单相气问题研究，并不能解决页岩气中严重的产水问题，但由于鄂尔多斯盆地陆相页岩气的压力较低，造成压裂液返排率较低，在整个生产过程中都伴有压裂液的产出，因此裂缝中存在气水两相流动，裂缝中气和水的质量守恒定律为

其中

气体和水的压缩系数分别为c_g和c_w，忽略裂缝中毛管压力，令C_ft＝S_fgc_g+S_fwc_w，得到裂缝***控制方程

整个研究区域外边界封闭，内边界定压。

(2)通过有限元方法，建立吸附气浓度、自由气浓度、裂缝气浓度和裂缝水饱和度有限元数值方程；

采用加权余量法，考虑页岩气藏外边界封闭，由(7)、(11)和式(18)得到吸附气浓度、自由气浓度、裂缝气浓度的加权余量方程分别为

基质区域采用三角型单元离散，裂缝区域采用一维线性单元离散，式(20)至式(22)进行分步积分，并应用Green原理，考虑本质边界条件，对时间采用向后差分格式，得到吸附气浓度、自由气浓度、裂缝气浓度和裂缝水饱和度的有限单元方程，将单元方程矩阵组合得到整体有限元方程，如下

其中

K_s ^e＝∫φ_f/Δt·Φ_fΦ_f ^Tdl (37)

(3)对方程求解后得到任意位置任意时刻的吸附气浓度、自由气浓度、裂缝气浓度和裂缝水饱和度的具体过程为：先对吸附气浓度、自由气浓度、裂缝气浓度总体有限元方程进行求解，获得任意节点上n+1时刻的吸附气浓度、自由气浓度以及裂缝气浓度，将n+1时刻的吸附气浓度、自由气浓度、裂缝气浓度带入总体裂缝水饱和度有限元控制方程中可求得任意节点上n+1时刻的裂缝含水饱和度。

对于待模拟的一口页岩气井，根据研究区域的大小和通过微地震资料取得的裂缝走向、方位、半长和开度等资料建立相应的地质模型。将地质模型中的基质区域划分为三角形网格，裂缝区域降维处理后划分为线型网格，得到基质三角形的单元数、单元节点数和单元节点坐标；裂缝线型单元数、单元节点数、单元节点坐标。以上这些参数是计算式(27)-式(38)所需的必要参数。按照单元节点编号和整体节点编号的关系将单元系数矩阵组合到整体方程的系数矩阵中得到总体有限元数值方程。如果基质模拟区域离散为m个节点，裂缝区域离散为f个节点，任意时间为n，将单元方程组装为整体有限元数值方程，这样便形成2m+f个方程和2m+f个未知数的整体矩阵方程式，具有组合后的整体有限元方程数值方程见下式(39)。通过对方程的求解，可得到任意节点上n+1时刻的吸附气浓度、自由气浓度、裂缝气浓度。之后将n+1时刻的吸附气浓度、自由气浓度、裂缝气浓度带入裂缝含水饱和度方程式可求得任意节点上n+1时刻的裂缝含水饱和度

公式中，k_a为吸附速率，s^-1；ρ_m为基质中气体浓度，kg/m³；θ为覆盖度，无因次；q_m为吸附气浓度，kg/m³；q_{m_max}为最大吸附气浓度，kg/m³；k_d为解吸附速率，s^-1；β_net为单位时间内净解吸附速率，kg/(m³·s)；J_ad为吸附气表面扩散质量通量，kg/(m²·s)；J_m为考虑吸附气表面扩散的自由气质量通量，kg/(m²·s)；δ_m-f为表基质和裂缝之间的关系，当基质节点和裂缝节点相邻，δ_m-f＝1，当基质节点和裂缝节点不相邻，δ_m-f＝0；ρ_fg为裂缝中气体浓度，kg/m³；φ_f为裂缝孔隙度，无因次；S_fg为裂缝气饱和度，无因次；J _fg和J _fw为裂缝气相和水相质量通量，kg/(m²·s)；k_f为裂缝渗透率，m²；k_fg为裂缝气相相对渗透率，无因次；k_frw裂缝水相相对渗透率，无因次；q_g和q_w为生产气量和水量，kg/(m³·s)；c_g和c_w分别为气体和水的压缩系数；δ_qm、δ_ρm和δρ_fg分别为吸附气浓度增量、自由气浓度增量和裂缝气浓度增量，kg/m³；Ω为积分区域。

本发明的技术效果在于：

本发明在常规离散裂缝模型基础上有四方面的改进，使其适应于页岩气的数值模拟：

一是在原有离散裂缝模型的控制方程中考虑了吸附气的影响。符合页岩气藏中存在吸附气的特征，而本发明模型考虑吸附气的影响，生产后期吸附气对产气量补充，产量下降较慢；二是将原有离散裂缝模型中的单相气拓展为气液两相，符合页岩气长期排液的现状，本发明模型初始阶段只产水，不产气，随着生产时间增长，产水量减少，逐渐开始产气。

三是吸附气的表征是通过引入了吸附气控制方程来实现的，传统的吸附气通过在自由气控制方程中加入源相作为吸附气的表征，对于这种情况，确定的兰格缪尔压力和兰格缪尔体积下对产量的影响是固定不变的，且吸附气和自由气在不同时间和不同位置各自的变化不可量化表征。但本发明模型可以计算任意位置任意时间吸附气的浓度变化。解吸速率和吸附速率共同影响着吸附气的净解吸附速率，从而影响页岩气的产量。

四是对页岩气的模拟，传统方法只考虑吸附气和自由气的脱附耦合，但忽略了吸附气在传输过程中可作为自由气传输的边界。本发明模型考虑表面扩散引起自由气变化时对产量的影响，其对产量的影响较大。

附图说明

图1为基质和裂缝网格划分示意图。

图2为模型对比示意图。

图3为自由气与吸附气浓度变化率示意图。

图4为净解吸附速率与裂缝的关系示意图。

图5为表面扩散传输引起自由气变化对产量的影响示意图。

图6为水平井模型网格划分示意图。

图7为YP井产量拟合示意图。

图8为单相气模型和气液两相模型相对误差示意图。

图9为YP井压力分布示意图。

具体实施方式

第一步：

选取一口页岩气井。建立边长为100m的正六边型地质模型，模型左上角和右下角对角线方向有一条长200m的裂缝，裂缝宽度0.3m，裂缝中心有一口页岩气井，如图1。基质划分为三角单元，单元数为54，节点数为37；裂缝划分为一维线性单元，单元数6，节点数7。模型计算所用各属性参数见表1。

第二步：

先求第一时间步长中各单元系数矩阵的计算通过表中参数计算式(27)-式(38)。

其余所有各式计算方法相同，不再重复。

第三步：

通过局部节点与整体节点的关系将系数以上所计算的系数矩阵进行组装，将每个单元的K_m1 ^e组装到整体系数矩阵中后得到37*37阶的矩阵。如表2、表3。同样将K_ρ1 ^e，f₁ ^e等所有单元矩阵都经过组装后得到相应的整体矩阵，这些矩阵也正是式(39)中的相关矩阵。通过求解式(39)的矩阵方程得到第一时间步长时的吸附气浓度、自由气浓度、裂缝气浓度，最后求得第一时间步长的裂缝气饱和度。如表4。

第四步：

以第一时间步中的各项计算结果为已知数，继续计算第二时间步的各参数(重复第二步和第三步的计算)，直至得到时间步长为2000天的吸附气浓度、自由气浓度、裂缝气浓度和裂缝气饱和度，如表5。

在以上步骤完成的基础上，可讨论本发明模型的计算结果和常规离散裂缝模型的计算结果，具体如下：本发明在常规离散裂缝模型基础上有四方面的改进，使其适应于页岩气的数值模拟：

一是在原有控制方程中考虑了吸附气的影响。符合页岩气藏中存在吸附气的特征，图2可以看出，Karimi单相气模型未考虑吸附的影响，因此在开发后期页岩气产量大幅度下降，而本发明模型考虑吸附气的影响，生产后期吸附气对产气量补充，产量下降较慢(这是由于吸附气和自由气之间的脱附耦合造成的)。

二是将原有模型中的单相气拓展为气液两相，符合页岩气长期排液的现状，可以看出，Karimi模型裂缝中只存在单相气，因此初期产气量最大，后期递减快，本发明模型初始阶段只产水，不产气，随着生产时间增长，产水量减少，逐渐开始产气。这是由于裂缝中气液两相流降低气体的流动性，因此裂缝中水(压裂液)的存在很大程度上会降低产气量；初始阶段裂缝中含水饱和度很高，含气饱和度低，水的流动能力较强，因此只产水，不产气，大量的压裂液排出，含水饱和度降低后，气逐渐开始流动，气井开始产气，如图2。

三是吸附气的表征是通过引入了吸附气控制方程(非稳态吸附)来实现的，传统的吸附气通过在自由气控制方程中加入源项作为吸附气的表征，对于这种情况，确定的兰格缪尔压力和兰格缪尔体积下对产量的影响是固定不变的，且吸附气和自由气在不同时间和不同位置各自的变化不可量化表征。但本发明模型在不同的解吸速率下有着不同的产量，解吸速率和吸附速率共同影响着吸附气的净解吸附速率，从而影响页岩气的产量，图2。

图3为自由气和吸附气浓度变化率。其中，ρ_m/ρ_{m_in}、q_m/q_{m_in}分别代表目前自由气浓度和初始自由气浓度比值、目前吸附气浓度和初始吸附气浓度比值。可以看出，离裂缝越远，ρ_m/ρ_{m_in}、q_m/q_{m_in}越大，且变化趋势逐渐平缓，这说明远离裂缝基质区域自由气和吸附气对产气量贡献较少。开发越久，ρ_m/ρ_{m_in}、q_m/q_{m_in}越小，越多的自由气和吸附气被采出。同样可以看出，对于同一区域，ρ_m/ρ_{m_in}小于q_m/q_{m_in}，且开发时间越久，二者之间的差值越大，这说明吸附气的解吸滞后于自由气浓度的变化，开发时间越长，滞后度越大。

图4为裂缝对净解吸附速率的影响，当页岩气井生产10天时，裂缝附近处基质的净解吸附速率较高，而远离裂缝处基质的净解吸附速率较低，这是因为，开发初期裂缝附近基质中自由气的浓度较低，净解吸附速率较大，而远离裂缝基质中的自由气浓度较高，净解吸附速率较小。当页岩气井生产1年时，曲线下降变平缓，裂缝附近基质的净解吸附速率在减小，而远离裂缝基质的的净解吸附速率在增大，二者差别在减小。当页岩气井生产10年时，曲线变化更平缓。以上分析说明，开发时间越久，离裂缝附近基质和远离裂缝基质处的净解吸附速率差别越小。

四是对页岩气的模拟，传统方法只考虑吸附气和自由气的脱附耦合，但忽略了吸附气在传输过程中可作为自由气传输的边界。本发明模型很好的解决了这一问题。图5为考虑表面扩散引起自由气变化时对产量的影响，计算所用表面扩散系数值为1×10^-4m²/s，图6中粉色线条为不考虑表面扩散引起自由气的变化，可以看出，其对页岩气产量的影响百分比最高为10^-4数量级，但如果考虑表面扩散对自由气的影响，其对产量的影响明显变大，最高可达1％。(这是由于吸附气和自由气之间的传输耦合造成的)

应用实例

水平井是目前页岩气开采的最佳方案，鄂尔多斯盆地中生界YP井实施大型压裂，共计8段，入地总液量17000m³，加砂总量727m³，累计返排液4200m³，返排率仅为25％左右，日平均产气3400m³，该井原始地层压力10MPa，压裂液氯根平均为7000mg/L，生产过程中的产出水氯根为7000mg/L左右，显示地层中不含可动水。建立长2000m，宽1200m的地质模型，模型单元数3469，节点数1750，8条裂缝，裂缝单元数160，裂缝节点数168。裂缝宽度、长度和方位按实际微地震监测结果给出，基质和裂缝网格划分见图6。不考虑水平井筒的压力降，水平井产量由每条裂缝产量累加得到。具体拟合参数见表1。具体计算过程在模型验证中已经详细演示，在此不重复累述。

图7为YP井实际生产数据与模拟数据对比。投产初期，裂缝中90％以上的流体为压裂液，大量的压裂液开始返排，产液量最高可达100m³/d，这时裂缝中的气体不可动，产气量为0m³/d；随着排液的进行，裂缝压力下降，基质中的自由气开始向裂缝中运移，裂缝含水饱和度也相应开始降低；排液三个月后，裂缝中的气饱和度达到一定程度，气相相对渗透率升高，开始产气，此时产液量迅速下降。产气量快速上升，并稳中有升，生产500天仍然比较稳定，没有下降的趋势，与海相页岩气井相比产量较低，但稳产程度较好。从模拟情况来看，如果采用单相流进行模拟，在同样的参数条件下，初期产气量较高，但产量下降较快，不符合陆相页岩的生产特征。本文建立的裂缝气水两相离散裂缝模型所模拟的产气量和产液量趋势都与实际生产数据较吻合，裂缝中的气液两相传输抑制了气体快速产出，使其稳定生产，能很好的反映真实页岩气的生产动态。

图8为裂缝单相气模型和裂缝气水两相模型计算结果与实际生产数据的偏差程度，其中单相气模型和气水两相模型计算结果与实际生产数据的相对误差的计算式为别为(q_sf)-q)/q和(q_mf)-q)/q，其中q_sf为裂缝单相气模型产气量，q_mf为裂缝气水两相模型产气量，q为实际产气量。可以看出，对裂缝单相气模型，平均误差为135％左右，平均拟合误差在一倍以上，如果用此模型模拟第二类页岩气井，会与实际不相符。而对于裂缝气水两相流模型，平均误差仅为1％左右，可以较准确的模拟第二类页岩气井。

图9为采用裂缝气水两相流模型模拟YP井时井周地层压力分布。图中黑色线条为井筒，图中的a、b、c和d分别为生产1月、1年、5年和10年时的地层压力分布。在生产1个月时，图a，整个模拟区域的压力基本不变，压力降低小于0.1MPa。在生产1年时，图b，离裂缝较近的区域压力稍有降低，而两条裂缝中间的压力基本没有降低；从左到右第二条和第三条裂缝的开度相对于其他裂缝较大，因此该两条裂缝周围稍远的区域也受到压力的波及，压力降低0.5MPa左右，而对于其他裂缝开度较小的裂缝，其周围区域的压力波及相对较小，压力降低0.2MPa左右；由于主裂缝的导流能力较强，因此裂缝长度方向上压力的变化不是很大。生产5年时，图c，压力波及范围开始增大，压力降低最多1.2MPa；第三条至第八条裂缝周围，由于裂缝间距相对较小，压力已波及到裂缝中间区域，两条裂缝中间区域压力降低0.5MPa左右；而第一条至第三条裂缝周围，压力还未完全波及到裂缝中间；由于第一条裂缝和第二条裂缝方位的影响，在裂缝北面两条裂缝间距较大区域，两条裂缝中间区域压力降低小于0.1MPa，在裂缝南面两条裂缝间距较小区域，两条裂缝中间压力降低小于0.3MPa。生产10年时，图d，整个压裂区域基本都受到压力波及，压力降低最多1.8MPa。

表1模型验证计算参数

表2第一时间步长中K_m1 ^e组合后得到的整体矩阵中1-20列计算结果

表3第一时间步长中K_m1 ^e组合后得到的整体矩阵中21--37列计算结果

表4第一时间步长中吸附气浓度、自由气浓度、裂缝气浓度和裂缝水饱和度计算结果

表5步长2000天时吸附气浓度、自由气浓度、裂缝气浓度和裂缝水饱和度计算结果

表6实例应用计算参数

Claims (4)

1.一种考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法，模型假设如下：(1)考虑基质孔隙表面吸附气的传输，基质孔隙中存在单相自由气；(2)人工裂缝中存在气水两相流体，不考虑天然裂缝的影响；(3)气相为单组份，基质中气相微可压缩，裂缝中气水均微可压缩，基质和裂缝不可压缩，粘度恒定；(4)吸附气和自由气不能瞬时达到平衡；(5)气藏在生产过程中温度保持不变；(6)不考虑构造对气藏的影响，基质均质、等厚、各向同性；其特征在于：包括如下步骤：

第一步骤：根据吸附气、自由气质量守恒方程，分别获得吸附气浓度和自由气浓度控制方程；根据裂缝中气相和水相质量守恒方程，获得裂缝中气浓度和水饱和度控制方程；

其中：根据吸附气质量守恒方程，获得吸附气浓度控制方程的具体实现步骤如下：气体的吸附过程和解吸过程是同时进行的，自由气分子吸附在孔隙表面的同时，一部分吸附在孔隙表面的吸附气分子在脱离孔隙表面；单位时间内吸附在单位固体体积上的吸附量通过兰格缪尔吸附模型表示

β_ad＝k_a(1-θ)ρ_m (1)

其中

单位时间内在单位固体体积上解吸量用Henry模型来表示

β_de＝k_dq_m (3)

单位时间内的净解吸附速率为：

β_net＝k_a(1-θ)ρ_m-k_dq_m (4)

β_net表示单位时间内的净解吸附速率，得到吸附气质量守恒方程：

其中

J_ad＝-D_s·▽q_m (6)

将式(4)和式(6)带入式(5)中，得到吸附气浓度控制方程

其中：根据自由气质量守恒方程，获得自由气浓度控制方程的具体实现步骤如下：

β_net表示单位时间内的净解吸附速率，得到自由气质量守恒方程：

页岩气藏中，自由气和吸附气接触，吸附气在吸附气浓度差下沿着孔隙壁面传输，因此自由气在与吸附气接触的边界处的传输速度不为0，根据毛细管层流物理模型得到考虑吸附气影响的自由气质量通量J_m，如下式，括号中第二项为吸附气的影响

将式(9)和式(10)带入式(8)得到自由气浓度控制方程

第二步骤：通过有限元方法，建立吸附气浓度、自由气浓度、裂缝气浓度和裂缝水饱和度有限元数值方程；

第三步骤：对上述有限元数值方程求解后得到任意位置任意时刻的吸附气浓度、自由气浓度、裂缝气浓度和裂缝水饱和度。

2.据权利要求1所述的一种考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法，其特征在于：所述的第一步骤中根据裂缝中气相和水相质量守恒方程，获得裂缝***数学模型的具体实现步骤如下：

离散裂缝模型大多应用在常规油气的数值模拟预测中，在页岩气藏中应用很少，即使有应用也是在特殊情况下只是针对单相气问题研究，并不能解决页岩气中严重的产水问题；但由于鄂尔多斯盆地陆相页岩气的压力较低，造成压裂液返排率较低，在整个生产过程中都伴有压裂液的产出，因此裂缝中存在气水两相流动，裂缝中气和水的质量守恒定律为

其中

气体和水的压缩系数分别为c_g和c_w，忽略裂缝中毛管压力，令C_ft＝S_fgc_g+S_fwc_w，得到裂缝***控制方程

整个研究区域外边界封闭，内边界定压。

3.据权利要求2所述的一种考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法，其特征在于：所述的第二步骤中建立有限元数值方程的具体过程为：

采用加权余量法，考虑页岩气藏外边界封闭，由式(7)、式(11)和式(18)得到吸附气浓度、自由气浓度、裂缝气浓度的加权余量方程分别为

基质区域采用三角型单元离散，裂缝区域采用一维线性单元离散；对式(20)至式(22)进行分步积分，并应用Green原理，考虑本质边界条件，对时间采用向后差分格式，得到吸附气浓度、自由气浓度、裂缝气浓度和裂缝水饱和度的有限元数值方程，将单元方程系数矩阵组合，得到吸附气浓度、自由气浓度、裂缝气浓度和裂缝水饱和度的总体有限元数值方程，

其中

K_s ^e＝∫φ_f/Δt·Φ_fΦ_f ^Tdl (37)

4.据权利要求3所述的一种考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法，其特征在于：所述的第三步骤对方程求解后得到任意位置任意时刻的吸附气浓度、自由气浓度、裂缝气浓度和裂缝水饱和度的具体过程为：先对吸附气浓度、自由气浓度、裂缝气浓度总体有限元方程进行求解，获得任意节点上n+1时刻的吸附气浓度、自由气浓度以及裂缝气浓度，将n+1时刻的吸附气浓度、自由气浓度、裂缝气浓度带入总体裂缝水饱和度有限元控制方程中可求得任意节点上n+1时刻的裂缝含水饱和度；

对于待模拟的一口页岩气井，根据研究区域的大小和通过微地震资料取得的裂缝走向、方位、半长和开度建立相应的地质模型；将地质模型中的基质区域划分为三角形网格，裂缝区域降维处理后划分为线型网格，得到基质三角形的单元数、单元节点数和单元节点坐标；裂缝线型单元数、单元节点数、单元节点坐标；以上这些参数是计算式(27)-式(38)所需的必要参数；按照单元节点编号和整体节点编号的关系将单元系数矩阵组合到整体方程的系数矩阵中得到总体有限元数值方程；如果基质模拟区域离散为m个节点，裂缝区域离散为f个节点，任意时间为n，将单元方程组装为整体有限元数值方程，这样便形成2m+f个方程和2m+f个未知数的整体矩阵方程式，具有组合后的整体有限元数值方程见下式(39)；通过对方程的求解，可得到任意节点上n+1时刻的吸附气浓度、自由气浓度、裂缝气浓度；之后将n+1时刻的吸附气浓度、自由气浓度、裂缝气浓度带入裂缝含水饱和度方程式可求得任意节点上n+1时刻的裂缝含水饱和度

公式中，k_a为吸附速率，s^-1；ρ_m为基质中气体浓度，kg/m³；θ为覆盖度，无因次；q_m为吸附气浓度，kg/m³；q_{m_max}为最大吸附气浓度，kg/m³；k_d为解吸附速率，s^-1；β_net为单位时间内净解吸附速率，kg/(m³·s)；J_ad为吸附气表面扩散质量通量，kg/(m²·s)；J_m为考虑吸附气表面扩散的自由气质量通量，kg/(m²·s)；δ_m-f为表基质和裂缝之间的关系，当基质节点和裂缝节点相邻，δ_m-f＝1，当基质节点和裂缝节点不相邻，δ_m-f＝0；ρ_fg为裂缝中气体浓度，kg/m³；φ_f为裂缝孔隙度，无因次；S_fg为裂缝气饱和度，无因次；J_fg和J_fw为裂缝气相和水相质量通量，kg/(m²·s)；k_f为裂缝渗透率，m²；k_fg为裂缝气相相对渗透率，无因次；k_frw裂缝水相相对渗透率，无因次；q_g和q_w为生产气量和水量，kg/(m³·s)；c_g和c_w分别为气体和水的压缩系数；δ_qm、δ_ρm和δρ_fg分别为吸附气浓度增量、自由气浓度增量和裂缝气浓度增量，kg/m³；Ω为积分区域。

Images