CN116861753B - 一种基于模拟有限差分法的油水两相流线模拟新方法 - Google Patents
一种基于模拟有限差分法的油水两相流线模拟新方法 Download PDFInfo
- Publication number
- CN116861753B CN116861753B CN202310946821.3A CN202310946821A CN116861753B CN 116861753 B CN116861753 B CN 116861753B CN 202310946821 A CN202310946821 A CN 202310946821A CN 116861753 B CN116861753 B CN 116861753B
- Authority
- CN
- China
- Prior art keywords
- phase
- water
- streamline
- calculation
- simulation method
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 124
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 title claims abstract description 104
- 238000004088 simulation Methods 0.000 title claims abstract description 51
- 238000004364 calculation method Methods 0.000 claims abstract description 56
- 230000035699 permeability Effects 0.000 claims abstract description 21
- 239000000126 substance Substances 0.000 claims abstract description 8
- 239000000243 solution Substances 0.000 claims description 9
- 230000006835 compression Effects 0.000 claims description 6
- 238000007906 compression Methods 0.000 claims description 6
- 239000012088 reference solution Substances 0.000 claims description 5
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 230000035515 penetration Effects 0.000 claims description 3
- 239000011435 rock Substances 0.000 claims description 3
- 239000000109 continuous material Substances 0.000 claims description 2
- 230000001788 irregular Effects 0.000 abstract description 5
- 238000004519 manufacturing process Methods 0.000 description 18
- 238000002347 injection Methods 0.000 description 8
- 239000007924 injection Substances 0.000 description 8
- 238000010586 diagram Methods 0.000 description 6
- 239000011159 matrix material Substances 0.000 description 5
- 239000006185 dispersion Substances 0.000 description 3
- 238000005325 percolation Methods 0.000 description 3
- 238000003745 diagnosis Methods 0.000 description 2
- 230000009977 dual effect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000011144 upstream manufacturing Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- CLOMYZFHNHFSIQ-UHFFFAOYSA-N clonixin Chemical compound CC1=C(Cl)C=CC=C1NC1=NC=CC=C1C(O)=O CLOMYZFHNHFSIQ-UHFFFAOYSA-N 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000011438 discrete method Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000012886 linear function Methods 0.000 description 1
- 239000007788 liquid Substances 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 230000000630 rising effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
- Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- Fluid Mechanics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于模拟有限差分法的油水两相流线模拟新方法，包括如下步骤：步骤1、基于MFD对各相连续物质的压力值进行计算；步骤2、基于流线对含水饱和度进行计算；步骤3、计算在各向异性、复杂几何情况时相较于传统FV和MFD的计算精度和计算效率。本发明采用上述的一种基于模拟有限差分法的油水两相流线模拟新方法，基于MFD的本发明能够在全渗透率张量和不规则网格情况下取得高精度的压力计算而更准确的估计计算域的渗流速度分布；相较于经典的FV和MFD，在相近网格数量情况下，本发明能够显著降低饱和度计算的耗散误差；随着网格数量的增加，本发明在计算效率上的优势也会更加显著，可应用于复杂油藏计算域并快速获取流线分布。
Description
技术领域
本发明涉及流线模拟技术领域，尤其是涉及一种基于模拟有限差分法的油水两相流线模拟新方法。
背景技术
油藏数值模拟是对地下油藏多相多组分渗流研究的重要手段，其将油藏渗流方程组进行离散求解。随着几十年的发展，基于满足局部物质守恒、紧致、单调的两点流量估计格式的有限体积离散方法最常见于各类商业或学术油藏数值模拟器中。然而，对于具有全张量渗透率并带有复杂几何的油藏模型，传统基于有限体积(差分)法的两点流量格式将难以给出高精度的油藏压力及渗流速度分布，而MFD被证实能在网格分布不规则、网格质量差的情况下取得比多点流量格式及混合有限元法更高的计算精度，因此开发基于MFD的流线模拟方法是十分自然且必要的，此外，在复杂网格情况下，传统的流线追踪方法会明显增加流线模拟工作流的复杂度。
发明内容
本发明的目的是提供一种基于模拟有限差分法的油水两相流线模拟新方法，基于MFD的本发明能够在全渗透率张量和不规则网格情况下取得高精度的压力计算而更准确的估计计算域的渗流速度分布；相较于经典的FV和MFD，在相近网格数量情况下，本发明能够显著降低饱和度计算的耗散误差；随着网格数量的增加，本发明在计算效率上的优势也会更加显著，可应用于复杂油藏计算域并快速获取流线分布。
为实现上述目的，本发明提供了一种基于模拟有限差分法的油水两相流线模拟新方法，包括如下步骤：
步骤1、基于MFD对各相连续物质的压力值进行计算；
步骤2、基于流线对含水饱和度进行计算；
步骤3、计算在各向异性、复杂几何情况时相较于传统FV和MFD的计算精度和计算效率。
优选的，所述步骤1中各相连续物质守恒方程如下：
式中a＝o(油相)或w(水相)，v是渗流速度，Sa是相a的饱和度，qa是相a的源汇项，t是时间，φ＝φ(p)是孔隙度，α是用于单位转换的常数；
渗透速度满足达西公式：
v＝-λaK▽pa，a＝o或w.....................................(2)
式中K是渗透率张量，λa＝kra/μa是相a的流度，其中kra＝kra(Sa)是相a的相对渗透率，μa是相a的粘度，pa是相a的压力，单位是MPa；
忽略毛管力，此时p＝pw＝po，将各相物质守恒方程叠加起来得到：
式中，λ＝krw/μw+kro/μo是总流度，Ct＝Cr+φ(SwCw+SoCo)是综合压缩系数，Cr、Cw和Co分别是岩石、水相和油相的压缩系数。
优选的，所述步骤2包括以下步骤：
步骤2.1、流线追踪方法；
步骤2.2、基于DG计算流线上饱和度。
优选的，所述步骤2.2中计算流线上含水饱和度分布的具体步骤是：
步骤2.2.1、计算
步骤2.2.2、根据计算/>
步骤2.2.3、利用限通器对解进行重构，计算出n+1时间步的含水饱和度分布：
式中L(·)是限通器算子。
优选的，所述步骤3中分别建立二维模型和三维模型，并通过比较三维模型的结果与二维模型的参考解验证本方法在三维四面体网格情况下的计算表现。
因此，本发明采用上述一种基于模拟有限差分法的油水两相流线模拟新方法，具有以下效果：
(1)本发明提出了基于MFD的流线模拟新方法，该方法利用MFD计算带有全渗透率张量的压力方程，基于简易的流线追踪方法，采用间断伽辽金法计算流线上的含水饱和度分布；
(2)相较于各类经典的流线模拟方法，基于MFD的本发明能够在全渗透率张量和不规则网格情况下取得高精度的压力计算而更准确得估计计算域的渗流速度分布。该方法构建了简单实用的适用于二维三角网格或三维四面体网格的流线追踪方法，能应用于复杂油藏计算域并快速获取流线分布；
(3)FV和MFD即使显著增加网格数量也难以达到与本发明相近的饱和度计算精度。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
图1是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例矩形各向异性油藏的网格及流量连续性条件示意图；
图2是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例某三角网格内的流线追踪示意图；
图3是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例矩形各向异性油藏中油藏计算域及网格离散示意图；
图4是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例矩形各向异性油藏基于本发明计算出的飞行时间分布；
图5是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例矩形各向异性油藏情况1中不同方法计算得到的100天和210天时含水饱和度分布对比图；
图6是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例矩形各向异性油藏中情况1不同方法计算得到的含水率曲线对比图；
图7是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例矩形各向异性油藏中情况2不同方法计算得到的100天和210天时含水饱和度分布对比图；
图8是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例矩形各向异性油藏中情况2不同方法计算得到的含水率曲线对比图；
图9是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例复杂几何油藏的油藏计算域及网格离散示意图；
图10是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例复杂几何油藏基于本发明计算出的飞行时间分布；
图11是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例复杂几何油藏中不同方法计算得到的100天和210天时含水饱和度分布对比图；
图12是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例复杂几何油藏中不同方法计算得到的含水率曲线对比图；
图13是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例基于本发明的复杂几何油藏模型流动诊断结果图；
图14是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例三维四面体应用油藏计算域及网格离散示意图；
图15是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例三维四面体应用中基于本发明计算出的飞行时间分布；
图16是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例三维四面体应用中流线上含水饱和度分布的计算结果对比图；
图17是是本发明一种基于模拟有限差分法的油水两相流线模拟新方法实施例三维四面体应用中生产井含水率曲线对比图。
具体实施方式
以下通过附图和实施例对本发明的技术方案作进一步说明。
除非另外定义，本发明使用的技术术语或者科学术语应当为本发明所属领域内具有一般技能的人士所理解的通常意义。
实施例一
如图所示，本发明提供了一种基于模拟有限差分法的油水两相流线模拟新方法，包括如下步骤：
步骤1、基于MFD对各相连续物质的压力值进行计算；
各相连续物质守恒方程如下：
式中a＝o(油相)或w(水相)，v是渗流速度，Sa是相a的饱和度，qa是相a的源汇项，t是时间，φ＝φ(p)是孔隙度，α是用于单位转换的常数；
渗透速度满足达西公式：
v＝-λaK▽pa，a＝o或w.....................................(2)
式中K是渗透率张量，λa＝kra/μa是相a的流度，其中kra＝kra(Sa)是相a的相对渗透率，μa是相a的粘度，pa是相a的压力，单位是MPa；
忽略毛管力，此时p＝pw＝po，将各相物质守恒方程叠加起来得到：
式中，λ＝krw/μw+kro/μo是总流度，Ct＝Cr+φ(SwCw+SoCo)是综合压缩系数，Cr、Cw和Co分别是岩石、水相和油相的压缩系数。
利用生成方法成熟且能适应于复杂几何的三角网格或四面体网格来离散计算域。为了方便叙述，以二维情况下的三角网格为例进行阐述，图1展示了计算域中的第i个网格，记该网格所在区域为Ωi、网格边为Aiβ(β＝1，2，3)、各边中心的位矢为xik(k＝1，2，3)、各边的单位外法向量为niβ。
对式(3)的Ωi积分，得到：
式中，a＝o或w。
对式(4)的左侧第一项利用散度定理，可以得到：
式中，Fiβ是该网格各边上的流量。
MFD中Fiβ的近似表达式为：
式中，λiβ是采用单点上游权格式的总流度，矩阵Ti是对称正定的第i个网格的传导率矩阵。表示矩阵Ti的β行η列。
MFD推导出矩阵Ti的计算表达式为：
Ti＝Ti1+Ti2，
式中，
Qi＝orth(AiXi)，Ai是由该网格各边边长构成的对角阵，Ki是该网格的渗透率张量。
根据式(6)和式(7)，得到左侧第一项近似为：
式中，
对时间偏导采用前向欧拉格式，则右侧近似为：
综合式(6)和式(7)，离散格式是：
此外，如图1所示，第i个网格和第j个网格间存在公共边，即η，s.t.Aiβ＝Ajη＝Am，则该边要求满足流量连续性条件，即要求分别从网格i、网格j得到的该边上的外法向流量代数和为0，具体表达式为
假设整个油藏计算域有nc个网格，有ne条边，其中有nie条内部边，压力自由度包括各网格中心压力值及各边中心压力值，共nc+ne个，全局线性方程组包括nc个方程、nie个方程以及由边界条件给出的ne-nie个方程，形成封闭的线性系统计算得到各压力值。
步骤2、基于流线对含水饱和度进行计算；
步骤2.1、流线追踪方法；
将给出一个简单的适用于三角网格或四面体网格的流线追踪方法，为了便于叙述，以三角网格为例进行阐述。
如图2所示，记该三角形网格的各顶点为1、2、3，各顶点坐标为x1＝(x1,y1)、x2＝(x2,y2)、x3＝(x3,y3)，各边为A、B、C，各边单位外法向量为nA、nB、nC，各边外法向渗流速度为VA、VB、VC，该网格流线入口坐标为xinlet＝(xinlet,yinlet)。
在步骤1中计算出压力分布后，利用式可以计算出各边上的外法向流量及渗流速度VA、VB、VC，假设该三角网格内的渗流速度处处相同，记作V＝(Vx,Vy)T，则有：
V·nA＝VA，V·nB＝VB，V·nC＝VC...................................(13)
利用最小二乘法计算得到V＝(Vx,Vy)T的近似值，具体表达式为：
V＝(GTG)-1(GTVe).................................................(14)
式中，G＝(nA,nB,nC)T，Ve＝(VA,VB,VC)T。
此时，由于该三角网格内的渗流速度处处相同，该网格内的流线轨迹将会成为一条直线，这使得该网格内的流线追踪及飞行时间计算变得简单。可以看到式的最小二乘解可能会使方程存在部分残差，但由于方程的方程个数3仅略大于未知数个数2，在三维四面体网格情况下，类似得，方程个数4(4个面)也仅略大于未知数个数3(渗流速度有x、y、z三个分量)，此时方程的残差将不会十分显著。但当网格变为四边形网格及六面体网格时，由于此时方程个数会较大于渗流速度分量数，使得网格内渗流速度处处相同的假设不再能取得很精确的流线追踪结果，此时，则可以使用假设渗流速度在网格内沿每个方向线性分布的经典Pollock追踪方法。
当得到V后，流线在该网格内的参数方程是：
假设该流线从A边离开该网格，边A的参数方程是：
式中，Δt是流线方程中的参数，是该流线在进入该网格Δt后的坐标，xinlet是网格流线入口坐标。
式中，Δr是A边方程中的参数，是A边上Δr对应的点坐标。
流线离开该网格处需满足：
WP＝U，Δtexit＞0................................................(17)
式中，W＝(V,x1-x2)，V是式(14)中的速度，x1和x2分别是三角形顶点1和顶点2的坐标，U＝x1-xinlet，xinlet是网格流线入口坐标。
则可以分别计算出流线在该网格内的飞行时间Δτexit(即流线从进入到离开该网格所需的真实时间，即按照真实速度而非渗流速度计算得到的时间)和流线离开该网格处的坐标xexit为：
Δτexit＝φΔtexit，xexit＝x1+Δrexit(x2-x1)............................(18)
式中，φ是该网格的孔隙度。
如上所分析的，本步骤流线追踪方法可以自然得推广到四面体网格，只是将式换成四面体网格面的参数方程即可。
步骤2.2、基于DG计算流线上饱和度。
当忽略毛管力、重力及油水相压缩性，流线上含水饱和度方程是：
式中，fw是含水率，v是沿着流线的油水相的总渗流速度，l是流线上的距离。
流线模拟通过引入飞行时间(TOF)将式改写为各条流线通用的方程：
一阶迎风差分格式一般被用来数值求解方程，在计算时，局部时间步长Δts与飞行时间步长Δτ的选取与计算压力的背景网格尺寸无关。当渗流速度分布稳定而无需更新流线后，对方程的计算将适用于所有流线，且各条流线的追踪和其上饱和度的计算均是相互独立、可高度并行的，因此，流线模拟的效率理论上远高于传统有限体积(差分)法，尤其是在网格尺寸较小而CFL条件对时间步限制较强的时候。
本发明为了能在相同飞行时间步长Δτ时取得比一阶迎风差分格式更高精度的流线上含水饱和度分布，本发明在流线上采用间断伽辽金法(DG)对方程进行计算，采用标准的有限元形函数，在每个Δτ表示的一维单元K上，假设单元左侧为上游，有
式中，上标r和l分别表示某点处右侧、左侧的取值，在DG中，同一个点处有左侧及右侧两个值，且该两个值可能不相等。下标K,1表示单元K左侧端点，下标K,1表示单元K右侧端点。
将式简记作： 和/>分别是单元K左侧及右侧的含水饱和度值，l(·,·)是对应的算子。采用简化的龙格库塔格式，并使用MUSCL格式的限通器对解进行重构，则利用DG计算流线上含水饱和度分布的具体步骤是：
步骤2.2.1、计算
步骤2.2.2、根据计算/>
步骤2.2.3、利用限通器对解进行重构，计算出n+1时间步的含水饱和度分布：
式中L(·)是限通器算子。
在计算得到各流线上含水饱和度分布后，则利用式计算出各网格平均含水饱和度：
式中，下标i表示第i个网格，下标j表示经过网格i的第j条流线，ns表示经过网格i的流线总数，Δτij和分别表示第j条流线在网格i内的飞行时间和平均含水饱和度。
步骤3、计算在各向异性、复杂几何情况时相较于传统FV和MFD的计算精度和计算效率。
分别建立二维模型和三维模型，并通过比较三维模型的结果与二维模型的参考解验证本方法在三维四面体网格情况下的计算表现。
矩形各向异性油藏：
左下角一口注水井，右上角一口生产井，储层各向异性渗透率张量见式，包含两种情况：其一是FV-TPFA能够处理的渗透率张量，其二是全渗透率张量。相渗数据及其他物性参数分别见表1和表2，可以看到，本例采用了特殊的线性相渗数据，结合表2中的粘度数据，含水率则是含水饱和度的线性函数。
从而，式(20)成为一个线性双曲守恒律方程，该黎曼问题的解表现出间断的水驱前缘，且在间断水驱前缘的前后应该分别是0.15和0.85(当含水率是含水饱和度的非线性函数时，会有一个油水两相区，但在线性情况下，则不会存在，这可以在Buckley-Leverett解析解中发现)。
因此，水驱前缘条带的宽度体现了数值耗散误差的大小，即条带宽度越小，数值耗散误差越小。该认识为对比不同方法对含水饱和度分布的计算误差提供了便利。
表1矩形各向异性油藏相渗数据
表2矩形各向异性油藏物性参数取值
图3中的(b)展示了本发明所使用的三角网格(网格数2408)，图3中的(b)和(c)是MFD所使用的三角网格(网格数分别是2408和38860)，图3中的(d)和(e)分别展示了FV所使用的笛卡尔网格(网格数分别是50×50＝2500和200×200＝40000)。图4展示了在两种渗透率张量情况下基于本发明计算出的300条流线上的飞行时间分布。
对于第一种情况，图5对比了FV、MFD和本发明计算得到的100天和210天时的含水饱和度分布，图6对比了不同方法计算出的生产井含水率曲线。可以看到：将各种方法计算结果中水驱前缘条带的中线看作计算出的水驱前缘，不同方法得到的水驱前缘是十分相近的，说明了各种方法的计算结果基本是可靠的。
在网格数相近的情况下，本发明的耗散误差远小于FV和MFD的耗散误差，将FV和MFD使用的网格数均扩大16倍左右到38860和40000时，FV和MFD的结果逐渐逼近于本发明计算结果，但耗散误差仍然高于网格数仅为2408时本发明的耗散误差。
这一点也反映在含水率曲线的对比中，因为理论上在该相渗数据情况下生产井的含水率在见水后会直接上升到1，而本发明得到的含水率曲线最接近这个特征，且随着网格数的增加，MFD和FV的含水率曲线也逼近于本发明得到的结果，也表明FV和MFD预测出的见水时间与真实见水时间会存在明显偏差。此外，还能看到在网格数相近时，MFD的计算精度会高于FV，这是已发现的结果，也揭示了本发明在MFD基础上更加显著得提高了计算精度。
表3对比了模拟500天生产动态时不同方法的CPU耗时(MFD和FV非线性求解器中最大时间步长设置为1天)，可以看到：在网格数相近时(2408和2500)，本发明CPU耗时均小于FV和MFD的CPU耗时，且远小于网格数40000时FV的CPU耗时。因此，本发明相较于FV、MFD具有计算精度和计算效率上的双重优势。
表3矩形各向异性油藏情况1中不同方法的CPU耗时对比
对于第二种情况，传统的FV-TPFA难以处理，图7和图8分别对比了本发明与MFD计算出的含水饱和度分布，表4对比了模拟600天生产动态时不同方法的CPU耗时，与第一种情况下的分析类似，也能发现本发明相较于MFD具有计算精度和计算效率上的双重优势。
式中，K为一维单元，mD为毫达西。
表4矩形各向异性油藏情况2中不同方法的CPU耗时对比
复杂几何油藏：
储层厚度为8米，有两口注水井和三口生产井。本例仍然采用与矩形各向异性油藏相同的相渗数据，因此可以根据前缘条带中线的位置和前缘条带的宽度来判断计算结果的合理性并估计相应的耗散误差大小，其他物性参数见表5。
表5复杂几何油藏计算物性参数取值
图9中的(b)展示了本发明及MFD所使用的三角网格，图9中的(c)和(d)分别展示了FV所使用的不同数量的PEBI网格(分别是1375、7537)。图10展示了本发明计算得到的飞行时间分布。图11对比了不同方法计算得到的100天和210天时的含水饱和度分布，图12对比了生产井的含水率变化曲线。从中可以看到：不同方法计算得到的水驱前缘条带中线差不多是相同的，说明了计算结果基本是可靠的，而本发明得到的水驱前缘条带最窄，说明了本发明对含水饱和度分布计算的数值耗散误差最小，从而对生产井见水时间的预测最准确，而FV和MFD由于耗散误差较大的原因预测的生产井见水时间过早。对于生产井pro#2，由于inj#2到pro#2的距离明显大于inj#1到pro#2的距离，因此pro#2的含水率曲线理论上应该体现出两个水驱前缘，即在含水率曲线中应存在两个陡峭上升的片段，可以看到，只有本发明做到了这一点，进一步说明了本发明的高精度。表6对比了复杂几何油藏中不同方法的CPU耗时，同样可以看出本发明在计算效率上的明显优势。
表6复杂几何油藏不同方法的CPU耗时对比
在此基础上，图13展示了基于本发明对该油藏模型的流动诊断结果，分别利用被不同颜色标记的流线绘制了便于现场工程师直观分析的各生产井的控制区域分布、各注水井的控制区域分布、各注采对之间的驱替区域分布、各生产井的实时卸油体积、各注水井的实时波及体积以及各注采井间的平面劈分系数，展示了本发明在油田现场应用的潜力。
三维四面体应用：
如图14中的(a)所示，本例中长方体油藏尺寸为200米×100米×10米，在左下角(0,0)和右上角(200,100)处分别有一口射穿整个储层的注水井和生产井，分别定注水速度80方/天和采液速度80方/天，其他物性参数与复杂几何油藏相同。本例采用图14中的(a)中的四面体网格离散该三维油藏计算域，并利用图14中的(b)中厚度为10米的二维模型计算结果(矩形各向异性油藏和复杂几何油藏已验证本发明在二维三角网格中应用的高精度)作为参考解来检验本发明在三维四面体网格情况下的计算精度。图15展示了本发明计算得到的飞行时间分布。图16对比了100天、200天、300天和400天时三维模型和二维模型分别得到的流线上含水饱和度分布。图17对比了三维模型和二维模型分别得到的生产井含水率曲线。可以看到：三维模型的结果与二维模型参考解十分相近，验证了本发明在三维四面体网格情况下也能取得良好的计算表现。
因此，本发明采用上述一种基于模拟有限差分法的油水两相流线模拟新方法，具有以下效果：
(1)本发明提出了基于MFD的流线模拟新方法，该方法利用MFD计算带有全渗透率张量的压力方程，基于简易的流线追踪方法，采用间断伽辽金法计算流线上的含水饱和度分布；
(2)相较于各类经典的流线模拟方法，基于MFD的本发明能够在全渗透率张量和不规则网格情况下取得高精度的压力计算而更准确得估计计算域的渗流速度分布。该方法构建了简单实用的适用于二维三角网格或三维四面体网格的流线追踪方法，能应用于复杂油藏计算域并快速获取流线分布；
(3)FV和MFD即使显著增加网格数量也难以达到与本发明相近的饱和度计算精度。
最后应说明的是：以上实施例仅用以说明本发明的技术方案而非对其进行限制，尽管参照较佳实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对本发明的技术方案进行修改或者等同替换，而这些修改或者等同替换亦不能使修改后的技术方案脱离本发明技术方案的精神和范围。
Claims (3)
1.一种基于模拟有限差分法的油水两相流线模拟新方法，其特征在于：包括如下步骤：
步骤1、基于MFD对各相连续物质的压力值进行计算；
步骤2、基于流线对含水饱和度进行计算；
步骤2.1、流线追踪方法；
步骤2.2、基于DG计算流线上饱和度；
步骤2.2.1、计算
步骤2.2.2、根据计算/>
步骤2.2.3、利用限通器对解进行重构，计算出n+1时间步的含水饱和度分布：
式中L(·)是限通器算子；
步骤3、计算在各向异性、复杂几何情况时相较于传统FV和MFD的计算精度和计算效率。
2.根据权利要求1所述的一种基于模拟有限差分法的油水两相流线模拟新方法，其特征在于：所述步骤1中各相连续物质守恒方程如下：
式中a＝o或w，其中，o为油相，w为水相，v是渗透速度，Sa是相a的饱和度，qa是相a的源汇项，t是时间，φ＝φ(p)是孔隙度，α是用于单位转换的常数；
渗透速度满足达西公式：a＝o或w.....................(2)
式中K是渗透率张量，λa＝kra/μa是相a的流度，其中kra是相a的相对渗透率，μa是相a的粘度，pa是相a的压力，单位是MPa；
忽略毛管力，此时p＝pw＝po，将各相物质守恒方程叠加起来得到：
式中，λ＝krw/μw+kro/μo是总流度，Ct＝Cr+φ(SwCw+SoCo)是综合压缩系数，Cr、Cw和Co分别是岩石、水相和油相的压缩系数。
3.根据权利要求1所述的一种基于模拟有限差分法的油水两相流线模拟新方法，其特征在于：所述步骤3中分别建立二维模型和三维模型，并通过比较三维模型的结果与二维模型的参考解验证本方法在三维四面体网格情况下的计算表现。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310946821.3A CN116861753B (zh) | 2023-07-28 | 2023-07-28 | 一种基于模拟有限差分法的油水两相流线模拟新方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310946821.3A CN116861753B (zh) | 2023-07-28 | 2023-07-28 | 一种基于模拟有限差分法的油水两相流线模拟新方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116861753A CN116861753A (zh) | 2023-10-10 |
| CN116861753B true CN116861753B (zh) | 2024-05-03 |
Family
ID=88226946
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310946821.3A Active CN116861753B (zh) | 2023-07-28 | 2023-07-28 | 一种基于模拟有限差分法的油水两相流线模拟新方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116861753B (zh) |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2008276357A (ja) * | 2007-04-26 | 2008-11-13 | Toyota Motor Corp | 気液二相流れ評価シミュレーション装置及び気液二相流れ評価シミュレーションプログラム |
| CN110083882A (zh) * | 2019-04-04 | 2019-08-02 | 河海大学 | 一种晃动水槽对波浪形态影响的模拟方法 |
| CN111104766A (zh) * | 2019-12-16 | 2020-05-05 | 中国石油大学(华东) | 基于离散裂缝模型的油水两相非达西渗流数值模拟方法 |
| CN111144008A (zh) * | 2019-12-27 | 2020-05-12 | 西安石油大学 | 一种基于高阶速度场拟合的流线追踪方法 |
| CN111210522A (zh) * | 2020-01-14 | 2020-05-29 | 西南石油大学 | 利用fem在三维非结构网格流场内追踪流线分布的方法 |
| CN112685970A (zh) * | 2020-12-25 | 2021-04-20 | 中国石油大学(华东) | 一种水驱油藏流动单元渗流界面定量表征方法及系统 |
| AU2021104861A4 (en) * | 2020-08-03 | 2021-09-30 | Southwest Petroleum University | Simulation method of unsteady-state gas-water two-phase seepage flow in gas reservoir based on pore-fracture dual media |
| CN113836695A (zh) * | 2021-08-23 | 2021-12-24 | 长江大学 | 一种基于无网格连接元的油藏数值模拟方法 |
| CN114021502A (zh) * | 2021-11-10 | 2022-02-08 | 长江大学 | 基于迎风gfdm的多孔介质油水两相流计算方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2016178934A1 (en) * | 2015-05-01 | 2016-11-10 | Schlumberger Technology Corporation | Multiphase flow in porous media |
| US20230125944A1 (en) * | 2021-10-25 | 2023-04-27 | Qatar Foundation For Education, Science And Community Development | Oil and gas reservoir simulator |
-
2023
- 2023-07-28 CN CN202310946821.3A patent/CN116861753B/zh active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2008276357A (ja) * | 2007-04-26 | 2008-11-13 | Toyota Motor Corp | 気液二相流れ評価シミュレーション装置及び気液二相流れ評価シミュレーションプログラム |
| CN110083882A (zh) * | 2019-04-04 | 2019-08-02 | 河海大学 | 一种晃动水槽对波浪形态影响的模拟方法 |
| CN111104766A (zh) * | 2019-12-16 | 2020-05-05 | 中国石油大学(华东) | 基于离散裂缝模型的油水两相非达西渗流数值模拟方法 |
| CN111144008A (zh) * | 2019-12-27 | 2020-05-12 | 西安石油大学 | 一种基于高阶速度场拟合的流线追踪方法 |
| CN111210522A (zh) * | 2020-01-14 | 2020-05-29 | 西南石油大学 | 利用fem在三维非结构网格流场内追踪流线分布的方法 |
| AU2021104861A4 (en) * | 2020-08-03 | 2021-09-30 | Southwest Petroleum University | Simulation method of unsteady-state gas-water two-phase seepage flow in gas reservoir based on pore-fracture dual media |
| CN112685970A (zh) * | 2020-12-25 | 2021-04-20 | 中国石油大学(华东) | 一种水驱油藏流动单元渗流界面定量表征方法及系统 |
| CN113836695A (zh) * | 2021-08-23 | 2021-12-24 | 长江大学 | 一种基于无网格连接元的油藏数值模拟方法 |
| CN114021502A (zh) * | 2021-11-10 | 2022-02-08 | 长江大学 | 基于迎风gfdm的多孔介质油水两相流计算方法 |
Non-Patent Citations (2)
| Title |
|---|
| A Two-Phase Flow Simulation of Discrete-Fractured Media using Mimetic Finite Difference Method;Zhaoqin Huang et al;Commun. Comput. Phys.;第799-816页 * |
| Zhaoqin Huang et al.A Two-Phase Flow Simulation of Discrete-Fractured Media using Mimetic Finite Difference Method. Commun. Comput. Phys..2014,第799-816页. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116861753A (zh) | 2023-10-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104899380B (zh) | 一种基于蒙特卡洛模拟的边坡稳定可靠度敏感性分析方法 | |
| CN104331621B (zh) | 一种风资源计算方法 | |
| Epshteyn et al. | Fully implicit discontinuous finite element methods for two-phase flow | |
| He et al. | Reduced-order flow modeling and geological parameterization for ensemble-based data assimilation | |
| CN106202746B (zh) | 模拟多孔介质中水流达西速度的Yeh-多尺度有限元方法 | |
| Hægland et al. | Improved streamlines and time-of-flight for streamline simulation on irregular grids | |
| Köppel et al. | A Lagrange multiplier method for a discrete fracture model for flow in porous media | |
| CN102915406A (zh) | 径向流条件下油水相对渗透率曲线的计算方法 | |
| CN104112057A (zh) | 一种大尺度裂缝性油藏数值模拟方法 | |
| CN105868508A (zh) | 一种基于气测录井信息的产能定量预测方法 | |
| CN102339325A (zh) | 一种分析离散裂缝性油藏流体流动的方法 | |
| CN113836695B (zh) | 一种基于无网格连接元的油藏数值模拟方法 | |
| CN110083853B (zh) | 模拟地下水流运动的有限体积Yeh多尺度有限元法 | |
| CN112347678B (zh) | 一种同时模拟地下水流和达西速度的新型多尺度有限元法 | |
| CN116861753B (zh) | 一种基于模拟有限差分法的油水两相流线模拟新方法 | |
| CN114330156A (zh) | 一种基于中心圆柱矩形槽比尺效应的渠道测流法 | |
| Dawson | A continuous/discontinuous Galerkin framework for modeling coupled subsurface and surface water flow | |
| Dong et al. | Adaptive moving grid methods for two-phase flow in porous media | |
| Ersland et al. | Numerical methods for flow in a porous medium with internal boundaries | |
| CN107341336A (zh) | 一种贮箱产品几何精度一致性评价方法 | |
| CN106596335B (zh) | 一种评价无填充粗糙单裂隙中污染物运移特性的方法 | |
| Lamine et al. | Higher order multidimensional upwind convection schemes for flow in porous media on structured and unstructured quadrilateral grids | |
| CN108664678B (zh) | 一种产量预测方法 | |
| CN117421939B (zh) | 一种基于轨迹分段线性化的页岩油裂缝系统模拟代理方法 | |
| Strub et al. | Comparison of two data assimilation algorithms for shallow water flows |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |