Sequential Local Mesh Refinement Solver with Separate Temporal and Spatial Adaptivity for Non-linear Two-phase Flow Problems

TL;DR

This paper presents a sequential local mesh refinement method for non-linear two-phase flow problems that enhances convergence and computational efficiency by adaptively refining the mesh in space and time.

Contribution

The proposed framework adaptively refines the mesh in space and time for non-linear flow problems, improving convergence and reducing computational cost without restricting the entire system to small time steps.

Findings

Achieved approximately 25 times speedup in solution time.

Confirmed accuracy of the adaptive algorithm compared to uniform fine discretization.

Effectively prevents convergence failure in large-scale models.

Abstract

Convergence failure and slow convergence rates are among the biggest challenges with solving the system of non-linear equations numerically. Although mitigated, such issues still linger when using strictly small time steps and unconditionally stable fully implicit schemes. The price that comes with restricting time steps to small scales is the enormous computational load, especially in large-scale models. To address this problem, we introduce a sequential local mesh refinement framework to optimize convergence rate and prevent convergence failure, while not restricting the whole system to small time steps, thus improving computational efficiency. We test the algorithm with the non-linear two-phase flow model. Starting by solving the problem on the coarsest space-time mesh, the algorithm refines the domain in time at water saturation front to prevent convergence failure. It then deploys fine spatial grid in regions with large saturation variation to assure solution accuracy. After each refinement, the solution from the previous mesh is used to estimate the initial guess of unknowns on the current mesh for faster convergence. Numerical results are presented to confirm accuracy of our algorithm as compared to the uniformly fine time step and fine spatial discretization solution. We observe approximately 25 times speedup in the solution time by using our algorithm.

Figures

Equations

$$\frac{\partial (\varphi {\rho }_{\alpha }{s}_{\alpha })}{\partial t}+\nabla \cdot {\mathbf{u}}_{\alpha }=q_{\alpha }in\Omega \times J$$

$${\mathbf{u}}_{\alpha }=-K{\rho }_{\alpha }\frac{k_{r\alpha }}{{\mu }_{\alpha }}(\nabla p_{\alpha }-{\rho }_{\alpha }\mathbf{g})in\Omega \times J$$

$${\mathbf{u}}_{\alpha }\cdot \mathbf{\nu }=0on\partial \Omega \times J$$

