Hybrid Multiscale Finite Volume method for two-phase flow in porous media

Journal of Computational Physics 250 (2013) 293–307

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

Hybrid Multiscale Finite Volume method for two-phase flow in porous media Pavel Tomin ⇑, Ivan Lunati Faculty of Geosciences and Environment, University of Lausanne, Switzerland

a r t i c l e

i n f o

Article history: Received 13 February 2013 Received in revised form 17 April 2013 Accepted 3 May 2013 Available online 24 May 2013 Keywords: Multi-physics modeling Hybrid multiscale methods Multiscale Finite Volume method Volume of Fluid method

a b s t r a c t We present a novel hybrid (or multiphysics) algorithm, which couples pore-scale and Darcy descriptions of two-phase flow in porous media. The flow at the pore-scale is described by the Navier–Stokes equations, and the Volume of Fluid (VOF) method is used to model the evolution of the fluid–fluid interface. An extension of the Multiscale Finite Volume (MsFV) method is employed to construct the Darcy-scale problem. First, a set of local interpolators for pressure and velocity is constructed by solving the Navier–Stokes equations; then, a coarse mass-conservation problem is constructed by averaging the pore-scale velocity over the cells of a coarse grid, which act as control volumes; finally, a conservative pore-scale velocity field is reconstructed and used to advect the fluid–fluid interface. The method relies on the localization assumptions used to compute the interpolators (which are quite straightforward extensions of the standard MsFV) and on the postulate that the coarse-scale fluxes are proportional to the coarse-pressure differences. By numerical simulations of two-phase problems, we demonstrate that these assumptions provide hybrid solutions that are in good agreement with reference pore-scale solutions and are able to model the transition from stable to unstable flow regimes. Our hybrid method can naturally take advantage of several adaptive strategies and allows considering pore-scale fluxes only in some regions, while Darcy fluxes are used in the rest of the domain. Moreover, since the method relies on the assumption that the relationship between coarse-scale fluxes and pressure differences is local, it can be used as a numerical tool to investigate the limits of validity of Darcy’s law and to understand the link between pore-scale quantities and their corresponding Darcy-scale variables. Ó 2013 Elsevier Inc. All rights reserved.

1. Introduction Flow through porous media is usually described by Darcy’s law, which relates the macroscopic velocity (volumetric flux density integrated over a large number of pores) to the pressure gradient. This empirical law [8] is used in place of a macroscopic-momentum conservation equation and relies on several hypothesis: slow flow and negligible inertia; short relaxation time; and spatial-scale separation. Although Darcy’s law has proven appropriate to many practical problems, it breaks down for flow regimes in which the interaction between nonlinearity of the flow and disorder of the porous medium lead to the emergence of long correlation structures. In these cases, the macroscopic solution becomes dependent on the pore-scale details, which cannot be neglected. A way to deal with this situation, is to construct a hybrid (or multiphysics) algorithm that couples a coarse-scale Darcy description with a local pore-scale description. This has been recently illustrated for single-phase reactive transport [3] and a ⇑ Corresponding author. Tel.: +41 216924412. E-mail addresses: [email protected] (P. Tomin), [email protected] (I. Lunati). 0021-9991/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.05.019

294

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307

hybrid algorithm was developed to divide the domain into regions where either a pore-scale description or a Darcy description is used [4]. In this approach, no global coarse-scale problem was introduced and the consistency between the different regions was achieved simply by an iterative algorithm to obtain continuity of fluxes across the boundary. A similar iterative coupling algorithm has been proposed [2] for single phase flow problems under the assumption of periodic isotropic porous media in application to non-linear flow phenomena. In this paper, we follow a different strategy and devise a general framework that allows constructing a global coarse-scale problem that couples a set of local numerical pore-scale solutions. We consider the simultaneous flow of two phases, in which the presence of surface forces and the different properties of the two fluids yield a two-way coupling between the flow and the transport problem. At the pore-scale, the flow is described by the Navier–Stokes equation of motion supplemented by the divergence-free condition for the velocity. The evolution of the fluid–fluid interface is modeled with the Volume of Fluid (VOF) method [12,17,27], which has been extensively examined to be able to deal with large viscosity contrast and complex geometries [10], as well as contact angle dynamics [1,23]. Starting from these governing equations at the pore-scale, we employ the Multiscale Finite Volume (MsFV) method [13,22,24] as a general framework to construct the coarse-scale mass conservation equations. A Darcy description of the flow at the coarse scale is obtained by postulating that the coarse pressures are the only macroscopic degrees of freedom, which allows dropping the coarse momentum-conservation equation. A set of local interpolators (basis and correction functions [19,21,15,22]) are constructed for pressure and velocity and used to build a Darcy-scale mass-balance equation. These interpolators are localized solutions of the Navier–Stokes equations and surface tension effects are accounted by the correction functions. Notice that the approach followed here is very different from previous applications of the MsFV method to solve the Navier–Stokes equations [6], where the MsFV procedure was simply proposed as an efficient algorithm for solving the pressure Poisson’s equation arising in numerical methods in a predictor–corrector solver. The paper is organized as follows. The pore-scale governing equations and the VOF method are described in Section 2, whereas the details of the hybrid algorithm are presented in Section 3. The results of numerical simulation performed for different flow regimes are presented and discussed in Section 4. Here, we use direct pore-scale simulation results [10] as references to check the quality of hybrid solutions. Finally, concluding remarks are given in Section 5.

