一种用于多孔介质气体流动数值模拟中传导率精确计算的积分方法

《Frontiers in Earth Science》：An integral method for accurate calculation of the transmissibility in the numerical simulation of gas flow in porous media

　　

准确评估传导率（transmissibility）是可靠进行多孔介质中气体流动数值模拟的核心，因为气体密度和粘度等性质具有压力依赖性。在商业油藏数值模拟器中，块间传导率采用上游平均（upstream averaging）计算，而在含井网格块中则在网格块压力处进行评估。该方法未能捕捉流动的物理本质，因为气体性质从一个网格块到另一个网格块逐渐变化，且在含井网格块内变化显著。本研究提出了一种新颖的基于积分的公式，其中压力依赖性质在压力范围内进行积分，为传导率项提供精确值，无需调用假设或数学平均。所提出的方法完全基于物理、稳健且易于在有限体积油藏数值模拟器中实现。由于其通用推导，该方法适用于天然气储层、涉及多相流的CCUS（Carbon Capture, Utilization and Storage）应用中的CO2注入与封存以及地下氢储存（Underground Hydrogen Storage, UHS）。该方法通过避免无依据的数学平均，有助于更精确地获取压力和产量值。将所提公式的结果与商业软件获得的结果进行比较表明，当网格块相对较大且绝对渗透率较低时，有显著改进。

研究背景与问题提出

在天然气储层、CO₂地质封存及地下氢储存（UHS）的数值模拟中，多孔介质气体流动的准确表征至关重要。在有限体积法（Finite-Volume）公式中，传导率（transmissibility）控制着相邻网格块之间及储层与井之间的质量传递，其准确性决定了数值解的精度。由于气体密度（ρ）和粘度（μ_g）随压力呈强非线性变化，传导率成为依赖压力的变量，必须谨慎处理以保持气体流动的物理本质。目前大多数商业油藏模拟器（如ECLIPSE或CMG）在计算块间传导率时采用上游平均（upstream averaging），在含井网格块中则直接在网格块压力点进行评估。然而，自然界中压力沿流线连续变化，而非在网格块间离散突变，且含井网格块内压降极为陡峭。在单一压力点评估气体流动性（mobility）歪曲了真实物理过程，尤其在低渗透地层或大网格块条件下误差显著。尽管早期研究建立了真实气体拟压力（real-gas pseudo-pressure）理论强调沿压力积分的必要性，后续模拟实践及MPFA等离散格式仍延续了单点流动性采样的做法，未在传导率项中显式整合气体粘度和密度对压力的显式依赖，导致在非常规储层近裂缝区域及CO₂临界态附近存在物理缺口。因此，开展一项完全基于物理、不依赖数学平均的传导率计算方法研究具有重要的理论与应用价值。该研究发表于《Frontiers in Earth Science》。

主要关键技术方法

研究人员采用质量平衡方法与达西定律（Darcy’s law）推导一维及三维离散方程，基于IMPES（Implicit Pressure Explicit Saturation）公式框架开展研究。核心方法是将原本在网格边界或含井网格块内采用的单点评估流动性项1/(B_gμ_g)（其中B_g为气体体积系数）重构为压力域内的积分形式，通过对两相邻网格块压力差区间或网格块至井筒压力区间的1/(B_gμ_g)进行定积分并除以压力差来获得平均传导率，以此替代上游加权平均。研究人员利用McCain Jr 1972的样本数据拟合μ_g、B_g及1/(B_gμ_g)的压力相关性方程并积分得到拟压力类函数mp，开发内部代码（in-house code）实施该积分方法，并与CMG、Eclipse商业软件在相同油藏属性（变渗透率0.1 mD与0.0001 mD、5000 ft长度、1000网格块等，见表2）下进行对比模拟。

研究结果

2 Mathematical development

研究人员从质量守恒出发推导了有限单元内的流动方程，指出离散化导致网格间压力突变但与自然连续变化不符。为解决此问题，研究人员在两网格中心间引入虚拟点，将达西定律适用的微分形式推广至有限压差，推导出块间界面处的1/(B_gμ_g)_i+1/2,j应等于压力从p_i+1,jⁿ⁺¹到p_i,jⁿ⁺¹对1/(B_gμ_g)的积分除以压力差（Equation 7）。同理，对于含井网格块流向井筒的流动（Equation 4），研究人员提出在网格块压力与井底流压间进行积分以获取精确的传导率项，因该区域内压降梯度最大且单点评估误差最显著。对于质量平衡方程右侧的B_g与c_g（气体压缩系数），研究人员建议在当前时间步计算并配合早期生产小时间步长以减小误差。

2.1 Inter-block transmissibility calculation

研究人员应用Equation 7计算各时间步网格块边界的1/(B_gμ_g)，指出其与上游压力点取值的差异即模拟误差来源。在网格加密条件下两者差异极小，但在大网格下积分法提供更真实的边界属性。

2.2 Transmissibility calculation in the well equation

这是研究的核心部分。研究人员图解说明了含井网格块内1/(B_gμ_g)随压力的曲线，指出积分法获得的传导率（曲线下面积平均）与在网格块压力单点计算值存在显著差异，此差异在生产早期、低渗储层及大尺寸立方体网格块中尤为突出。迭代更新压力或小时间步无法解决物理本质的缺失，唯有积分法能考虑网格块内渐变压降。

2.3 RHS of the mass balance equation (i.e., equation 3)

研究人员讨论了累积项中的压力依赖参数处理，确认在单时间步内网格压力变化不大时可安全使用当前时间步值，早期生产采用小时间步以避免大压降。

3 Results and discussion

研究人员基于表1气体属性数据与表2系统参数进行模拟对比。Case 1（k=0.1 mD）显示三种工具（CMG、Eclipse、积分法内部代码）仅在生产极早期有差异，随时间推移含井网格压力逼近井压，差异消失，一年累计误差可忽略，表明研究早期动态需积分法。Case 2（k=0.0001 mD）极低渗下含井网格压降极慢，CMG/Eclipse与积分法结果差异显著且持续长久（前100天及至27年仍不同），证实在低渗条件下必须采用基于物理的积分法计算1/(μ_gB_g)。研究还发现网格数量多时间积分数值接近上游法，径向小网格近井区积分法必要性降低，但在大立方网格中至关重要。

讨论与结论翻译

自然界中，压力从一个网格块到另一个网格块逐渐变化，在含井网格块中从网格块边缘到井筒逐渐变化。为解释此行为，提出一种简单的积分方法来计算传导率项中的1/(B_gμ_g)，提供了两网格块间或含井网格块内压力范围内的平均值。若系统中网格块数量多导致相邻网格块压力非常接近，则无需用此方法获取块间传导率。但在井方程中，尤其当含井网格块较大时，必须使用此方法。当该网格块较小或使用径向离散（中心径向网格块极小）时，1/(B_gμ_g)可在井压处安全计算。当渗透率较低时，更有必要采用此积分方法。采用更小时间步或迭代程序无助于获得正确结果，因为这些属于数值手段而非物理本质。所提方法完全基于物理、稳健且易于在有限体积油藏模拟器中实现。由于其通用推导，该方法适用于广泛的产气主导储层系统，包括天然气储层、涉及多相流的CCUS（Carbon Capture, Utilization and Storage）应用中的CO₂注入与封存以及地下氢储存（Underground Hydrogen Storage, UHS）。该方法为气体流动模拟提供了统一可靠的框架。