低渗透油藏CO2驱油数值模拟方法与应用

鞠斌山，于金彪，吕广忠，曹伟东

（1.中国地质大学（北京）能源学院，北京 100083；2.中国石化胜利油田分公司勘探开发研究院，山东东营 257015）

近年来CO2驱油在理论研究和矿场实践方面取得了较快发展［1-3］，前人研究表明CO2驱油可较大幅度地提高采收率［4-8］，相对于非混相驱，CO2混相驱具有更高的驱油效率［9］。对于CO2驱油提高采收率研究领域，前人研究主要集中在提高采收率机理［10-12］、混相压力预测方法［13-14］和室内实验［15-19］。对于实际油藏内CO2驱后驱油前缘推进、浓度分布和驱油指标，主要采用商业油藏模拟器开展研究［20-21］，但是商业油藏模拟器将CO2驱油处理为等温渗流，未考虑CO2在油层中的扩散作用。因此，对于CO2驱油数值模拟技术仍然有待于深入研究。此外，对非均质油藏CO2驱油过程中的“指进”现象、采油井开关机时机和顺序优化等关键问题也缺乏系统的研究。

正理庄油田樊142 块经过天然能量开发后，平均油层压力下降到不足20 MPa，油井产能大幅度降低，为了改善该区块的开发效果，于2013 年7 月开始选取樊142-7-斜4 注采井组的7 口井进行CO2驱油矿场试验。为了实现混相，将注入井周围的6 口井关闭，经过2 a 的CO2注入，生产井压力监测和数值模拟结果表明，该井组90%的油层区域达到混相驱条件（最小混相压力为31.0 MPa）。但是由于油层非均质性和注采井间压力梯度差异，出现了CO2驱油指进现象，如果所有生产井同时开井生产，势必造成个别井气窜严重，降低波及系数和驱油效果，影响采收率。因此非常有必要进行开井顺序优化，达到提高驱油效果的目的。为了实现最佳开井顺序，实现最大的驱油效果的目标，建立多组分CO2驱油数学模型，编制多组分气驱油藏数值模拟软件。利用该软件跟踪模拟该井组CO2驱油注入过程，根据CO2在储层中的分布特征和注采井对应关系，制定6个开发方案，并应用所研制软件，进行模拟与预测，得出最佳的开井顺序，预测CO2驱油指标。

1 CO2驱油多组分渗流模拟数学模型

1.1 假设条件

数学模型假设条件包括：①CO2驱油为多组分渗流，其组成为CO2，H2O 和n个原油拟组分。②考虑CO2在岩石孔隙表面上的吸附，可以更加准确地描述CO2在油层孔隙中起到驱油作用的比例。③考虑CO2在原油中的溶解作用，以便定量预测溶解量对原油物性的影响。④考虑油藏中流体动态相变，通过最小混相压力和溶解度参数判定驱油过程中的相变及相数。⑤考虑温度、原油组成对动态最小混相压力的影响，通过动态最小混相压力判别混相驱油还是非混相驱油，提高驱油相态的预测精度。⑥考虑非等温渗流，以便更加准确地预测注入井附近储层温度变化对渗流过程的影响。⑦在焖井或停注阶段，CO2扩散作用有助于扩大波及系数，因此需考虑CO2在油层中的扩散作用。

1.2 模型建立

根据质量和能量守恒定律，多相渗流的质量和能量守恒的数学描述通式为［22］：

1.2.1 质量项

质量项的表达式为：

其中：

1.2.2 质量交换速度

质量交换速度的表达式为［23］：

考虑温度和原油性质的最小混相压力关系式为：

2 数值离散与求解

2.1 差分离散

物质和能量守恒方程的数值离散形式为：

（8）式中，各相的质量流量可用达西定律进行描述：

其中：

（10）式中的γij定义为：

2.2 数值求解

采用全隐式差分离散格式，通常采用牛顿迭代方法求解，将（6）式写成余量［23］的形式为：

N个节点有（Nc+1）×N个非线性方程，每个节点有Nc+1个变量，需要用牛顿迭代方法求解Nc+1个方程［23］。

3 樊142-7-斜4 井组注CO2压力恢复过程数值模拟

采用所建立的CO2驱油数学模型和数值离散求解方法，在Windows 操作系统下，利用Visual Fortran语言编制的组分油藏数值模拟软件，对正理庄油田樊142-7-X4井组进行数值模拟。

3.1 樊142块概况

正理庄油田樊142 块油层的平均孔隙度为12.5%，平均渗透率为1.89 mD。地面原油密度为0.877 1 kg/L，地层原油密度为0.746 3 kg/L；地面原油黏度为20.4 mPa·s，地层原油黏度为1.18 mPa·s。原始油层压力为44.87 MPa，含油饱和度为0.72。樊142-7-斜4井组为CO2驱油矿场试验井组，包括1口注气井（樊142-7-斜4井），6口采油井，井位见图1。该井组于2013 年7 月至2015 年8 月注入CO2，日注入量达10～20 t/d。

