A priori error analysis for transient problems using Enhanced Velocity approach in the discrete-time setting

TL;DR

This paper provides theoretical error estimates for transient subsurface flow problems using the Enhanced Velocity Mixed FEM combined with backward Euler and Crank-Nicolson methods, aiding in optimal time and space discretization.

Contribution

It introduces a priori error analysis in discrete-time for transient problems using Enhanced Velocity Mixed FEM, extending previous adaptive mesh refinement applications.

Findings

Validated backward Euler error estimates through numerical tests

Error estimates assist in choosing appropriate time step and mesh size

Enhanced Velocity scheme is effective for heterogeneous porous media

Abstract

Time discretization along with space discretization is important in the numerical simulation of subsurface flow applications for long run. In this paper, we derive theoretical convergence error estimates in discrete-time setting for transient problems with the Dirichlet boundary condition. Enhanced Velocity Mixed FEM as domain decomposition method is used in the space discretization and the backward Euler method and the Crank-Nicolson method are considered in the discrete-time setting. Enhanced Velocity scheme was used in the adaptive mesh refinement dealing with heterogeneous porous media [1, 2] for single phase flow and transport and demonstrated as mass conservative and efficient method. Numerical tests validating the backward Euler theory are presented. This error estimates are useful in the determining of time step size and the space discretization size. References. [1] Yerlan Amanbek, Gurpreet Singh, Mary F Wheeler, and Hans van Duijn. Adaptive numerical homogenization for upscaling single phase flow and transport. ICES Report,12:17, 2017. [2] Gurpreet Singh, Yerlan Amanbek, and Mary F Wheeler. Adaptive homogenization for upscaling heterogeneous porous medium. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, 2017.

Tables

Table 1. Table 1: Accuracy results of pressure and velocity for various levels.

Level	$h$	$H$	$\Delta t$	$erro{r}_p$	$erro{r}_{𝐮}$
1	1/52	1/26	1/50	6.33e-04	1.51e-01
2	1/100	1/50	1/100	3.32e-04	1.02e-01
3	1/120	1/60	1/120	2.79e-04	9.15e-02
4	1/152	1/76	1/150	2.26e-04	7.97e-02

Table 2. Table 2: Accuracy results of pressure and velocity for various levels.

Level	$h$	$H$	$\Delta t$	$erro{r}_p$	$erro{r}_{𝐮}$
1	1/100	1/50	1/7	7.28e-01	8.64e-01
2	1/128	1/64	1/8	6.38e-01	7.71e-01
3	1/164	1/82	1/9	5.69e-01	6.87e-01
4	1/200	1/100	1/10	5.13e-01	6.24e-01

Equations

$$\frac{\partial }{\partial t}\left(\varphi \rho \right)+\nabla \cdot (\rho \mathbf{u})=q\text{in}\Omega \times J$$

$$\varphi \frac{\partial \rho }{\partial t}=\varphi \frac{\partial \rho }{\partial p}\frac{\partial p}{\partial t}=\varphi C_1\frac{\partial p}{\partial t}$$

$$\text{u}=-\frac{\mathbf{K}}{\mu }\left(\nabla p-\rho \text{g}\right),$$

$$p=g\text{on}\partial \Omega \times J.$$

$$p(x,0)=p^0(x).$$

$$=-\mathbf{K}\nabla p\text{ in}\Omega \times J,$$

$$\frac{\partial p}{\partial t}+\nabla \cdot \text{u}$$

$$p$$

$$k_{min}{\xi }^T\xi \le {\xi }^T\mathbf{K}(x)\xi \le k_{max}{\xi }^T\xi \forall x\in \Omega \forall \xi \in {\mathbb{R}}^d,d=1,2,3.$$

$$\left(K^{-1}\text{u},\mathbf{v}\right)$$

$$\left(\frac{\partial p}{\partial t},w\right)+\left(\nabla \cdot \text{u},w\right)$$

$$\displaystyle\left(p,w\right)\biggr{\rvert}_{t=0}$$

$$\left({\mathbf{K}}^{-1}{\text{u}}_h,\mathbf{v}\right)$$

$$\left(\frac{\partial p_h}{\partial t},w\right)+\left(\nabla \cdot {\text{u}}_h,w\right)$$

