CN114021502A - 基于迎风gfdm的多孔介质油水两相流计算方法 - Google Patents

Abstract

本发明公开了一种基于迎风GFDM的多孔介质油水两相流计算方法，包括以下步骤：S1、建立多孔介质油水两相流模型；S2、基于迎风GFDM对所述多孔介质油水两相流模型进行离散，得到所述多孔介质油水两相流模型的全隐式离散格式；S3、利用基于牛顿迭代和自动微分的非线性求解器对所述多孔介质油水两相流模型的全隐式离散格式进行求解，实现压力和含水饱和度的高精度无网格计算。本发明提供一种基于迎风GFDM的多孔介质油水两相流计算方法，为形成不依赖网格的多孔介质流动通用型无网格数值模拟器提供了参考基础。

Description

基于迎风GFDM的多孔介质油水两相流计算方法

技术领域

本发明涉及油水两相流高精度计算领域。更具体地说，本发明涉及一种基于迎风GFDM的多孔介质油水两相流计算方法。

背景技术

广义有限差分法(GFDM)是近些年兴起的一种无网格法，该方法基于子区域内多元函数泰勒级数展开和加权最小二乘近似，将控制方程中未知量的各阶偏导数表示为子区域内相邻节点函数值的线性组合，克服了传统有限差分法对网格的依赖性。目前，该方法在国内外发展迅速，被广泛地应用于求解各种科学和工程问题，包括浅水方程，高阶偏微分方程，瞬态热传导分析，应力分析，水波相互作用，逆热源问题，地震波传播问题，耦合热弹性问题，地下水流动随机分析，以及一些典型的微分方程问题(伯格方程，非线性对流差分方程，时间分数扩散方程)。总得来说，广义有限差分法只需要在计算域内布置一组节点即可实现控制方程的准确求解,可以节省了有限元、有限差分和边界元法中对计算域复杂几何特征费时费力的网格划分及某些无网格法中所需的大量数值积分,并且Benito等提出了一种自适应广义有限差分法，实现了根据精度要求局部范围的自动配点。一方面，这些研究成果表明广义有限差分法是一种高效、高精度的无网格数值建模方法，另一方面，既然GFDM比传统有限差分或其他网格类方法具有理论上的优势，因此GFDM的潜力(例如基于GFDM构建通用性的物理场模拟器)尚未被完全发掘。

对于物理场问题中典型的流动问题，往往需要对方程中的某些物理量进行迎风处理，在无网格方法中，目前一般采用改造影响域的方式来实现上游权，包括向上游方向移动中心节点位置的迎风影响域及将上游节点更多得囊括在中心节点影响域的偏影响域。然而，由于实际流场可能十分复杂，改造影响域的做法难以形成稳定的迎风效果，且计算精度难以得到有效保证。Rao等认为广义有限差分法虽是一种无网格法，但也具有有限差分方法的特征，因此可以在广义有限差分方法中直接借鉴有限差分方法中的迎风格式提出了一种无需改造节点影响域的迎风GFDM方法来处理单相流温耦合问题，并基于该迎风GFDM实现了单相流温耦合的高精度建模。

Rao的工作暗示了迎风GFDM在多孔介质流动问题中的巨大潜力，不混溶两相流是多孔介质流动的一个基本场景，例如地下油藏中的油水两相流。两相流方程中相对渗透率相是一个典型的需要取迎风格式的物理量，为了进一步验证迎风GFDM在多孔介质流动问题中应用的有效性。

发明内容

本发明的目的是提供一种基于迎风GFDM的多孔介质油水两相流计算方法，为形成不依赖网格的多孔介质流动通用型无网格数值模拟器提供了参考基础。

为了实现根据本发明的这些目的和其它优点，提供了一种基于迎风GFDM的多孔介质油水两相流计算方法，包括以下步骤：

S1、建立多孔介质油水两相流模型；

S2、基于迎风GFDM对所述多孔介质油水两相流模型进行离散，得到所述多孔介质油水两相流模型的全隐式离散格式；

S3、利用基于牛顿迭代和自动微分的非线性求解器对所述多孔介质油水两相流模型的全隐式离散格式进行求解，实现压力和含水饱和度的高精度无网格计算。

