基于渗流-温度双场耦合的油藏型储气库数值模拟

何海燕,刘先山,耿少阳,孙军昌,孙彦春,贾倩

（1.中国石油冀东油田公司，河北 唐山 063200；2.成都理工大学能源学院，四川 成都 610059；3.东北石油大学环渤海能源研究院，秦皇岛 河北 066004）

国家能源局再次强调要持续大力推进天然气“产供储销”体系建设，储气库作为其中重要一环备受关注。“十四五”规划预计2025 年建成350×108m3调峰及储备能力，但目前中国储气库建设优质资源缺乏，因此，油藏被纳入建库选址范畴。

中国大陆经历多期次构造运动，其东部形成一系列复杂断块，中西部形成褶皱构造，导致建库油藏断块多、非均质性强、边底水及人工注水侵入地层后流体关系复杂[1-2]。因此，目前无论是已建库的京58 储气库，还是正在建库的冀东油田南堡1-29储气库，都是复杂断块油藏型储气库。复杂断块油藏改建储气库后，多周期高速注采过程均为油气水三相流动，存在注采周期短、气体流速高、压力波及范围小等特征。注气周期注入的冷气会扰动储层温度场，油气水的黏度、气油比等高压物性参数受温度影响十分严重。在油藏衰竭开发阶段，有较多学者针对注气、水、压裂液等对储层产生的温度场扰动开展了研究。王增林等[3]使用Fluent 软件模拟均质储层条件下使用不同管柱的油藏温度场变化。郑少婧等[4]通过实验探究了储气库交替注采工况下储层渗透率温度敏感性，基于实验结果建立了考虑渗透率温度敏感性的气井产能方程。郭肖等[5]应用热应力理论，推导了渗透率随温度变化的理论模型，并采用变围压、变内压应力敏感实验进行了验证。前人通过实验和理论证实了储层中温度场被扰动后，会对岩石渗透率、弹性模量、泊松比及流体黏度、体积系数等高压物性产生较大影响。此外，不同于气藏衰竭开发阶段中流体的低速流动只需克服黏滞阻力产生的压降，储气库运行阶段采气速度是气藏衰竭开发阶段的20~30 倍，井筒附近天然气高速流动产生的惯性力附加压降已不容忽视。EL-ZEHAIRY 等[6]基于XCT 数据研究了多孔介质微观非均质性对惯性流的影响，研发了孔隙网络模型（PNM）来模拟多孔介质中非达西流动，与均匀多孔介质相比，非均质多孔介质中由于连接较差，孔喉中存在更多停滞区域，减少了流体流动的有效面积，样品内从停滞区到吼道连接处速度分布由低到高，更容易观察到流体惯性效应。WANG 等[7]学者指出毛细管力和边界层效应是油藏中非达西流动产生的主要原因，提出了一种油藏低速注水过程中非达西流体动力学表征方法，通过该方法，给定毛细管力和边界层厚度即可计算出单个毛细管内流体的流速。NIE 等[8]学者也研究了流体在未固结介质（中等强度支撑剂填充介质）与固结介质（露头、岩心等）中流动时的非达西紊流因子测定方法及表征公式。

传统油藏工程方法和常规数值模拟将储层视为恒温[9-11]，岩心实验虽能模拟流体高速非达西效应及温度变化对储层岩石和流体高压物性的影响，但受限于尺度，仅能代表储层中一个点。然而，在复杂断块油藏中，同一断块内储层存在非均质性，不同断块中储层具有不同温压及流体系统，在改建储气库后更是存在短周期、高气体流速以及井控范围小等特征。因此，传统油藏工程方法、岩心实验和常规数值模拟难以同时刻画复杂断块油藏型储气库交替注采工况、温度场扰动及高速非达西效应附加压降，导致储气库多周期运行过程调峰、单井注采能力等生产动态指标预测精度低，最终造成新钻井数量、投资预测等误差大。因此，结合高速非达西实验及流体黏温实验，重点研究复杂断块油藏型储气库周期注采过程储层温度交替变化及高速非达西附加压力损失对储层流体渗流及生产动态的影响。研究成果可以为复杂断块油藏建库方案设计提供理论指导，为储气库安全高效运行奠定基础。

1 工区概况

储气库由带气顶弱边底水的饱和油藏改建而成，储层平均净毛比为0.35，平均孔隙度为0.08，平均渗透率为19×10-3µm2。储层自上而下有3 个不同的流体系统，其中气层中束缚水饱和度为0.28，另外2套油层中束缚水饱和度分别为0.36 和0.4，油层之下为含水层。工区被断层分为7个断块，各断块气油水界面不统一，流体分布（图1）较为复杂。

