Recovery of the Interface Velocity for the Incompressible Flow in Enhanced Velocity Mixed Finite Element Method

TL;DR

This paper improves the accuracy of interface velocity approximation in incompressible flow simulations using an enhanced velocity mixed finite element method and a post-processed pressure reconstruction, leading to better convergence.

Contribution

It introduces a novel velocity reconstruction method at interfaces using post-processed pressure to enhance approximation accuracy in mixed finite element methods.

Findings

Improved convergence rate at interfaces.

Enhanced velocity approximation accuracy.

Numerical validation confirms effectiveness.

Abstract

The velocity, coupling term in the flow and transport problems, is important in the accurate numerical simulation or in the posteriori error analysis for adaptive mesh refinement. We consider Enhanced Velocity Mixed Finite Element Method for the incompressible Darcy flow. In this paper, our aim to study the improvement of velocity at interface to achieve the better approximation of velocity between subdomains. We propose the reconstruction of velocity at interface by using the post-processed pressure. Numerical results at the interface show improvement on convergence rate.

Tables

Table 1. Table 1: Convergence test 1, velocity and recovered velocity error using the post-processed pressure at interface.

$n$	$e_{{𝐮}_h,\Gamma }$		$e_{{\widetilde{𝐮}}_h,\Gamma }$
	error	order	error	order
8	1.47e-01	—–	3.55e-01	—–
16	7.70e-02	0.93	1.12e-01	1.67
32	3.94e-02	0.97	3.73e-02	1.58
48	2.65e-02	0.98	2.09e-02	1.43

Table 2. Table 2: Convergence test 2: velocity and recovered velocity error using the post-processed pressure at interface.

$n$	$e_{{𝐮}_h,\Gamma }$		$e_{{\widetilde{𝐮}}_h,\Gamma }$
	error	order	error	order
8	1.78e-01	—–	3.78e-01	—–
16	8.89e-02	1.00	1.00e-01	1.91
32	4.43e-02	1.00	2.87e-02	1.81
48	2.96e-02	1.00	1.56e-02	1.51

Equations

$${\Vert \mathbf{u}-{\mathbf{u}}_h\Vert }_{\Omega }\le C\left({\Vert p\Vert }_{1,\infty ,{\Omega }^{\ast }}+{\Vert \mathbf{u}\Vert }_{1,\Omega }{h}^{1/2}\right)h^{1/2}$$

$${\Vert \mathbf{u}-{\mathbf{u}}_h\Vert }_{{\Omega }^{{}^{\prime }}}\le C_{\epsilon }\left({\Vert p\Vert }_{1,\infty ,{\Omega }^{\ast }}+{\Vert \mathbf{u}\Vert }_{1,\Omega }\right)h^{r-\epsilon }$$

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

$$\nabla \cdot \text{u}$$

$$p$$

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

$$\left(\nabla \cdot \text{u},w\right)$$

$$\Omega ={\left(\bigcup _{i=1}^{N_b}{\bar{\Omega }}_i\right)}^o,{\Gamma }_{i,j}=\partial {\Omega }_i\bigcap \partial {\Omega }_j,\Gamma ={\left(\bigcup _{i,j=1}^{N_b}{\bar{\Gamma }}_{i,j}\right)}^o,{\Gamma }_i={\Omega }_i\bigcap \Gamma =\partial {\Omega }_j\setminus \partial \Omega .$$

$${\mathbf{V}}_h(T)=\{\mathbf{v}=(v_1,v_2)\text{or}\mathbf{v}=(v_1,v_2,v_3):v_l={\alpha }_l+{\beta }_l{x}_l:{\alpha }_l,{\beta }_l\in \mathbb{R};l=1,..d\},$$

$$W_h(T)=\{w=\text{constant}\}.$$