$$\displaystyle\left(p_{h},w\right)\biggr{\rvert}_{t=0}$$

$$\left({\mathbf{K}}^{-1}(\text{u}-{\text{u}}_h),\mathbf{v}\right)-\left(p-p_h,\nabla \cdot \mathbf{v}\right)$$

$$\left(\frac{\partial }{\partial t}(p-p_h),w\right)+\left(\nabla \cdot (\text{u}-{\text{u}}_h),w\right)$$

$$E_p^I=p-\hat{p},$$

$$E_{\mathbf{u}}^I=\mathbf{u}-{\Pi }^{\ast }\mathbf{u},$$

$$\left(p-\hat{p},w\right)=\left(E_p^I,w\right)=0\forall w\in W_h.$$

$$⟨{\Pi }^{\ast }\mathbf{q}\cdot \nu ,1{⟩}_e=⟨\mathbf{q}\cdot \nu ,1{⟩}_e$$

$$\left(\nabla \cdot ({\Pi }^{\ast }\mathbf{q}-\mathbf{q}),w\right)=0\forall w\in W_h$$

$$\left(\nabla \cdot \left({\Pi }^{\ast }\mathbf{u}-\mathbf{u}\right),w\right)=\left(\nabla \cdot E_{\mathbf{u}}^I,w\right)=0\forall w\in W_h$$

$$\left(\frac{\partial }{\partial t}{E}_p^I,w\right)=0\forall w\in W_h$$

$$\Vert E_p^I\Vert \le C{\Vert p\Vert }_r{h}^{r}0\le r\le 1,$$

$$\Vert E_{\mathbf{u}}^I\Vert \le C{\Vert \mathbf{u}\Vert }_{1}h.$$

$$ab\le \frac{1}{2\epsilon }{a}^2+\frac{\epsilon }{2}{b}^2$$

$$\Vert \nabla \cdot {\mathbf{u}}_h\Vert \le C{h}^{-1}\Vert {\mathbf{u}}_h\Vert$$

$$f^n=f(x,t_n),0\le n\le N,$$

$$f^{n,\theta }=\frac{1}{2}(1+\theta )f^{n+1}+\frac{1}{2}(1-\theta )f^n,0\le n\le N-1.$$

$${\Vert f\Vert }_{l^{\infty }(L^2)}=\underset{0\le n\le N}{\max }{\Vert f^n\Vert }_{L^2}$$

$${\Vert f\Vert }_{l^2(L^p)}={\left(\sum _{n=0}^{N-1}{\Vert f^{n,\theta }\Vert }_{L^p}^2\Delta t\right)}^{\frac{1}{2}}.$$

$$\begin{split}f^{n+1}=f\biggr{\rvert}_{t=t^{n,\theta}}+\frac{1}{2}(1-\theta)\Delta t\frac{\partial f}{\partial t}\biggr{\rvert}_{t=t^{n,\theta}}+\frac{1}{8}(1-\theta)^{2}(\Delta t)^{2}\frac{\partial^{2}f}{\partial t^{2}}\biggr{\rvert}_{t=t^{n,\theta}}+O(\Delta t^{3})\end{split}$$

Full text

Keywords.  a priori error analysis, enhanced velocity, mixed finite element method, error estimates, Darcy flow.

1 Introduction

The most subsurface flow equations are dynamic and time-dependent problems. In decision making process, numerical simulation of flow plays vital role in many engineering applications such as oil and gas production evaluation, CO2 sequestration and contaminate transport problems. It is natural to deal with non-matching multiblock grids in the reservoir simulation since subsurface parameters such as permeability or porosity can vary over subdomains substantially. The accuracy of simulation can depend on discretization method of space and time variables. For space discretization, we are concerned with a well-established domain decomposition method, i.e. Enhanced Velocity Mixed Finite Element Method (EVMFEM), which provides similar accuracy as the Multiscale Mortar Mixed FEM [3]. EVMFEM is a mass conservative and an efficient domain decomposition method. By using this method, several applications such as single, two-phase flow, bio-remediation simulation and others were considered in [3, 4]. Recently, an adaptive mesh refinement strategy, which is based on Enhanced Velocity scheme, has been proposed in the numerical simulations of flow and transport through heterogeneous porous media [1, 2, 5, 6]. Such a novel approach demonstrates the efficiency and accuracy of simulation in the heterogeneous porous media by allowing to capture an important features of flow and transport problems.