图1 不同断块中油气水3相初始分布Fig.1 Initial distribution of oil,gas and water in different fault blocks

2 储气库渗流-温度耦合数学模型

为描述复杂断块油藏型储气库储层中低温天然气注入后流体高压物性变化，以及温度场扰动对渗流场的影响，基于多孔介质中渗流-温度耦合作用机理，建立了渗流-温度耦合数学模型。模型基本假设条件如下：①枯竭油藏型储气库中构造由7个断块组成，采用笛卡尔网格剖分；②储层中流体由油气水三相组成，密度、黏度是温度和压力的函数；③天然气在储层中高速流动时，除了受孔喉壁面及流体间黏滞阻力影响外，还受高速惯性力产生的附加压降影响。

描述了储气库多孔介质中油气水三相流动的控制方程，主要包含质量守恒方程和能量守恒方程。每一组方程都有3个相的子集，加在一起就组成了完整控制方程。质量守恒方程和能量守恒方程适用于所有相。因此，控制方程可以简化为单相形式。

多孔介质中流体质量守恒方程：

式中：φ为储层孔隙度；ρm为流体密度，单位kg/m3，m表示油气水三相；t为时间，单位s；v→m为流体速度，单位m/s；qm为流体产量，单位m3/s。

考虑高速非达西效应的流动方程采用经典的Forchheimer方程：

式中：p为储层压力，单位MPa；μm为流体黏度，单位mPa·s；k为渗透率，单位10-3µm2；βm为地层流体体积系数，单位m3/m3。

地层流体密度是压力与温度的函数，表达式为：

式中：pi为参考压力，单位MPa；Ti为参考温度，单位K；ρi为参考密度，单位kg/m3；pm为储层中油气水的压力，单位MPa；Tm为储层中油气水的温度，单位K；ρm为储层压力pm和储层温度Tm对应的油气水的密度，单位kg/m3；ρi为参考压力pi和参考温度Ti对应的油气水的密度，单位kg/m3；CL为流体弹性压缩系数；αL为流体热膨胀系数。

式中：μm为储层压力pm和储层温度Tm对应的油、气、水的黏度，单位mPa·s；μi为参考压力pi和参考温度Ti对应的油、气、水的黏度，单位mPa·s；γm与ηm均为储层中流体的黏度方程系数。

流体饱和度方程：

式中：nm为流量数量；Sm为储层流体饱和度。

将毛细管压力定义为单位面积的力，可以表示为：

式中：pc为毛细管压力，单位MPa；r为毛细管半径，单位m；σ为表面张力，单位N/m；θ为接触角，（°）。

毛细管压力与重力作用过程，向上与向下的力是平衡的，可以得到：

式中：g为重力加速度，单位m/s2；h为润湿相在毛细管中上升的高度，单位m。

油气间毛管压力可以表示为：

式中：pcgo为油气间的毛细管压力，单位MPa；pg为气相压力，单位MPa；po为油相压力，单位MPa。

饱和度和毛细管压力之间关系：

式中：pcow为油水间的毛细管压力，单位MPa；pw为油相压力，单位MPa；Sw为含水饱和度。

通常使用LEVERETT[12]提出的J函数：

为了模拟多相流动，需要确定储层内初始饱和度场分布。如果我们知道油水界面的位置，就可以通过结合式（7）和式（10）来确定地层中饱和度随深度的分布：

初始条件：

式中：z为纵向距离，单位m；z0为参考垂向距离，单位m；p0为参考垂向距离z0对应的储层压力，单位MPa。

为了表征流体与流体、流体与岩石间的热量交换关系，引入能量守恒方程：

式中：hm为流体的焓，单位J；λm为流度系数，单位10-3µm2/（mPa·s），其中，λm=k·kr,m/μm，kr,m为储层流体相对渗透率，单位10-3µm2；（ρU）eff为单位质量的有效内能，单位J；Λeff为有效导热系数，单位W/（m·K）。

式中：Um为流体的内能，单位J；Ur为岩石的内能，单位J；ρr为岩石密度，单位kg/m3。

式中：Λm为流体导热系数，单位W/（m·K）；Λr为岩石导热系数，单位W/（m·K）。

基于MATLAB 软件中油藏数值模拟工具箱（Matlab Reservoir Simulation Toolbox，简称MRST），非线性方程离散采用有限体积法（Finite Volume Method，简称FVM），在空间上采用两点通量近似（Two-Point-Flux-Approximation，简称TPFA）有限体积格式，在时间上采用后向（隐式）欧拉格式对方程进行耦合离散求解[13]，井模型采用Peaceman 模型[14]进行离散求解。