$$\displaystyle\mathbf{V}_{h,i}=\{\mathbf{v}\in H({\rm div};\Omega_{i}):\mathbf{v}\biggr{\rvert}_{T}\in\mathbf{V}_{h}(T),\forall T\in\mathcal{T}_{h,i}\}\qquad i\in\{1,...n\}$$

$${\mathbf{V}}_h=\underset{i=1}{\overset{n}{⨁}}{\mathbf{V}}_{h,i}.$$

$$\displaystyle\mathbf{v}_{T_{k}}\cdot\nu=\begin{cases*}1,\qquad{\rm on}\;e_{k}\\ 0,\qquad\rm other\;edges\end{cases*}$$

$${\mathbf{V}}_h^{\ast }=\underset{i=1}{\overset{n}{⨁}}{\mathbf{V}}_{h,i}^0⨁{\mathbf{V}}_h^{\Gamma }\cap H(\mathrm{div};\Omega ).$$

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

$$\left(\nabla \cdot {\text{u}}_h,w\right)$$

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

$$\left(\nabla \cdot {\text{u}}_h,w\right)$$

$$\displaystyle\left(\mathbf{v},\mathbf{q}\right)_{M,T}=\begin{cases*}\left(v_{x},q_{y}\right)_{T_{x},M_{y}}+\left(v_{y},q_{y}\right)_{M_{x},T_{y}}\quad&\text{if } d =2,\\ \left(v_{x},q_{y}\right)_{T_{x},M_{y},M_{z}}+\left(v_{y},q_{y}\right)_{M_{x},T_{y},M_{z}}+\left(v_{z},q_{z}\right)_{M_{x},M_{y},T_{z}}\quad&\text{if } d=3.\end{cases*}$$

$${\left({\mathbf{K}}^{-1}{\mathbf{u}}_h,\mathbf{v}\right)}_{{\Omega }_i,M,T}-{\left(p_h,\nabla \cdot \mathbf{v}\right)}_{{\Omega }_i}$$

$${\left(\nabla \cdot {\mathbf{u}}_h,w\right)}_{{\Omega }_i}$$

$$\sum _{i=1}^{N_b}\{{\left({\mathbf{K}}^{-1}{\mathbf{u}}_h,{\mathbf{v}}^{EV}\right)}_{{\Omega }_i,M,T}-{\left(p_h,\nabla \cdot {\mathbf{v}}^{EV}\right)}_{{\Omega }_i}\}$$

$${\mathbf{V}}_h^{\ast }=\left({\mathbf{V}}_{h,1}^0\oplus {\mathbf{V}}_{h,2}^0\oplus {\mathbf{V}}_h^{\Gamma }\right)$$

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

$${\mathbf{u}}_h^{\Gamma }=A_1{p}_L+A_2{p}_R$$

$${\mathcal{P}}_h^i:{\mathbf{V}}_h^{\ast }(\Gamma )\to {\mathbf{V}}_{h,i}({\Gamma }_i)for\psi \in L^2(\Gamma ),⟨\left(\psi -{\mathcal{P}}_h^i\psi \right)\cdot {\mathbf{\nu }}_i,\mathbf{v}\cdot {\mathbf{\nu }}_i{⟩}_{\Gamma }=0\forall \mathbf{v}\in {\mathbf{V}}_{h,i}.$$

$${\mathbf{u}}_{h,i}^{\Gamma }={\mathcal{P}}_h^i{\mathbf{u}}_h^{\Gamma },i=LorR$$

$${\left({\mathbf{K}}^{-1}{\tilde{\mathbf{u}}}_h,\mathbf{v}\right)}_{M,T,{\Omega }_i}$$

$${\text{u}}_{h,L}^{\Gamma }=A_1^L{p}_L+A_2^L{p}_L^e$$

$${\mathcal{P}}_h^L{\text{u}}_h^{\Gamma }={\mathcal{P}}_h^L{A}_1{p}_L+{\mathcal{P}}_h^L{A}_2{p}_R$$

