基于聚类分析的降阶模型在低渗透油藏数值模拟中的应用

张艳玉，李威威，陈会娟

（中国石油大学（华东）石油工程学院，山东青岛 266580）

文章编号：1001⁃246X（2015）05⁃0603⁃07

基于聚类分析的降阶模型在低渗透油藏数值模拟中的应用

张艳玉，李威威，陈会娟

（中国石油大学（华东）石油工程学院，山东青岛 266580）

基于全隐式法求解考虑低渗透油藏流动过程中启动压力梯度影响的二维油水两相流动数学模型；利用本征正交分解法从获取的快照中提取基函数，将原模型投影到基函数展成的低维子空间中构造相应的降阶模型；考虑到依据等时间间隔获取的快照的缺陷，利用CVT⁃Lloyd算法对快照进行聚类分析.实例计算表明：通过本征正交分解法构造降阶模型能很好地对原模型进行近似；通过基于聚类分析处理快照后，一方面减少数据冗余，生成空间分布更均匀的快照集合从而提高降阶模型的精度，另一方面过少的聚类数可能导致快照信息的缺失，影响最后模型的精度.

低渗透油藏；聚类分析；降阶模型；本征正交分解

0 引言

根据国内外大量的研究，流体在低渗透介质中的流动不符合达西渗流规律，存在着启动压力梯度［1－3］.在油藏数值模拟中由于启动压力梯度的引入，使得方程的非线性增强.在利用差分法求解方程的过程中，目前通常采用全隐式法来提高计算的稳定性，但该方法计算量大，求解所需时间较长［4］，为进一步的生产制度优化和历史拟合造成一定困难.为了解决这个问题，考虑利用模型降阶的思想用以加快模拟速度.降阶模型是一类有效近似大型系统的方法，近几年已被成功地应用于许多高新工业技术应用领域.目前本征正交分解法（proper orthogonal decomposition，POD）是构造非线性系统降阶模型最流行的方法之一.该方法的基本思想是首先通过实验或者数值模拟获取快照并从中提取基函数，而后通过投影法将原始系统的控制方程投影到优选出的少数主要基函数上从而得到相应的降阶模型［5－8］.然而快照一般采用等时间间隔法获取，这样一来可能出现的问题是快照在空间分布可能不均匀，这意味着有些空间部位的快照数相对集中而其它部位则分布稀疏.2006年，Burkardt［9］引入CVT（Centroidal Voronoi Tesselation）用于构造降阶模型的基函数.运用该方法有两点好处，首先它提供了一个空间上分布更加统一的快照，因此可以有效地减少基函数的个数.其次进行CVT聚类分析后的矩阵在进行奇异值分解（Singular Value Decomposition，SVD）分解时效率更高.2009年，Cardoso［10］将聚类分析用于常规油水两相渗流的油藏数值模拟，指出通过对快照进行聚类分析可以有效删除数据的冗余，降低降阶模型构造的难度.

本文建立考虑启动压力梯度的低渗透油藏二维油水两相数学模型，并利用全隐式有限差分法对数学模型进行求解，过程中通过本征正交分解（POD）法构造原系统的降阶模型加快计算速度，最后通过实例分析计算研究基于POD法的降阶模型在低渗透油藏数值模拟中的计算效率以及精度.

1 低渗透油藏数值模拟数学模型

对于低渗透油藏二维油水两相渗流，根据质量守恒方程建立模型［1］：

辅助方程

式中，pcow为油水两相毛管力，0.1MPa.

初始条件

外边界条件

式中，n为油藏外边界G的外法线方向.

内边界条件

2 数学模型的求解及降阶模型的构造

结合式（1）－（5），基于五点差分格式全隐式法进行求解.在构造降阶模型的过程中，首先以一定生产制度运行全隐式法求解的模型，以固定时间间隔记录下模拟过程中每一时刻下各个网格的压力和饱和度（即快照）［11］，再利用CVT⁃Lloyd算法对快照进行聚类分析，而后结合本征正交分解法的思想求出基函数，将基函数和全隐式法求解过程相结合，即可得到相应的降阶模型，具体过程如下.