优选的是，所述的基于迎风GFDM的多孔介质油水两相流计算方法中，所述多孔介质油水两相流模型包括油相物质守恒方程、水相物质守恒方程和辅助方程，其中多孔介质中油水两相流互不相溶；

其中，所述油相物质守恒方程为：

上式中：k是渗透率，mD；k_ro＝k_ro(S_w)是与含水饱和度S_w有关的油相相对渗透率；μ_o是油相粘度，mPa·s；p_o是油相压力，MPa；q_o是油相源汇项，1/天；t是时间，天；S_o是油相饱和度；φ＝φ(p)是与压力有关的储层孔隙度；

所述水相物质守恒方程为：

上式中：k_rw＝k_rw(S_w)是与含水饱和度有关的水相相对渗透率；μ_w是水相粘度，mPa·s；p_w是水相压力，MPa；q_o是水相源汇项，1/天；

所述辅助方程为：

S_o＝1-S_w；

p_w＝p_o-P_c(S_w)；

φ＝φ(p_o)＝φ₀+C_r(p_o-p_o ⁰)；

k_rw＝k_rw(S_w)；

k_ro＝k_ro(S_w)；

(3)

上式中：p_o为油相压力；S_w为含水饱和度；C_t＝C_t(S_w)是综合压缩系数，1/MPa；p_o ⁰是储层原始油相压力，φ₀是p_o ⁰对应的孔隙度。

优选的是，所述的基于迎风GFDM的多孔介质油水两相流计算方法中，S2中得到多孔介质油水两相流计算方法如下：

2.1、基于GFDM的油相物质守恒方程的离散格式为：

上式中：P_c,(i,j) ^t+Δt＝P_c(S_w,(i,j) ^t+Δt)；P_c,i ^t+Δt＝P_c(S_w,i ^t+Δt)；φ_i ^t+Δt＝φ(p_o,i ^t+Δt)；φ_i ^t＝φ(p_o,i ^t)；

2.2、基于GFDM的水相物质守恒方程的离散格式为：

上式中：P_c,(i,j) ^t+Δt＝P_c(S_w,(i,j) ^t+Δt)；P_c,i ^t+Δt＝P_c(S_w,i ^t+Δt)；

2.3、油相的相对渗透率的迎风格式如下：

水相的相对渗透率的迎风格式如下：

本发明的基于迎风GFDM的多孔介质油水两相流计算方法，首次给出了基于迎风GFDM的多孔介质两相流计算方法，实现了对多孔介质两相流的无网格高精度计算。其中提出了GFDM处理多相流问题时无需改造节点影响域的相渗的迎风格式，并通过计算检验了该迎风格式的高性能，从而给出了需要取迎风格式的物理量是基础变量函数时该物理量的高性能迎风格式，进一步完善了迎风GFDM的基础理论。相较于传统网格类方法，本发明的基于迎风GFDM的多孔介质油水两相流计算方法利用无网格节点刻画计算域，避免了计算域高质量网格生成的困难，在处理具有复杂几何特征的计算域时具有更灵活、更简单的优点。

本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。

附图说明

图1为本发明的实施例1中在y＝40m截线上的油相压力及含水饱和度计算结果对比图；

图2为本发明的实施例1中压力及含水饱和度分布计算结果对比图；

图3为本发明的实施例3中计算域示意图；

图4为本发明的实施例3中不规则布点时的计算结果示意图；

图5为本发明的实施例3中不规则边界形状时的计算结果示意图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

需要说明的是，在本发明的描述中，术语“横向”、“纵向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，并不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

<实施例1>

一种基于迎风GFDM的多孔介质油水两相流计算方法，包括以下步骤：

S1、建立多孔介质油水两相流模型；

所述多孔介质油水两相流模型包括油相物质守恒方程、水相物质守恒方程和辅助方程，其中多孔介质中油水两相流互不相溶；

其中，所述油相物质守恒方程为：