$$A_2^L{p}_L^e={\mathcal{P}}_h^L{A}_2{p}_R$$

$$p_L^e={\left(A_2^L\right)}^{-1}{\mathcal{P}}_h^L{A}_2{p}_R.$$

$$p_R^e={\left(A_1^R\right)}^{-1}{\mathcal{P}}_h^R{A}_1{p}_L$$

$$⟨{\lambda }_{h,T},{\mathbf{v}}_h\cdot {\mathbf{n}}_T{⟩}_e:={\left({\mathbf{K}}^{-1}{\mathbf{u}}_h,{\mathrm{v}}_{\mathrm{h}}\right)}_T-{\left(p_h,\nabla \cdot {\mathbf{v}}_h\right)}_T\forall {\mathbf{v}}_h\in {\hat{\mathbf{V}}}_h\left(T\right)$$

Full text

Keywords.  domain decomposition, velocity improvement, enhanced velocity, mixed fem.

1 Introduction

The numerical reservoir simulations have been utilized in many subsurface applications such as groundwater remediation, reservoir well evaluation, and contaminate transport problems. For such applications, it is common to deal with the flow and transport problem. The main component or coupling term of the flow and transport systems is the velocity and its accuracy the mostly achieved by employing classical mixed finite element system. Due to the heterogeneity of porous media multiphysics problems could be categorized systematically in which one physical phenomena influences within a subdomain and another physical phenomena dominates within another subdomain. Such solutions are coupled through continuity of normal flux at interface, shared region between differently discretized subdomains. To deal with these problems there are the well-known methods such as Multiscale Mortar and Enhanced Velocity schemes that are established in various applications. Recently, a novel adaptive method was studied in subsurface applications [1, 2, 3, 4] using Enhanced Velocity scheme. The main idea is here to utilize the EVMFEM as domain decomposition method to couple different discretized subdomains with more accurate upscaled subsurface parameters.

In the simulation of flow with adaptivity, the results obtained in [5] suggest that pressure values could be interpolated using neighboring elements values to approximate auxiliary pressure values within provided elements. Selection of interpolants is based on convex combinations of vertical and horizontal oriented pressure values. In related reference [6], it was studied that the interface error of solution between subdomains for different numerical methods including Mortar Multiscale MixedFEM which provided better approximation for second-order elliptic problems. One of reasons is the iterative procedure in the mortar scheme that is a key in coupling two subdomains physics. According to author in [6] mortar scheme is general method in coupling for practical multiphysics problems. On the other hand, the efficient Enhanced Velocity scheme has not been investigated from the point of view of the improvement solution including velocity at interface in the previous studies.

The challenge here is to construct the velocity approximation of EVMFEM and specifically at interface to have a better velocity between subdomains that leads accurate approximation in the flow and transport problems. In [7], a priori error analysis states that the global error is

[TABLE]

and away from the interface $\Gamma$ the velocity error convergence rate is better, since

[TABLE]

where $\varepsilon>0$, $r=1$ if $d=2$ and $r=5/6$ if $d=3$, and $\Omega^{{}^{\prime}}_{i}$ is compactly contained in $\Omega_{i}$, $\Omega^{{}^{\prime}}=\bigcup_{i=1}^{N_{b}}\Omega^{{}^{\prime}}_{i}$. This implies that the discrete velocity should be approximated more precise near interface region $\Omega^{*}$. On the question of pressure approximation, the convergence rate of pressure approximation is $\mathcal{O}(h^{1})$, if $d=2$, and $\mathcal{O}(h^{5/6})$, if $d=3$ [7, 8]. If one compare the error of velocity (1) and pressure approximation these results indicate that the velocity convergence rate is not strong as pressure in $\Omega$. Similar a priori error result was shown in [3] for transient problems. Nevertheless, there are still problems including the velocity approximation at the interface to be addressed.

