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

Accurate prediction of the macroscopic flow parameters needed to describe flow in porous media relies on a good knowledge of flow field distribution at a much smaller scale—in the pore spaces. The extent of the inertial effect in the pore spaces cannot be underestimated yet is often ignored in large-scale simulations of fluid flow. We present a multiscale method for solving Oseen’s approximation of incompressible flow in the pore spaces amid non-periodic grain patterns. The method is based on the multiscale finite element method [MsFEM Hou and Wu in J Comput Phys 134:169–189, 1997)] and is built in the vein of Crouzeix and Raviart elements (Crouzeix and Raviart in Math Model Numer Anal 7:33–75, 1973). Simulations of inertial flow in highly non-periodic settings are conducted and presented. Convergence studies in terms of numerical errors relative to the reference solution are given to demonstrate the accuracy of our method. The weakly enforced continuity across coarse element edges is shown to maintain accurate solutions in the vicinity of the grains without the need for any oversampling methods. The penalisation method is employed to allow a complicated grain pattern to be modelled using a simple Cartesian mesh. This work is a stepping stone towards solving the more complicated Navier–Stokes equations with a nonlinear inertial term.

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 ≈ 0.001, three orders of magnitude smaller than what suggested in the literature (Re ≈ 1) based on studies of homogeneous media, such as bead packs, which are poor representations of the heterogeneous reservoir rocks of practical interest.
The flow of viscous fluid is adequately governed by Navier-Stokes equation where u and p are velocity and pressure, respectively, μ is the dynamic viscosity, and ρ is the density. Typically the nonlinear, inertial term ρ(u · ∇u) would lead to nonlinear algebraic equations which require iterations of the linearised equation to solve. Our ultimate goal is to study the application of MsFEM on solving (1) in both fine and coarse-scales. However, we find it reasonable to first apply our method on the Oseen's approximation of Navier-Stokes equation where the inertial term is substituted by Oseen term ρ(U · ∇u), where U is a known velocity vector. Proudman and Pearson (1957) point out that Oseen's solution provides a uniformly valid zeroth order approximation to the Navier-Stokes equation at finite Reynolds number. Oseen system includes important inertial factors, while retaining the solution linearity, and thus simplicity (Oseen equation is analytically tractable). As a result, Oseen system has been utilised as a testbed for validating different computational approaches (e.g. stabilisation of finite element discretisations (Abraham et al. 2004), weak Galerkin finite element method (Liu et al. 2016), or artificial boundary conditions on truncated computational domains (Chun-xiong 2002) for it represents an alternative model to the Stokes' in which the nonlinear term is eliminated entirely.
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 singlephase 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 twophase 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.

Problem Formulation
We define Ω, a two-dimensional domain consisting of a grain domain Ω grain , and a pore domain Ω pore perforated by grains, see Fig. 1. Then, let ε denote the diameter of the smallest grain. The steady-state Oseen's problem is to find velocity u which is the solution to: The boundary condition of Eq.
(2) is given by: where f is a source function, and g is a function fixed at the boundary ∂Ω. In this paper, we consider only a no-slip condition on the grain boundaries: u = 0, on ∂Ω grain .

Fig. 2
An illustration of the discretised domain τ H , and ω E which is the support space for E . In all applications in this paper, ε/ h ≥ 5 is maintained

Discretisation
We discretise Ω into a two-dimensional homogeneous Cartesian coarse mesh τ H (see Fig. 2

Crouzeix-Raviart Functional Spaces
The functional spaces for velocity V H , and for pressure M H are given below: The key here is to maintain the continuity of (only) the average of v across an edge E: ] is the jump of the value v across E. We wish to retain the advantage of our approach which has been successfully applied on Advection-Diffusion, and Stokes problems, namely: the weak imposing of continuity across element boundaries allows adaptive boundary conditions which relaxes the sensitivity of our method to random arrangements of grains, without the need of applying the more cumbersome oversampling methods.

Construction of a Crouzeix-Raviart Basis
We also define supp( E,i ) ⊂ ω E , the ensemble of two quadrangles T k , k = 1, 2 in τ H which share an edge E. Hence find, in each of these quadrangles, E,i and π E,i , which are the solutions to: is the set of all the edges of the quadrangle T k . To solve Eq. (6), we use Q1-Q1 finite element spaces in which both velocity and pressure degrees of freedom are defined on the same set of grid points. This arrangement is chosen due to the ease of programming and the computational efficiency. A stabilisation method is however necessary when this approach is considered. In a homogeneous Cartesian coordinate with fine element width of h, a stable solution can be achieved by perturbing the condition ∇ · E,i = 0 with a pressure Laplacian term (see Brezzi and Fortin 1991). In a weak form, Eq. (6) reduces to finding E,i ∈ H 1 (T k ∩ Ω pore ), π E,i ∈ L 2 0 (T k ∩ Ω pore ), and the Lagrange multipliers λ F , Here, v h and q h occupy the same finite element spaces as E,i and π E,i , respectively. θ is the stabilisation parameter which we set as 0.01 for all our simulations (see Brezzi and Pitkäranta 1984).

Crouzeix-Raviart MsFEM Coarse-Scale Solution
Here, we describe the coarse-scale solution of Eq. (2). By discretising p into p H and u into u H , we can solve Eq.
(2) in τ H and rewrite it in a weak form as: The solution to problem (2) can then be approximated as linear combination of multiscale basis functions E, j = E e j , j = 1, 2 with {e 1 , e 2 } being the canonical basis of R 2 , such that (Fig. 3):

Penalisation Method
Often Ω 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 Ω pore , we solve ( Fig. 4): in which In our simulations, the chosen fine-scale element width h always satisfies ε/ h ≥ 5. The penalisation coefficient σ κ then forces the solution u to vanish inside the obstacles. Other variants of penalisation methods are studied in Angot et al. (1999).

Reference Solution
We use a Q1-Q1 finite element method to compute the reference solutions, as we do for computing , in the global fine mesh τ h . In a weak form, the solution to Eq. (2) are u h and p h such that (Fig. 5) where v h , and p h occupy Q1-Q1 finite element spaces. We use the stabilisation parameter θ = 0.01 throughout this paper (Elman et al. 2005). Note that other stable or stabilised elements can be used to provide the reference solution.

Boundary Condition
Note that the boundary condition u = g on ∂Ω is not included in Eq. (4). This is possible since we approximate the boundary condition in Eq. (6) only in a weak sense, i.e Together with Eq. (12), we apply on the boundary ∂Ω (Fig. 6): This approach has been applied successfully in Muljadi et al. (2015) for the Stokes equation and is a modification with respect to previous work (Bris et al. 2013(Bris et al. , 2014 where the boundary conditions were strongly incorporated in the definition of V H . Our approach therefore gives more flexibility when implementing non zero g, i.e every basis functions E,i including those on boundary ∂Ω can be computed in the same fashion-according to Eq. (6).

Numerical Results
We consider a channel domain Ω = [−2 ≤ x ≤ 2, −1 ≤ y ≤ 1] containing a porous medium spanning from x = −1 to x = 1. We then assign ρ = 1, μ = 0.001, and f = 0. At the inlet, the theoretical incompressible Poiseuille solution (parabolic velocity profile) is applied for all cases, i.e u = 1 − y 2 , 0 on x = −2, whereas the Neumann boundary condition ∂u/∂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). 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 ε = 0.06. Pattern (b) consists of 900 grains with ε = 0.0067 (Fig. 10).
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 × 1280 quadrangles. This ensures the ratio ε/ h ≥ 5.

Poiseuille Flow
We apply a vector field U = (0.002, 0) corresponding to Re ≈ 4 where Re = ρε|U| μ (in the absence of grains, we assume that ε 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).

Pattern (a)
Here, we test our method in solving flow pass a porous body with a random pattern of grains. We consider pattern (a) where each grain is a rectangle with width ε = 0.06, Re ≈ 0.12. In Figs. 4 and 5, the contours of velocity components u x and u y computed on a number coarse meshes are given. The results are compared to the reference solution. The contours computed using Crouzeix-Raviart MsFEM at 256×128 are identical to the reference solution; however, even at 32 × 16, the flow features already resemble those of the reference solution, and at 64 × 32 without any appreciable difference. Norms of error relative to the reference solution in L 1 , L 2 , and H 1 spaces are given in Table 2 showing a converging behaviour.

Pattern (b)
Here, we simulate flow pass pattern (b) which contains finer (ε = 0.0067) and much more grains. We use the same kinds of coarse meshes, and vector field U as in the previous tests, which corresponds to Re ≈ 0.013. This is a more challenging test than the previous ones due to much larger H/ε ratios-ranging from 18.74 to 2.34. This means each coarse elements in the vicinity of the grains has higher chances of being occupied by more than one grain, and therefore suffers more oscillations.
In Figs. 6 and 7 the contours of u x and u y are given. As in the previous tests, the results are compared to the reference solution. The contours of u x and u y computed using Crouzeix-

Setting with Dense Grains
Here, we simulate flow pass a porous body with dense grains. For the same grain size with that in pattern (b), ε = 0.0067, we now lay 3600 grains randomly. The results computed with Crouzeix-Raviart MsFEM and the reference solution are given in Fig. 10. At 100 times coarser mesh, the result obtained using Crouzeix-Raviart MsFEM has no noticable difference compared with the reference solution.

Heterogeneous Velocity Field U
To further test the feasibility of our method, we apply a heterogeneous velocity field In Fig. 8, the contours of u x through pattern (b) computed using Crouzeix-Raviart MsFEM on 128 × 64 coarse elements; and the reference solution are compared. The flow features are noticeably different than the previous results especially in regions away from the porous body. Indeed at the farthest from grains, the flow experiences Re ≈ 80 where the inertial effect manifests the most. In Table 4, the convergence study is given for a range of coarse meshes. Again we can see that our method gives good qualitative and quantitative agreement with the reference solution. In Fig. 9, the contours of the magnitude of velocity computed using Crouzeix-Raviart MsFEM on 128 × 64 coarse elements are displayed along with their streamlines and compared with the reference solution. We notice that the Crouzeix-Raviart MsFEM gives an excellent agreement in terms of flow pattern with the reference solution.
The effects of inertia away from the porous bodies can be instantly noticed when comparing the flow structures in Figs. 6 and 8 (computed on the same domain and grains structure, but with different U). To help us understand the extent of these effects in a close vicinity of grains, we compare the fine-scale basis functions with, and without Oseen term. In Figs. 11, we plot the multiscale basis function computed in the framed, and zoomed patch that lies in a computational domain with pattern (b). Figure 11a shows the basis functions E,1 computed with U = 0-the Oseen problem therefore reduces to a Stokes problem. Figure 11b shows the basis functions E,1 computed with U according to Eq. (18). First, we note the effects of Oseen term on the shapes of the basis functions, which our method successfully captures. Second, in both cases the Crouzeix-Raviart basis functions succesfully accommodate the coincidence between a grain and an edge E and maintain E E,i = e i .

Concluding Remarks
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 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 f = 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.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.