2.1 全隐式求解

将方程（1）左右两端乘以单元体体积VB＝ΔxiΔyjh并引入线性差分算子，简化方程为

计算过程中流度和启动压力梯度按“上游权”选择［4］，综合辅助方程（2），水相差分方程表示为

其中，l为每个时间步内的迭代次数；t为时间，s；ϕ0为地层初始孔隙度，无因次；CR为岩石压缩系数，0.1MPa－1；

产量可以表示为［12－13］

其中，re为Peaceman等效半径，cm；rw为井眼半径，cm；pBH为井底流压，0.1MPa；pe为外边界压力，0.1 MPa；Ql为地面条件下产量，cm3·s－1；S为表皮系数；

对定井底流压条件进行全隐式求解时，油水产量可以分别表示为

其中，WI＝（2πhKe）/［ln（re/rw）＋S］，λl＝Krl/（Blμl）．

考虑Nc个网格的油水两相流动，全隐式法求解的具体过程参考文献［4］，可得

其中，J为雅可比矩阵；δx＝［δpo1，δSw1，δpo2，δSw2，…，δpoNc－1，δSwNc－1，δpoNc，δSwNc］T，

g ＝ ［go1，gw1，go2，gw2，…，goNc－1，gwNc－1，goNc，gwNc］T．

2.2 基于本征正交分解法的降阶模型构造

本征正交分解法的基本思想就是将高维空间模型投影到低维空间模型之中.具体来说，首先通过实验或者计算实例得到的快照构造相关矩阵，利用求得的相关矩阵的特征向量以及对应特征值大小优选出少数正交基函数来刻画原模型.而后通过将原始模型投影到POD基函数上得到降阶模型，从而减少计算过程中的未知数数量［5，10，14］.

2.2.1 获取快照

由于压力和饱和度物理性质差异很大，因此将二者分开处理，但处理的过程类似.考虑Nc个网格中油水两相流动，定义状态向量：x＝［po1，…，poNc－1，poNc］T或［Sw1，…，SwNc－1，SwNc］T（以下统一用x表示）．按照

2.2.2 快照的聚类分析

根据以上步骤，将所得质心点集合作为新的快照集合参与下一步运算.

2.2.3 基函数的建立

针对处理后快照形成的矩阵，构造相应的相关矩阵C＝XTX．依据文献［4］所提出的方法通过直接对X进行奇异值分解可以得到计算相关矩阵C的特征值λ和特征向量Ψ．

基函数可表示为快照的线性组合Φ ＝XΨ．将相关矩阵特征值依次递减排列，即：λ1≥λ2≥λ3≥λ4….根据文献［17－18］所述筛选主要的特征基函数，可得压力保留的基函数个数设为lp，饱和度的基函数个数设为ls．将上一部筛选出的特征值所对应的基函数构成新矩阵Φl，它包含l（l＝lp＋ls≪2Nc）列，这样就将Φ降阶至Φl．

2.2.4 将POD基函数和原模型求解过程结合得到降阶模型

将降阶后的状态向量定义为z，全阶模型下的油水状态向量x就可以表示为

将（12）式带入（11），并将等式两端同时乘以Φl，可得牛顿迭代法的降阶表达形式：

3 实例计算

利用某油田区块实际地质资料建立概念模型对前文所述方法进行验证.其中基础数据为：油藏顶部深度1 450m，油层厚度为20m，孔隙度0.3，平均渗透率为0.01μm2，初始含油饱和度0.7，油水启动压力梯度为0.0009MPa·m－1.采用五点井网，网格数为Nc＝25×25×1，注水井定井底流压为14MPa，生产井定井底流压为9MPa，地层原始压力为12MPa.

3.1 降阶模型的建立