In this paper, we introduce the way to improve velocity accuracy at interface in the Enhanced Velocity MFEM for incompressible flow using the post-processed pressure from [9]. This improvement is important in flow coupled with transport problems and it also can be a good candidate for a recovery-based error estimate evaluation. In a recent work [4], a posteriori error analysis was shown for the incompressible flow problems without recovery of velocity.

The remaining part of the paper proceeds as follows. Section 2 of this paper will describe model formulation with different view of EVMFEM. In Section 3, the proposed numerical method will be discussed. Section 4 shows numerical results. Section 5 summarizes the results of this work and draws conclusions.

2 Model formulation

We start by giving the model formulation for the incompressible single-phase flow. For the convenience of reader we repeat the relevant material of domain decomposition method, discrete formulation with Enhanced Velocity from [7]. We next describe the proposed different view of Enhanced Velocity Discrete Scheme with projection operator.

2.1 Governing equations of the incompressible flow

We consider the incompressible single-phase flow model for pressure $p$ and the Darcy velocity u:

[TABLE]

where $\Omega\in\mathbb{R}^{d}(d=2$ or $3$) is multiblock domain, $f\in L^{2}(\Omega)$ 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$ $k_{min}\xi^{T}\xi\leq\xi^{T}\mathbf{K}(x)\xi\leq k_{max}\xi^{T}\xi\qquad\forall x\in\Omega\quad\forall\xi\in\mathbb{R}^{d}$, under the Dirichlet boundary condition.

A weak variational form of the fluid flow problem $(\ref{eq:a})-(\ref{eq:c})$ is to find a pair $\textbf{u}\in\mathbf{V}$, $p\in W$

[TABLE]

where ${\bm{\nu}}$ is the outward unit normal to $\partial\Omega$, V is $H({\rm div};\Omega)=\{\mathbf{v}\in\left(L^{2}(\Omega)\right)^{d}:\nabla\cdot\mathbf{v}\in L^{2}(\Omega)\}$ and equipped with the norm $\left\lVert\mathbf{v}\right\rVert_{V}=\left(\left\lVert\mathbf{v}\right\rVert^{2}+\left\lVert\nabla\cdot\mathbf{v}\right\rVert^{2}\right)^{\frac{1}{2}}$ and the pressure the space is $W=L^{2}(\Omega)$ and the corresponding norm $\left\lVert w\right\rVert_{W}=\left\lVert w\right\rVert.$.

Discrete formulation

Let $\Omega$ be decomposed into non-overlapping small subdomains, see Fig. 1. We consider

[TABLE]

This implies that the domain is divided into $N_{b}$ subdomains, the interface between $i^{th}$ and $j^{th}$ subdomains($i\neq j$), the interior subdomain interface for $i^{th}$ subdomain and union of all such interfaces, respectively. Let $\mathcal{T}_{h,i}$ be a conforming, quasi-uniform and rectangular partition of $\Omega_{i}$, $1\leq i\leq N_{b}$, with maximal element diameter $h_{i}$. We then set $\mathcal{T}_{h}=\cup_{i=1}^{n}\mathcal{T}_{h,i}$ and denote $h$ the maximal element diameter in $\mathcal{T}_{h}$; note that $\mathcal{T}_{h}$ can be nonmatching as neighboring meshes $\mathcal{T}_{h,i}$ and $\mathcal{T}_{h,j}$ need not match on $\Gamma_{i,j}$. We assume that all mesh families are shape-regular.

In Enhanced Velocity scheme setting, the velocity basis functions are based on the traditional Raviart-Thomas spaces of lowest order on rectangles for $d=2$ and bricks for $d=3$. The $RT_{0}$ spaces are defined for any element $T\in\mathcal{T}_{h}$ by the following spaces:

[TABLE]