2. Pore-scale governing equations and the VOF method We consider the flow of two incompressible, immiscible fluids in a complex pore geometry. For simplicity (but without loss of generality) the porous medium is idealized as a horizontal two-dimensional domain filled with impermeable rigid circular obstacles (Fig. 1), and gravity is therefore neglected.

Fig. 1. Pore geometry.

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307

295

We use a whole-domain formulation, which considers the two phases as a single fluid with space dependent viscosity, l, and density, q (see, e.g., [27]). Then, velocity, v , and pressure, p, are governed by the Navier–Stokes equations,

@ qv þ r  ðqvv Þ ¼ r  r þ f @t

ð1Þ

sa ;

supplemented by the divergence-free condition for the velocity,

r  v ¼ 0:

ð2Þ 



Here r ¼ pI þ l rv þ rv is the stress tensor (with I the identity tensor), and f sa the surface force between two fluids. In the following we will consider flow at low Reynolds number, Re ¼ qUa=l < 101 (a is the mean pore throat size), and we will neglect the non-linear term in (1), which leads to a momentum equation of the form T

@ qv ¼ rp þ r  lðrv þ rv T Þ þ f @t

sa :

ð3Þ

Notice that for generality we keep the partial time derivative in the equation above; this allows to consider problems in which temporal scales other than the advective-time scale are present. This is important for two-phase flow processes, which are essentially time-dependent due to the evolution of the fluids distribution and require accounting for timederivative term in momentum equation (see, e.g., [9]). Moreover, an algorithm capable of handling an equation of the form of Eq. (3) can be readily extended to solve full Navier–Stokes equations with explicit treatment of the inertia term. In the whole-domain formulation the distribution of the fluids is described by the color function, a, which represents the volumetric fluid fraction and is advected with velocity v , i.e.,

@a þ r  ðav Þ ¼ 0: @t

ð4Þ

Here we assume a to be equal to 1 in the regions occupied by the nonwetting fluid, nw, and to 0 where only the wetting fluid, w, is present; the interface between the fluids is represented by the surface of discontinuity of a. The VOF method is a rather straightforward discretization of the equations above, which requires two main ingredients: a numerical technique to limit the numerical diffusion while solving (4), and an appropriate model for the surface tension. The simulations presented in the following are performed with OpenFOAM [25], which controls numerical diffusion by introducing an artificial compression term in (4) [26]. The surface tension acting at the fluid–fluid interface is modeled by the continuum surface force method [7], which approximates the surface force, f sa , in Eq. (3) with a volume (body) force, f sv , i.e.,

f

sa

f

sv

¼ cjra;

ð5Þ

where c is surface tension, and

j ¼ r  n ¼ r 

ra jraj

ð6Þ

is the interface curvature. On the solid boundary, @ Xs , zero velocity is imposed, and wetting is modeled by assigning the boundary condition

nj@ Xs ¼ ns cos heq þ t s sin heq ;

ð7Þ

where ns and t s are the unit vectors normal and tangential to solids boundary, respectively; and heq is the equilibrium contact angle. Finally, we have the set of Eqs. (2)–(4), which are coupled through the velocity, v , and the color function, a, because the density and the viscosity are functions of a, i.e.,

qðaÞ ¼ aqnw þ ð1  aÞqw ; lðaÞ ¼ alnw þ ð1  aÞlw ;

ð8Þ