由于从基础模型中获取快照的质量直接影响到最后所生成的基函数能否精确反映原模型物理量的变化，此处采用随机改变生产井井底流压的生产制度运行程序获取快照.如图1所示，对四口生产井（井底流压相同）每隔1小时随机改变其井底流压（但保证其在8.5MPa到9.5MPa之间）并保持该压力生产一小时，而注入井保持井底流压为14MPa.

图1 生产井井底压力、产水量和产油量随时间变化Fig.1 Evolution of bottom hole pressure，water production rate and oil production rate

在上述生产制度下，以等时间间隔获取600张快照，依据CVT⁃Lloyd算法分别取聚类数为600、500、200、100、50、20和10对快照进行预处理，最后利用上文提出的本征正交分解法构造基函数，最后建立原模型对应的降阶模型.对于建立的降阶模型，采用基础模型中使用的生产制度（由于井底流压随机改变，因此每种情况具体值都不同），比较原模型和降阶模型计算精度.

3.2 结果分析

由图2可知，整体上经过快照预处理所构造的降阶模型对原模型都有较好近似.随着聚类数的减少，降阶模型和原模型在相同生产制度下的偏差逐渐变大.为了进一步定量比较原模型和降阶模型的区别，定义降阶模型相对原模型的误差计算公式

其中：ns为快照数；为第i快照时间处原模型产量，cm3·s－1；为第i快照时间处降阶模型产量，cm3· s－1；为原模型的平均日产量，cm3·s－1.

图2 聚类数为500、200、100、50、20和10时生产井井底压力和产油量随时间变化Fig.2 Evolution of bottom hole pressure and oil production rate as number of clusters is set at 500，200，100，50，20 and 10

结合图2和图3可知，一方面随着聚类数进一步减小，模型的误差不断增加.另一方面随着聚类数的增加，模型的误差也略微增大.产生此现象的原因主要是在获取快照的过程中，采用间隔固定的时间步长获取，这样获取的在时间上均匀分布的快照可能在空间上不均匀分布.当聚类数不断减小时，原本不同类型的快照合并导致信息缺失，最终由此构造的降阶模型无法有效近似原模型.而当聚类数偏大时，原本相似的快照可能被聚并成多类，那么最终生成的基函数就会倚重某些快照，造成误差增大.相比于直接利用快照构造降阶模型，对快照进行一定大小的聚类一方面可以降低数据的冗余，防止生成的降阶模型偏重于某些状态空间，提高模型的精度，另一方面可以减小相关矩阵维数，从而降低降阶模型构造的难度.

图3 降阶模型误差随着聚类数变化Fig.3 Errors of reduced⁃ordermodel at various numbers of clusters

4 结论

结合全隐式法求解考虑启动压力梯度的二维油水两相渗流模型的基本过程和基于本征正交分解构造降阶模型的基本思想，研究了利用聚类分析预处理快照对所构造出的降阶模型的影响.研究表明：

1）基于本征正交分解法基本思想所建立的降阶模型对原模型有很好的近似，在一定程度上减小了计算量.

2）基于CVT⁃Lloyd聚类分析处理快照可以有效减小数据冗余从而提高模型精度，但过少的聚类数使得信息缺失，影响降阶模型对原模型的近似.

［1］ Huang Yanzhang.Flow mechanism in low permeability oil layers［M］.Beijing：Petroleum Industry Press，1998：93－97.

［2］ Song Fuquan，Liu Ciqun.Analysis of two⁃phase fluid flow in low permeability reservoir with the threshold pressure gradient ［J］.Journal of the University of Petroleum，China，1999，23（3）：47－56.

［3］ Li Daopin.Development of low permeability reservoir［M］.Beijing：Petroleum Industry Press，1999：25－30.

［4］ Han Dakuang，Chen Qinlei，Yan Cunzhang.Fundamentals of reservoir simulation［M］.Beijing：Petroleum Industry Press，2001.