The pressure finite element approximation space on $\Omega$ is taken to be as W_{h}(\Omega)=\{w\in L^{2}(\Omega):w\biggr{\rvert}_{E}\in W_{h}(T),\forall T\in\mathcal{T}_{h}\}. In addition, a vector function in $\mathbf{V}_{h}$ can be determined uniquely by its normal components $\mathbf{v}\cdot\nu$ at midpoints of edges (in 2D) or face (in 3D) of $T$. The degrees of freedom of $\mathbf{v}\in\mathbf{V}_{h}(T)$ were created by these normal components. The degree of freedom for a pressure function $p\in W_{h}(T)$ is at center of $T$ and piecewise constant inside of $T$.

Let us formulate $RT_{0}$ space on each subdomain $\Omega_{i}$ for partition $\mathcal{T}_{h}$

[TABLE]

and then

[TABLE]

Although the normal components of vectors in $\mathbf{V}_{h}$ are continuous between elements within each subdomains, the reader may see $\mathbf{V}_{h}$ is not a subspace of $H({\rm div};\Omega)$, because the normal components of the velocity vector may not match on subdomain interface $\Gamma$.

Let us define $\mathcal{T}_{h,i,j}$ as the intersection of the traces of $\mathcal{T}_{h,i}$ and $\mathcal{T}_{h,j}$, and let $\mathcal{T}^{\Gamma}_{h}=\bigcup_{1\leq i\leq j\leq N_{b}}\mathcal{T}_{h,i,j}$. We require that $\mathcal{T}_{h,i}$ and $\mathcal{T}_{h,j}$ need to align with the coordinate axes. Fluxes are constructed to match on each element $e\in\mathcal{T}^{\Gamma}_{h}$. We consider any element $T\in\mathcal{T}_{h,i}$ that shares at least one edge with the interface $\Gamma$, i.e., $T\cap\Gamma_{i,j}\neq\emptyset$, where $1\leq i,j\leq N_{b}$ and $i\neq j$. Then newly defined interface grid introduces a partition of the edge of $T$. This partition may be extended into the element $T$ as shown in Fig. 2.

Such partitioning helps to construct fine-scale velocities that is in $H(\textbf{ div},\Omega)$. So we represent a basis function $\mathbf{v}_{T_{k}}$ in the $\mathbf{V}_{h}(T_{k})$ space ($RT_{0}$) for given $T_{k}$ with the following way:

[TABLE]

i.e. a normal component $\mathbf{v}_{T_{k}}\cdot\nu$ equal to one on $e_{k}$ and zero on all other edges(faces) of $T_{k}$. Let $\mathbf{V}^{\Gamma}_{h}$ be span of all such basis functions defined on all sub-elements induced the interface discretization $\mathcal{T}_{h,i,j}$. Thus, the enhanced velocity space $\mathbf{V}^{*}_{h}$ is taken to be as

[TABLE]

where $\mathbf{V}^{0}_{h,i}=\{\mathbf{v}\in\mathbf{V}_{h,i}:\mathbf{v}\cdot\nu=0\text{ on }\Gamma_{i}\}$ is the subspace of $\mathbf{V}_{h,i}$. The finer grid velocity allows to velocity approximation on the interface and then form the $H({\rm div},\Omega)$ conforming velocity space. Some difficulties arise, however, in analysis of method and implementation of robust linear solver for such modification of $RT_{0}$ velocity space at all elements, which are adjacent to the interface $\Gamma$. We now formulate the discrete variational form of equations $(\ref{eq:a})-(\ref{eq:c})$ as: Find $\mathbf{u}_{h}\in\mathbf{V}^{*}_{h}$ and $p_{h}\in W_{h}$ such that

[TABLE]

2.2 A different view of the EVMFEM in the Discrete Variational Formulation

We consider the discrete variational form that is given in (8)-(9). Find $\mathbf{u}_{h}\in\mathbf{V}^{*}_{h}$ and $p_{h}\in W_{h}$ such that

[TABLE]

We exploit the approximation inner product and for $\mathbf{v},\mathbf{q}\in\mathbb{R}^{d}$

