Multiscale Method for Oseen Problem in Porous Media with Non-periodic Grain Patterns

Modelling of flow through porous bodies is a topic of high importance in various fields of engineering, chemical, biological or geological applications. One of the most significant challenges persisting in virtually all these areas is the disparity between the spatial scales at which flow and transport can be understood; and the scales at which practical model predictions are needed (Scheibe et al. 2015). This disparity in scales forces a trade-off between building models which suffice for practical application, and models that solve the problem ab-initio but which may not be able to cope with large-scale problems adequately. To place this in context, in geological media, X-ray techniques now allow three-dimensional images to be acquired routinely (Blunt et al. 2013). The pore spaces of these rocks are typically of order microns across. However, for practical applications in oil recovery, carbon dioxide storage and contaminant transport, flow over 100s m to km needs to be predicted. This enormous range of scales precludes the use of a direct method that resolves pore-scale flow while determining reservoir-scale behaviour. Instead, techniques that can approximate the flow over distances much larger than the pore scale are needed. A number of multiscale simulation paradigms have been developed to bridge first-principles and empirical methods, and provide a link between micro, and macroscale models. An additional problem is that in many applications, such as flow in fractured rock and near-well bore flows, the nonlinear, or inertial term in the Navier–Stokes equation are significant. This means that at the large-scale, the application of the linear Darcy-law for flow is inaccurate.

Our choice of a particular multiscale method is based on the following. First, we consider direct Navier–Stokes simulation on the pore geometry as the holy-grail of microscale simulation for it is considered to be the most complex, and highly resolved spatially (although Navier–Stokes itself can be seen as an upscaled representation of molecular-scale interactions, with effective parameters such as viscosity and density). Such a microscale model strikes a balance between the appropriate level of complexity with current technological advances. For example, recent developments in both computational algorithms, and increases in computer power, coupled with the availability of pore-scale images, have enabled the routine prediction of permeability, with direct Navier–Stokes calculations on samples containing up to a billion voxels (Blunt et al. 2013; Mostaghimi et al. 2012). Second, we are interested in a method capable of resolving the microscale model directly over the domain of interest without losing any degrees of freedom—which rules out other multiscale methods which borrow their philosophy from homogenisation theory (e.g formal upscaling with closure approximation). Several multiresolution solvers are designed for this purpose, i.e to provide computationally efficient ways of obtaining a complete solution on the fine grid, for example: multigrid solvers and preconditioners (Wesseling 1992), multiscale finite element methods (MsFEM; Hou and Wu 1997; Aarnes et al. 2005; Jenny et al. 2003), and multiscale mimetic methods (Lipnikov et al. 2011). We choose to develop an adaptation of MsFEM dedicated for solving flow in a pore domain left void by non-periodic grain patterns, which is a representation of all natural pore structures.

The challenge in applying MsFEM in a non-periodic setting is to avoid an intersection between a coarse element boundary and a grain. On the other hand, the overall performance of MsFEM rely on the accuracy of the multiscale basis function which is very sensitive the treatment of subgrid boundary condition. The application of oversampling methods (Efendiev et al. 2013; Chu et al. 2008; Henning and Peterseim 2013) was intended to circumvent this problem by broadening the domain in which basis functions are sampled. While the methods perform satisfactorily in the context of perforated media (Bris et al. 2014; Chung et al. 2015), nevertheless it necessitates an ad hoc parameterisation and results in a larger computational problem. Another alternative is to adopt a nonconforming finite element method and impose only a weak continuity between coarse element boundaries and therefore allowing the coarse element boundaries to adapt to random patterns of grains. In our previous works, the nonconforming Crouzeix–Raviart element has been adopted successfully for solving advection–diffusion and Stokes equations (Bris et al. 2013; Degond et al. 2015; Muljadi et al. 2015).

Creeping or Stokes flow is often assumed in porous media. This ceases to apply, as mentioned above, for example, near propped fractures, or boreholes in reservoirs where inertial forces becomes dominant. Even in the absence of fractures, Muljadi et al. (2016b) studied the non-Darcy flow behaviour in porous media with different pore heterogeneities and found that the cessation of the Darcy relationship in Estaillades limestone already takes place at Re \(\approx 0.001\), three orders of magnitude smaller than what suggested in the literature (Re \(\approx 1\)) based on studies of homogeneous media, such as bead packs, which are poor representations of the heterogeneous reservoir rocks of practical interest.