$$\{\begin{array}{l}\begin{array}{r}{p}_{\alpha }=p_{\alpha }^0\\ s_{\alpha }=s_{\alpha }^0\end{array}\end{array}at\Omega \times \{t=0\}$$

$${\rho }_{\alpha }={\rho }_{\alpha ,ref}\cdot e^{c_{f,\alpha }(p_{\alpha }-p_{\alpha ,ref})}$$

$$k_{r\alpha }=f(s_{\alpha })$$

$$p_c=g(s_{\alpha })=p_{nw}-p_w$$

$$\sum _{\alpha }{s}_{\alpha }=1$$

$$\bm{V}=H(div;\Omega)=\bigg{\{}\bm{v}\in\big{(}L^{2}(\Omega)\big{)}^{d}:\nabla\cdot\bm{v}\in L^{2}(\Omega)\bigg{\}},\quad W=L^{2}(\Omega),$$

$$\bm{V}_{h,i}^{n,m}=\Bigg{\{}\bm{v}\in L^{2}\Big{(}I_{m};H(div;\Omega_{i})\Big{)}:\bm{v}(\cdot,\bm{x})\Big{|}_{F}\in\bm{V}_{h}(F)\ and\ \bm{v}(t,\cdot)\Big{|}_{T}=\displaystyle\sum_{a=1}^{l}\bm{v}_{a}t^{a}\ with\ \bm{v}_{a}\in\bm{V}_{h}(F),\forall\ E\in\mathcal{T}_{h,i}^{n,m}\Bigg{\}},$$

$$W_{h,i}^{n,m}=\Bigg{\{}w\in L^{2}\Big{(}I_{m};L^{2}(\Omega_{i})\Big{)}:w(\cdot,\bm{x})\Big{|}_{F}\in W_{h}(F)\ and\ w(t,\cdot)\Big{|}_{T}=\displaystyle\sum_{a=1}^{l}w_{a}t^{a}\ with\ w_{a}\in W_{h}(F),\forall\ E\in\mathcal{T}_{h,i}^{n,m}\Bigg{\}}.$$

$$\bm{V}_{h}(F)=\bigg{\{}\bm{v}=(v_{1},v_{2})\text{ or }\bm{v}=(v_{1},v_{2},v_{3}):v_{l}=\alpha_{l}+\beta_{l}x_{l};\;\alpha_{l},\beta_{l}\in\mathbb{R},l=1,\cdots,d\bigg{\}}$$

$$W_{h}(F)=\bigg{\{}w\text{ is a constant in }F\bigg{\}}$$

$$\bm{V}^{n,*}_{h}=\Big{(}\oplus_{i,m}\bm{V}_{h,i}^{n,m,0}\Big{)}\oplus\bm{V}^{\Gamma}_{h}$$

$$\displaystyle f_{\tau}(t,\cdot)\Big{|}_{(\tau_{m},\tau_{m+1}]\times\Omega_{i}}=\frac{t-\tau_{m}}{\tau_{m+1}-\tau_{m}}f(\tau_{m+1}^{-},\cdot)\Big{|}_{\Omega_{i}}+\frac{\tau_{m+1}-t}{\tau_{m+1}-\tau_{m}}f(\tau_{m}^{-},\cdot)\Big{|}_{\Omega_{i}},$$

$${\int }_{{\tau }_m}^{{\tau }_{m+1}}{\partial }_t{f}_{\tau }(t,\cdot )=f({\tau }_{m+1},\cdot )-f({\tau }_m,\cdot )$$

$$\int_{J_{n}}\int_{\Omega}\partial_{t}\Big{(}\phi C_{\alpha,h,\tau}^{n}\Big{)}w+\int_{J_{n}}\int_{\Omega}\Big{(}\nabla\cdot\bm{u}_{up,\alpha,h}^{n}\Big{)}w=\int_{J_{n}}\int_{\Omega}q_{\alpha}w\quad\forall w\in W_{h}^{n}$$

$${\int }_{J_n}{\int }_{\Omega }{K}^{-1}{\tilde{\mathbf{u}}}_{\alpha ,h}^n\cdot \mathbf{v}={\int }_{J_n}{\int }_{\Omega }{p}_{\alpha ,h}^n\nabla \cdot \mathbf{v}\forall \mathbf{v}\in {\mathbf{V}}_h^{n,\ast }$$

$${\int }_{J_n}{\int }_{\Omega }{\mathbf{u}}_{\alpha ,h}^n\cdot \mathbf{v}={\int }_{J_n}{\int }_{\Omega }{\lambda }_{\alpha }{\tilde{\mathbf{u}}}_{\alpha ,h}^n\cdot \mathbf{v}\forall \mathbf{v}\in {\mathbf{V}}_h^{n,\ast }$$

$${\lambda }_{\alpha }=\frac{k_{r\alpha }{\rho }_{\alpha }}{{\mu }_{\alpha }}$$

$${\int }_{J_n}{\int }_{\Omega }{\mathbf{u}}_{up,\alpha ,h}^n\cdot \mathbf{v}={\int }_{J_n}{\int }_{\Omega }{\lambda }_{\alpha }^{\ast }{\tilde{\mathbf{u}}}_{\alpha ,h}^n\cdot \mathbf{v}\forall \mathbf{v}\in {\mathbf{V}}_h^{n,\ast }$$

$$\eta_{t,r,\alpha,i}^{n,m}=\big{|}\tau_{m+1}-\tau_{m}\big{|}\Bigg{(}\int_{E_{i}^{m}}\Big{|}\partial_{t}\Big{(}\phi C_{\alpha,h,\tau}^{n}\Big{)}+\nabla\cdot\bm{u}_{up,\alpha,h,\tau}^{n}-q_{\alpha}\Big{|}^{2}\Bigg{)}^{\frac{1}{2}}$$

$$\eta_{s,r,\alpha,i}^{n,m}=\big{|}\Omega_{i}\big{|}\Bigg{(}\int_{E_{i}^{m}}\Big{|}\partial_{t}\Big{(}\phi C_{\alpha,h,\tau}^{n}\Big{)}+\nabla\cdot\bm{u}_{up,\alpha,h,\tau}^{n}-q_{\alpha}\Big{|}^{2}\Bigg{)}^{\frac{1}{2}}$$

$$\eta_{t,f,\alpha,i}^{n,m}=\Bigg{(}\int_{E_{i}^{m}}K^{-1}\Big{|}\bm{u}_{\alpha,h}^{n}-\bm{u}_{\alpha,h,\tau}^{n}\Big{|}^{2}\Bigg{)}^{\frac{1}{2}}$$

$$\eta_{s,f,\alpha,i}^{n,m}=\Bigg{(}\int_{E_{i}^{m}}K^{-1}\Big{|}\bm{u}_{up,\alpha,h}^{n}-\bm{u}_{\alpha,h}^{n}\Big{|}^{2}\Bigg{)}^{\frac{1}{2}}$$

$$\eta_{t,p,\alpha,i}^{n,m}=\Bigg{(}\int_{E_{i}^{m}}K^{-1}\Big{|}\tilde{\bm{u}}_{\alpha,h}^{n}-\tilde{\bm{u}}_{\alpha,h,\tau}^{n}\Big{|}^{2}\Bigg{)}^{\frac{1}{2}}$$

$$\eta_{s,p,\alpha,i}^{n,m}=\Bigg{(}\big{|}\Omega_{i}\big{|}^{2}\int_{E_{i}^{m}}\Big{|}\nabla\times\Big{(}K^{-1}\tilde{\bm{u}}_{\alpha,h}^{n}\Big{)}\Big{|}^{2}+\sum_{e\subset\overline{E_{i}^{m}}}|e|\int_{e}\Big{[}\Big{(}K^{-1}\tilde{\bm{u}}_{\alpha,h}^{n}\Big{)}\times\bm{n}_{e}\Big{]}^{2}\Bigg{)}^{\frac{1}{2}}$$

$$X^{n}=L^{2}\bigg{(}J_{n};H^{1}(\Omega)\bigg{)}\bigcap H^{1}\bigg{(}J_{n};L^{2}(\Omega)\bigg{)},$$

$$\big{\|}\psi\big{\|}_{X^{n}}=\int_{J_{n}}\sum_{\Omega_{i}\in\Omega}\bigg{(}\big{\|}\psi\big{\|}_{L^{2}(\Omega_{i})}^{2}+\big{\|}K^{\frac{1}{2}}\nabla\psi\big{\|}_{L^{2}(\Omega_{i})}^{2}+\big{\|}\psi_{t}\big{\|}_{L^{2}(\Omega_{i})}^{2}\bigg{)}.$$

$$\left|\left|\left|(s_{\alpha }^n-s_{\alpha ,h}^n,p_{\alpha }^n-p_{\alpha ,h}^n,{\tilde{\mathbf{u}}}_{\alpha }^n-{\tilde{\mathbf{u}}}_{\alpha ,h}^n)\right|\right|\right|:=N_{\alpha }^n+N_{\alpha ,p}^n$$

Full text

Keywords.

Space-time domain decomposition, Mixed finite element method, Sequential local mesh refinement, Iterative solver, Non-linear system

1   Introduction

Complex multi-phase flow and reactive transport in subsurface porous media is mathematically modeled by a system of non-linear equations. Due to the significant non-linearity, solving such system with Newton’s method usually suffers from convergence issues even when applying strictly small time steps and using unconditionally stable fully implicit schemes. This problem becomes much more severe in large-scale models. The enormous number of unknowns makes each Newton iteration computationally exhaustive. Therefore by reducing the size of the model using multiscale techniques and optimizing the convergence rate of Newton’s method, we can achieve orders of magnitude greater computational efficiency.

Adaptive homogenization [3, 28] reduces the number of unknowns in the model by replacing fine mesh with coarse mesh in regions where non-linearity and variable (eg. saturation) variation is negligible. However, fine and coarse discretization in space requires different time scales for stable numerical solution. Forcing the coarse mesh to accommodate the fine mesh by taking fine time steps fails to reduce the number of unknowns in time. Space-time domain decomposition addresses this issue by allowing different time scales for different spatial grid. Several space-time domain decomposition approaches have been proposed in the past. [19, 20] introduced space-time finite element method for elastodynamics with discontinuous Galerkin (DG) in time. The method has also been applied to other types of problems such as diffusion with different time discretization schemes [5, 6, 22, 23].

The aforementioned literatures applied space-time domain decomposition method to mostly mechanics problems. On the other hand, prior work regarding flow mainly focused on linear single phase flow and transport problems where flow is naturally decoupled from the advection-diffusion component transport [17, 18]. [31] first validates the space-time approach for non-linear coupled multiphase flow and transport problems on static grid. [31] enforces strong continuity of fluxes at non-matching space-time interfaces with enhanced velocity introduced in [32, 34, 1, 2]. It also constructs and solves a monolithic system to avoid computational overheads associated with iterative solution schemes introduced in [17], that require subdomain problems to be solved iteratively until weak continuity of fluxes is satisfied at interfaces. [30] upgrades the method by allowing adaptive mesh refinement, thus improving computational efficiency while maintaining accuracy as compared to the uniformly fine scale solution. The upgraded algorithm uses initial residual as a cheap error estimator to search for regions that need refinement. A number of other error estimators are mentioned in [26, 33]. As shown in Fig.1, the normalized non-linear residual becomes the largest in regions with the highest non-linearity (water saturation front) and thus consumes most computational resources. Refining such regions in time ensures Newton convergence while refining in space maintains solution accuracy.

The adaptive local mesh refinement approach demonstrated by [30] allows only one level of refinement in both space and time, thus restricting the largest coarse time step allowed for stable numerical convergence. [25] extended such approach by allowing more refinement levels, similar to the algorithm introduced in [8], which however only considers spatial adaptivity. When solving problems on each coarse space-time domain, regions with large non-linear residual and saturation variation are sequentially refined to the finest resolution to ensure solution convergence and accuracy. After each refinement, before solving the problem on the new mesh, the initial guess for the unknowns are populated by the solution on the previous mesh using linear projection. The initial guess provided in such manner is naturally closer to the true solution. Therefore, the non-linear solver convergence is not only guaranteed, but also accelerated. Although achieving 5 times speedup on solution time with iterative linear solver, [25] relied on isotropic space-time refinement which produces a significant number of unnecessary elements. Regions with large saturation variation behind the front are forced to take redundant fine time steps. Preventing such over-refinement will further improve computational performance. Another problem associated with isotropic refinement scheme is that, the error indicator used to pinpoint refinement location combines both temporal and spatial saturation variations. Error indicator calculated in this fashion sometimes misleads the refinement process, especially in channelized permeability field, thus damaging the solution accuracy and hindering numerical convergence. In this work, we extend the method demonstrated in [25] by separating temporal and spatial adaptivity to further improve computational performance and solution accuracy.

In this paper, we restrict ourselves to non-linear two-phase flow problems in subsurface porous media. We intend to approach more complicated non-linear models such as three-phase black oil model in the near future. The rest of the paper begins by describing the governing equations for two phase flow, the functional spaces for space-time domain decomposition, the enhanced velocity weak variational formulation and its fully discrete form in Section 2. Then we will discuss the error estimator analysis used for the refinement process in Section 3. Following the math ground work, we present the algorithm for the sequential local mesh refinement solver in Section 4, including the mechanism to search for refinement regions and separately adapt the mesh in temporal and spatial dimensions. Afterwards, we demonstrate results from two numerical experiments using the proposed algorithm in Section 5.

2   Two phase flow formulation

We consider the following well-known two-phase, slightly compressible flow in porous medium model, with oil and water phase mass conservation, constitutive equations, boundary and initial conditions.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$\phi$ is porosity and $K$ is permeability tensor. $\rho_{\alpha}$, $s_{\alpha}$, $\bm{u}_{\alpha}$ and $q_{\alpha}$ are density, saturation, velocity and source/sink, respectively for each phase. The phases are slightly compressible and the phase densities are calculated by Eqn.(2.5),

[TABLE]

with $c_{f,\alpha}$ being the fluid compressibility and $\rho_{\alpha,ref}$ being the reference density at reference pressure $p_{\alpha,ref}$. In the constitutive equation (2.2) given by Darcy’s law, $k_{r\alpha}$, $\mu_{\alpha}$ and $p_{\alpha}$ are the relative permeability, viscosity and pressure for each phase. Relative permeability is a function of saturation. Pressure differs between wetting phase and non-wetting phase in the presence of capillary pressure, which is also a function of saturation.

[TABLE]

[TABLE]

The saturation of all phases obeys the constrain (2.8).

[TABLE]

The boundary and initial conditions are given by Eqn.(2.3) and (2.4). $J=(0,T]$ is the time domain of interest while $\Omega$ is the spatial domain.

Now we will give a brief introduction of enhanced velocity formulation in space-time domain. Let $J=(0,T]$ be partitioned in to a number of coarse time intervals $\{t_{n}\}_{n=1}^{N}$ where $0=t_{1}<t_{2}<\cdots<t_{N}=T$. $J_{n}=(t_{n},t_{n+1}]$ is the nth partition of the time domain of interest. Consider $J_{n}\times\Omega$ as an union of some non-overlapping subdomains \big{\{}I_{m}\times\Omega_{i}\big{\}}, namely J_{n}\times\Omega=\cup_{i,m}\big{(}I_{m}\times\Omega_{i}\big{)}, where $I_{m}=(\tau_{m},\tau_{m+1}]$ is a sub-interval of $J_{n}=(t_{n},t_{n+1}]$ and $\Omega_{i}$ is a subdomain of $\Omega$. The interfaces of the subdomains are defined as \Gamma_{i,j}^{m,s}=\partial\Big{(}I_{m}\times\Omega_{i}\Big{)}\cap\partial\Big{(}I_{s}\times\Omega_{j}\Big{)}, $\Gamma=\cup_{m,s,i,j}\Gamma_{i,j}^{m,s}$ and \Gamma_{i}^{m}=\Gamma\cap\partial\Big{(}I_{m}\times\Omega_{i}\Big{)}. We use space-time enhanced velocity method similar as [34] to discretize the system. The functional spaces for mixed weak formulation are

[TABLE]

with finite dimensional subspaces as $\bm{V}_{h}$ and $W_{h}$. Let $\mathcal{T}_{h,i}^{n,m}$ be a rectangular partition of $I_{m}\times\Omega_{i}$, $E=T\times F$ be an space-time element in such partition and $\mathcal{T}_{h}^{n}=\cup_{i,m}\mathcal{T}_{h,i}^{n,m}$. Define velocity and pressure/saturation spaces as

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

Functions in $\bm{V}_{h,i}^{n,m}$ and $W_{h,i}^{n,m}$ along time dimension are represented by polynomials with degrees up to $l$. As described in [31], following the discontinuous Galerkin (DG) discretization in time [4, 21], the $DG_{0}$ (polynomial of degree zero) scheme makes \bm{v}(t,\cdot)\big{|}_{T} and w(t,\cdot)\big{|}_{T} constant. Then we define the product spaces as $\bm{V}_{h}^{n}=\oplus_{i,m}\bm{V}_{h,i}^{n,m}$. We remark that $\bm{V}_{h}^{n}$ is not a subspace of $\bm{V}$. To obtain a finite element space containing basis functions with continuous normal flux, we need to modify the basis functions on the space-time interface $\Gamma_{i,j}^{m,s}$. Let $\mathcal{E}_{h,i,j}^{n,m,s}$ be the rectangular partition of $\Gamma_{i,j}^{m,s}$ obtained from the intersection of the traces of $\mathcal{T}_{h,i}^{n,m}$ and $\mathcal{T}_{h,j}^{n,s}$. For each $e\in\mathcal{E}_{h,i,j}^{n,m,s}$, we define a $RT_{0}$ basis function $\bm{v}_{e}$ with a normal component equal to one on $e$ , namely \bm{v}_{e}\big{|}_{e}\cdot\bm{\nu}=1. We then define the space $\bm{V}^{\Gamma}_{h}$ to be the span of all these basis function, $\bm{v}_{e}$. Then the space-time mixed finite element velocity space $\bm{V}^{n,*}_{h}$ is

[TABLE]

where $\bm{V}^{n,m,0}_{h,i}$ is the subspace of $\bm{V}_{h,i}^{n,m}$ with zero normal component on $\Gamma_{i}^{m}$. Similarly, the pressure/saturation space is $W_{h}^{n}=\oplus_{i,m}W_{h,i}^{n,m}$.

Now consider any function $f$ piecewise in time (eg. functions in $\bm{V}_{h}^{n,*}$ and $W_{h}^{n}$), define $f_{\tau}$ as the linear interpolation along time direction as

[TABLE]

and we have

[TABLE]

To simplify the notation, let C_{\alpha,h}^{n}=\rho_{\alpha}\Big{(}p_{\alpha,h}^{n}\Big{)}s_{\alpha,h}^{n} be the phase mass concentration. Then space-time enhanced velocity method formulates Eqn.(2.1) and (2.2) as: find $\bm{u}_{\alpha,h}^{n}\in\bm{V}_{h}^{n,*}$, $\tilde{\bm{u}}_{\alpha,h}^{n}\in\bm{V}_{h}^{n,*}$, $s_{\alpha,h}^{n}\in W_{h}^{n}$, $p_{\alpha,h}^{n}\in W_{h}^{n}$ such that

[TABLE]

[TABLE]

[TABLE]

The mobility ratio in (2.11) is defined as

[TABLE]

and $\bm{u}_{up,\alpha,h}^{n}$ is the upwind velocity calculated by

[TABLE]

The additional auxiliary phase fluxes $\tilde{\bm{u}}_{\alpha,h}^{n}$ is used to avoid inverting zero phase relative permeability [27]. Calculation of the upwind mobility ratio $\lambda_{\alpha}^{*}$ is referred to Eqn.(A.19). The variational form for specifically oil-water system and its fully discrete formulation is attached in A. The discrete formulation provides us a non-linear system of equations for pressure and saturation. We approximate such system in linear form and use Newton iteration to approach the true solution. Depending on the level of non-linearity and the closeness between the initial guess and the true solution, Newton’s method could take numerous iterations before achieving convergence. We will use the sequential local mesh refinement algorithm to accelerate the Newton convergence. Before introducing refinement algorithm, in the next section, we will first present analysis for error estimator used for searching refinement regions.

3   A Posteriori Error Estimate

In this section, we discuss the error estimate analysis as an extension to the work presented in [33]. In contrast to the previous work, our approach to calculate a posteriori error estimate does not rely on computationally expensive local reconstruction of fine scale solution from coarse scale solution. In this section, let $E_{i}^{m}=(\tau_{m},\tau_{m+1}]\times\Omega_{i}\in\mathcal{T}_{h}^{n}$ be a space-time element, we define local error estimators $\eta_{t,r,\alpha,i}^{n,m},\eta_{s,r,\alpha,i}^{n,m},\eta_{t,f,\alpha,i}^{n,m},\eta_{s,f,\alpha,i}^{n,m}$, $\eta_{t,p,\alpha,i}^{n,m}$, $\eta_{s,p,\alpha,i}^{n,m}$ as follow

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Eqn.(3.1)-(3.2) are residual estimators, Eqn.(3.3)-(3.4) are flux estimators and Eqn.(3.5)-(3.6) are nonconformity estimators. Eqn.(3.1) through (3.6) will provide upper bound for the error measure we are about to introduce. It is common to use energy norm as an error measurement for linear problems. However, it is much more complicated for nonlinear problems. Instead we use the dual norm of the residual, which is also widely applied, as our error measure. We denote

[TABLE]

and for any $\psi\in X^{n}$

[TABLE]

Let $s_{\alpha}^{n}$, $p_{\alpha}^{n}$, $\bm{u}_{\alpha}^{n}$ and $s_{\alpha,h}^{n}$, $p_{\alpha,h}^{n}$, $\bm{u}_{\alpha,h}^{n}$ be the exact and numerical saturation, pressure and velocities solutions. The error measure ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ is defined as

[TABLE]

where

[TABLE]

and

[TABLE]

$N^{n}_{\alpha}$ represents the dual norm of the residual and $N^{n}_{\alpha,p}$ measures the non-nonconformity of the numerical solutions.

Assumption 1**.**

We assume there exist a subspace M_{h}^{n}\subset L^{2}\Big{(}J_{n};H(curl,\Omega)\Big{)} such that

[TABLE]

Moreover, for all $v\in H(curl,\Omega)$, there exist $\Pi_{M^{n}_{h}}(v)\in M_{h}^{n}$ such that

[TABLE]

where $e$ is an edge of $E_{i}^{m}$ and $\Pi$ is the $L^{2}$ projection operator.

Next, we will present a posteriori error estimate for this error measure. In the following lemma, we estimate the dual norm of the residual of the mass balance equations.

Lemma 1**.**

Let $\eta^{n,m}_{t,r,\alpha,i},\eta^{n,m}_{s,r,\alpha,i},\eta^{n,m}_{t,f,\alpha,i}$ and $\eta^{n,m}_{s,f,\alpha,i}$ be the error indicators defined in Eqn.(3.1)-(3.4). There exist constants $C,C_{poin}>0$ such that

[TABLE]

Proof.

Since $s_{\alpha}^{n}$, $p_{\alpha}^{n}$ and $\bm{u}_{\alpha}^{n}$ are exact saturation, pressure and velocity, using the phase mass concentration formulation introduced in Section 2 to simplify the notation, for each $\psi\in X^{n}$ we have

[TABLE]

We split the term $\int_{J_{n}}\int_{\Omega}\bm{u}_{\alpha,h,\tau}^{n}\cdot\nabla\psi$ into two parts such that

[TABLE]

Next, we split $\int_{J_{n}}\int_{\Omega}\bm{u}_{\alpha,h}^{n}\cdot\nabla\psi$ into two parts as

[TABLE]

Therefore, by Eqn.(3.7), (3.8) and (3.9), the dual norm can be separated into three terms such that

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

Since $s_{\alpha,h}^{n}$, $p_{\alpha,h}^{n}$, $\bm{u}_{\alpha,h}^{n}$ are the numerical solutions of Eqn.(2.9) and (2.10) we have the following

[TABLE]

We take $w=\Pi_{W_{h}^{n}}\psi$ and by Poincaré inequality obtain the following bound for $I_{1}$

[TABLE]

Next for $I_{2}$, by Cauchy–Schwarz inequality we have

[TABLE]

Similarly, $I_{3}$ goes as

[TABLE]

By the definition of $N_{\alpha}$, the inequality for the dual norm of the residual is proved. ∎

In the following lemma, we will provide an upper bound estimate for the non-nonconformity error measure.

Lemma 2**.**

Assuming $K\in C^{1}(\Omega)$, there exist constant C such that

[TABLE]

Proof.

First, using Helmholtz decomposition, we have

[TABLE]

with $\phi_{0}\in H^{1}(\Omega)$, $\phi_{1}\in H(curl,\Omega)$ for $d=3$ and $\phi_{1}\in H^{1}(\Omega)$ for $d=2$. Since $\psi\in X^{n}$, $\nabla\psi\in L^{2}(J_{n}\times\Omega)$. Meanwhile $\nabla\phi_{0}\in L^{2}(\Omega)$ so $(\nabla\phi_{0})_{\tau}\in L^{2}(J_{n}\times\Omega)$. Thus we have

[TABLE]

By Eqn.(3.18) and $\nabla\phi_{0}\perp\nabla\times\phi_{1}$, we have

[TABLE]

Using Assumption 1, there exist a $\Pi_{M^{n}_{h}}\phi_{1}\in M^{n}_{h}$ such that for all $E_{i}^{m}=(\tau_{m},\tau_{m+1}]\times\Omega_{i}\in\mathcal{T}_{h}^{n}$,

[TABLE]

Then $\nabla\times\Pi_{M^{n}_{h}}\phi_{1}\in V_{h}^{n,*}$, by Eqn.(2.10) we have

[TABLE]

Therefore, we have

[TABLE]

Apply Eqn.(3.21) to Eqn.(3.20) we obtain

[TABLE]

Then the nonconformity error measure is bounded by

[TABLE]

The inequality of the nonconformity error measure is proved. ∎

The error estimator introduced in this section is used to search for refinement regions. In the next section, we will introduce our sequential local mesh refinement algorithm to reduce the size of the system and minimize the number of iterations required for convergence, while maintaining accuracy as compared to uniformly fine scale solution.

4   Solution algorithm

In this section we present the sequential local mesh refinement solver algorithm. The procedure starts by solving the problem at its coarsest resolution in space-time domain and then sequentially refines certain regions to its finest resolution. The coarsest time step is chosen such that the numerical convergence is guaranteed on the coarsest spatial grid. During the sequential refinement process, the solver first keeps the spatial mesh static at its coarsest level and searches for regions to refine in time. Once the last level of temporal refinement is implemented, the temporal discretization is finalized and the solver refines the mesh in space until reaching the finest resolution. Then the spatial grid is restored to the coarsest resolution, the solver marches forward in time with the coarsest time step and the whole process reiterates. The complete algorithm is illustrated in Fig.2. We always start from the coarsest mesh and refines into deeper levels due to the tree data structure inherited from [25, 16, 13]. The tree structure is represented by a group of pointers linked to each other. Allowing both refinement and agglomeration requires inserting and removing pointers in the middle of the tree and then re-associating hanging pointers. The toll caused by such complex operation will counteract the computational efficiency improvement. By constraining the operation to solely refinement, we are only required to evolve the tree by adding new levels on the bottom, which has a much smaller operation count.

Fig.3 demonstrates a sample semi-structured grid generated with the algorithm stated above. Here the $z$ axis represents time in a 2-D spatial problem. Separating temporal and spatial refinement makes the mesh construction more flexible. As observed from the plot, subdomains can have temporal refinement, spatial refinement or both. This flexibility reduces the total number of elements by two to three times as compared to the isotropic refinement scheme implemented in [25], thus improving computational performance. Note that the solver always refine spatially to the finest resolution for cells with well contained, for accurate estimate of production rate and bottom-hole flowing pressure. Adding temporal refinements for these cells depends on whether the saturation front is sweeping through the well or not.

We use the error estimators defined in Section 3 to scan regions for refinement that in return diminishes the upper bound of the error measure. First we study the spatial and temporal residual estimators defined by Eqn.(3.2) and (3.1), similar to [30]. Provided by linear projection of the solution on the previous mesh, the initial guess of unknowns after each refinement procedure is naturally close to the true solution. Consequently, the residual estimators is only useful for indicating the main refinement region on the coarsest space-time resolution. Large estimator values appears only sporadically on all other resolutions as demonstrated by Fig.4. The sporadic appearance infers strong heterogeneity of the underlying petrophysical properties.

The most important estimators to ensure solution convergence and accuracy are the flux estimators. Let us first review the temporal one represented by Eqn.(3.3). We expand the original formulation as the following.

[TABLE]

Then the output is mainly controlled by the temporal flux difference term \bm{u}_{\alpha,h}\big{|}_{\tau_{m+1}}-\bm{u}_{\alpha,h}\big{|}_{\tau_{m}}, which we can further expand to

[TABLE]

The second term in the above equation is effectively zero in slightly compressible flow since density variation caused by pressure is negligible and pressure gradient stays fairly constant in time. Moreover, temporal refinement is not necessary if large estimator output is caused by the leading constant \frac{K}{\mu_{\alpha}}\rho_{\alpha}\big{(}p_{\alpha,h}|_{\tau_{m+1}}\big{)}\nabla p_{\alpha,h}|_{\tau_{m+1}} in the first term (eg. regions around the well with large pressure gradient), since pressure solution is smooth in time and does not trigger any convergence issues. We need to apply temporal refinement in regions with large $\eta_{t,f,\alpha,i}^{n,m}$ caused specifically by significant change in relative permeability k_{r\alpha}\big{(}s_{\alpha,h}|_{\tau_{m+1}}\big{)}-k_{r\alpha}\big{(}s_{\alpha,h}|_{\tau_{m}}\big{)}. Therefore we calculate the temporal water saturation gradient

[TABLE]

and applied refinement exclusively to regions with both $\eta_{t,f,\alpha,i}^{n,m}$ and $\varepsilon_{t,i}^{n,m}$ values exceeding the threshold.

Similarly, the spatial flux estimator $\eta_{s,f,\alpha,i}^{n,m}$ in Eqn.(3.4) can be expanded to

[TABLE]

with the output mainly controlled by the flux spatial difference $\bm{u}_{up,\alpha,h}-\bm{u}_{\alpha,h}$. We can also expand this term as the following.

[TABLE]

Surely, we need to refine regions with significant change in relative permeability k_{r\alpha}\big{(}s_{up,\alpha,h}\big{)}-k_{r\alpha}\big{(}s_{\alpha,h}\big{)} to accurately represent the features of the reservoir. Furthermore regions with large estimator output caused by the leading constant $\frac{K}{\mu_{\alpha}}\rho_{\alpha}\nabla p_{\alpha,h}$ also need special care. Such regions are characterized by rapid mass flow and refining them facilitates convergence. Therefore we calculate the spatial water saturation gradient.

[TABLE]

and apply refinement to regions with either $\eta_{s,f,\alpha,i}^{n,m}$ or $\varepsilon_{s,i}^{n,m}$ values exceeding the threshold. Please note that we are taking some extra steps when calculating spatial saturation gradient by looking at the previous time step values. This mechanism ensures a more accurate exposure of features in the system, especially in channelized permeability distributions. Fig.5 and 6 show flux estimator and saturation gradient at each refinement level.

Instead of setting a subjective threshold, we outline the regions for refinement by distribution percentile. We first define $[0.01,1]$ as the analysis range of flux error estimator and saturation gradient. Values below $0.01$ are considered too small and thus neglected. The threshold is determined by distribution. The cumulative distribution function of flux error estimator and saturation gradient at each refinement level is plotted in Fig.7 against sample data collected during simulation. As illustrated by the graphs, the data for both variables generally follow log-normal distribution trend. During temporal adaptation, we refine cells with $50\%$ largest values in both flux estimator and saturation gradient. Therefore, we use the log-mean, which covers approximately $50\%$ of the analysis range, as the threshold. We notice that during temporal refinement, the cumulative distribution functions of both variables are better described by normal distribution. However, we still choose the log-mean as the threshold since it leads to a slight over-refinement in time and thus better guarantees Newton convergence. During spatial adaptation, we refine cells with either $50\%$ largest values in the saturation gradient or $10\%$ largest values in flux estimator. Therefore, the thresholds for the two variables are the log-mean and one standard deviation above the log-mean, corresponding to their respective distribution.

Finally for the non-conformity estimators, $\eta_{t,p,\alpha,i}^{n,m}$ is naturally diminished during temporal refinement using $\eta_{t,f,\alpha,i}^{n,m}$ since the two estimators have similar formulations. $\eta_{s,p,\alpha,i}^{n,m}$ represents the tangential gradient of flux on non-conformal grid interfaces. To reduce this term, we apply mesh smoothing algorithm. Adjacent grid cells cannot be more than one refinement level apart within the hierarchical tree structure in both spatial and temporal dimensions. Such algorithm also ensures a smooth transition from fine grid into coarse grid and thus facilitates convergence.

5   Numerical results

In this section we will show results from two numerical experiments on 2-D two-phase flow model. Both experiments use the same fluid data from the SPE10 dataset [7]. The oil and water reference densities in Eqn.(2.5) are taken to be $53\ lb/ft^{3}$ and $64\ lb/ft^{3}$ and compressibilities are $1\times 10^{-4}\ psi^{-1}$ and $3\times 10^{-6}\ psi^{-1}$ respectively. We use Brooks’s Corey model for both relative permeability and capillary pressure. The equations for relative permeability are

[TABLE]

The endpoint values are $s_{or}=s_{wirr}=0.2$ and $k_{ro}^{0}=k_{rw}^{0}=1.0$ while the model exponents are $n_{w}=n_{o}=2$. The equation for capillary pressure is

[TABLE]

with $P_{en,cow}=10\ psi$ and $l_{cow}=0.2$. Fig.8 visualizes the relative permeability and capillary pressure curve. The two experiments uses gaussian-like and channelized permeability and porosity distributions from SPE 10 dataset [7] layer 20 and 52, respectively. The reservoir size is $56ft\times 216ft\times 1ft$. We place a water rate specified injection well at the bottom left corner and a pressure specified production well at the upper right corner. The water injection rate is $1\ ft^{3}/day$ and production pressure is $1000\ psi$. Furthermore, the initial pressure and water saturation are set to be $1000\ psi$ and $0.2$.

5.1 Gaussian-like Permeability Distribution

The gaussian-like permeability field comes from SPE 10 dataset layer 20. The fine scale petrophysical data are shown in Fig.9, assuming isotropic permeability. We allow three refinement levels in both space and time in our experiment. Although the framework allows different refinement ratios between levels, for the sake of simplicity we set the same ratio, a factor of 2, between all levels. We use the numerical homogenization technique introduced in [3] to upscale the fine scale permeability to different coarse levels. This calculation only needs to be performed once at the beginning of the experiment. The homogenized permeability distribution in $X$ and $Y$ directions, which does not manifest high anisotropy, is illustrated in Fig.10. The porosity is upscaled simply by weighted volumetric average and therefore is not visualized. The computational domain is $56\ ft\times 216\ ft\times 1\ ft\times 1000\ days$ with coarsest and finest element size of $8\ ft\times 8\ ft\times 1\ ft\times 10\ days$ and $1\ ft\times 1\ ft\times 1\ ft\times 1.25\ days$.

The adaptive water saturation profile with its mesh as compared to fine scale solution at 100 and 500 days are plotted in Fig.11. We observe the finest mesh stays concentrated at the water front to correctly capture the dramatic changes in saturation. In this region, mass transfer is not dominated by either oil or water phase and thus contributes the most non-linearity and requires temporal refinement for stable Newton convergence. Elements behind the water front is gradually coarsened due to the decreased saturation variation. Overall, the saturation profile provided by the sequential refinement solver looks similar to the fine scale solution. Fig.12 shows the production rates and cumulative recoveries of the two solutions, which are nearly identical. The oil rate from the sequential refinement solver appear to be slightly smoother at the early time and has tiny oscillation around water break though, which is caused by the coarse mesh.

The program execution time is presented in Fig.13. The total execution time consists of system setup which constructs the linear system, solving the linear system and data handle which mainly involves copying and pasting data from the current to the previous time step. Since the experiment problem size is still small, we use both direct and iterative solver to resolve the linear system. The semi-structured space-time mesh results in highly non-symmetric matrices and therefore we use GMRES with ILU preconditioner as our iterative solver. We observe 8 and 4 times speedup on system setup and data handle using direct solver. These two types of operations are highly dependent upon the number of time steps taken and total number of refinement levels. Hence, the speedup scales linearly with the total temporal refinement ratio and the same runtime reduction behavior is observed when using iterative solver. The speedup on solving the linear system best represents the computational performance improvement. Since our problem size is small, the efficiency gain is not substantial when using direct solver. On the contrary, we observe 25 times speed up on solving the linear system when using iterative solver. Additional techniques on solving non-symmetric linear systems iteratively, such as relaxing linear solver tolerance using forcing function [12, 24] and applying specialized preconditioners [14, 10] for Krylov-based method, may be utilized for additional acceleration. Note that as we move towards more complex models such as 3-D black oil, the solution to the corresponding linear system is only accessible through iterative methods. Thus we should expect significant improvement on computational efficiency once we approach those types of problems.

5.2 Channelized Permeability Distribution

The channelized permeability field comes from SPE 10 dataset layer 52. The fine scale petrophysical data are shown in Fig.14. Unlike the gaussian case, the permeability and porosity here is highly structured and some structures are subjected to destruction during the homogenization process. We allow two refinement levels in both space and time for this experiment. The refinement ratio is also set to a factor of 2 between all levels. During numerical homogenization, we impose oversampling technique introduced by [11] and [9] to preserve channel connectivity as much as possible. The homogenized permeability distribution in $X$ and $Y$ direction is illustrated in Fig.15. On the contrary to the gaussian case, the upscaled channel permeability is highly anisotropic. The computational domain is $56\ ft\times 216\ ft\times 1\ ft\times 1000\ days$ with coarsest and finest element size of $8\ ft\times 8\ ft\times 1\ ft\times 5\ days$ and $1\ ft\times 1\ ft\times 1\ ft\times 1.25\ days$.

The adaptive water saturation profile with its mesh as compared to the fine scale solution at 200 and 400 days are plotted in Fig.16. The overall production profile also resembles each other between the two solutions. Here, the fine mesh not only concentrates at the water front, but also outlines the channel structure. The channel boundary is characterized by dramatic contrast of permeability, thus resulting in steep water saturation gradient. The refinement algorithm detects these features and deploys mesh with appropriate size accordingly. Many low permeability spots inside the main high permeability channel are also accurately identified and represented. Fig.17 shows the production rates and cumulative recoveries of the two solutions. The adaptive and fine scale rates also look similar, however with obvious discrepancies. The rates from sequential refinement solver looks smoother than the fine scale solution at early time. It also suffers from slightly early water breakthrough. The oil and water cumulative production from the two solutions nearly overlap.

We also approach the solution by both direct and iterative method. The program execution time is shown in Fig.18. The speedup on system setup and data handle also scales linearly with total temporal refinement ratio. The solution time reduction by direct solver remains low. However, we still observe a 25 times speedup using iterative solver, even when the coarse time step is halved as compared to the gaussian case. The substantial improvement is caused by two main reasons. First of all, the flow and transport is constrained within the channel structure behind the water front, making the saturation variation effectively zero in other part of the reservoir. Consequently, the number of grid cells required to represent the channel structure and saturation front is relatively small, causing the adaptive solution easier to acquire. Secondly, the fine scale system consists of dramatic permeability contrast, resulting in the related linear system to have eigenvalues close to zero. Solving such linear system with Krylov-based iterative methods requires many iterations, making the fine scale solution harder to obtain.

6   Conclusions

We have introduced an algorithm that constructs adaptive mesh using error estimators to solve non-linear two-phase flow problems with reduced execution time. The procedure sequentially refines the mesh from coarsest to finest resolution in large non-linearity regions, with temporal and spatial adaptivity separated to accurately expose features in the system while ensuring numerical convergence. After each refinement, the initial guess for the new mesh is generated by the solution on the previous mesh through linear projection, which accelerates convergence rate. Results from two numerical experiments are demonstrated. Rates and cumulative production from both experiments resembles well between the adaptive and fine scale solution. The water saturation profiles also look similar. We observe approximately 25 times speedup in solution time for the gaussian-like and the channelized permeability field. The channel case suffers from a slightly early water breakthrough, which could be mitigated through loosening refinement criterion. However, doing so will counteract the computational efficiency improvement. With the promising results from two-phase flow problems, we plan to test our algorithm on more complex models such as 3-D three-phase black oil system in the near future.

Acknowledgement

The author would like to thank the Center for Subsurface Modeling for supporting this research. We would also like to extend our gratitude to Gurpreet Singh and Uwe Köcher for the discussion related to this work.

A   Fully Discrete Formulation

Consider the oil-water system, the variational form of Eqn.(2.1) through (2.4) is: find $\bm{u}_{\alpha,h}^{n}\in\bm{V}_{h}^{n,*}$, $\tilde{\bm{u}}_{\alpha,h}^{n}\in\bm{V}_{h}^{n,*}$, $s_{w,h}^{n}\in W_{h}^{n}$, $p_{o,h}^{n}\in W_{h}^{n}$ such that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

for all $\bm{v}\in\bm{V}_{h}^{n,*}$ and $w\in W_{h}^{n}$. The conversion between auxiliary and actual phase flux is referred to Eq.(2.11). The oil saturation and water pressure are eliminated by the saturation constrain and the capillary pressure relation (assume oil phase being the non-wetting phase).

For the fully discrete formulation, we will start by stating the basis functions in $RT_{0}\times DG_{0}$ discretization scheme. In spatial dimensions, the pressure and saturation are piecewise constants while velocity is piecewise linear. Meanwhile all variables are piecewise constants in temporal dimension as stated in Section 2. Let $E_{i}^{m}=(\tau_{m},\tau_{m+1}]\times\Omega_{i}$ be a space-time element, we have

[TABLE]

[TABLE]

The solution to Eqn.(A.1) through (2.11) can be written in discrete form using basis functions as

[TABLE]

We remove the superscript $n$ and subscript $h$ in the above solution variables for this section since we need to use $n$ to pair basis functions. While keeping the solution in discrete form, we now substitute the testing functions in the variational forms of mass conservation and constitutive equation with $w_{j}^{n}$ and $\bm{\varphi}_{j+\frac{1}{2}}^{n}$. For the first term in Eqn.(A.3) and (A.4) we obtain

[TABLE]

Here, \Big{|}e_{j+\frac{1}{2}}^{n}\Big{|} is an edge of a space-time element. Since the framework uses backward Euler scheme in time to avoid Courant-Fredricks-Levy condition, we have the construction

[TABLE]

The second term in Eqn.(A.3) and (A.4) can be reformulated as

[TABLE]

When non-matching grid caused by different time scales at $(j+\frac{1}{2})^{-}$ and $(j+\frac{1}{2})^{+}$ is encountered, assume the ratio between coarse and fine time step to be $\frac{\delta t_{c}}{\delta t_{f}}=\ell$, then for each $0\leq k\leq\ell-1$

[TABLE]

The variational form of capillary pressure term can be revised in similar way as Eqn.(A.10) and (A.11). Now we evaluate the mass conservation equation. The first term in Eqn.(A.2) becomes

[TABLE]

In fine time scales, Eqn.(A.12) can be altered as follow.

[TABLE]

The second term is calculated as

[TABLE]

The approach to handle non-matching grid is a little different for this term. Assume the fine time partition stays on $(j+\frac{1}{2})^{-}$ side, then on fine time elements we have

[TABLE]

while for the coarse time element we have

[TABLE]

Eqn.(A.13) and (A.16) will cause the accumulation and transmissibility matrix to have extra temporal bands forming in the lower triangle, making the corresponding linear system non-symmetric. The oil phase mass conservation equation is similar. Combining the equations for both phases will provide the expression for the total mass conservation equation. The two sides of Eqn.(2.11) is estimated as

[TABLE]

[TABLE]

The $\lambda_{\alpha,j+\frac{1}{2}}^{*,n}$ is the upwind mobility for stable numerical solution and is defined as

[TABLE]

The matrix corresponding to the above discrete formulation has sparsity pattern of three, five or seven non-zero diagonals, depending on the spatial dimension of the problem, with one extra temporal diagonal in the lower triangle. Forming such matrix in block format is referred to [31, 15, 29].