[TABLE]

where $T_{(\cdot)}$ and $M_{(\cdot)}$ denote the the trapezoidal and midpoint quadrature rules in each coordinate direction respectively, see [10]. In particularly, we take $\mathbf{v}=\mathbf{K}^{-1}\mathbf{u}_{h}$ and $\mathbf{q}=\mathbf{v}$.

It is easily proven that the finite variational form (10)-(11) is equivalent to finding $\mathbf{u}_{h}\in\mathbf{V}^{*}_{h}$, $p_{h}\in W_{h}$, $1\leq i\leq N_{b}$, such that

[TABLE]

We note that similar the discrete variational formulation was proposed in [11] with conjugate gradient method. We want to share the idea for small number of discretization elements that can be applied for a large number of elements. Thus, we consider two subdomains,i.e., $\Omega=\bar{\Omega}_{1}\cup\bar{\Omega}_{2}$ and $\Gamma$ is the interface. Then

[TABLE]

Consider equations

[TABLE]

These allow us to express $\mathbf{u}^{\Gamma}_{h}$ in terms of the one-element layers along $\Gamma$, it is shown in Fig. 3:

[TABLE]

Now we consider each subdomain separately with ghost layers. We define $L_{2}$-projection of Enhanced Velocity space at interface $\Gamma_{i,j}$ to each subdomain space $\partial\Omega_{i}\cap\Gamma_{i,j}$ such that $\mathcal{P}^{i}_{h}:\mathbf{V}^{*}_{h}\rightarrow\mathbf{V}_{h,i}$.

[TABLE]

We denote

[TABLE]

In subdomain $\Omega_{i}$, we define $p^{e}_{i}$ in the following way

[TABLE]

where $\Omega^{*}$ is union of all elements $T$ that shares edge (2D) or face (3D) with $\Gamma_{i}$ and $p^{e}_{i}$ ghost layers pressure values, and $\widetilde{\mathbf{u}}_{h}=\mathcal{P}^{i}_{h}(\mathbf{u}_{h})$. Such ghost layers are depicted in the Fig. 4. Then, for $i=L$, we have

[TABLE]

We compare equation (18) and the projected to $\Omega_{i}$ pressure equation (16):

[TABLE]

Since $\mathbf{u}^{\Gamma}_{h,L}=\mathcal{P}^{L}_{h}\textbf{u}^{\Gamma}_{h}$, $A^{L}_{1}=\mathcal{P}^{L}_{h}A_{1}$ , we have the following

[TABLE]

$A^{L}_{2}$ is non-singular and diagonal matrix, since $\mathbf{K}$ is SPD.

[TABLE]

Similarly, we can obtain

[TABLE]

In non-linear problems including slightly compressible flow or multiphase flow in heterogeneous porous media, this approach could be applied analogously by taking into account ghost layers values arising from $e_{i}$ in each Newton iteration. So during Block Jacobi iteration variables $p^{e,k-1}_{L}$, $p^{e,k-1}_{R}$ is computed by utilizing given $p^{k-1}_{L}$, $p^{k-1}_{R}$ and then solve decoupled subdomain problems with Dirichlet boundary conditions $p^{e,k-1}_{i}$, $i=L,R$ to find $u^{k},p^{k}$.

3 Methods

We use the postprocessing procedure associated to pressure and velocity. We first apply locally postprocessing algorithm for given pressure $p_{h}$ and velocity $\mathbf{u}_{h}$ which was previously proposed in [9] and then Oswald interpolation operator [12, 13, 14, 15] to have better pressure values. At the interface, we use two-point flux computation method in order to have better approximation of pressure. As a result, the Enhanced Velocity scheme solution of flux can be improved by using a post-processed pressure. The key idea is illustrated in Fig. 5 for resulting approximation of EV scheme that is shown in Fig. 3.