［5］ Jiang Yaolin.Methods of reduced⁃ordermodeling［M］.Beijing：Science Press，2010.

［6］ Quarteroni A，Veneziani A.Analysis of a geometricalmultiscalemodel based on the coupling of PDE’s and ODE’s for blood flow simulations［J］.Multiscale Modeling and Simulation，2003，1（2）：173－195.

［7］ Vermeulen P T M，Heemink A W.Inverse modeling of groundwater flow using model reduction［J］.Water Resources Research，2005，41（6）：1－13.

［8］ Yang Y J，Shen K Y.Nonlinear heat⁃transfermacromodeling for MEMS thermal devices［J］.Journal of Micromechanics and Microengineering，2005，15（2）：408－418.

［9］ Burkardt J，Gunzburger M，Lee H.POD and CVT⁃based reduced⁃order modeling of Navier⁃Stokes flows［J］.Computer Methods in Applied Mechanics and Engineering，2006，196（1）：337－355.

［10］ Cardoso M A，Durlofsky L J.Use of reduced⁃ordermodeling procedures for production optimization［J］.SPE Journal，2010，15（2）：426－435.

［11］ Sirovich L.Turbulence and the dynamics of coherent structures part I⁃III［J］.Quarterly of Applied Mathematics，1987，45 （3）：561－590.

［12］ Zhang Bo.Research on numerical simulation of low permeability［D］.Qingdao：China University of Petroleum（East China），2011.

［13］ Wang Jianzhong，Yao Jun，Wang Wei.Numerical simulation of oil⁃water two⁃phase flow in low permeability reservoirsconsidering well⁃grid influence［J］.Journal of Southwest Petroleum University，2013，35（5）：81－86.

［14］ Kunisch K，Volkwein S.Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics［J］. SIAM Journal on Numerical analysis，2002，40（2）：492－515.

［15］ Li Chengtao，Xiao Yiqing，Ou Jinping.Study on themodel reduction for flexible structure based on clustering algorithm and DOFs concentration［J］.Chinese Journal of Computational Mechanics，2012，29（2）：236－241.

［16］ Hao Pan.A connected area weighted centroidal Voronoi tessellation model for image segmentation［D］.Hangzhou：Zhejiang University，2010.

［17］ Wu Xuehong，Tao Wenquan，Lv Yanli，Gong Yi.Reduced order model for fast computation of incompressible fluid flow problem［J］.Proceedings of the CSEE，2010，30（26）：69－73.

［18］ Wang Zhongxia.Reduced ordermodels based on PODmethod for wind turbine blade design［D］.Beijing：Chinese Academy of Sciences，2011

Reduced⁃order M odel Based on Cluster Analysis in Numerical Simulation of Low Permeability Reservoir

ZHANG Yanyu，LIWeiwei，CHEN Huijuan
（College of Petroleum Engineering，China University ofPetroleum，Qingdao 266580，China）

A mathematical model for numerical simulation of 2D water⁃oil phase reservoir concerning start⁃up pressure gradient is solved with full implicit finite difference scheme.A reduced⁃ordermodel is established by projecting originalmodel in low dimension space formed by basis functions，which is produced from collected snapshots；CVT⁃Lloyd algorithm is employed in cluster analysis to tackle with flaws from collecting snapshots of same time interval.Computational results show that reduced⁃order model established based on POD is able to approximate the originalmodel effectively.With cluster analysis，data reduction is eliminated，which could create uniform snapshots in space resulting in enhanced accuracy of reduce⁃order model.On the other hand，extremely limited clustering groupsmay lead to loss of information and thus impairing accuracy of reduced⁃ordermodel.

low permeability reservoir；cluster analysis；reduced⁃ordermodel；proper orthogonal decomposition；

TE34

A

2014－08－07；

2014－12－28

张艳玉（1963－），女，博士，教授，从事油藏工程的教学和研究工作，E⁃mail：yyzhang＠upc.edu.cn

Received date： 2014－08－07；Revised date： 2014－12－28