3.2 樊142-7-斜4井组数值模拟模型

图1 樊142-7-斜4井组油层顶部深度等值线及井位Fig.1 Depths isograms of oil formation top and well locations of F142-7-X4 well group

目的层自下而上划分为4 个模拟层，第1—4 层命名为12层、13层、14层和38层，运用Petrol软件进行地质建模，x，y方向网格尺寸均为30 m，z方向网格厚度不同，其平均值为0.92 m；网格数为66×35×4共9 240个。樊142-6-2井、樊141-1井及樊142-6-3 井附近有效厚度较大，樊142-8-斜4 井有效厚度最小。

3.3 注CO2压力恢复阶段数值模拟结果

由于CO2混相驱效果明显好于非混相驱，为了达到CO2混相驱的目的，采取只注不采的方式，提高油层压力，以便油层压力高于最小混相压力。在该阶段，最关注的问题是压力分布特征是否具备混相驱油条件与前缘推进特征。平面压力分布（图2）表明：压力恢复阶段末期，樊142-8-3 井与樊142-7-3井附近区域受储层连通性差等因素影响，油层压力低于最小混相压力（31 MPa），未能达到混相条件，其余4 口井所在油藏的压力均高于31 MPa，达到混相驱油条件。CO2浓度前缘出现明显指进现象，驱油前缘向樊142-8-3 井、樊142-7-3 井与樊142-8-斜4 井方向推进速度较慢。CO2浓度前缘已经到达樊142-6-2 井，接近樊142-6-3 井和樊141-1-1 井，如果同时开井或开井顺序不当，则会在这3 口井造成气窜，降低整体波及系数和驱油效果。

4 6 口井采油方案设计及开井组合顺序优化

4.1 方案设计

为了获得最佳开井时机和开井顺序，根据6 口井的油层压力高低（表1），设计6 种开井顺序方案（表2），对于油层压力高于最小混相压力31 MPa的，设定井底流压为31 MPa，以保证混相驱，各井井底流压设定和压力恢复阶段末各井对应的油层压力见表1。模拟生产5 a，对不同预测方案进行数值模拟计算，得到各网格压力、CO2浓度及各生产开发指标等数据。根据综合前缘推进特征与CO2驱油效果评价对比，优选出最佳开井顺序方案。

图2 樊142-7-斜4井组压力恢复阶段末CO2浓度前缘与压力前缘分布Fig.2 Fronts of CO2concentration and pressure at the end of pressure building stage for F142-7-X4 well group

表1 各井井底流压与压力恢复阶段末油层压力Table1 Pressures of wellbore flow and formation pressure of each well at the end of pressure building stage

4.2 注采5 a后CO2浓度及前缘对比

开发至第1 800 d 时，方案二的CO2浓度前缘波及范围明显比其他方案的更大（图3）。结合CO2浓度分布及前缘推进特征可知，在相同开井数量下，优先开低压井（指井所在的油层压力低）比先开高压井CO2浓度前缘波及面积大，开发效果好。其中方案二与方案四波及面积分别排第1和2，可以作为推荐候选方案，最优方案的确定还需从各方案的驱油效果分析上做出最终评价。

4.3 驱油效果评价

为了对6个开井顺序方案的驱油效果进行评价和优选，对6 个开井顺序方案和1 个同时开井生产方案进行数值模拟，所有方案的CO2日注入量为15 t/d，连续注入5 a。将模拟计算得到的累积产油量、采出程度、CO2埋存率和换油率4个指标进行综合对比分析，确定最佳的开井生产方案。

累积产油量 在CO2注入速度和累积注入量相同的情况下，方案二的累积产油量最高，5 a 累积产油量为49 900×104t，与全部同时开井相比增产油量为5 407 t，其次是方案四（表3）。

表2 6种开井顺序方案Table2 Open sequence schemes of six wells

采出程度 由各方案的采出程度随时间变化（图4）可见，在第1 a 结束时方案六采出程度最高；然而1 a 后，各方案的采出程度排序发生了明显变化，方案二和方案四在第600 d 时，采出程度均超过方案六，从第3 a开始方案二采出程度始终保持最高（图4）。

图3 生产5 a后第14层CO2浓度前缘Fig.3 Fronts of CO2concentration of Layer 14th after 5 years

表3 累积产油量对比Table3 Comparison of cumulative oil productions

CO2埋存率 埋存率可间接地反映驱替效率。在CO2突破到生产井之间，埋存率为1，意味着注入的全部CO2占据油藏的孔隙空间，这部分空间的原油被驱替走；突破气窜之后，部分CO2从生产井采出，采出的这部分CO2没有留存在油藏内，埋存率降低，也意味着气窜程度加剧，因此埋存率越高越好。由所有方案的埋存率随时间变化（图5）可见，在第1 a内，埋存率为0.68～1.0，埋存率较高，随着生产时间延长，气窜导致埋存率下降。其中方案二在600 d后的埋存率最高，表明注入CO2利用率和驱油效率高，与该方案的累积产油量和采收率最高的结论吻合。

