0 引言
产量不稳定分析方法（RTA）可用于生产井产量预测、可采储量评估和增产措施优化[4-5]。常规储集层多相流不稳定压力分析方法[6-7]为非常规储集层多相流产量不稳定分析提供了重要理论基础。Perrine[8]最先在1956年用总流度替换单相流度以考虑多相流动。Fetkovich[9]在边界控制流方程中提出用两相拟压力模拟两相流动。Fevang等[10]提出凝析气藏拟压力方程积分的劈分法，这种方法在常规凝析气藏中应用广泛，但是在致密凝析气藏中会产生极大的误差。现有的凝析气井产能评价方法大都基于储集层稳态渗流假设[11]或采用数值模型进行生产动态分析[12]，只考虑径向流或边界控制流，并不适用于非常规油气开发中常见的长期线性流，而且在非常规油气开发RTA中考虑多相流动尤为困难。近年来，非常规油气RTA技术迅速发展[13]，Hamdi等[14]提出了一种半解析方法用于分析致密油三相瞬态线性流条件下的产量数据，并使用矿场实例数据演示了RTA的工作流程，但并未考虑边界控制流对油井产能的影响。Qanbari等[15]基于动态泄油面积的概念提出了一种富液致密油/页岩油RTA方法，用来分析不稳定线性流产量数据，但该方法针对早期生产数据，并不适用于气井长期生产数据分析。Behmanesh等提出了致密凝析气藏不稳定线性流和边界控制流的RTA分析方法，其研究表明线性流系统中饱和度和压力的关系不能简单地用传统油藏开发理论表征[16⇓-18]，但并未考虑全流态生产数据分析问题。Li等[19]提出了一个半解析的RTA模型用于分析具有非均质性和多相流特征的挥发油藏生产数据，但只研究线性流态，并未考虑边界控制流。吴永辉等[20]建立了页岩凝析气三线性流产能预测模型，但并没有考虑裂缝间边界控制流对产能的影响。Zhang等[21]提出了一种考虑多相渗流的早期生产数据分析方法，计算了流体饱和度和压力的关系，但并不适用于气井中晚期生产数据分析。Luo等提出了一种用于低渗透挥发性油藏生产数据分析的RTA方法[22-23]，强调全流态分析并提出了挥发油藏中流体饱和度和压力的经验关系式，但并不适用于致密凝析气藏的拟压力计算。
本文考虑凝析气藏相态变化、储集层中油气两相线性渗流和边界控制流，提出了致密凝析气井全路径油相饱和度与压力关系的预测方法，建立了致密凝析气井多相流不稳定产量预测模型。最终通过与数值模拟对比以及矿场实例应用，验证模型的可靠性和实用性。
1 致密凝析气井渗流模型
1.1 物理模型
为模拟致密凝析气藏多级压裂水平井动态生产过程，建立了油气两相渗流物理模型（见图1 ）。基本假设如下：①储集层均质等厚，渗透率和孔隙度均匀分布，流体流动遵循达西定律。②压裂后各簇水力裂缝的几何形态和供气能力均相同，裂缝完全穿透储集层，每个裂缝区域均对称，忽略改造区以外的流体供给。③水力裂缝具有无限导流能力，生产初期裂缝间线性流为主要流态。④储集层中的凝析油气均为微可压缩流体，气体满足真实气体状态方程，高压物性（PVT）参数随压力变化。
1.2 数学模型
2 凝析气藏含油饱和度和压力关系计算
So-p关系受多种因素的影响，如生产过程中的压力变化、地层条件、井的生产制度（定压或定产）、流态分布（线性流或径向流）以及相态变化等。稳态流动模型、物质平衡模型、等组分膨胀实验（CCE）和等容衰竭实验（CVD）等方法均可以确定饱和度和压力关系，但不同方法适用于不同气井生产阶段。将上述方法与数值模拟得到的So-p关系进行对比分析，寻找凝析气井最优的So-p关系路径组合，从而获取准确的拟压力计算结果。所用凝析气样品组分组成如表1 所示。
表1 凝析气组分表 |
| 凝析气组分 | 占比/% | 凝析气组分 | 占比/% |
|---|---|---|---|
| CO2 | 8.05 | nC4 | 0.26 |
| N2 | 0.48 | iC5 | 0.51 |
| C1 | 58.85 | nC5 | 1.18 |
| C2 | 6.57 | C6 | 0.76 |
| C3 | 13.35 | C7+ | 8.96 |
| iC4 | 1.03 |
2.1 常用So-p关系研究方法
稳态流动模型：两相稳态流动期间，近井筒区域油气同时流动，整个区域的生产气油比（GOR）保持恒定，该阶段So-p关系可采用稳态流动模型描述。从相渗曲线上可以获取相对渗透率和饱和度的关系，从PVT实验中可以获取流体物性参数和平均地层压力的关系，将这两种关系组合到一起即可得到稳态流动模型的So-p关系，可以表示为[22]：
CCE实验用来模拟流体在恒定温度下的衰竭过程，测试期间总成分保持不变，是一种用来测试PVT参数的方法。
CVD实验也是PVT参数测试方法，在实验中，每次容器压力降低后都移除固定体积的气体，该过程接近于气藏的衰竭过程。
2.2 So-p关系数值模拟研究
为优选最佳So-p关系计算方法，假设平面对称的裂缝几何形状相同，每个裂缝区域的性质均相同，构建具有M条相同双翼裂缝的多级压裂水平井数值模型（见图2 ）。为了减少模拟运行的时间，选择了单个裂缝中1/4裂缝区域作为单元模型，该井的产能为模拟单元区域产能的4M倍。本文使用tNavigator数值模拟软件进行模拟，裂缝半长70 m，裂缝间距100 m，地层孔隙度10%，基质渗透率0.01×10−3 μm2，储集层有效厚度20 m。模型PVT参数根据表1 所示凝析气组分，通过PR（Peng-RobinSon）状态方程经过闪蒸计算模拟得到。
储集层的相对渗透率数据使用广义Brooks-Corey型幂律方程近似，对于油气系统有：
设凝析气藏的露点压力等于原始地层压力，绘制从原始地层压力下降到井底压力过程中的So-p关系曲线。采用数值模拟软件生成CCE、CVD模拟实验结果，基于数值模拟预测生产数据采用物质平衡模型和稳态流动模型进行解析计算，4种方法得到的So-p关系如图3 所示，可见不同方法计算得到的So-p曲线差异明显，说明选择合适的计算方法非常重要。
由图3 可见，稳态流动模型计算结果明显偏离实际情况，该方法不予考虑。CCE模拟实验、CVD模拟实验数据的变化趋势与凝析气藏的温度和凝析气组分相关，不同的温度和组分会生成不同的So-p关系曲线。物质平衡模型描述了地下流体不断被采出的过程，反应了储集层平均压力与平均饱和度的关系，预测的趋势大致与模拟实验数据一致，但预测结果往往大于模拟实验结果。
物质平衡模型的计算需要PVT数据和生产动态数据，从而限制了该模型的适用性，因此该模型不能用于原始地层压力到井底流压范围的拟压力计算，不能预测整个生产阶段的So-p关系。CCE、CVD实验过程与实际储集层开采过程不完全相符，会导致So-p的预测误差，故需优选最佳的So-p计算方法。
2.3 So-p关系最优组合路径
将物质平衡模型、CCE实验、CVD实验与数值模拟得到的So-p关系进行对比分析，寻找凝析气井最优的So-p关系组合路径。设计15个模拟方案（见表2 ），采用上述方法获取致密凝析气井压力-温度（简称“p-T”）相图与So-p关系曲线。
表2 模拟方案基础数据表 |
| 模拟方案编号 | 模拟条件 | 温度/℃ | 原始地层 压力/MPa | 井底流压/ MPa |
|---|---|---|---|---|
| 1 | 恒定井底 流压 | 130 | 25.15 | 5.00 |
| 2 | 150 | 24.81 | 5.00 | |
| 3 | 180 | 23.49 | 5.00 | |
| 4 | 210 | 21.12 | 5.00 | |
| 5 | 240 | 17.50 | 5.00 | |
| 6 | 恒定井筒处 凝析油比例 | 130 | 25.15 | 4.14 |
| 7 | 140 | 25.04 | 4.96 | |
| 8 | 160 | 24.48 | 6.97 | |
| 9 | 180 | 23.49 | 9.72 | |
| 10 | 200 | 22.04 | 16.18 | |
| 11 | 恒定温度 | 180 | 23.49 | 20.00 |
| 12 | 180 | 23.49 | 15.00 | |
| 13 | 180 | 23.49 | 12.50 | |
| 14 | 180 | 23.49 | 10.00 | |
| 15 | 180 | 23.49 | 5.00 |
设井底流压为5 MPa，设计恒定井底流压模拟方案1—方案5（见图4 ）。方案1原始储集层温度、压力条件离临界点最近，方案5原始储集层条件离临界点最远。方案1靠近临界点，凝析油含量较高，物质平衡模型与数值模拟结果匹配程度较高。方案5远离临界点，凝析油含量较低，当压力较高时物质平衡模型与模拟结果匹配程度较高，当压力较低时物质平衡模型与模拟结果出现偏差，但偏离程度较小。CVD模拟实验与模拟结果偏差均较大，但在靠近井底地层压力较低地带CCE模拟实验与模拟结果具有较高的一致性。综合5次模拟结果可知，地层压力与原始地层压力比较接近时，适合采用物质平衡模型预测So-p关系；近井区适合采用CCE实验预测So-p关系。
随着生产的进行，凝析油会在裂缝处大量聚集，设计恒定凝析油比例方案6—方案10，模拟方案均从原始地层压力下降至90%等气量线，凝析油含量随着地层压力的下降而逐渐增加，随后出现反凝析现象，凝析油蒸发，油相含量逐渐下降，气相含量逐渐提升，最终体系内油相比例稳定在10%（见图5a ）。对比方案6—方案10的So-p关系发现，物质平衡模型与模拟结果在整个压力变化范围内具有很好的一致性（见图5b 、图5c 、图5d ），凝析油比例恒定时，适合采用物质平衡模型预测So-p关系。
通过上述15个模拟方案的对比分析可知：物质平衡模型可以较准确地表征生产初期或离井底较远地层压力较高区域的So-p关系，而CCE实验数据则可较好预测生产后期或离井底较近地层压力较低区域的So-p关系。故实际应用中可利用实际生产井前期的生产动态数据，根据物质平衡模型计算生产初期或离井底较远地层压力较高区域的So-p关系，当物质平衡模型计算结果出现偏离时，利用CCE实验数据计算生产后期或离井底较近地层压力较低区域的So-p关系，然后根据两段计算所得数据利用分段三次Hermite插值多项式求取中间未知部分的So-p关系，形成完整的全路径So-p关系。图7 为方案2凝析气藏全路径So-p关系预测结果，可以看到，预测结果与模拟结果除局部存在一些小偏差外，整体匹配程度较高，该方法可以表征生产过程中的全路径So-p关系。
3 不稳定产量预测方法及应用
将本文提出的致密凝析气井不稳定产量预测方法与产能数值模拟结果对比验证，并将新方法应用于现场实例中，验证新方法的可靠性和实用性。
3.1 方法验证
验证方案模型原始地层压力等于露点压力，储集层参数采用模拟方案3中的基础数据，PVT实验数据采用tNavigator软件模拟生成。方案3全路径So-p关系预测结果见图8a ，随后根据相渗曲线可求得相对渗透率和含油饱和度的关系，最后预测凝析油气产量并绘制产量变化曲线。为验证方法的可靠性，采用模拟方案3中的PVT数据和储集层参数进行了产能数值模拟计算（见图8 ）。对比本文方法预测结果和数值模拟结果产量曲线可以看到，两种方法得到的结果吻合度非常高，说明新方法可靠，可用于预测致密凝析气井多相流不稳定产量。同时从图8b 、图8c 中还可以看到，在早期瞬态线性流动期间，产量曲线具有一段−1/2斜率的直线，随着生产的进行，压力波传播到位于两个相邻裂缝中间的干扰点时，产量快速下降，曲线表现出明显的边界控制流特征。
3.2 应用实例
渤海油田致密凝析气藏多级压裂水平井A井已投产约1.5年，原始地层压力49.40 MPa，露点压力40.07 MPa，地层孔隙度11.2%，储集层有效厚度56 m，基质渗透率0.015×10−3 μm2，岩石压缩系数1.26×10−4 MPa−1，共压裂25段，各段间距约60 m，裂缝半长约164 m。流体物性参数和油气相渗曲线均来自实验室测试报告。全路径So-p关系曲线采用本文方法进行计算，根据本文计算流程，预测并绘制了凝析油气产量与物质平衡时间的关系曲线[26]（见图9 ）。由图可以看出，物质平衡时间大于90 d时，实测数据与预测曲线吻合良好，物质平衡时间大于1 600 d时，压力波将传播到位于两个相邻裂缝中间的干扰点，此时产量曲线下降速度加快，表现出边界控制流的特征。实例验证说明新方法具有较好解决实际问题的能力，可以预测致密凝析气井生产中、晚期产量，评估单井可采储量。
4 结论
凝析油饱和度和压力关系曲线是计算拟压力的关键，在生产初期或离井底较远地层压力较高区域，可利用实际生产井前期的生产动态数据，采用物质平衡模型进行计算；生产后期或离井底较近地层压力较低区域，可采用CCE实验数据进行计算；生产中期或地层压力位于中间水平时，则可根据前两段计算所得数据插值求得，形成完整的全路径油相饱和度和压力关系曲线。
经模拟对比与矿场实例的验证，本文方法具有较好的可靠性和实用性，可用于致密凝析气井生产中、晚期产量预测和单井可采储量的评估。
符号注释：
Bg——气相体积系数，m3/m3；Bo——油相体积系数，m3/m3；c——广义压缩系数，MPa−1；h——储集层厚度，m；G——天然气原始地质储量，m3；Gp——累计天然气产量，m3；GOR——生产气油比，m3/m3；K——储集层平均渗透率，10−3 μm2；Krg——气相相对渗透率，无因次；Kro——油相相对渗透率，无因次；Krg,max——最大气相相对渗透率，无因次；L——裂缝间距的一半，m；M——裂缝条数；n——求和变量，无因次；ng——油气系统中气相相对渗透率曲线的指数，无因次；nog——油气系统中油相相对渗透率曲线的指数，无因次；N——石油原始地质储量，m3；Np——累计原油产量，m3；p——压力，MPa；pi——原始地层压力，MPa；ppg——气相拟压力，MPa/（mPa·s）；ppg,i——原始地层拟压力，MPa/（mPa·s）；ppg,wf——井底拟压力，MPa/（mPa·s）；ppo——油相拟压力，MPa/（mPa·s）；pref——参考压力，MPa；pwf——井底压力，MPa；qgsc——凝析气产量，m3/d；Rs——溶解气油比，m3/m3；Rv——凝析油气比，m3/m3；Sg——凝析气饱和度，%；Sgc——临界含气饱和度，%；So——凝析油饱和度，%；Srog——油气系统的残余油饱和度，%；Swc——临界含水饱和度，%； ——平均含油饱和度，%；t——生产时间，d；T——温度，℃；x——x方向距离，m；xf——裂缝半长，m；ε——单位转换系数，取值0.086 4，无因次；λgsc——气相组分流度，1/（mPa·s）；μg——天然气黏度，mPa·s；μo——原油黏度，mPa·s；ϕ——孔隙度，%。