A little attention has been given to discrete time setting analysis. Theoretical convergence analysis of EVMFEM has been shown in [4] for slightly compressible single phase flow for general continuous in time approximations. However, we could not find the numerical error tests that compare with analytical solution. A few authors have begun to implement time domain decomposition method to be flexible on selection of time step size [7, 8, 9]. The key idea is to extend the EVMFEM in space to time discretization by constructing a monolithic system without subdomain iteration.

In this paper, we are concerned with the solution of time dependent problem that is discretized by the backward Euler method or the Crank-Nicolson method combined with EVMFEM, which uses the lowest order Raviart-Thomas spaces on non-matching multiple subdomains. In particularly, we focus on deriving a priori error estimates for the transient subsurface problems in discrete-time setting. This gives asymptotive behavior of the numerical error for a given mesh size, time step size and others. This analsys allows us to conclude about the convergence and the guarantee of stability of the numerical method. The reader is referred to [10, 12] for more different time and space discretization in the transient problems. We also provide numerical tests of the error estimate.

This paper is organized as follows. In the next section, we describe the slightly compressible flow model formulation as well as give a strong and weak formulation using Enhanced Velocity space. Section 3 is devoted to the error analysis with preliminary projections, definitions and discrete formulation. This analysis carried out for backward Euler scheme and Crank-Nicolson schemes simultaneously in time discretization settings. Numerical results are presented in Section 4; Section 5 concludes the paper.

2 Model Formulation

We describe in this section the slightly compressible flow formulate with initial and boundary conditions. Next, the transient problem is presented in the strong and weak formulations.

2.1 Slightly Compressible Flow Formulation

Our focus is a single phase and slightly compressible fluid in heterogeneous porous media. The classical mass conservation equation is defined by

[TABLE]

where $\Omega\in\mathbb{R}^{d}(d=1,2$ or $3$), $J=(0,T]$, $d$ is the number of spatial dimensions, $q$ is the source/sink term, $\phi$ is the porosity, $\rho$ is the phase density, and u is the phase velocity. We remark that a Peacemann correction is used for modeling source/sink terms [11].

In slightly compressible fluid, the phase density is given by $\rho=\rho_{ref}e^{C_{f}(p-p_{ref})}$, where where, $C_{f}$ is the fluid compressibility, and $\rho_{ref}$ is the reference density at reference pressure $p_{ref}$. Using Taylor series expansion we obtain $\rho\approx\rho_{ref}(1+C_{f}(p-p_{ref}))$. Then it follows that

[TABLE]

for invariant-in-time $\phi$ and for $C_{1}=C_{f}\rho_{ref}$. The phase velocity u is defined by Darcy’s law as,

[TABLE]

where, $\mu$ is the viscosity, $\mathbf{K}$ is the permeability (absolute permeability) tensor, $\rho$ is the density of the fluid and g is the gravity vector. Although more general global boundary conditions can also be treated, we restrict ourselves to the following,

[TABLE]

Additionally, the initial condition is given by,

[TABLE]

From now on our analysis focus on a transient (parabolic) problems which might be involved to various applications problems.

2.2 Transient problem with EVMFEM

We start with the strong formulation of the transient (slightly compressible) flow problems governing single phase flow model for pressure $p$ and the velocity u, which is also case of slightly compressible single phase flow model:

[TABLE]

where $\Omega\in\mathbb{R}^{d}(d=2$ or $3$) is multiblock domain, $J=[0,T]$ and $\mathbf{K}$ is a symmetric, uniformly positive definite tensor representing the permeability divided by the viscosity with $L^{\infty}(\Omega)$ components, for some $0<k_{min}<k_{max}<\infty$

[TABLE]

The Dirichlet boundary condition is considered for convenience. A weak solution of parabolic Eqns. (6) - (9) is a pair $\{\mathbf{u}_{h},p_{h}\}:J\rightarrow\mathbf{V}^{*}_{h}\times W_{h}$, i.e. Enhanced Velocity Mixed Finite Element approximation

[TABLE]