上式中：k是渗透率，mD；k_ro＝k_ro(S_w)是与含水饱和度S_w有关的油相相对渗透率；μ_o是油相粘度，mPa·s；p_o是油相压力，MPa；q_o是油相源汇项，1/天；t是时间，天；S_o是油相饱和度；φ＝φ(p)是与压力有关的储层孔隙度；

所述水相物质守恒方程为：

上式中：k_rw＝k_rw(S_w)是与含水饱和度有关的水相相对渗透率；μ_w是水相粘度，mPa·s；p_w是水相压力，MPa；q_o是水相源汇项，1/天；

所述辅助方程为：

S_o＝1-S_w；

p_w＝p_o-P_c(S_w)；

φ＝φ(p_o)＝φ₀+C_r(p_o-p_o ⁰)；

k_rw＝k_rw(S_w)；

k_ro＝k_ro(S_w)；

(3)

上式中：p_o为油相压力；S_w为含水饱和度；C_t＝C_t(S_w)是综合压缩系数，1/MPa；p_o ⁰是储层原始油相压力，φ₀是p_o ⁰对应的孔隙度。

S2、基于迎风GFDM对所述多孔介质油水两相流模型进行离散，得到所述多孔介质油水两相流模型的全隐式离散格式；

具体包括以下步骤：

2.1、基于GFDM的油相物质守恒方程的离散格式为：

上式中：P_c,(i,j) ^t+Δt＝P_c(S_w,(i,j) ^t+Δt)；P_c,i ^t+Δt＝P_c(S_w,i ^t+Δt)；φ_i ^t+Δt＝φ(p_o,i ^t+Δt)；φ_i ^t＝φ(p_o,i ^t)；

2.2、基于GFDM的水相物质守恒方程的离散格式为：

上式中：P_c,(i,j) ^t+Δt＝P_c(S_w,(i,j) ^t+Δt)；P_c,i ^t+Δt＝P_c(S_w,i ^t+Δt)；

2.3、油相的相对渗透率的迎风格式如下：

水相的相对渗透率的迎风格式如下：

GFDM是近年来迅速发展起来的一种无网格配点方法。在局部节点影响域中，该方法利用泰勒展开和移动最小二乘法，得到了未知变量在节点处的各阶偏导数的近似格式。该近似格式是节点影响域中每个节点的未知变量值的线性组合，类似于FDM中的差分表达式，因此被称为GFDM。

对于计算域内某点M₀＝(x₀,y₀)，设该节点影响域内包含着另外n个节点，记作{M₁,M₂,M₃,…,M_N}，其中M_i＝(x_i,y_i)。将这n个节点处的场变量值{u(M_i),i＝1,…n}在节点M₀＝(x₀,y₀)处进行Taylor展开，得到：

其中，x_i＝x₀-x_i,Δy_i＝y₀-y_i。

可以知道，M_i＝(x_i,y_i)离中心节点M₀＝(x₀,y₀)越近，则式(4)中去掉余项后的表达式更应该成立，因此定义式(5)中带权重的误差函数B(D_u)，其中，M_i距离M₀越近，则权重应该越大。

其中，

ω_j＝ω(Δx_j,Δy_j)为权函数，Pan和Li[16]论证了不同的权函数类型对计算结果的影响不大，而一般选取式(6)中的四次样条函数作为式(5)中的权重函数：

其中，

是节点M_j到节点M₀的欧氏距离，r_m是节点M₀影响域半径。

基于移动最小二乘，要求B(D_u)取极小值，则此时B(D_u)关于自变量D_u中各分量的偏导等于零：

将式(11)整理成线性方程组如下：

AD_u＝b

式中，A＝L^TωL，b＝L^TωU，L＝(L₁ ^T,L₂ ^T,…,L_n ^T)^T，

ω＝diag(ω₁ ²,ω₂ ²,…,ω_n ²)，U＝(u₁-u₀,u₂-u₀,…,u_n-u₀)^T。

根据上式，可以推到得到：

D_u＝(u_x0,u_y0,u_xx0,u_yy0,u_xy0)^T＝A^-1b＝A^-1L^TωU＝EU

Where E＝A^-1L^Tω。

从而得到中心节点M₀处一阶即二阶偏导的近似表达式是：