The velocity at the edge or face is computed by using pressure values between subdomains $\Omega_{i}$ and $\Omega_{j}$. To be specific, $p_{h}\in\Omega^{*}$ is required in the original velocity for constructing in Enhanced Velocity MFEM. However, the post-processed pressure leads to the improved velocity and the visual representation is in Fig. 5. In case of multiscale setting, it is important to be able to approximate better pressure values nearby the interface. The recovery of velocity computation requires three steps

Compute locally $\tilde{p}_{h}$ from given $\left(p_{h},\mathbf{u}_{h}\right)$ 2. 2.

Obtain $s_{h}$ by using Oswald operator 3. 3.

Compute the velocity at the interface using the two-point flux scheme for $s_{h}$

We describe construction of $\tilde{p}_{h}$ and then $s_{h}$ below.

Construction of $\tilde{p}_{h}$.

In the Enhance Velocity setting, we may identify $\widehat{\mathbf{V}}_{h}$ be spaces omitting interface constraints $\mathbf{V}^{\Gamma}$, so $\widehat{\mathbf{V}}_{h,i}:=\bigoplus_{i=1}^{n}\mathbf{V}_{h,i}(T)$ and then $\widehat{\mathbf{V}}_{h}:=\bigoplus_{i=1}^{n}\widehat{\mathbf{V}}_{h,i}$. Let $\mathbf{u}_{h}$, $p_{h}$ be the solution of equations (8)-(9). Initially, Lagrange multipliers can be computed in each element. In other words, we define $\lambda_{h,T}\in\Lambda_{h}$, which is piecewise constant polynomials at edge or face,

[TABLE]

where the element $T\in\mathcal{T}_{h}$ and its side $e$. We employ the $L^{2}$ projected velocity from the interface, which has a finer enhanced velocity approximation, to the edge or face of subdomain element and the formulation is provided in Subsection 2.2. We denote polynomial space $\widetilde{W}_{h}$ in the following manner

[TABLE]

where $\mathbb{Q}_{m}$ is standard polynomial space that is defined in [9, 12, 16]. We next set the post-processed $\tilde{p}_{h}$ which is proposed in [9] and the construction is performed with the following properties, for each $T\in\mathcal{T}_{h}$

[TABLE]

Construction of $s_{h}$.

We propose to construct $s_{h}$ in each subdomain $\Omega_{i}$ that has the conforming mesh in order to be a computational efficient. Construction of $s_{h}$ involves the averaging operator $\mathcal{I}_{\rm av}:\mathbb{Q}_{k}(\mathcal{T}_{h})\rightarrow\mathbb{Q}_{k}(\mathcal{T}_{h})\cap H_{0}^{1}(\Omega_{i})$. For definition of $\mathbb{Q}_{m}$ we refer reader to [16]. The operator is called Oswald operator and appeared in [12, 13, 14, 15] and the analysis can be found in [17, 14]. It is interesting to note that the mapping of the gradient of pressure through Oswald operator also considered in [18]. For given $\varphi_{h}\in\mathbb{Q}_{m}(\mathcal{T}_{h})$, we regard the values of $\mathcal{I}_{\rm av}(\varphi_{h})$ as being defined at a Lagrange node $V\in\Omega$ by averaging $\varphi_{h}$ values associated this node,

[TABLE]

where $|A|$ is cardinality of sets $A$ and $\mathcal{T}_{h}$ is all collection of $T\in\mathcal{T}_{h}$ for fixed $V$. One can see that $\mathcal{T}_{h}(V)=\varphi(V)$ at those nodes that are inside of given $T\in\mathcal{T}_{h}$. We set the value of $\mathcal{I}_{\rm av}(\varphi_{h})$ is zero at boundary nodes. Now in our setting we define recovered pressure $s_{h}$ for the locally post-processed $\tilde{p}_{h}$ as follows.

[TABLE]

3.1 Implementation steps