In addition, there is an initial condition

[TABLE]

We formulate the variational problem in semi-discrete space as: Find $\{\mathbf{u}_{h},p_{h}\}:J\rightarrow\mathbf{V}^{*}_{h}\times W_{h}$ such that

[TABLE]

In addition, there is an initial condition

[TABLE]

Subtracting Eqns. (11) - (12) from Eqns. (14) - (15) yields

[TABLE]

3 Error estimates

In this section, we start with preliminaries including projections operators and some notations. Using them we present discrete formulations for analysis. Next, we derive auxiliary error estimates and a priori error estimate theorems.

3.1 Projections

We shall write a projection and define auxiliary error of pressure and velocity as follows:

[TABLE]

Note that $\mathbf{u}-\mathbf{u}_{h}=E^{I}_{\mathbf{u}}+E^{A}_{\mathbf{u}}$, $p-p_{h}=E^{I}_{p}+E^{A}_{p}$. We used $\hat{p}$ the $L^{2}$-projection of $p$ that is defined as

[TABLE]

We know from original work [3] that the projection operator $\Pi^{*}$ was introduced and was utilized for a priori error analysis of elliptic problems. For convenience of the reader, we repeat the relevant and brief definition. Thus, we denote by $\Pi^{*}$ the projection operator that maps $(H^{1}(\Omega))^{d}$ onto $\mathbf{V}^{*}_{h}$ that defined locally for any element $T\in\mathcal{T}_{h}$ and any $\mathbf{q}\in(H^{1}(T))^{d}$ such that for all $\mathbf{q}\in(H^{1}(T))^{d}$

[TABLE]

where $e$ is either any edge in 2D (or face in 3D) of $T$ not lying on $\Gamma$ or an edge in 2D (or face in 3D) of a sub-element, $T_{k}$. Such projection is developed prior to conducting error analysis for a priori estimate. As can be seen in Figure 1, $T_{k}$ has a common edge with the interface grid $\mathcal{T}^{\Gamma}$. According to divergence theorem, we have

[TABLE]

For $\mathbf{u}\in H^{1}(\Omega)$

[TABLE]

Lemma.

[TABLE]

Proof.

$0=\frac{d}{dt}\left(E^{I}_{p},w\right)=\left(\frac{\partial}{\partial t}E^{I}_{p},w\right)+\cancelto{0}{\left(E^{I}_{p},w_{t}\right)}$, by 21, since $w_{t}\in W_{h}$. ∎

Useful inequalities of projections , see [3]:

[TABLE]

[TABLE]

Recall Young’s Inequality: for $a,b\geq 0$

[TABLE]

The Inverse Inequality can be given as

[TABLE]

In this inequality, we have been working under the assumption that $\mathcal{T}_{h,i}$ quasi-uniform rectangular partition of $\Omega_{i}$.

3.2 Definitions

In this section, we make analysis of discrete in time error estimates. Firstly, some definitions are made: for $\Delta t=\frac{T}{N}$, $N$ is a positive integer, $t_{n}=n\Delta t$ and for given $\theta\in[0,1]$,

[TABLE]

[TABLE]

Let’s make also the following definitions:

[TABLE]

We note that the difference can be expressed $t^{n,\theta}-t^{n}=\frac{1}{2}(1+\theta)\Delta t$. Next, using the Taylor series expansion about $t=t^{n,\theta}$, for any sufficiently smooth function $f(t)$, we obtain:

[TABLE]

[TABLE]

After multiplying the first equation by $\frac{1}{2}(1+\theta)$ and the second equation by $\frac{1}{2}(1-\theta)$ and then summing them, we obtain

[TABLE]

Note that if $\theta=1$ then f^{n,\theta}=f\biggr{\rvert}_{t=t^{n,\theta}}+O(\Delta t^{3}). In addition, we can get second order approximation of $\Delta t$, details in [13] : $p(\mathbf{x},t^{n,\theta})\approx p^{n,\theta}\text{ and }u(\mathbf{x},t^{n,\theta})\approx u^{n,\theta}.$ According to Taylor series expansion [13], we obtain

[TABLE]

where $\rho^{p,n,\theta}$ depends on time-derivatives of $p$ and $\Delta t$

[TABLE]