式中，e_1j是矩阵E的元素。

因此，GFDM可根据中心节点影响域内的各节点坐标就可以得到中心节点处一阶导和二阶导的估计表达式，并且可迭代得到中心节点处更高阶导数的估计表达式，从而能够在仅对计算域配点的情况下实现对物理场控制微分方程的离散，而无需复杂的网格剖分。对于具有复杂几何的计算域，无网格GFDM相较于基于网格剖分的有限差分法(也包括有限元、有限体积法等)具有更好的适应性。

虽然GFDM是一种无网格方法，但也具有FDM的特点。因此，GFDM可以直接借鉴FDM中的迎风格式，Rao et al.的工作也证明了迎风GFDM可以在不修改节点影响域的情况下实现单相流温耦合问题中稳定的迎风效果。

S3、利用基于牛顿迭代和自动微分的非线性求解器对所述多孔介质油水两相流模型的全隐式离散格式进行求解，实现压力和含水饱和度的高精度无网格计算。

<实施例2>

本实施例选取了一个规则的矩形计算域([0m,200m]×[0m,80m])，令左右边界(分别记作Γ₁、Γ₂)均为第一类边界条件、上下边界(分别记作Γ₃、Γ₄)为封闭边界条件，源汇项为0，表1展示了包括相渗表达式、计算域初值的相关物性参数取值，则具体方程组见方程(13)。

p₀＝10,S_w＝0.8,

本实施例将小空间步长及小时间步长(Δx＝1m,Δy＝1m,Δt＝0.005d)得到的FDM计算结果作为参考解，图1对比了Δx＝4m,Δy＝4m情况下的FDM及选取节点影响域半径

的迎风GFDM在y＝40m的一维截线上的计算结果,可以看到迎风GFDM基本能取得较高的精度，尤其是在压力计算上，只是在具有间断的水驱前缘处对含水饱和度的计算精度相对较低，但也是能令人满意的。

然后，在迎风GFDM中选取不同的影响域半径

图2直观对比了参考解、FDM解及不同节点影响域半径的含水饱和度及压力计算结果，表2定量得对比了相对计算误差(RE)，可以看到：各种方法对满足扩散方程的压力计算相对误差都较小，对于含水饱和度，影响域半径为

时的迎风GFDM计算精度与迎风FDM计算精度相当，随着影响域半径的增大，迎风GFDM对含水饱和度的计算误差略有下降，这一点从图2中的含水饱和度场图中也能看出，随着影响域半径的增大，含水饱和度前缘间断处的耗散误差增大，但计算误差仍然没有高于FDM很多。对比结果说明了迎风GFDM差不多能够达到同等节点/网格密度下迎风FDM的计算效果，然而本文迎风GFDM是一种无网格方法，能够有效适用于具有复杂几何特征的计算域，而porous flow发生的实际地下储层往往具有复杂的几何特征(断层、复杂边界、裂缝等)，也说明了本文迎风GFDM在porous flow中的巨大应用潜力。

考虑到形成通用的数值模拟器往往需要基于控制方程组的全隐式格式，因此本文采用的是基于牛顿迭代和自动微分的全隐式非线性求解器。表3对比了不同节点影响域情况下的牛顿迭代步数，可以看到不同影响域半径情况下的牛顿迭代步数相同，即影响域的增加没有增加系统的非线性程度，说明在实际计算时可以灵活得选用不同的节点影响域。总得来说，本实施例的结果为基于迎风GFDM构建强鲁棒性的多孔介质流动数值模拟器提供了依据。

表1物理参数

Parameter	Value	Parameter	Value
				Permeability k	100mD	Initial water saturation	0.2
Initial porosityφ<sub>0</sub>	0.3	Irreducible water saturation S<sub>wc</sub>	0.2
				Rock compressibility C<sub>r</sub>	0MPa<sup>-1</sup>	Residual Oil Saturation S<sub>or</sub>	0.2
Oil viscosity	5mPa·s	Tolerated max time step size	20days
				Water viscosity	1mPa·s	Initial time step size	0.01day
Initial pressure	10MPa	Tolerated error	10<sup>-6</sup>