where qnw and lnw (resp. qw and lw ) are the constant density and viscosity of the nonwetting (resp. wetting) fluid. Notice that for zero surface tension, c ¼ 0, and equal properties of the fluids, qnw ¼ qw ¼ q and lnw ¼ lw ¼ l, we obtain the set of equations describing the transport of an ideal tracer. This allows us to describe the hybrid MsFV method starting from the simpler single-phase case, and then extending the result to the two-phase case. 3. The Hybrid MsFV method The hybrid multiscale method is based on three main steps: the solution of a set of local problems to obtain the pore-scale velocity and pressure, which are used to interpolate the coarse-scale degrees of freedom; the construction of a coarse Darcy problem, which is solved to couple the local problems; and the solution of a fluid-advection problem, which can be adaptively solved at the pore or at the Darcy scale. As in the standard MsFV method [13,22], we employ an auxiliary coarse grid together with its dual (Fig. 2). First, a set of ~ e , by solving localized problems. In contrast with standard application of MsFV, interpolators is defined in each dual cell, X

296

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307

Fig. 2. Coarse (solid) and dual (dashed) grids. The coarse-scale degrees of freedom are defined at the coarse cell centers (squares), whose connection forms the dual grid. Shown is the case of a 2  2 dual grid.

here pore-scale velocity and pressure gradient are no longer related through local Darcy’s law, and velocity interpolators, wej , must be explicitly defined together with scalar interpolators for the pressure, /ej . (Here, the subscript j indicates the node, at ~ e .) the center of the coarse cell Xj , for which the velocity or the pressure is interpolated in the dual cell X In principle, different coarse degrees of freedom are required for pressure and velocity. For a hybrid algorithm, however, we are interested in coupling pore-scale Navier–Stokes solution with coarse-scale Darcy equations; this means that coarse velocities depend only on coarse pressure gradients. Therefore, no independent degrees of freedom for velocity are required, and velocity and pressure interpolators share the same degrees of freedom: the coarse scale pressure values, pj . This is our assumption of Darcy flow at the coarse scale and allows constructing only a scalar coarse mass-balance equation (no coarse momentum-conservation equation is necessary). The coarse problem is obtained by requiring mass conservation over coarse cells, which serve as control volumes. The two steps above define an upscaling procedure with numerical closure provided by the definition of the interpolators. The pore-scale pressure and velocity distributions can be obtained by a simple interpolation. However, since the localized problems were solved independently, the interpolated velocity is not conservative across dual cell boundaries and would introduce significant mass-balance errors if transport problem is solved at the fine scale (this has been clearly demonstrated for Darcy flow, see, e.g., [13]). To avoid this issue, a velocity reconstruction (downscaling) step can be applied to obtain an approximate, but fully mass conservative fine-scale velocity field. Analogously to the standard MsFV method [20], we solve a local problem in each coarse cell with prescribed velocity (flux) boundary conditions, which are consistent with the coarsescale fluxes. This reconstructed conservative velocity field is then used in the transport equation. In the following we illustrate the details of these three steps. For clarity of the exposition, we first describe the Hybrid Multiscale Finite Volume (HMsFV) method for single-phase flow. 3.1. Single-phase flow In case of steady single-phase flow, Eqs. (2) and (3) reduce to the Stokes equations:



r  v ¼ 0; rp þ lr2 v ¼ 0:

ð9Þ

~ e , we use the decomposition In each dual cell, X e

pjX~ e ¼

Nc X pj /ej ;

e

v jX~

e

¼

j¼1

Nc X pj wej ;

ð10Þ

j¼1

where N ec is the number of coarse nodes on the dual-cell boundary (for a cartesian grid, N ec ¼ 2d , where d is the number of dimensions); and the interpolators (or basis functions), /ej and wej , are solutions of the local problems