图4 采出程度对比Fig.4 Comparison of oil recoveries

图5 CO2埋存率对比Fig.5 Comparison of CO2sequestration rates

累积换油率 换油率从侧面反映CO2驱油的置换率，换油率随着时间是变化的（图6），在生产720 d后，方案二累积换油率最高，其次是方案四。

4.4 最佳开井顺序优选

最佳开井顺序需要综合考虑2 个方面：①在前缘推进特征方面，同时满足CO2波及范围大，而且地层压力较高的方案为优选标准。从以上分析得知，方案二和方案四2 个方案较好，可作为候选方案。②在驱油效果方面，兼顾采出程度、累积产油量、埋存率、增油量和累积换油率较高（表3，图4—图6），方案二为推荐的最佳方案。

图6 换油率对比Fig.6 Comparison of oil replacement rates of CO2

4.5 矿场实施效果

根据坝砂3 口井的试采结果，注CO2后开井求产期单井自喷日产油量为5.6～7.6 t/d，平均单井日产油量为6.5 t/d；周期注采阶段平均单井日产油量为5.8 t/d，产能远大于注气前泵抽的1.8 t/d，较试验前提升2.8～5.1 倍，平均增产3.5 倍，自喷累积产油量为961 t。为了更加客观地反映注入CO2后油井增产能力，采用采油指数进行效果评价（表4）。注CO2后开井，3 口井采油指数明显增大，平均为原来的4.76 倍，其中樊141-1 井的采油指数为原来的6.93倍，表明注CO2后，采油能力大幅度提高。

5 结论

利用多相多组分渗流理论，建立了CO2驱油渗流模型，进行了全隐式离散，并采用牛顿迭代方法实现了求解。

樊142 块樊142-7-斜4 注采井组由于油层非均质性和注采井间压力梯度差异导致了CO2在油层注入过程中出现了明显的指进现象，预测结果表明CO2首先在樊142-6-2井突破，其次是樊141-1井和樊142-6-3井。

数值模拟结果表明：生产井开井顺序对采出程度、CO2埋存率、换油率和CO2波及程度影响显著，根据CO2驱油前缘和浓度分布特征进行开井时机和顺序优化有助于提高CO2驱油效果。

表4 注CO2前后油井的采油指数对比Table4 Comparison of production indexes of oil production wells before and after CO2injection

通过6个开井顺序方案和全部开井方案的模拟结果对比与分析，处于低压区的樊142-8-3 井和樊142-7-3 井先开井，生产1 a 后其余4 口井再开井生产为最佳方案。生产5 a 情况下最佳方案比同时开井方案的采收率高1.34%，增产油量为5 407 t。

符号解释

Ak——某单元体内组分k（k≤Nc）的每单位体积质量项，kg/m3；或某单元体内每单位体积的热量，kJ/m3；上标k——组分数，当k=1，2，3，…，Nc时，（1）式为每个组分（共Nc个组分）的质量守恒方程；当k=Nc+1 时，（1）式为能量守恒方程；Nc——组分总数；t——时间，s；Fk——质量交换速度，kg/s；或热量交换速度，kJ/s；qk——某单元体内的每单位体积的质量源汇项，kg（/s·m3）；或某单元体内的每单位体积的热量源汇项，kJ（/s·m3）；ϕ——孔隙度；β——相标志，其为g，w 和o 时分别表示气相、水相和油相；ρβ——相密度，kg/m3；Sβ——相饱和度；——组分k在相中的摩尔分数；——组分k在单位体积岩石孔隙表面的吸附量，kg/m3；下标S——固体岩石；ρS——岩石颗粒密度，kg/m3；——β相分配系数，m3/kg；vβ——渗流速度，m/s——综合扩散系数，m2/s；pmm——最小混相压力，MPa；——C7+相对分子质量；TRC——油藏温度，℃；fRF——C2—C4摩尔组成，%；——第i节点第k组分在n+1时刻单位体积的质量，kg/m3；——第i节点第k组分在n时刻单位体积的质量，kg/m3；Vi——第i节点体积，m3；下标j——节点个数；ηi——第i节点的外表面；——第i节点第k组分在n+1 时刻的流入速度，kg/s；——第i节点第k组分在n+1 时刻的源汇项，kg/s；——组分k的对流项，kg/s——组分k的扩散项，kg/s；Fβ,ij——β相的对流项，kg/s；nij——相邻两节点i和j之间的单位向量；Aij——相邻两节点i和j之间的接触面积，m2；——β相的流度系数，（kg/（m3·Pa·s）；γij——导流系数，m3；Ψβj——第j节点β相的压强，Pa；Ψβi——第i节点β相的压强，Pa；Krβ——β相的相对渗透率；μβ——β相的黏度，Pa·s——第i和j节点的平均渗透率，m2；Di——第i节点中心到接触面的垂直距离，m；Dj——第j节点中心到接触面的垂直距离，m——余量，kg/s。