so $\left\lVert\rho^{p,n,\theta}\right\rVert=\mathcal{O}(\Delta t^{2-\theta})$.

3.3 Discrete formulation

We formulate variational form in semi-discrete space as: Find $\{\mathbf{u}_{h},p_{h}\}:J\rightarrow\mathbf{V}^{*}_{h}\times W_{h}$ such that

[TABLE]

In addition, there is an initial condition

[TABLE]

where $l_{1}$ and $l_{2}$ are bounded linear functionals, i.e.

[TABLE]

With this definitions, Eqns. (34) - (35) become as: Find $\{\mathbf{u}^{n,\theta}_{h},p^{n,\theta}_{h}\}\in\mathbf{V}^{*}_{h}\times W_{h}$, $n=1,2,...,N-1$, such that

[TABLE]

Note that if $\theta=1$ then the time discretization is the backward Euler(Implicit) method, and if $\theta=0$ then the Crank-Nicolson scheme.

We consider true solution $\mathbf{u}\in L^{2}\left(J,\mathbf{V}\right)$ and $p_{h}\in H^{1}\left(J,W\right)$ of Eqns. (11) and (12) at time $t=t^{n,\theta}$ in the continuous in time with spatially discrete scheme. We used Eqn. (32) and additional remark related to the Taylor series expansion in order to obtain the following equations with at least order of $\mathcal{O}(\Delta t)$:

[TABLE]

We approximate the vector integrals type $\left(\mathbf{v},\mathbf{q}\right)_{T,M}$ by trapezoidal-midpoint quadrature rules and $\left(\mathbf{K}^{-1}\mathbf{q},\mathbf{v}\right)_{T}$ by trapezoidal quadrature rules respectively. In [14], the equivalence between finite volume methods and the mixed finite element method was established for special quadrature rule for $\mathbf{K}$ diagonal tensor and using the lowest-order Raviart Thomas spaces on rectangles. We emphasize that the EVMFEM with special quadrature and velocity elimination in the discrete system can be reduced to well-known a cell-centered finite difference method.

For each time step we use the Newton method to solve the system, in case of the slightly compressible flow: the fluid compressibility $C_{f}$ term brings us to a nonlinear system. That is why consideration of it would be beneficial for nonlinear problems in the future.

3.4 Analysis

We first derive the bounds of auxiliary error terms.

Theorem 1** (Auxiliary error estimate).**

For the velocity $\mathbf{u_{h}}$ and pressure $p_{h}$ of the mixed method spaces $\textbf{V}^{*}_{h}\times W_{h}$ satisfying equations (37) - (38), assume $\Delta t$ is sufficiently small and positive, $\mathbf{K}$ is uniformly positive definite and sufficient regularity of true solution in equations (6)-(9). Then, there exist a constant $C$ such that

[TABLE]

where $C=C(T,\mathbf{K},\mathbf{u},p)$ and

[TABLE]

Proof.

Subtracting Eqns. ($\ref{eq:3_8})-(\ref{eq:3_9}$) from Eqns. ($\ref{eq:3_6p})-(\ref{eq:3_7p}$) respectively yields

[TABLE]

Take $\mathbf{v}=\Pi^{*}\mathbf{u}^{n,\theta}-\mathbf{u}^{n,\theta}_{h}=E_{\mathbf{u}}^{A\;n,\theta}$ and $w=E_{p}^{A\;n,\theta}$ in (43) and (42) respectively.

[TABLE]

[TABLE]

After adding them, we can rewrite as

[TABLE]

[TABLE]

[TABLE]

It follows immediately that

[TABLE]

By using 45, multiply by $2\Delta t$ and sum from [math] to $N-1$ in Equation 44.

[TABLE]

[TABLE]

[TABLE]

We use the Holder inequality and the Young inequality to get

[TABLE]

We should note that

[TABLE]

since \nabla\cdot E^{A}_{u}\biggr{\rvert}_{\Omega\backslash\Omega^{*}}\in W_{h} and the property 21.

[TABLE]

Remark 1*.*

We used the following properties:

[TABLE]

and

[TABLE]

Next, we know that

[TABLE]

Therefore,

[TABLE]

We can multiply by 2 and make $\varepsilon$ small enough in order to have LHS with positive coefficients. Later take minimum and divide both side of inequality.