(

r  wej ¼ 0; r/ej þ lr2 wej ¼ 0;

~ e , with the following boundary conditions: which are solved in each dual cell, X

ð11Þ

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307

(1) (2) (3) (4)

297

~ e; unit pressure for node j and zero for all other nodes of X ~ e; linear pressure variation between the coarse nodes on @ X ~ e (which is the same as for original MsFV method); zero normal gradient of the velocity on @ X zero velocity and zero pressure gradient on solid boundary @ Xs .

Conditions 1–3 define the localization assumptions, whereas conditions in 4 are the standard no-slip boundary condition assigned on solid walls. From no-flow requirement in case of constant pressure, we have the constraints e

e

Nc X /j ¼ 1;

Nc X wj ¼ 0;

j¼1

j¼1

ð12Þ

which allow solving only N ec  1 independent problems for the interpolators. To formulate coarse-scale mass-balance equation for cell Xi we use the velocity interpolators to calculate the coarse-scale transmissibilities, T ij , by simple integration of the velocity along the segments of control volume boundary, i.e.,

T ij ¼

XZ

~e @ Xi \X

e

wej  g ds;

ð13Þ

~ e , and @ Xi \ X ~ e is the boundary of the coarse cell relative where wej is the velocity interpolator of node j defined in the dual X to node i contained in the same dual. Using Eq. (13), we obtain coarse-scale mass-balance equations of the form

X T ij ðpi  pj Þ ¼ 0;

i 2 ½1; Nc ;

ð14Þ

j–i

where N c is the number of coarse nodes, and summation is taken over all neighboring nodes of node i, which results in a 3d point stencil discretization of the coarse problem. Next, we proceed with time-dependent Stokes equations and apply first-order fully-implicit time differencing scheme, i.e.,

(

r  v ¼ 0; q vDtv ¼ rp þ lr2 v ;

ð15Þ

 is the velocity at previous time step. where v  =Dt. The most effective The second equation in (15) can be considered as an inhomogeneous equation with source term v way to deal with a source term is to add local correction functions and leave the basis functions unchanged [19,21,15,22]. Therefore we introduce two local correction functions, /e and we for pressure and velocity, respectively, which leads to the new decomposition e

pjX~ e ¼

Nc X pj /ej þ /e ;

e

v jX~

e

¼

j¼1

Nc X pj wej þ we :

ð16Þ

j¼1

Inserting Eq. (16) into Eq. (15) and equating to zero the coefficients of the coarse scale degrees of freedom, pj , we obtain the following splitting between the basis and correction functions:

(

r  wej ¼ 0;

ð17Þ

we

q Djt ¼ r/ej þ lr2 wej ; and

(

r  we ¼ 0;

ð18Þ

e

q wDt v ¼ r/e þ lr2 we :

The boundary conditions are the same as in the steady flow case, except for the pressure correction function for which on ~ e we have @X

/e j@ X~ e ¼ 0:

ð19Þ

The splitting above is similar to the one employed to model compressible Darcy flow [11] and allows keeping interpolators independent of degrees of freedom. Once the velocity correction function, we , is calculated, the induced fluxes across the coarse-cell boundary, @ Xi , are added to the right-hand side of Eq. (14), and we obtain the new coarse problem

X XZ T ij ðpi  pj Þ ¼  j–i

e

~e @ Xi \X

we  g ds:

ð20Þ

298

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307

Note that only the correction functions need to be updated at each time step, while the interpolators are computed only once. 3.2. Two-phase flow Keeping in mind that density and viscosity (8) now depend on the fluid distribution through the color function a, the approach described above can be naturally extended to two-phase flow. For a given distribution of the color function a, the basis functions are the local solutions of the system

8 < r  wej ¼ 0; e

: qðaÞ wj ¼ r/ej þ r  Dt

h



i

lðaÞ rwej þ rweT ; j

ð21Þ

whereas the correction functions are obtained by solving

(

r  we ¼ 0; e







qðaÞ wDt v ¼ r/e þ r  lðaÞ rwe þ rweT þ cjra; 

ð22Þ

where the surface tension term in Eq. (3), f sa  cjra, has been included in the correction functions. The boundary conditions described in the previous section must be supplemented by Eq. (7) as boundary condition for the normal to the fluid– fluid interface at the solids. The coarse problem has again the form of Eq. (20). The transmissibilities defined by Eq. (13) represent an extension of the classic Darcy-scale total mobility in the sense that they do not simply depend on the average value of the saturation, but on the actual spatial distribution of the two fluids. In principle interpolators and transmissibilities should be recomputed every time that the fluid–fluid interface moves and the fluid distribution changes. Since frequent updates might considerably increase the computational costs, it is practically important to devise adaptivity techniques to balance efficiency and accuracy of the method. A variety of strategies have been developed for Darcy-scale applications [14,18,16] and can be naturally incorporated in the present hybrid framework without any significant changes of latter. For Darcy-scale problems, MsFV basis functions are recomputed based on the local fine-scale changes of total mobility [14]. For pore-scale basis functions, an analogous criterion would clearly be too restrictive and integral measures should be adopted. Here, for simplicity, we compare the integrated absolute changes of color function field with a predefined threshold (here we choose a ¼ 103 ), and we update the interpolators if

1 ~ jXe j

Z ~e X

jan  an1 jdV < a ;

ð23Þ

where superscripts n and n  1 denote the new and the old time steps, respectively. 3.3. Mass-conservative velocity field reconstruction and transport The two steps described above (computation of interpolators and solution of coarse-scale problem) constitute a numerical upscaling procedure. To solve the problem with pore-scale resolution, however, we need a pore-scale conservative velocity field that can be used in Eq. (4) to advect the color function. As discussed above, the approximate velocity obtained by simple interpolation, Eq. (16), cannot be used because it is non-conservative across dual cell boundaries due to the localization assumptions. Therefore, we reconstruct the massconservative velocity field by solving local pore-scale problems in the coarse cells, Xi , with fixed velocity boundary conditions obtained by interpolation. Mass conservation is ensured by the consistency between the integral pore-scale mass-fluxes along the boundary @ Xi and the corresponding coarse-scale fluxes; solvability of the local momentum-balance equations is guaranteed by the fact that obstacles can act as momentum sinks. The calculated velocity is then locally used to solve the transport problem in each cell; this localization allows to solve pore-scale transport only in regions of interest (e.g. front regions), whereas in the rest of the domain equations can be solved only for coarse scale (similarly to what previously done for Darcy-scale problem [16]). For the two-phase case, we compute Eq. (4) locally with a global time step defined by the stability restrictions [26] and fixed values for the color function on @ Xi . After a new distribution of fluids is obtained, we close the algorithm loop and proceed to the next time step, starting from construction and solution of the coarse problem. For tracer tracking studies we reconstruct mass-conservative velocity only once and then solve transport equation globally. 4. Numerical simulations In this section, the quality of the HMsFV solution is investigated by comparison with the reference pore-scale solution obtained by solving the same problem on the fine grid. First, we consider a single-phase problem in which an ideal tracer is transported in a steady velocity field; then, we simulate drainage for two pairs of fluids varying the viscosity ratio,

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307

299

Table 1 Pore geometry properties. Domain size Mean obstacle size Mean pore throat size Number of obstacles Number of cells Permeability Porosity

45 mm  56.8 mm 0.6 mm 0.4 mm 2 507 711 720 108 m2 70%

M ¼ llnw , from favorable to unfavorable to illustrate the capabilities of the HMsFV algorithm to handle stable and unstable w flow regimes.

4.1. Pore-geometry and coarse grids We consider two-dimensional fluid flow in idealized pore geometry (Fig. 1), whose properties are summarized in Table 1; and we assign no-slip boundary conditions at the left and right boundaries, as well as on the solid obstacles. For the HMsFV method, several coarse- and dual-grid sizes are considered. We start defining a regular cartesian dual grid which divides the domain into N  N dual cells for which the basis functions (interpolators) are computed as described in Section 3.1. An example of pressure and velocity interpolators is shown in Fig. 3. The coarse degrees of freedom are assigned at the node of the dual cells and the coarse grid is then constructed by connecting the dual centers, which results in ðN þ 1Þ  ðN þ 1Þ coarse cells. An example of the resulting coarse and dual grids is given in Fig. 2 for the case of a 2  2 dual grid. Notice that coarse degrees of freedom are also assigned at the domain boundaries.

Fig. 3. Example of pressure (/i ) and velocity (wi ) interpolators. Blue (resp. red) corresponds to low (resp. high) values. For the velocity, colors indicate the absolute value. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

300

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307

4.2. Single-phase flow A fluid (viscosity

l ¼ 0:1 kg m1 s1 and density q ¼ 1000 kg m3) is injected at the top boundary with inlet velocity

v in ¼ 104 m s1, while the outlet pressure at the bottom boundary is fixed to pout ¼ 0. As the corresponding Reynolds num-

ber is Re  104 , inertia can be neglected and the use of Stokes equations, Eq. (9), is justified. The pressure and velocity fields obtained with a direct pore-scale simulations and with the HMsFV method are compared in Fig. 4 for a 4  4 dual grid (the HMsFV conservative velocity is shown). The HMsFV results are in good agreement with the reference; in particular, velocity field details such as high-flow channels are well reproduced also across coarse-cell boundaries. To further evaluate the performance of the HMsFV method with respect to coarse grid of different sizes, we simulate tracer injection and compare the concentration distribution as a function of time, which can be considered as an integral

Fig. 4. Comparison of pressure (top) and magnitude of the velocity (bottom): (a) direct fine-scale solution, (b) HMsFV results (coarse grid is 4  4). Hybrid pressure field is obtained by simple interpolation. Blue (resp. red) corresponds to low (resp. high) values. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

301

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307 Table 2 Number of dual cells, number of coarse cells, dual-cell sizes and scale ratios for different coarse grids (H is dual-cell size and a is the mean pore-throat size). Dual-grid size

Coarse-grid size

H, mm

a=H

22 44 88 16  16

33 55 99 17  17

22.5 11.25 5.625 2.8125

0.018 0.036 0.071 0.142

Fig. 5. Tracer concentration at t ¼ 300 s after injection: (a) fine-scale solution, (b)–(e) HMsFV solutions for different coarse-grid sizes.

measure of the velocity error. The tracer is injected at the top boundary and then transported in accordance to the computed stationary velocity field. We consider purely advective tracer transport and numerical dispersion is controlled by adding an artificial compression term acting only at the front [26]. The results for fine-scale velocity field and velocities computed using different coarse grids (see Table 2) are presented in Fig. 5. (Notice the good qualitative agreement with experiments reported in the literature for similar geometries, e.g., [5]). For the 2  2 dual grid, the local concentration is in excellent agreement with the pore-scale reference solutions. However, for the 16  16 dual grid, local deviation from the reference can be clearly observed (Fig. 6). This deterioration of the solution is expected, because our Darcy assumption at the coarse scale requires separation between the ‘‘scale of observation’’ (the dual-cell size, H) and the pore-size (the mean pore-throat size, a). When this condition is satisfied, the interaction with many obstacles contained in a dual is the dominating momentum-transfer mechanism, and the fluid–fluid viscous exchange at the dual boundary is negligible. As the momentum transfer to the obstacles is proportional to the mean velocity,

302

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307

Fig. 6. Difference between concentration fields (a) and (e) from Fig. 5. The hybrid coarse grid has 16  16 cells, and the ratio between pore and cell sizes is a=H  0:142. Blue (resp. red) corresponds to negative (resp. positive) values. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7. Horizontally averaged tracer concentration at different times for reference finescale solution and hybrid solutions with different coarse grid sizes. HMsFV concentration profiles are in a very good agreement with the fine-scale reference for all coarse grid sizes.

303

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307

Fig. 8. HMsFV errors of pressure (top) and velocity (bottom) as a functions of coarse-grid size.

Table 3 Fluids properties. Property

Nonwetting fluid

Wetting fluid

Density, kg m3 Viscosity, kg m1 s1

1000 0.15 or 0.01

800 0.05

Viscosity ratio Interfacial tension, kg s2 Contact angle, deg.

3 or 0.2 0.064 30

this results in Darcy’s law relating mean velocity and pressure drop. For a 2  2 dual grid, we have the scale ratio a=H  102 , whereas for a 16  16 grid it is a=H  101 , which explains the less accurate solution of the latter. Despite the local fluctuation, however, the average concentration profiles (Fig. 7) are in excellent agreement for all coarse grids. Finally, we study the convergence behavior of the HMsFV method as a function of coarse grid size. Fig. 8 shows the L1 - and the L2 -norm of the HMsFV pressure and velocity errors (for the velocity we consider only dominating y-component). As also indicated by the tracer simulation, the errors of hybrid solution increase with coarse grid size (Fig. 8) due to the weaking of the scale separation hypothesis necessary for our assumption of a Darcy-like flow at the coarse scale. This is a quite different situation than in the standard MsFV method (applied to a fine-scale Darcy problem) where a larger coarse grid size yields to a more accurate solution due to the increased number of coarse degrees of freedom [13]. For Stokes flow, the results of hybrid algorithm are (locally) reliable when scale ratio is small: a=H < 0:1. 4.3. Two-phase flow Next, we consider two drainage experiments in the same pore geometry. Varying the viscosity of the nonwetting phase,

lnw , we examine the HMsFV method for stable and unstable displacement regimes. The properties of the fluids are summa-

rized in Table 3. The nonwetting fluid is injected at the top boundary with constant inlet velocity v in ¼ 102 m s1, whereas the pressure at the bottom is pout ¼ 0. The pore-scale Reynolds number is Re  102 , and capillary number is Ca ¼ lw a2 v in =ðkrÞ  0:1 (where k is the absolute permeability, which is measured numerically from a single-phase simulation in the geometry under consideration). We start from the stable case which corresponds to the viscosity ln ¼ 0:15 kg m1 s1 and viscosity ratio M ¼ 3. The phase distributions for t ¼ 3:05 s are shown in Fig. 9 both for the pore-scale reference and for the HMsFV solution. The

304

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307

Fig. 9. Phases distributions at t ¼ 3:05 s after injection for viscosity ratio M ¼ 3.

fluid–fluid interface is stable in both cases, but the hybrid solution shows a slightly larger amount of trapped wetting phase behind the invading front. From visual observation of the fluid distribution it appears that trapping is slightly enhanced in the proximity of coarse cell boundaries. This is caused by the perturbation of the vorticity field that is introduced when the velocity is forced to be conservative [16]. This amplified vorticity yields a rougher front in this region, which makes trapping more likely. To quantitatively compare the results, we plot the horizontally averaged saturation profile in Fig. 10 and the temporal evolution of the fluid–fluid interface (area is normalized by the initial value) in Fig. 11. The main characteristics of the flow are captured, although, as effects of lager trapping, the hybrid solution shows slightly larger interfacial area and tip velocity. Next, we consider the unstable regime and specify the viscosity of the nonwetting fluid equal to ln ¼ 0:01 kg m1 s1, which correspond to the viscosity ratio M ¼ 0:02. The transitions to viscous instability is very well captured by the hybrid approach: the solution exhibit branched fingers and large clusters of trapped wetting fluid that appear qualitatively similar to those visible in the pore-scale reference solution (Fig. 12). Notice that exact agreement between the solution is impossible in unstable regimes because the hybrid method employs an approximate velocity field; since small differences are amplified with time, the exact position of fingers and trapped clusters will differ and only their statistics can be capture. To evaluate this we consider again the horizontally averaged saturation profile (Fig. 13) and the temporal evolution of the fluid–fluid interface (Fig. 14). These two quantities are in excellent agreement with the reference. This shows the ability of the HMsFV method to correctly describe the transition from stable displacement to viscous unstable regime.

5. Conclusions We have presented a novel hybrid multiscale method for multiphase flow in porous media. Local Navier–Stokes solutions are computed to obtain the pore-scale velocity field and the VOF method is used to model the evolution of the fluid–fluid interface; a coarse mass-conservation problem, which couples the local solutions, is obtained by averaging the pore-scale velocities over the cells of a coarse grid, which act as control volumes. The HMsFV method relies on an extension of the standard MsFV method, but, in addition to pressure interpolators, it requires specific interpolators (basis and correction functions) for the velocity, which at the pore-scale is not simply proportional to the gradient of the pressure. In addition to the localization assumptions, which are quite straightforward extensions of the standard MsFV (i.e., linear pressure variation along the boundaries and zero normal gradient of the velocity), we have postulated that the coarse-scale fluxes are proportional to the coarse-pressure differences. This yields a coarse massconservation equation with fluxes that obey a generalized Darcy equation and transmissibilities that depend on the actual fluid distribution rather than simply on the averaged saturation.

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307

305

Fig. 10. Average nonwetting phase saturation profiles at different times for viscosity ratio M ¼ 3.

Fig. 11. Time evolution of the normalized fluid–fluid interfacial area for viscosity ratio M ¼ 3.

We have demonstrated that these assumptions provide hybrid solutions that are in good agreement with reference porescale solutions for single- and two-phase flow. In particular, the HMsFV is able to model the transition from viscous stable to viscous unstable flow, capturing the statistics of the morphology of the invading front when this spans several coarse cells. Notice that, as the methods employs an approximate velocity, a deterministic reproduction of viscous fingers is impossible, because small differences grow with time due to the unstable nature of the flow. Although in the simulations presented here we have computed a pore-scale velocity field in the entire domain, this is in general not necessary and our framework allows considering pore-scale fluxes only in some regions, while Darcy fluxes (based simply on the coarse-scale velocity) are used in the rest of the domain. The methods can applied to problems in which local flow conditions yield unstable flow regimes that require a finer description of the flow; in this case appropriate adaptive criteria can be introduced to identify the unstable-flow regions (similar to what proposed by Künze and Lunati [16] for Darcy-scale gravity fingering).

306

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307

Fig. 12. Phases distributions at t ¼ 1:9 s after injection for viscosity ratio M ¼ 0:2.

Fig. 13. Average nonwetting phase saturation profiles at different times for viscosity ratio M ¼ 0:2.

Finally, we remark that the coarse-scale description of the flow is based on volume averaging and on the assumption of scale separation, which is necessary to assume that the relationship between coarse-scale fluxes and pressure differences is local and does not require a true coarse-scale momentum-balance equation. If the latter assumption is not satisfied, the quality of the HMsFV solution might deteriorate. We have, therefore, a numerical tool to investigate the limits of validity of the Darcy equation and to understand the link between pore-scale quantities and their corresponding Darcy-scale variables.

P. Tomin, I. Lunati / Journal of Computational Physics 250 (2013) 293–307

307

Fig. 14. Time evolution of the normalized interfacial area for viscosity ratio M ¼ 0:2.

Acknowledgments This project is supported by the Swiss National Science Foundation Grant No. PP00P2-123419/1, and by the University of Lausanne. The authors wish to thank Andrea Ferrari and Rouven Künze for many constructive discussions. References [1] S. Afkhami, S. Zaleski, M. Bussmann, A mesh-dependent model for applying dynamic contact angles to VOF simulations, Journal of Computational Physics 228 (15) (2009) 5370–5389. [2] S. Alyaev, E. Keilegavlen, J.M. Nordbotten, Analysis of control volume heterogeneous multiscale methods for single phase flow in porous media, Multiscale Modeling and Simulation, 2013 (submitted for publication). [3] I. Battiato, D.M. Tartakovsky, A.M. Tartakovsky, T. Scheibe, On breakdown of macroscopic models of mixing-controlled heterogeneous reactions in porous media, Advances in Water Resources 32 (11) (2009) 1664–1673. [4] I. Battiato, D.M. Tartakovsky, A.M. Tartakovsky, T. Scheibe, Hybrid models of reactive transport in porous and fractured media, Advances in Water Resources 34 (9) (2011) 1140–1150. [5] A. Birovljev, K.J. Måløy, J. Feder, T. Jøssand, Scaling structure of tracer dispersion fronts in porous media, Physical Review E 49 (6) (1994) 5431–5437. [6] G. Bonfigli, P. Jenny, An efficient multi-scale Poisson solver for the incompressible Navier–Stokes equations with immersed boundaries, Journal of Computational Physics 228 (12) (2009) 4568–4587. [7] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, Journal of Computational Physics 100 (2) (1992) 335–354. [8] H. Darcy, Les Fontaines Publiques de la Ville de Dijon, Dalmont, Paris, 1856. [9] B. Eckhardt, J. Buehrle, Time-dependent effects in high viscosity fluid dynamics, The European Physical Journal Special Topics 157 (1) (2008) 135–148. [10] A. Ferrari, I. Lunati, Direct numerical simulations of interface dynamics to link capillary pressure and total surface energy, Advances in Water Resources 57 (2013) 19–31. [11] H. Hajibeygi, P. Jenny, Multiscale finite-volume method for parabolic problems arising from compressible multiphase flow in porous media, Journal of Computational Physics 228 (14) (2009) 5129–5147. [12] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, Journal of Computational Physics 39 (1) (1981) 201–225. [13] P. Jenny, S.H. Lee, H.A. Tchelepi, Multi-scale finite-volume method for elliptic problems in subsurface flow simulation, Journal of Computational Physics 187 (1) (2003) 47–67. [14] P. Jenny, S.H. Lee, H.A. Tchelepi, Adaptive multiscale finite-volume method for multi-phase flow and transport in porous media, Multiscale Modeling and Simulation 3 (1) (2005) 50–64. [15] P. Jenny, I. Lunati, Modeling complex wells with the multi-scale finite-volume method, Journal of Computational Physics 228 (3) (2009) 687–702. [16] R. Künze, I. Lunati, An adaptive multiscale method for density-driven instabilities, Journal of Computational Physics 231 (17) (2012) 5557–5570. [17] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, G. Zanetti, Modelling merging and fragmentation in multiphase flows with SURFER, Journal of Computational Physics 113 (1) (1994) 134–147. [18] S.H. Lee, H. Zhou, H.A. Tchelepi, Adaptive multiscale finite-volume method for nonlinear multiphase transport in heterogeneous formations, Journal of Computational Physics 228 (24) (2009) 9036–9058. [19] I. Lunati, P. Jenny, The multiscale finite volume method: a flexible tool to model physically complex flow in porous media, in: Proceedings of ECMOR X, Amsterdam, The Netherlands, September 4–7, 2006. [20] I. Lunati, P. Jenny, Multiscale finite-volume method for compressible multiphase flow in porous media, Journal of Computational Physics 216 (2) (2006) 616–636. [21] I. Lunati, P. Jenny, Multiscale finite-volume method for density-driven flow in porous media, Computational Geosciences 12 (3) (2008) 337–350. [22] I. Lunati, S.H. Lee, An operator formulation of the multiscale finite-volume method with correction function, Multiscale Modeling and Simulation 8 (1) (2009) 96–109. [23] I. Lunati, D. Or, Gravity-driven slug motion in capillary tubes, Physics of Fluids 21 (052003) (2009). [24] I. Lunati, M. Tyagi, S.H. Lee, An iterative multiscale finite volume algorithm converging to the exact solution, Journal of Computational Physics 230 (5) (2011) 1849–1864. [25] OpenFOAM, The Open Source CFD Toolbox: User Guide. OpenFOAM Foundation, version 2.1.1, 2012. [26] H. Rusche, Computational fluid dynamics of dispersed two-phase flows at high phase fractions, Ph.D. Thesis, Imperial College London, 2002. [27] R. Scardovelli, S. Zalesky, Direct numerical simulation of free-surface and interfacial flow, Annual Review of Fluid Mechanics 31 (1999) 567–603.