In the context of multiscale methods for flow in porous media, Abdulle and Bud (2015) presented a multiscale method based on the coupling of an effective Darcy equation on a macroscopic mesh with unknown permeabilities recovered from micro finite element calculations for Stokes problems. Alyaev et al. (2014) presented a heterogeneous multiscale method that utilises fine-scale information directly to solve problems for general single-phase flow on the Darcy scale. Their work focuses on the nonlinear flow regime, but does not employ Oseen’s approximation. Bonfigli and Jenny (2009) presented a multi-scale pressure solver for the incompressible Navier–Stokes equations with immersed boundaries. For two-phase flows in porous media, Tomin and Lunati (2013) presented hybrid Multiscale Finite Volume (MsFV) method, which couples pore-scale and Darcy descriptions of two-phase flow in porous media. The finite volume discretisation ensures conservation of mass at the coarse-scale, and a conservative pore-scale velocity field. The above-mentioned methods apply Darcy modelling at the coarse-scale whereas in this paper, Oseen equations are solved at both fine, and coarse-scales.

To avoid having to work with complicated boundary-fitted or even unstructured meshes, we employ the penalisation method (Angot et al. 1999) when modelling non-periodic grain patterns. Here, we simply force the solution to vanish within the grain boundaries. Consequentially, this approach allows the modelling of a complicated grains pattern on a simple Cartesian mesh.

This paper is organised as follows. The formulation of the problem is given in Sect. 2. The construction of Crouzeix–Raviart elements is presented in Sect. 3. In Sect. 4, the application of penalisation method on our problem is described. Then, the description of the computation of the reference solutions is given in Sect. 5 followed by some remarks on the treatment of the boundary condition in Sect. 6. Numerical tests are presented and the results discussed in Sect. 7 followed by some concluding remarks.

Open image in new window

Here, we explain the application of our method starting from the definition of the coarse and fine meshes. We then introduce the functional spaces for our multiscale basis functions and describe the construction of these functions within each coarse elements.

3.4 Crouzeix–Raviart MsFEM Coarse-Scale Solution

Here, we describe the coarse-scale solution of Eq. (2). By discretising p into \(p_H\) and \(\mathbf {u}\) into \(\mathbf {u}_H\), we can solve Eq. (2) in \(\tau _H\) and rewrite it in a weak form as:

$$\begin{aligned} a(\mathbf {u}_H,\mathbf {v}_H) + c(\mathbf {u}_H,\mathbf {v}_H) + b(\mathbf {v}_H,p_H)= & {} (\mathbf {v}_H,\mathbf {f}) ,\quad \forall \mathbf {v}_H \in V_H\nonumber \\ b(\mathbf {u}_H,q_H)= & {} \quad 0,\quad \forall q_H \in M_H \end{aligned}$$

(8)

where

$$\begin{aligned} a(\mathbf {u},\mathbf {v})= & {} \int _{\varOmega _\text {pore}} \mu \nabla \mathbf {u}: \nabla \mathbf {v}\, d{\varOmega _\text {pore}} \end{aligned}$$

(9)

$$\begin{aligned} c(\mathbf {u},\mathbf {v}) = c(\mathbf {U};\mathbf {u},\mathbf {v})= & {} \int _{\varOmega _\text {pore}}\rho (\mathbf {U}\cdot \nabla \mathbf {u})\cdot \mathbf {v}\, d{\varOmega _\text {pore}} \end{aligned}$$

(10)

$$\begin{aligned} b(\mathbf {v},q)= & {} -\int _{\varOmega _\text {pore}} q \nabla \cdot \mathbf {v}\, d{\varOmega _\text {pore}}. \end{aligned}$$

(11)

The solution to problem (2) can then be approximated as linear combination of multiscale basis functions \(\mathbf {\Phi }_{E,j}={\Phi }_{E}\mathbf {e}_j\), \(j = 1,2\) with \(\{\mathbf {e}_{1},\mathbf {e}_{2}\}\) being the canonical basis of \(\mathcal {R}^{2}\), such that (Fig. 3):

$$\begin{aligned} \mathbf {u}_H(x,y)=\sum _{E\in \mathcal {E}_H,\,j=1,2}{u}_{E,j}\mathbf {\Phi }_{E,j}(x,y). \end{aligned}$$

(12)

Open image in new window

Often \(\varOmega _\text {pore}\) is a complicated pore structure in which solving Eq. (2) may require a complicated boundary-fitted, or even an unstructured mesh. In order to confine our computations in a homogeneous Cartesian mesh, we employ the penalisation method (Angot et al. 1999). Henceforth, instead of solving Eq. (2) directly in \(\varOmega _\text {pore}\), we solve (Fig. 4):

$$\begin{aligned} \mu _\kappa \nabla ^2 \mathbf {u}+\rho _\kappa (\mathbf {U}\cdot \nabla \mathbf {u}) + \sigma _\kappa \mathbf {u} - \nabla p= & {} \mathbf {f}_\kappa \text {, in } \varOmega \nonumber \\ \nabla \cdot \mathbf {u}= & {} 0 \text {, in }\varOmega \nonumber \\ \mathbf {u}= & {} \mathbf {g} \text {, on } \partial \varOmega \end{aligned}$$

(13)

in which