表2计算相对误差表

表3不同节点影响域半径时的牛顿迭代步数

<实施例3>

实施例2论证了迎风GFDM能够取得足够精度的计算结果，并且节点影响域的变化对计算精度的影响有限。因此，相较于依赖高质量网格剖分的有限差分、有限元或有限体积方法等，无网格的迎风GFDM能够利用灵活的布点实现对计算域复杂几何的刻画，从而将模拟计算的前处理难度大大降低。

如图3(a)所示，本实施例针对实施例2，采用更灵活的布点方式，而不再是规则的笛卡尔式的布点，图4展示了计算得到的饱和度及压力分布，可以看到该计算结果与图2(b)中的结果基本是匹配的，说明了本发明迎风GFDM可在任意布点情况下实现有效的计算。然后，本实施例再给出图3(b)所示的一个边界形状不规则的计算域及布点，初边值条件和控制方程与3.1节相同，仅是边界Γ₁和Γ₂的形状与例一中不同，图5展示了计算得到的看着符合物理的含水饱和度与压力分布，说明了本实施例迎风GFDM具有复杂几何的计算域情况下能够通过灵活布点实现对two-phase porous flow的有效计算。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的实施例。

Claims (3)

1.一种基于迎风GFDM的多孔介质油水两相流计算方法，其特征在于，包括以下步骤：

S1、建立多孔介质油水两相流模型；

S2、基于迎风GFDM对所述多孔介质油水两相流模型进行离散，得到所述多孔介质油水两相流模型的全隐式离散格式；

S3、利用基于牛顿迭代和自动微分的非线性求解器对所述多孔介质油水两相流模型的全隐式离散格式进行求解，实现压力和含水饱和度的高精度无网格计算。

2.如权利要求1所述的基于迎风GFDM的多孔介质油水两相流计算方法，其特征在于，所述多孔介质油水两相流模型包括油相物质守恒方程、水相物质守恒方程和辅助方程，其中多孔介质中油水两相流互不相溶；

其中，所述油相物质守恒方程为：

上式中：k是渗透率，mD；k_ro＝k_ro(S_w)是与含水饱和度S_w有关的油相相对渗透率；μ_o是油相粘度，mPa·s；p_o是油相压力，MPa；q_o是油相源汇项，1/天；t是时间，天；S_o是油相饱和度；φ＝φ(p)是与压力有关的储层孔隙度；

所述水相物质守恒方程为：

上式中：k_rw＝k_rw(S_w)是与含水饱和度有关的水相相对渗透率；μ_w是水相粘度，mPa·s；p_w是水相压力，MPa；q_o是水相源汇项，1/天；

所述辅助方程为：

S_o＝1-S_w；

p_w＝p_o-P_c(S_w)；

φ＝φ(p_o)＝φ₀+C_r(p_o-p_o ⁰)；

k_rw＝k_rw(S_w)；

k_ro＝k_ro(S_w)；

(3)

上式中：p_o为油相压力；S_w为含水饱和度；C_t＝C_t(S_w)是综合压缩系数，1/MPa；p_o ⁰是储层原始油相压力，φ₀是p_o ⁰对应的孔隙度。

3.如权利要求1所述的基于迎风GFDM的多孔介质油水两相流计算方法，其特征在于，S2中得到多孔介质油水两相流计算方法如下：

2.1、基于GFDM的油相物质守恒方程的离散格式为：

上式中：P_c,(i,j) ^t+Δt＝P_c(S_w,(i,j) ^t+Δt)；P_c,i ^t+Δt＝P_c(S_w,i ^t+Δt)；φ_i ^t+Δt＝φ(p_o,i ^t+Δt)；φ_i ^t＝φ(p_o,i ^t)；

2.2、基于GFDM的水相物质守恒方程的离散格式为：

上式中：P_c,(i,j) ^t+Δt＝P_c(S_w,(i,j) ^t+Δt)；P_c,i ^t+Δt＝P_c(S_w,i ^t+Δt)；

2.3、油相的相对渗透率的迎风格式如下：

水相的相对渗透率的迎风格式如下：

Images