3 流体高速非达西及高压物性实验

3.1 高速非达西实验

结合Forchheimer 方程（式16）可拟合紊流因子β及渗透率k[15-17]，为获取紊流因子设计了此次驱替实验。实验主要采用智能驱替模拟系统，选取N 储气库代表性岩心，分别设计采用0.50～0.95 MPa 共10组驱替压力，测试不同驱替压力下驱完2 mL 气体的气体流速（测试体积/测试时间）。实验气体为氮气，测试氮气体积为2 mL，实验岩心长度为5.12 cm，直径为2.53 cm，驱替环压为10 MPa，实验结果见表1。

表1 紊流因子测试实验数据Table 1 Experimental data for turbulence factor

式中：p为压力，单位MPa；X为流体流动方向；μ为气体黏度，单位mPa·s；v为流体流速，单位cm/s；k为渗透率，单位10-3µm2；β为紊流因子，单位108/m；ρ为流体密度，单位g/cm3。

在此次驱替实验条件下，Forchheimer 方程可以表达为式（17）。以MA（p12-p22）/(2zRTμlρpQp)为纵坐标，ρpQp/（μA）为横坐标，可以拟合出一条直线，该直线截距为1/k，斜率为β，即紊流因子。

式中：M为气体分子质量，单位g/mol；A为实验样品横截面积，单位cm2；p1为实验样品入口压力，单位MPa；p2为实验样品出口压力，单位MPa；z为气体压缩因子；R为通用气体常数，R=8.314 472 m3·Pa/（K·mol）；T为温度，单位K；l为实验样品长度，单位cm；ρp为泵中流体密度，单位g/cm3；Qp为泵中流体流量，单位cm3/h。

令式（17）中MA（p12-p22）/(2zRTμlρpQp)为y，ρpQp/（μA）为x，式（17）可以表达为式（18）。

通过紊流因子岩心实验，以MA（p12-p22）/（2zRTμlρpQp）为纵坐标，ρpQp/μA为横坐标，拟合绘制直线（图2）的斜率（即紊流因子β）为1.63×108/m。

图2 紊流因子岩心实验拟合Fig.2 Turbulence factor regression curve base on experiment data

3.2 流体黏温关系

实验室一般使用黏度计或旋转式流变仪来测量原油的黏度[18]。从N 油藏储气库储层中获取原油样品，采用自动密度黏度测定仪，测定剪切速率为60 s-1时，不同温度（10.13～90.04 ℃）条件下原油黏度，实验结果见图3。气水黏度计算通常采用4 种方法：Lohrentz-Brey-Clark（LBC 方法）、PFCT 方法、SUPERTRAPP 方法、Vesovic-Wakeham（VW 方法）。研究表明SUPERTRAPP 方法误差较小[19-21]。因此，储层压力（22.6 MPa）条件下的油气水黏度与温度关系基于SuperTrapp软件数据包，结合SuperTrapp模型计算获得（图3）。

图3 油气水黏温关系（压力：22.6MPa）Fig.3 Viscosity-temperature relationship curves of oil,gas and water

4 复杂断块油藏储气库数值模拟

结合前期油气藏数值模拟、油气藏改建储气库数值模拟方法[22]，建立了考虑储气库多轮注采过程冷气注入扰动储层温度场及高速非达西效应等机理模拟的油藏型储气库数值模拟方法和技术流程（图4）。依据流体高压物性随温度变化实验及高速紊流实验，获取流体高压物性参数随温度的函数关系及高速非达西紊流因子，分别用于储气库注气周期冷气注入过程温度场扰动数值模拟以及对流动方程（式2）中惯性力产生的附加压降进行校正。数值模拟采用笛卡尔网格剖分，网格数约为94万，X、Y、Z方向平均步长分别为40 m×40 m×4 m。

图4 油藏型储气库数值模拟技术流程Fig.4 Workflow of numerical simulation for UGS rebuilt from oil reservoir

4.1 生产动态历史拟合

建模阶段属性模型是通过地质统计学数据分析、插值生成的，存在较强的不确定性。因此，需结合衰竭开发阶段生产动态监测资料来反演井间储层物性参数。首先对区块油、气等进行拟合（图5a、图5b），从而保证区块物质平衡，然后检查单井瞬时产量拟合情况，在此基础上依次开展单井静压、流压拟合（图5c、图5d）。

4.2 冷气注入对储气库生产动态的影响