$$\begin{aligned} \mu _\kappa ;\rho _\kappa = \left\{ \begin{array}{cc} \frac{1}{h} &{} \text{ in } \varOmega _\text {grain} \\ \mu ;\rho &{} \text{ in } \varOmega _\text {pore} \end{array} \right. , \sigma _\kappa = \left\{ \begin{array}{cc} \frac{1}{h^3} &{} \text{ in } \varOmega _\text {grain} \\ 0 &{} \text{ in } \varOmega _\text {pore} \end{array} \right. , \mathbf {f}_\kappa = \left\{ \begin{array}{cc} 0 &{} \text{ in } \varOmega _\text {grain} \\ \mathbf {f} &{} \text{ in } \varOmega _\text {pore} \end{array} \right. . \end{aligned}$$

(14)

Open image in new window

In our simulations, the chosen fine-scale element width h always satisfies \(\varepsilon /h\ge 5\). The penalisation coefficient \(\sigma _\kappa \) then forces the solution \(\mathbf {u}\) to vanish inside the obstacles. Other variants of penalisation methods are studied in Angot et al. (1999).

We consider a channel domain \(\varOmega = [-2\le x \le 2, -1 \le y \le 1]\) containing a porous medium spanning from \(x=-1\) to \(x=1\). We then assign \(\rho =1\), \(\mu = 0.001\), and \(\mathbf {f}=0\). At the inlet, the theoretical incompressible Poiseuille solution (parabolic velocity profile) is applied for all cases, i.e \(\mathbf {u}=\left( 1-y^2,0\right) \) on \(x=-2\), whereas the Neumann boundary condition \(\partial u/\partial n = 0\) is assumed at the outlet, \(x=2\). No–slip boundary conditions at the top and bottom walls are applied (Figs. 7, 8, 9).

Open image in new window

Open image in new window

Open image in new window

First we apply our method on simple Poiseuille flow without any porous bodies. Then, two grain patterns are included as shown in see Fig. 3. From here on we refer to them as pattern (a) depicted in Fig. 3a, and pattern (b) depicted in Fig. 3b. Pattern (a) consists of 100 randomly placed grains, each with width \(\varepsilon =0.06\). Pattern (b) consists of 900 grains with \(\varepsilon =0.0067\) (Fig. 10).

Open image in new window

We compare all our results to the reference solutions. When computing the reference solutions, we employ Q1-Q1 finite element method on a fine mesh consisting of \(2560\times 1280\) quadrangles. This ensures the ratio \(\varepsilon /h\ge 5\).

7.1 Poiseuille Flow

We apply a vector field \(\mathbf {U}=\left( 0.002,0\right) \) corresponding to \({Re}\approx 4\) where \({Re}=\frac{\rho \varepsilon |\mathbf {U}|}{\mu }\) (in the absence of grains, we assume that \(\varepsilon \) is the channel diameter). In Table 1, the norms of error of the Crouzeix–Raviart MsFEM solutions relative to the reference solution on a number of coarse meshes are given, showing a convincingly converging trend (Fig. 11).

Open image in new window

The Crouzeix–Raviart MsFEM has been developed and tested for solving Oseen’s approximation for incompressible flow around solid bodies. The method performs very well in the presence of non-periodic grain formations. The weakly enforced continuity across coarse element edges ensures accurate basis function solutions without any oversampling methods. The basis functions are shown to successfully capture the effects of homogeneous and inhomogeneous vector field \(\mathbf {U}\). The penalisation method is seamlessly incorporated into our method allowing an extensive utilisation of simple Cartesian mesh.

This method is developed as a stepping stone towards solving then more complicated Navier–Stokes equation. Although only two-dimensional cases are considered, the extension of this work on three-dimensions is straightforward. Similarly, the method can be applied for inhomogeneous Oseen’s problem with \(\mathbf {f}\ne 0\). The reconstruction of fine-scale pressure is not the focus of the current work although it is possible (see Muljadi et al. 2015). The calculations of MsFEM basis functions within a coarse element are done independent of the neighbouring elements which makes it suitable for the application of parallel programming. The fine-scale discretisation of our method can be reformulated in a conservative fashion, which will improve accuracy. For example, following the proposition in Multiscale Finite Volume Method (MsFVM; Tomin and Lunati 2013; Jenny et al. 2003), or mass conservative Generalised Multiscale Finite Element Method (GMsFEM; Presho and Galvis 2016); this is one plausible future endeavour.

This work will lay the groundwork for applying various types of basis function enrichment. For example, bubble functions have been efficiently, and effectively integrated into the Crouzeix–Raviart MsFEM formulation (Degond et al. 2015) to solve advection–diffusion problems in a setting with very dense grains. A similar idea can be applied to the current method, and is the subject of our future work.

For practical applications, the method is a promising development towards simulations capable of handling a wide range of spatial scales, while accommodating nonlinear effects.