[TABLE]

We are thus apply the discrete Gronwall lemma, for sufficiently small $\Delta t$, to obtain:

[TABLE]

where

[TABLE]

We used the fact that $\sum^{N-1}_{n=0}\Delta tg_{n}\leq CT\sum^{N-1}_{n=0}g_{n}$ , $|\Omega^{*}|\leq Ch$ and property that is given in eqn. (33). This finishes the proof of theorem. ∎

The auxiliary error estimates theorem allows us to conclude the following theorem:

Theorem 2** (Error estimate).**

Assume the same conditions as in the previous theorem. Then,

[TABLE]

where $C=C(T,\mathbf{K},\mathbf{u},p)$ and

[TABLE]

Proof.

By applying triangle inequality,the Interpolation Error Inequalities and Theorem 1 results we obtain:

[TABLE]

∎

4 Numerical Examples

In this section, we conduct numerical experiment to verify the numerical accuracy of parabolic problem solution using EVMFEM in space and the backward Euler in time. Based on our a priori error analysis estimates we assume a sufficiently smooth analytical solution. In numerical examples, we set $\Omega=\left(0,1\right)\times\left(0,1\right)$, $\mathbf{K}_{i,j}=\delta_{i,j}$ and the domain $\Omega$ is divided into four subdomains $\Omega_{i}$; $\Omega_{1}$ and $\Omega_{4}$ have fine grids, $\Omega_{2}$ and $\Omega_{3}$ have coarse grids; such mesh discretization is illustrated in Figure 2.

4.1 Numerical example 1

We use the known solution

[TABLE]

and use it to compute the forcing $f$, the Dirichlet boundary data $g$, and the initial data $p_{0}$. We carry out several levels of uniform grid refinement in each subdomains. The time step and the element size are almost equal to each other, see Table 1. The simulation time interval is $(0;0.1)$, i.e. $T=0.1$, and we use the Backward Euler method to integrate with regard to time with uniform time step. We are interested in finding the exact error using a given true solution, so the pressure is true error and the velocity error is normalized error. On applying sufficient Newton iterations at each time step provided the residual is within the machine-precision tolerance, we obtain the numerical solution for evaluating of the error in specified norm. We compute $error_{p}$, corresponds to $\left\lVert p-p_{h}\right\rVert_{l^{\infty}(L^{2})}$, which is maximum of values among time steps that resulting for given time step a discrete pressure $L^{2}$-norm that associates only the function values at the cell-centers in space. Also, $error_{u}$ is defined as $\left\lVert\mathbf{u}-\mathbf{u}_{h}\right\rVert_{l^{2}(L^{2})}$ where in space a discrete $L^{2}$-norm that associates only the normal vector components at the midpoint edges and then normalized by $\left\lVert\mathbf{u}\right\rVert_{L^{2}}$ and $l^{2}$-norm in time. The convergence rate is illustrated in Figure 3.

4.2 Numerical example 2

We use manufactured the known solution

[TABLE]

and use it to compute the forcing $f$, the Dirichlet boundary data $g$, and the initial data $p_{0}$. The time step is equal to the root of the coarse mesh size. Thus first order convergence is expected from theoretical result. Such mesh discretization are depicted in the Figure 2. The simulation time interval is $(0,2)$, i.e. $T=2$, and we use the Backward Euler method to integrate with regard to time with uniform time step, see Table 2. The convergence rate is illustrated in Figure 4.

5 Conclusion

This research has provided a priori error analysis for transient problems or slightly compressible flow problems through the heterogeneous porous media using Enhanced Velocity scheme as the domain decomposition method in space that coupled with backward Euler or Crank-Nicolson method in the time setting. In these discretization settings, we obtained the first order convergence rate for the backward Euler method and the second order convergence rate for the Crank-Nicolson method. Numerical experiments are provided. The results suggest that this approaches could also be useful for the engineering subsurface applications including CO2 sequestration, etc. In our future research, we plan to concentrate on parareal algorithms to achieve efficiency in time discretization that allow to run simulation efficiently for the long time range.

Acknowledgements

First author thanks Drs. T. Arbogast and I. Yotov for some helpful discussions during analysis of method.