基于历史拟合的模型，选取主力断块5口储气库井，单井平均以11.2×104m3/d 注气200 d，20×104m3/d 采气120 d，注采平衡期均为15 d，设计考虑温度场及高速非达西效应影响的方案。其中模型的注入气温度为25 ℃，储层中深温度为87.8 ℃，储层初始温度梯度为3 ℃/hm，流体黏度随温度变化规律依据黏温实验数据（图3）赋值，紊流因子依据非达西实验结果赋值为1.63×108/m。研究将流体黏温、高速紊流实验结果敷设于模型里，把不考虑温度场及非达西效应的方案设为基础方案。

模拟结果表明，相较于基础方案，注冷气扰动温度场方案由于井控范围内储层温度下降，导致累产油量（图6b）下降，累产水量增加（图6c），而累产油下降幅度小于累产水增加幅度，使地层采出液量增多，地层压力下降（图6d）。考虑非达西效应影响时，一方面，除了黏滞阻力产生压降外，高速惯性力还多产生了一部分压降，因此，相同配产配注条件下，天然气注入后部分采不出，随着储气库多周期运行，注气末天然气储量及压力逐渐增加；另一方面，在定产气量生产条件下考虑非达西效应，需增大生产压差才能产出相同的天然气量，因此，多周期运行后油、水累产量增加。总的来说，不考虑冷气注入温度场扰动和高速非达西效应，将造成第三采气期末累产油、气量分别偏小3.29%、10.52%，累产水量偏大21.07%。

此外，模拟了4 组不同注入气温度（10、30、50、70℃）的储气库多周期运行方案。结果表明，温度场扰动对气体渗流的影响微弱，远小于高速非达西效应对气体渗流的影响。但储层温度降低后对原油渗流能力降低幅度影响较大（图7），不考虑温度场扰动的基础方案3周期内累产油26 201 m3，考虑注入气温度为10 ℃时，原油采出量降低4 399 m3，降幅为16.79%。

图7 注入气温度对原油累产的影响Fig.7 Effect of injection gas temperature on oil production cumulative

选取构造高部位QK4 井，设计考虑高速非达西与不考虑高速非达西效应各6项注采方案（注200 d，采120 d），方案中日采气速度分别为（10、20、40、60、80、100）×104m3，对应的日注气速度分别为（5.6、11.2、22.4、33.6、44.9、56.1）×104m3。结果表明，在当前物性及配产配注条件下，注气末井控温度范围随注气速度呈对数上升，低速阶段上升较快，高速阶段由于受储层物性及有限时率强注强采限制，上升速度变缓；考虑高速非达西效应后，低速注采阶段流体惯性力产生的附加压降较小，对流体渗流影响较小，注气末井控温度范围几乎重合；随着注采速度上升，井控温度范围逐渐增大，但注气速度上升到44.9×104m3/d 后，考虑高速非达西效应的方案注气末井控温度范围几乎恒定，未考虑高速非达西效应的方案仍在增加，这表明注采速度较高时，受非达西效应附加压降影响，该物性下单井注采能力已达到极限，而不考虑非达西效应影响的单井注采能力还有较大提升空间（图8）。因此，基于该极限，可以确定不同储层物性下受温度场扰动及高速惯性力影响的单井合理注采气能力。

图8 注气速度对注气末期井控温度范围的影响Fig.8 Effect of gas injection rate on well control temperature range at the end of gas injection cycle

5 结论

1）建立了渗流-温度双场耦合数学模型，并基于有限体积法（FVM），在空间上采用两点通量近似方案（TPFA），在时间上采用后向（隐式）欧拉格式对模型进行耦合离散求解。该模型考虑了储气库注气周期冷气注入扰动温度场以及高速注采过程中高速非达西效应对流体渗流的影响，更符合储气库的特殊工况。

2）基于实验回归了紊流因子与流体黏温关系，明确了油气黏度大幅上升时水黏度几乎不变，这一特性使得油气相对渗流能力下降时水相对渗流能力反而上升。基于实验数据及建立的模型，开展了储气库多周期运行数值模拟。模拟实例表明，温度场扰动、高速非达西效应分别是累产油、气量误差的主控因素。储层温度场扰动减少的产油量比高速非达西效应增加的产油量多。不考虑冷气注入温度场扰动和高速非达西效应，将造成第三采气期末的累产油、气量分别偏小3.29 %、10.52 %，累产水量偏大21.07%。

3）生产动态敏感性分析结果表明，冷气注入造成井控范围内储层温度下降，该范围随注气速度增加呈对数上升，低速阶段上升较快，高速阶段由于受储层物性及有限时率强注强采限制，上升速度变缓；高速惯性力附加压力损失使相同配产配注条件下天然气注入后部分采不出，随着储气库多周期运行，注气末天然气储量及压力逐渐增加，定天然气量生产所需生产压差增大，油水累产量随之增加。