For simplicity, we provide key steps of numerical implementation of post-processed pressure in two dimensional case. However, it can be extended for general cases. Based on piecewise pressure and velocity from the lowest order Raviart-Thomas spaces over rectangles our aim to reconstruct smoother pressure $s_{h}$. For given element $T\in\mathcal{T}_{h}(\Omega_{i})$, the main steps are

Evaluate $\lambda_{h,T}$ at edge $e_{j}$, $j=1,..4$ based on $(\mathbf{u}_{h},p_{h})$, 2. 2.

Compute $\tilde{p}_{h}$ from known $\lambda_{h,T}$, and $p_{h}$ by using (23), 3. 3.

Based on $\tilde{p}_{h}$ compute $s_{h}$ equation (27) at Lagrange nodes in $\Omega_{i}$.

Step 1 is standard computation of Lagrange multiplier for each element. In step 2, we are relying on higher order polynomial, in our case, it is Span$\{1,x,y,x^{2},y^{2}\}$. It is sufficient to store coefficients of polynomials. In step 3, we use

Span$\{1,x,y,x^{2},y^{2},xy,x^{2}y,xy^{2},x^{2}y^{2}\}$ and 9 Lagrange nodes of rectangle elements that are four rectangle nodes, four midpoints at edge and center of rectangle. This case each node requires to find neighboring elements values to compute coefficients of $s_{h}$.

4 Numerical Examples

In this section, numerical results are presented to demonstrate challenging problems of velocity approximation at the interface of non-matching multiblock grids. We have conducted tests for several examples and we concentrate our attention on the interface error for homogeneous and heterogeneous permeability coefficients, respectively. We set same domain $\Omega=\left(0,1\right)\times\left(0,1\right)$ for all tests and for some the ratio is $H/h=4$. Initial subdomains grids $\mathcal{T}_{h}$ are chosen in way that has a checkerboard pattern for subdomains. Example of such discretization is shown in Fig. 6. The discrete $L^{2}$ velocity error $e_{\mathbf{u}_{h},\Gamma}$ is based on the values of the normal component at the midpoint of the edges and is normalized by the analytical solution.

Numerical test 1

First example tests for uniform permeability, so $\mathbf{K}=\mathbf{I}$. We report the velocity error and the improved velocity error. We compute the source term and boundary conditions according to the analytical solution, which is taken as follows

[TABLE]

We set the ratio $H/h=4$ for the result that is shown in below table.

We list the error of velocity and recovered velocity in Tabel 1. The results shows the convergence rate improvement for recovered velocity approximation compare to the provided velocity approximation at the interface. We note that the improvement of convergence rate is significant from order $O(h^{0.95})$ to $O(h^{1.5})$.

Numerical test 2

We consider the a diagonal oscillating tensor coefficient as follows.

[TABLE]

We impose the source term $f$ and Dirichlet boundary condition according to the analytical solution

[TABLE]

We set the ratio $H/h=4$ for the result that is shown in Table 2.

From Table 2, we see a significant increase on the convergence rate for recovered velocity while the convergence rate of provided velocity stays $\mathcal{O}(h^{1.0})$. We observe that the numerical method is an effective way to improve velocity at the interface between subdomains.

5 Conclusion

The present study of velocity in Enhanced Velocity Mixed Finite Element Method was designed to investigate the effect of the post-processed pressure on velocity in the interface of subdomains. In this paper, the focus of attention is on the incompressible Darcy flow in the non-matching multiblock grid setting. Multiple numerical results demonstrate that the interface velocity approximation can be improved with using the post-processed pressure. These findings can contribute in several ways to our approximation of velocity and provide a good construction of velocity for a posteriori error analysis such as the recovery-based estimate.

6 Acknowledgments

First author would like to thank Drs. I. Yotov and T. Arbogast for discussions on formulation of the different view of EVMFEM. This research is supported by Faculty Development Grant, Nazarbayev University.