1 Introduction
Multiscale transport phenomena are fundamentally important for multiphase porous medium applications. Two scales are commonly of interest: the microscale and the macroscale. At the microscale, the locations of all phases are resolved explicitly. Complete knowledge of quantities such as the distribution of the solid phase, the resultant pore morphology and topology, fluid velocities, the location of boundaries between phases, and curves that form along the joint boundary formed at the common edge of three phases are fully resolved. The maximum length scales that can typically be studied for natural porous medium systems using microscale approaches are typically less than a centimetre (Blunt et al. Reference Blunt, Bijeljic, Dong, Gharbi, Iglauer, Mostaghimi, Paluszny and Pentland2013; Ovaysi, Wheeler & Balhoff Reference Ovaysi, Wheeler and Balhoff2014). At the macroscale, system properties at a ‘point’ are described in an averaged sense (Hubbert Reference Hubbert1956; Bear Reference Bear1972; Luckner, van Genuchten & Nielsen Reference Luckner, van Genuchten and Nielsen1989). Macroscale approaches are essential to describe transport processes that occur over length scales of metres to kilometres that are routinely of interest in subsurface systems. At these length scales, the system is represented as an overlapping continuum, using some set of averaged quantities such as porosity, volume fractions and specific interfacial areas. While the macroscale is the typical scale at which porous medium continuum models are posed, solved and applied, microscale modelling can be used to evaluate, close and validate macroscale models (Whitaker Reference Whitaker1986; Wood Reference Wood2009; Gray & Miller Reference Gray and Miller2014).
Multiscale modelling approaches can only be successful when the closure relations used to produce solvable macroscale models are based upon and consistent with the behaviour of microscale systems in some averaged sense; the microscale details of such systems must be understood, at least for some model cases, to connect the scales and realistically and rigorously close macroscale models. Thus microscale modelling is of considerable importance for advancing fundamental understanding and for the evolution of porous medium models that seek to rigorously connect disparate length scales. As an example, the thermodynamically constrained averaging theory (TCAT) is an approach for precisely describing macroscale transport phenomena in terms of specifically defined microscale quantities, making the connection across scales clear and transparent (e.g. Gray & Miller Reference Gray and Miller2005, Reference Gray and Miller2014; Miller & Gray Reference Miller and Gray2005; Jackson, Miller & Gray Reference Jackson, Miller and Gray2009). Within this context, sources of microscale information represent an essential model development tool, providing a resource to measure quantities that impact macroscopic behaviour, assess model assumptions and inform closure relationships. In particular, interface and common curve phenomena are thought to be important to fully characterize the behaviour of multiphase porous medium systems; existing macroscopic models do not track this behaviour explicitly (Miller et al. Reference Miller, Dawson, Farthing, Hou, Huang, Kees, Kelley and Langtangen2013).
Pore-scale models of multiphase flow in porous medium systems have been deployed widely as a way to obtain microscale information based on first principles while mimicking macroscopic behaviour (Celia, Reeves & Ferrand Reference Celia, Reeves and Ferrand1995; Reeves & Celia Reference Reeves and Celia1996; Blunt & Hilpert Reference Blunt and Hilpert2001; Pan, Hilpert & Miller Reference Pan, Hilpert and Miller2004; Schaap et al. Reference Schaap, Porter, Christensen and Wildenschild2007; Joekar-Niasar, Hassanizadeh & Leijnse Reference Joekar-Niasar, Hassanizadeh and Leijnse2008; Joekar-Niasar, van Dijke & Hassanizadeh Reference Joekar-Niasar, van Dijke and Hassanizadeh2012). Pore-scale studies have proved useful as a way to study limitations of macroscopic models for multiphase flow by directly simulating phenomena that are not incorporated in traditional macroscopic formulations (Li, Pan & Miller Reference Li, Pan and Miller2005; Porter, Schaap & Wildenschild Reference Porter, Schaap and Wildenschild2009; Porter et al. Reference Porter, Wildenschild, Grant and Gerhard2010). A topic of considerable interest is the relationship between transport phenomena in a multiphase porous medium system at the microscale and a corresponding macroscale representation of the same system (Joekar-Niasar et al. Reference Joekar-Niasar, van Dijke and Hassanizadeh2012; Blunt et al. Reference Blunt, Bijeljic, Dong, Gharbi, Iglauer, Mostaghimi, Paluszny and Pentland2013; Herring et al. Reference Herring, Harper, Andersson, Sheppard, Bay and Wildenschild2013). Many constitutive relationships have been posited in evolving models for two-fluid-phase flow in porous medium systems (Gray & Miller Reference Gray and Miller2014), but only a fraction of these have been studied. A prominent example that has been a topic of sustained study is the equilibrium relationship involving fluid saturations, pressures and interfacial areas (Reeves & Celia Reference Reeves and Celia1996; Held & Celia Reference Held and Celia2001a ,Reference Held and Celia b ; Dalla, Hilpert & Miller Reference Dalla, Hilpert and Miller2002; McClure et al. Reference McClure, Pan, Adalsteinsson, Gray, Miller, Miller, Farthing, Gray and Pinder2004; Niessner & Hassanizadeh Reference Niessner and Hassanizadeh2008; Porter et al. Reference Porter, Schaap and Wildenschild2009, Reference Porter, Wildenschild, Grant and Gerhard2010; Raeesi & Piri Reference Raeesi and Piri2009; Joekar-Niasar, Hassanizadeh & Dahle Reference Joekar-Niasar, Hassanizadeh and Dahle2010a ; Joekar-Niasar et al. Reference Joekar-Niasar, Prodanović, Wildenschild and Hassanizadeh2010b ; Porter & Wildenschild Reference Porter and Wildenschild2010; Ahrenholz et al. Reference Ahrenholz, Niessner, Helmig and Krafczyk2011; Landry, Karpyn & Piri Reference Landry, Karpyn and Piri2011; Joekar-Niasar & Hassanizadeh Reference Joekar-Niasar and Hassanizadeh2012). Beyond the determination of interfacial areas, the measurement of microscale interfacial curvatures has also been applied to porous medium systems (Armstrong, Porter & Wildenschild Reference Armstrong, Porter and Wildenschild2012). Approaches have also been advanced to study the contact angle in both simulated and experimental systems (Huang et al. Reference Huang, Thorne, Schaap and Sukop2007; Andrew, Bijeljic & Blunt Reference Andrew, Bijeljic and Blunt2014).
While computational approaches have been advanced to analyse static configurations, macroscale models for two-phase flow must be able to accurately model dynamic processes. Interfacial displacements occur for any cases that involve changing fluid saturations. Macroscale models to predict the rate of change in saturation are therefore of great interest, and an area of active research. Among the quantities of interest are the kinematic velocities of the interfaces and the common curve (Gray et al. Reference Gray, Dye, McClure, Pyrak-Nolte and Miller2015). Accurate calculation of these quantities is often hindered by the presence of spurious currents that arise in the vicinity of interfaces as a consequence of numerical instabilities. These instabilities are well known for many numerical methods, including lattice-Boltzmann methods (LBMs) (Hou et al. Reference Hou, Shan, Zuo, Doolen and Soll1997), as well as finite-element method (FEM) schemes and other methods that represent the interfacial region in a diffuse manner (Hysing Reference Hysing2012; Zahedi, Kronbichler & Kreiss Reference Zahedi, Kronbichler and Kreiss2012). Spurious currents have been observed in many multiphase LBMs and have been shown to arise due to an imbalance in the discretization of the phase pressure terms and the interfacial stresses (Lee & Fischer Reference Lee and Fischer2006; Lee & Liu Reference Lee and Liu2008; Lee Reference Lee2009). While a variety of strategies have been developed to reduce or eliminate these currents in the LBM (Lee & Fischer Reference Lee and Fischer2006; Lee Reference Lee2009) and using other approximation methods (Aulisa, Manservisi & Scardovelli Reference Aulisa, Manservisi and Scardovelli2006; Jafari, Shirani & Ashgriz Reference Jafari, Shirani and Ashgriz2007; Nourgaliev, Liou & Theofanous Reference Nourgaliev, Liou and Theofanous2008; Chang, Deng & Theofanous Reference Chang, Deng and Theofanous2013), higher-order approximations are typically required. The associated computational costs make these methods impracticable for simulations of multiphase flow in large complex geometries such as porous media. Methods that can address these errors at reduced computational cost are therefore desirable for these applications.
The overall goal of this work is to develop an improved numerical approach to model the dynamic behaviour of the fluid–fluid interface and common curve in two-fluid-phase flow, and to evaluate this method within a multiscale framework. The specific objectives of this work are as follows:
-
(1) to implement an efficient method to compute accurate interface and common curve velocities that is not sensitive to spurious currents;
-
(2) to use the approach to show that the static and dynamic contact angles can be accurately predicted based on a multiphase LBM; and
-
(3) to demonstrate that the interface and common curve velocities can be computed consistently for multiphase flow in complex geometries.
2 Approach
2.1 System description
Evolving models of multiphase flow explicitly contain measures of the pore-scale physics that originate from a first-principles description of microscale systems (e.g. Gray & Miller Reference Gray and Miller2014). An important use of microscale simulations is to characterize the macroscale state of the system. This information can be used to supplement laboratory experiments often used to determine the form and parameter values of a closure relation. For two-fluid-phase flow we consider simulations performed within a domain
${\it\Omega}$
. Three phases are present, including the wetting fluid phase (
$w$
), the non-wetting fluid phase (
$n$
) and the solid (
$s$
). Each of the phases occupies a three-dimensional subset of
${\it\Omega}$
, denoted as
${\it\Omega}_{w}$
,
${\it\Omega}_{n}$
and
${\it\Omega}_{s}$
, respectively. Since three phases are present, three interfaces are possible, each of which occupies a two-dimensional domain within
${\it\Omega}$
. These include the interface between the wetting and non-wetting fluids,
${\it\Omega}_{wn}$
, the interface between the wetting fluid and the solid,
${\it\Omega}_{ws}$
, and the interface between the non-wetting fluid and the solid,
${\it\Omega}_{ns}$
. Finally, a common curve can exist where all three phases meet, which is a one-dimensional subset of
${\it\Omega}$
denoted by
${\it\Omega}_{wns}$
. The complete set of entities for the two-fluid-phase model includes all phases and interfaces in addition to the common curve, which together form a set of entities with a corresponding index set
$\mathscr{I}=\{w,n,s,wn,ws,ns,wns\}=\mathscr{I}_{P}\cup \mathscr{I}_{I}\cup \mathscr{I}_{C}$
, where
$\mathscr{I}_{P}$
is the index set of phases
$\mathscr{I}_{P}=\{w,n,s\}$
,
$\mathscr{I}_{I}$
is the index set of interfaces
$\mathscr{I}_{I}=\{wn,ws,ns\}$
, and
$\mathscr{I}_{C}$
is the index set of common curves
$\mathscr{I}_{C}=\{wns\}$
. The fluid phase domains are
$\mathscr{D}_{\mathit{f}}=\{{\it\Omega}_{w},{\it\Omega}_{n}\}$
with the corresponding index set
$\mathscr{I}_{f}=\{w,n\}$
.
2.2 Microscale computational approach
Macroscale state variables can be expressed in terms of microscale information that includes the distribution of phases, and the pressure, density and velocity fields of the fluids. Any approach that yields all of the needed information with sufficient accuracy, precision and efficiency could be used to compute the macroscopic state. The microscale geometry of the solid was provided initially and was invariant with respect to time. The solid phase morphology and topology were specified analytically by constructing a signed distance function to represent the boundary of the solid phase.
An LBM was used to simulate the microscale state. A multiple-relaxation-time (MRT) ‘colour’ LBM method was used, including in situ analysis capabilities (McClure et al.
Reference McClure, Wang, Prins, Miller and Feng2014b
). This approach included the solution of a conservation-of-momentum equation and a conservation-of-mass equation for each of the two fluids. The colour LBM was chosen for several reasons. First, the method conserves exactly both local mass and local momentum. Second, the method is able to resolve immiscible fluid features without the non-physical shrinkage of bubbles that has been associated with some other LBMs (Zheng et al.
Reference Zheng, Lee, Guo and Rumschitzki2014). Finally, the method does not require higher-order approximations to improve numerical stability, which allows for computationally efficient implementations that are suitable for simulations of large three-dimensional flows in complex geometries. A three-dimensional 19-velocity-vector (D3Q19) set was used to construct an LBM that describes the momentum transport. The three-dimensional seven-velocity (D3Q7) set was used to describe the conservation of mass. Density values
${\it\rho}_{w}$
and
${\it\rho}_{n}$
were obtained by solving a separate pair of equations that recover the conservation of mass for each fluid phase, yielding density fields,
${\it\rho}_{w}$
and
${\it\rho}_{n}$
, for the wetting and non-wetting entities, respectively. The associated distribution functions were determined based on the local density, momentum and colour gradient (Latva-Kokko & Rothman Reference Latva-Kokko and Rothman2005). The momentum transport equation was coupled to the mass-conservation equation through the evolution of the non-dimensional density field
$$\begin{eqnarray}{\it\phi}=\frac{{\it\rho}_{n}-{\it\rho}_{w}}{{\it\rho}_{n}+{\it\rho}_{w}},\end{eqnarray}$$
where the densities of the two fluids are different.
The momentum transport equation was formulated using an MRT collision operator to approximate the Navier–Stokes equations with anisotropic contributions to the stress tensor due to interfacial forces (Ahrenholz et al.
Reference Ahrenholz, Toelke, Lehmann, Peters, Kaestner, Krafczyk and Durner2008). The momentum transport solution provides the microscale pressure fields,
$p_{w}$
and
$p_{n}$
, and the velocity fields,
$\boldsymbol{v}_{w}$
and
$\boldsymbol{v}_{n}$
. The relevant parameters for the momentum transport are the viscosities for each fluid,
${\it\mu}_{w}$
and
${\it\mu}_{n}$
, and the interfacial tension,
${\it\gamma}_{wn}$
, which were chosen explicitly. The units of all parameters are specified in terms of the lattice mass
$m$
, the lattice length
${\it\delta}x$
and the lattice time
${\it\delta}t$
. The viscosities are specified in units
$m({\it\delta}x{\it\delta}t)^{-1}$
, the interfacial tension in units
$m{\it\delta}t^{-2}$
and pressure in units
$m{\it\delta}x^{-1}{\it\delta}t^{-2}$
. The lattice quantities can be associated to physical units with the understanding that the speed of sound
$c_{s}={\it\delta}x/(\sqrt{3}{\it\delta}t)$
does not match the physical speed of sound, and the approach will only yield a valid approximation for the continuum mechanics when the Mach number is small.
Conditions at the solid boundary were imposed by enforcing a bounce-back rule for the momentum transport equation to set a no-slip condition (Ginzburg & d’Humiéres Reference Ginzburg and d’Humiéres2003; Lallemand & Luo Reference Lallemand and Luo2003; Pan, Luo & Miller Reference Pan, Luo and Miller2006). For two-fluid-phase flow, setting the equilibrium contact angle is essential to accurately model the system behaviour. A boundary condition was determined by setting a constant value
${\it\phi}(\boldsymbol{x}_{k})={\it\phi}_{s}$
for
$\boldsymbol{x}_{k}\in {\it\Omega}_{s}$
for the solid phase to yield the desired contact angle (McClure, Prins & Miller Reference McClure, Prins and Miller2014a
).
2.3 Approximations of microscale geometric and kinematic quantities
The mean curvature of the
$wn$
interface is one of the most important quantities to evaluate in order to advance macroscopic modelling approaches (McClure et al.
Reference McClure, Adalsteinsson, Wildenschild, Gray and Miller2006; Armstrong et al.
Reference Armstrong, Porter and Wildenschild2012). The mean curvature is defined as
$J_{w}=\boldsymbol{{\rm\nabla}}^{\prime }\boldsymbol{\cdot }\boldsymbol{n}_{w}$
, where
$\boldsymbol{n}_{w}$
is the outward normal vector to the wetting phase boundary. The normal vector to the wetting phase can be computed as the gradient of the phase density field evaluated on
${\it\Omega}_{wn}$
by
$$\begin{eqnarray}\boldsymbol{n}_{w}=\frac{\boldsymbol{{\rm\nabla}}{\it\phi}}{|\boldsymbol{{\rm\nabla}}{\it\phi}|}\Bigg|_{{\it\phi}=0}.\end{eqnarray}$$
The normal to the solid phase can be computed in a similar way using a signed distance function
${\it\psi}$
to represent the solid interface location, where the vector normal to the solid surface is obtained by computing
$$\begin{eqnarray}\boldsymbol{n}_{s}=\frac{\boldsymbol{{\rm\nabla}}{\it\psi}}{|\boldsymbol{{\rm\nabla}}{\it\psi}|}\Bigg|_{{\it\psi}=0}.\end{eqnarray}$$
Since the divergence of the normal vector is the curvature,
$J_{w}$
can be computed by relying on the information from the phase density field
${\it\phi}$
evaluated on the lattice (Sethian Reference Sethian1999) as
$$\begin{eqnarray}\displaystyle J_{w}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\frac{\boldsymbol{{\rm\nabla}}{\it\phi}}{|\boldsymbol{{\rm\nabla}}{\it\phi}|}\Bigg|_{{\it\phi}=0}. & & \displaystyle\end{eqnarray}$$
Finite-difference approximations were used for the spatial derivatives to evaluate the right-hand side of (2.4), and linear interpolation was used to estimate the interface curvature from the curvature of
${\it\phi}$
on the mesh.
Many multiphase implementations of the LBM are hindered by the presence of spurious currents in the vicinity of the interfacial region. These spurious currents arise due to a small disparity in the way the pressure and interfacial contributions to the stress tensor are discretized, and alternative formulations have been proposed to address the issue (Lee & Lin Reference Lee and Lin2005; Lee & Fischer Reference Lee and Fischer2006; Lee & Liu Reference Lee and Liu2008; Lee Reference Lee2009; Kuzmin & Mohamad Reference Kuzmin and Mohamad2010). However, the colour model implementation used in this work exhibits spurious currents and for this reason the interfacial velocity cannot be evaluated accurately by interpolating the values from the simulated velocity field. This model does, however, have significant computational advantages compared to the higher-order methods used to produce accurate interfacial velocities. Therefore, we compute the interfacial velocity by rearranging a level-set equation to provide the interface velocity in the direction normal to the interface, which we denote as
${\it\zeta}$
. The level-set equation (LSE) is (Sethian Reference Sethian1999)
$$\begin{eqnarray}\frac{\partial {\it\phi}}{\partial t}+{\it\zeta}|\boldsymbol{{\rm\nabla}}{\it\phi}|=0.\end{eqnarray}$$
Most typically, (2.5) is solved to obtain the phase density field
${\it\phi}$
and as a result update the position of the interface. In our case, the value of
${\it\phi}$
is already determined from the LBM, and we wish to determine
${\it\zeta}$
. We rearrange the LSE to provide
$$\begin{eqnarray}{\it\zeta}=-\frac{\partial {\it\phi}}{\partial t}\Big/|\boldsymbol{{\rm\nabla}}{\it\phi}|,\end{eqnarray}$$
which must then be evaluated numerically. Both
$\boldsymbol{{\rm\nabla}}{\it\phi}$
and the time derivative of
${\it\phi}$
are computed using second-order finite-difference approximations.
In our formulation, the normal velocity
${\it\zeta}$
at time
$t$
cannot be computed until time
$t+{\it\delta}t$
. However, since this information is not needed to advance the simulation, it has no impact on the viability of the analysis procedure. In order to advance the analysis, a short time history must be retained including
${\it\phi}(x,t-{\it\delta}t)$
,
${\it\phi}(x,t)$
and
${\it\phi}(x,t+{\it\delta}t)$
. Once the time derivative of
${\it\phi}$
has been computed at all lattice sites, the value of
${\it\zeta}$
can be approximated at points on the interface by interpolation. The kinematic interfacial velocity vector at points on the interface is then provided by
$$\begin{eqnarray}\boldsymbol{w}_{wn}={\it\zeta}\boldsymbol{n}_{w}.\end{eqnarray}$$
The contact angle is an indicator of the tendency of each fluid to wet the solid surface, which is commonly called the wettability. The ability to measure contact angles is valuable since the wettability of natural porous media often varies in both time and space, and it has a significant impact on phase connectivity and transport phenomena. At the microscale, the contact angle can be expressed in terms of the normal vectors to the boundaries of wetting phase and solid phase surfaces, defined in (2.2) and (2.3), respectively. Interpolating the normal vectors to points on the common curve allows the microscale contact angle to be computed as
$$\begin{eqnarray}\cos {\it\varphi}_{ws,wn}^{\,}=\boldsymbol{n}_{w}\boldsymbol{\cdot }\boldsymbol{n}_{s}.\end{eqnarray}$$
The averaging for the macroscale contacts is defined such that the Pythagorean identity holds at the macroscale. To accomplish this, the sine of the contact angle is also evaluated at the microscale,
$$\begin{eqnarray}\sin {\it\varphi}_{ws,wn}^{\,}=\sqrt{1-\cos ^{2}{\it\varphi}_{ws,wn}^{\,}}.\end{eqnarray}$$
The approach described to evaluate the kinematic velocity of the
$wn$
interface can be extended to consider the kinematic velocity of the common curve. The key modification needed to accomplish this is to restrict the colour gradient operator so that it includes only those components tangential to the solid surface,
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}^{\prime }{\it\phi}=\boldsymbol{{\rm\nabla}}{\it\phi}-(\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{n}_{s})\boldsymbol{n}_{s}.\end{eqnarray}$$
The normal speed of the common curve is then computed in an analogous fashion to (2.6),
$$\begin{eqnarray}{\it\zeta}^{\prime }=-\frac{\partial {\it\phi}}{\partial t}\Big/|\boldsymbol{{\rm\nabla}}^{\prime }{\it\phi}|.\end{eqnarray}$$
Using the fact that both
${\it\phi}$
and
${\it\psi}$
are constant on the common curve, it is clear that the normal to the common curve within the solid surface is
$$\begin{eqnarray}\boldsymbol{n}_{ws}=\frac{\boldsymbol{{\rm\nabla}}^{\prime }{\it\phi}}{|\boldsymbol{{\rm\nabla}}^{\prime }{\it\phi}|}.\end{eqnarray}$$
The kinematic velocity of the
$wns$
common curve is then
$$\begin{eqnarray}\boldsymbol{w}_{wns}={\it\zeta}^{\prime }\boldsymbol{n}_{ws}.\end{eqnarray}$$
2.4 Measures of macroscale state
Macroscale measures of the system state are computed by averaging of microscale quantities over the phases, interfaces and the common curve. The averages of interest can be computed with an averaging operator of the general form
$$\begin{eqnarray}\left\langle P_{i}\right\rangle _{{\it\Omega}_{{\it\alpha}},{\it\Omega}_{{\it\beta}},W}=\frac{\displaystyle \int _{{\it\Omega}_{{\it\alpha}}}WP_{i}\,\text{d}\mathfrak{r}}{\displaystyle \int _{{\it\Omega}_{{\it\beta}}}W\,\text{d}\mathfrak{r}},\end{eqnarray}$$
where
$P_{i}$
is the microscale quantity being averaged and
$W$
is a weighting function, which is typically either 1 or some grouping of microscale variables. If
$W$
is not explicitly noted in the specification of the averaging operator, then its value is unity. While the regions of integration have been denoted as domains in (2.14), these regions can also be specified as some subset of the boundary of
${\it\Omega}$
or some boundary of an entity within
${\it\Omega}$
, which will also be a lower-dimensional entity. For example, the domain of the wetting phase
${\it\Omega}_{w}$
has a boundary
${\it\Gamma}_{w}={\it\Gamma}_{we}\cup {\it\Gamma}_{wi}$
that consists of an external boundary where
${\it\Gamma}_{we}\subset {\it\Gamma}$
and an internal boundary where
${\it\Gamma}_{wi}={\it\Omega}_{ws}\cup {\it\Omega}_{wn}$
. When the region of integration can be specified with a domain of integration for an entity, the notation in (2.14) will be used. However, when it is necessary to integrate over an external boundary, this will be denoted explicitly.
Four types of averages will be used to denote macroscale measures of the state of the system, an intrinsic average
$$\begin{eqnarray}f^{{\it\alpha}}=\left\langle \,f_{{\it\alpha}}\right\rangle _{{\it\Omega}_{{\it\alpha}},{\it\Omega}_{{\it\alpha}}},\end{eqnarray}$$
a density-weighted average
$$\begin{eqnarray}f^{\overline{{\it\alpha}}}=\left\langle \,f_{{\it\alpha}}\right\rangle _{{\it\Omega}_{{\it\alpha}},{\it\Omega}_{{\it\alpha}},{\it\rho}_{{\it\alpha}}},\end{eqnarray}$$
and an average of a property of a higher-dimensional entity over a lower-dimensional entity (e.g. phase evaluated at an interface or common curve)
$$\begin{eqnarray}f_{{\it\alpha}}^{{\it\beta}}=\left\langle \,f_{{\it\alpha}}\right\rangle _{{\it\Omega}_{{\it\beta}},{\it\Omega}_{{\it\beta}}},\end{eqnarray}$$
where
$\text{dim}\,{\it\beta}<\text{dim}\,{\it\alpha}$
. Uniquely defined averages for some variables, which are denoted as a variable that is superscripted with a double-barred symbol, are explicitly defined.
Specific extent measure quantities such as the volume fractions, specific interfacial areas and specific common curve length are examples of specially defined averages, because the form does not fit any of the three standard defined averages. These measures can be formulated in terms of defined averages of the microscale indicator function defined as
$$\begin{eqnarray}{\it\Upsilon}_{{\it\alpha}}(\boldsymbol{x})=\left\{\begin{array}{@{}ll@{}}1, & \boldsymbol{x}\in {\it\Omega}_{{\it\alpha}},\\ 0, & \boldsymbol{x}\notin {\it\Omega}_{{\it\alpha}},\end{array}\right.\end{eqnarray}$$
where
$\boldsymbol{x}$
is a spatial position vector. Equation (2.18) can be used to formulate the specific entity measures as
$$\begin{eqnarray}{\it\epsilon}^{\overline{\overline{{\it\alpha}}}}=\langle {\it\Upsilon}_{{\it\alpha}}\rangle _{{\it\Omega},{\it\Omega}}\quad \text{for }{\it\alpha}\in \mathscr{I},\end{eqnarray}$$
where volume fractions correspond to the case in which
${\it\alpha}\in \mathscr{I}_{P}$
, specific interfacial areas correspond to the case in which
${\it\alpha}\in \mathscr{I}_{I}$
, and the specific common curve length corresponds to the case in which
${\it\alpha}\in \mathscr{I}_{C}$
.
The porosity can be formulated as
$$\begin{eqnarray}{\it\epsilon}^{\overline{\overline{}}}=\langle {\it\Upsilon}_{\mathit{f}}\rangle _{{\it\Omega},{\it\Omega}},\end{eqnarray}$$
where
$f$
is an index that refers to the pore space, which is occupied by a fluid phase. Fluid saturations are a common related set of fluid entity extent measures defined as
$$\begin{eqnarray}s^{\overline{\overline{{\it\alpha}}}}=\left\langle {\it\Upsilon}_{{\it\alpha}}\right\rangle _{{\it\Omega}_{\text{f}},{\it\Omega}_{\text{f}}}\quad \text{for }{\it\alpha}\in \mathscr{I}_{f}.\end{eqnarray}$$
Macroscale fluid pressures can be computed in a variety of ways that are relevant for multiphase flow in a porous medium system. An intrinsic average of the fluid pressure is expressed as
$$\begin{eqnarray}p_{}^{{\it\alpha}}=\left\langle p_{{\it\alpha}}\right\rangle _{{\it\Omega}_{{\it\alpha}},{\it\Omega}_{{\it\alpha}}}\quad \text{for }{\it\alpha}\in \mathscr{I}_{f},\end{eqnarray}$$
where
$p_{{\it\alpha}}$
is the microscale fluid pressure. Traditional laboratory and computational experiments to measure the state of a multiphase systems often control the fluid pressures and fluid flow at the boundaries of the domain, leading to a fluid pressure that is averaged over a segment of the boundary of the system given as
$$\begin{eqnarray}p_{{\it\alpha}}^{{\it\Gamma}_{i}\,}=\left\langle p_{{\it\alpha}}\right\rangle _{{\it\Gamma}_{i},{\it\Gamma}_{i}}\quad \text{for }{\it\alpha}\in \mathscr{I}_{f},\end{eqnarray}$$
where
${\it\Gamma}_{i}$
is a segment (e.g. one face) of the boundary of
${\it\Omega}$
. A third pressure of interest for two-fluid-phase systems is the interface-averaged pressure defined as
$$\begin{eqnarray}p_{{\it\alpha}}^{{\it\kappa}\,}=\left\langle p_{{\it\alpha}}\right\rangle _{{\it\Omega}_{{\it\kappa}},{\it\Omega}_{{\it\kappa}}}\quad \text{for }{\it\alpha}\in \mathscr{I}_{f},~{\it\kappa}\in \mathscr{I}_{c{\it\alpha}},\end{eqnarray}$$
where
$\mathscr{I}_{c{\it\alpha}}$
is the index set of entities connected to the
${\it\alpha}$
fluid within
${\it\Omega}$
, which must consist of an interface entity. For example, the
$w$
fluid meets the
$n$
fluid at the
$wn$
interface and the solid phase at the
$ws$
interface.
The macroscale pressures defined in (2.22)–(2.24) all arise as important measures in describing two-fluid-phase systems. In general, these pressures will not be equivalent. Thus care is needed in analysing the state of a two-fluid-phase system. It should also be noted that, in general, only the pressure defined by (2.23) is typically measured in traditional laboratory experiments, and this is often true even for state-of-the-art experiments that include high-resolution imaging. On the other hand, computational approaches provide a means to compute all of the defined pressures, yielding a basis to deduce a more complete understanding of the macroscale behaviour of the system than is possible when only the fluid pressures on the boundary of the domain can be controlled and observed.
Phase densities are computed as intrinsic averages by
$$\begin{eqnarray}{\it\rho}^{{\it\alpha}}=\left\langle {\it\rho}_{{\it\alpha}}\right\rangle _{{\it\Omega}_{{\it\alpha}},{\it\Omega}_{{\it\alpha}}}\quad \text{for }{\it\alpha}\in \mathscr{I}_{P},\end{eqnarray}$$
where the density is accessible through microscale computational results, which in turn will rely upon a microscale continuum equation of state.
Two types of velocity are of interest in evolving models of multiphase flow (e.g. Gray & Miller Reference Gray and Miller2014). The first is the density-weighted average velocity
$\boldsymbol{v}_{}^{\overline{{\it\alpha}}}$
for
${\it\alpha}\in \mathscr{I}$
. A second velocity of interest is the kinematic velocity of an interior boundary, which is the velocity of an interface or a common curve. The kinematic velocity is the velocity of the boundary of an entity, not the velocity of the material within the entity. To see the difference, consider steady-state co-current flow of two fluids through a pore with a rectangular cross-section. The wetting fluid will occupying the corners of the domain and a non-wetting fluid will occupy the central portion of the domain. Since the system is at steady state, the interface that forms the boundary between the two fluids is stationary, so the kinematic velocity vanishes. However, material may move through the pore without moving the interfaces, so the density-weighted interface velocity need not vanish. A boundary entity (interface or common curve) can have a velocity in the normal directions to the entity, and the boundary of the boundary entity itself can also change in the tangential direction. For example, the edge of an interface in a two-fluid-phase system is a common curve. This common curve can slide in a direction tangent to the interface.
We define a macroscale kinematic velocity of an interface as
$$\begin{eqnarray}\boldsymbol{w}^{\overline{\overline{{\it\alpha}}}}=\left\langle \boldsymbol{n}_{{\it\beta}}\boldsymbol{n}_{{\it\beta}}\,\boldsymbol{\cdot }\,\boldsymbol{v}_{{\it\alpha}}\right\rangle _{{\it\Omega}_{{\it\alpha}},{\it\Omega}_{{\it\alpha}}}=\left\langle \left(\unicode[STIX]{x1D644}-{\unicode[STIX]{x1D644}^{\prime }}_{{\it\alpha}}\right)\boldsymbol{\cdot }\boldsymbol{v}_{{\it\alpha}}\right\rangle _{{\it\Omega}_{{\it\alpha}},{\it\Omega}_{{\it\alpha}}}\quad \text{for }{\it\alpha}\in \mathscr{I}_{I},\end{eqnarray}$$
where
${\it\Omega}_{{\it\alpha}}=\bar{{\it\Omega}}_{{\it\beta}}\cap \bar{{\it\Omega}}_{{\it\gamma}}$
,
$\bar{{\it\Omega}}_{{\it\beta}}={\it\Omega}_{{\it\beta}}\cup {\it\Gamma}_{{\it\beta}}$
is a closed domain, which states that the
${\it\alpha}$
interface is defined as the internal boundary between the
${\it\beta}$
and
${\it\gamma}$
phases,
$\unicode[STIX]{x1D644}$
is the identity tensor, and
${\unicode[STIX]{x1D644}^{\prime }}_{{\it\alpha}}$
is the identity tensor in the surface. Thus
$\boldsymbol{w}^{\overline{\overline{{\it\alpha}}}}$
is the macroscale kinematic velocity of the
${\it\alpha}$
interface in the direction normal to the interface.
Similarly, the macroscale normal kinematic velocity of the common curve is
$$\begin{eqnarray}\boldsymbol{w}^{\overline{\overline{wns}}}=\left\langle \left(\unicode[STIX]{x1D644}-\boldsymbol{l}_{wns}\boldsymbol{l}_{wns}\right)\boldsymbol{\cdot }\boldsymbol{v}_{wns}\right\rangle _{{\it\Omega}_{wns},{\it\Omega}_{wns}}=\left\langle \left(\unicode[STIX]{x1D644}-{\unicode[STIX]{x1D644}^{\prime \prime }}_{wns}\right)\boldsymbol{\cdot }\boldsymbol{v}_{wns}\right\rangle _{{\it\Omega}_{wns},{\it\Omega}_{wns}},\end{eqnarray}$$
where
$\boldsymbol{l}_{wns}$
is the unit vector tangent to the common curve, and
${\unicode[STIX]{x1D644}^{\prime \prime }}_{wns}$
is the identity tensor in the common curve. Thus, the common curve can move in two directions that are normal to the curve. These directions of kinematic common curve movement are tangent and normal to the solid phase surface. If the solid phase is incompressible and fixed, the normal component of the kinematic velocity vanishes.
The curvature of the boundary of a phase
${\it\beta}$
is defined at the microscale as
$$\begin{eqnarray}J_{{\it\beta}}=\boldsymbol{{\rm\nabla}}^{\prime }\boldsymbol{\cdot }\boldsymbol{n}_{{\it\beta}}\quad \text{for }{\it\beta}\in \mathscr{I}_{P},~\boldsymbol{x}\in {\it\Gamma}_{{\it\beta}i},\end{eqnarray}$$
where
$\boldsymbol{{\rm\nabla}}^{\prime }\boldsymbol{\cdot }=(\unicode[STIX]{x1D644}-\boldsymbol{n}_{{\it\beta}}\boldsymbol{n}_{{\it\beta}})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }={\unicode[STIX]{x1D644}^{\prime }}_{{\it\alpha}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }$
is the microscale divergence operator restricted to an interface
${\it\alpha}$
. Note that the internal boundary of phase
${\it\beta}$
denoted
${\it\Gamma}_{{\it\beta}i}$
is formed at the transition between phases, which is an interface. Because interfaces are a subset of
$\Re ^{2}$
and do not have thickness, the curvature of the boundary of the phase is also the curvature of the interface that comprises the boundary. The macroscale interfacial curvature is specified by the microscale phase normal direction and the interface over which the average is computed, yielding the standard form
$$\begin{eqnarray}J_{{\it\beta}}^{\,{\it\alpha}}=\left\langle J_{{\it\beta}}\right\rangle _{{\it\Omega}_{{\it\alpha}},{\it\Omega}_{{\it\alpha}}}=-\left\langle J_{{\it\gamma}}\right\rangle _{{\it\Omega}_{{\it\alpha}},{\it\Omega}_{{\it\alpha}}}\quad \text{for }{\it\alpha}\in \mathscr{I}_{I},\end{eqnarray}$$
where
${\it\Omega}_{{\it\alpha}}=\bar{{\it\Omega}}_{{\it\beta}}\cap \bar{{\it\Omega}}_{{\it\gamma}}$
.
The capillary pressure of the
$wn$
interface is given at the microscale as
$$\begin{eqnarray}p_{wn}=-{\it\gamma}_{wn}J_{w},\end{eqnarray}$$
where
${\it\gamma}_{wn}$
is the interfacial tension of the
$wn$
interface. The capillary pressure is only approximately equal to the difference in the volume-averaged pressures of the fluids on each side of the interface when the system is in equilibrium. The macroscale capillary pressure is defined as
$$\begin{eqnarray}p^{wn}=-\left\langle {\it\gamma}_{wn}J_{w}\right\rangle _{{\it\Omega}_{wn},{\it\Omega}_{wn}},\end{eqnarray}$$
which, for the case of a constant interfacial tension at the microscale, allows for
$$\begin{eqnarray}p^{wn}=-{\it\gamma}^{wn}J_{w}^{wn}.\end{eqnarray}$$
For the case described by (2.32), a two-fluid-phase porous medium system at equilibrium has the property that
$$\begin{eqnarray}p^{wn}=p_{n}^{wn}-p_{w}^{wn}\approx p_{}^{n}-p_{}^{w}.\end{eqnarray}$$
The contact angle is a microscale common curve property, appearing under equilibrium conditions in both Young’s equation and the Young–Laplace equation. In a non-equilibrium system, the contact angle varies depending on the velocity of the common curve (Jiang, Oh & Slattery Reference Jiang, Oh and Slattery1979; Li & Slattery Reference Li and Slattery1991; Brochardwyart & DeGennes Reference Brochardwyart and DeGennes1992; Sagis & Slattery Reference Sagis and Slattery1995; Seppecher Reference Seppecher1996; Dhori & Slattery Reference Dhori and Slattery1997; Shikhmurzaev Reference Shikhmurzaev1997; Jacqmin Reference Jacqmin2000; Siebold et al.
Reference Siebold, Nardin, Schultz, Walliser and Oppliger2000; Pomeau Reference Pomeau2002; Slattery, Oh & Fu Reference Slattery, Oh and Fu2004; Blake Reference Blake2006). The microscale contact angle is expressed as
${\it\varphi}_{ws,wn}^{\,}$
and is the angle between the
$ws$
and the
$wn$
interfaces. A macroscopic average contact angle can be computed according to
$$\begin{eqnarray}\displaystyle \cos {\it\varphi}^{\overline{\overline{ws,wn}}}=\frac{\left\langle \cos {\it\varphi}_{ws,wn}^{\,}\right\rangle _{{\it\Omega}_{wns},{\it\Omega}_{wns}}}{\left(\left\langle \cos {\it\varphi}_{ws,wn}^{\,}\right\rangle _{{\it\Omega}_{wns},{\it\Omega}_{wns}}^{2}+\left\langle \sin {\it\varphi}_{ws,wn}^{\,}\right\rangle _{{\it\Omega}_{wns},{\it\Omega}_{wns}}^{2}\right)^{1/2}}, & & \displaystyle\end{eqnarray}$$
where
$\cos {\it\varphi}^{\overline{\overline{ws,wn}}}$
is defined to ensure that the identity
$$\begin{eqnarray}\cos ^{2}{\it\varphi}^{\overline{\overline{ws,wn}}}+\sin ^{2}{\it\varphi}^{\overline{\overline{ws,wn}}}=1\end{eqnarray}$$
is satisfied.
The state properties defined above require averaging of microscale quantities, which in turn requires a representation of the microscale state variables. Both high-resolution experimental and computational approaches can yield some of the state variables. Only computational approaches are currently able to yield all microscale state variables. Our LBM simulator is coupled to an in situ data analysis framework that is used to continuously analyse the simulation state as it progresses such that the temporal behaviour of averages can be studied at high resolution (McClure et al. Reference McClure, Wang, Prins, Miller and Feng2014b ). In this approach, the phase volumes, interfaces and the common curve are determined from image processing (McClure et al. Reference McClure, Adalsteinsson, Pan, Gray and Miller2007). The interface and common curve entities were constructed, and averages were determined by applying low-order interpolation and integration to approximate the averages listed in this section.
3 Results and discussion
3.1 Validation
Validation cases were designed to verify that the simulator and analysis framework recovers correct averages for both static and dynamic processes. Three cases are included: (1) a bubble test, which demonstrates that the interfacial tension is modelled correctly by the LBM; (2) simulation of a bubble trapped in a cylindrical capillary tube, which demonstrates that various equilibrium contact angles can be recovered; and (3) piston displacement in a cylindrical capillary tube, which demonstrates that the kinematic velocities of the interface and common curve are consistent with other flow measures and that the dynamic contact angle matches well-known behaviour based on our implementation of the LBM. Quantities are plotted in non-dimensional form when appropriate.
The bubble test is a standard validation for multiphase LBMs (Shan & Chen Reference Shan and Chen1994; Martys & Chen Reference Martys and Chen1996; Chen & Doolen Reference Chen and Doolen1998; Lee & Fischer Reference Lee and Fischer2006; Shan & Chen Reference Shan and Chen2007). It shows that the measured difference in the equilibrium phase pressures at the microscale is equivalent to the product of the interfacial tension and the curvature as given by Laplace’s law,
$$\begin{eqnarray}p_{n}-p_{w}=-{\it\gamma}_{wn}J_{w}.\end{eqnarray}$$
In the bubble test, a sequence of bubbles of one fluid are immersed in a second fluid. Simulations were performed using a lattice length
${\it\delta}x=2~{\rm\mu}\text{m}$
, dynamic viscosities
${\it\mu}_{w}={\it\mu}_{n}={\it\mu}=1\times 10^{-3}~\text{Pa}~\text{s}$
and densities
${\it\rho}_{w}={\it\rho}_{n}=1.0~\text{g}~\text{cm}^{-3}$
. The radius of curvature is determined at equilibrium based on the volume of each immersed bubble. Plotting the pressure difference as a function of the radius of curvature for several bubble sizes provides the interfacial tension as the slope of the line. Here, we compute the interfacial curvature directly from the microscale, so the test can also be used to confirm that the average phase pressure and interfacial curvature are consistent with the value of
${\it\gamma}^{wn}$
. At equilibrium, the macroscale phase pressures, interfacial tension and interfacial curvature are constant, which means that
$$\begin{eqnarray}p_{}^{n}-p_{}^{w}=-{\it\gamma}^{wn}J_{w}^{\,wn}.\end{eqnarray}$$
The microscale interfacial curvature is averaged over the fluid–fluid interface to obtain the macroscopic curvature
$J_{w}^{\,wn}$
. A range of bubble sizes were instantiated for five different interfacial tensions
${\it\gamma}^{wn}$
, with the results for each bubble test plotted in figure 1.
Figure 1. Bubble tests show the equilibrium relationship between the average phase pressure difference
$p_{}^{n}-p_{}^{w}$
, the curvature
$J_{w}^{\,wn}$
and the interfacial tension
${\it\gamma}^{wn}$
.
Similar to the bubble test, a trapped bubble of non-wetting phase within a cylindrical tube provides a way to numerically measure the equilibrium contact angle that results from a particular choice of simulation parameters (Huang et al.
Reference Huang, Thorne, Schaap and Sukop2007; McClure et al.
Reference McClure, Prins and Miller2014a
). Based on the equilibrium configuration, averages were computed for
$p_{}^{w}$
,
$p_{}^{n}$
,
${\it\varphi}^{\overline{\overline{ws,wn}}}$
and
$J_{w}^{\,wn}$
.
In a tube with radius
$R$
, the capillary pressure can be analytically determined based on the equilibrium contact angle
${\it\varphi}_{ws,wn}^{\,}$
according to
$$\begin{eqnarray}p_{n}-p_{w}=\frac{2{\it\gamma}_{wn}\cos {\it\varphi}_{ws,wn}^{\,}}{R}.\end{eqnarray}$$
In the absence of external forces and at equilibrium, the pressures
$p_{n}$
and
$p_{w}$
are constant, and the contact angle is invariant along the common curve. As a result, for an equilibrium configuration, it must be true that
$p_{}^{n}=p_{n}$
,
$p_{}^{w}=p_{w}$
and
${\it\varphi}^{\overline{\overline{ws,wn}}}={\it\varphi}_{ws,wn}^{\,}$
. The interfacial tension
${\it\gamma}^{wn}$
is independently determined based on the bubble test as described in the previous section. Inserting the averages into (3.3), we conclude that, at thermodynamic and mechanical equilibrium, the contact angle should satisfy
$$\begin{eqnarray}\cos {\it\varphi}^{\overline{\overline{ws,wn}}}=\frac{(p_{}^{n}-p_{}^{w})R}{2{\it\gamma}^{wn}}=\frac{p^{wn}R}{2{\it\gamma}^{wn}}.\end{eqnarray}$$
The contact angle was computed geometrically as described in the previous section. At equilibrium, the value should agree with the right-hand side of (3.4) and with the interfacial curvature. These three quantities are plotted in figure 2 based on simulations performed by inserting a bubble of non-wetting phase into a cylindrical tube with radius
$R=18{\it\delta}x$
. To change the contact angle, the value of the phase density field within the solid phase was chosen as
${\it\phi}_{s}=0.0$
, 0.15, 0.3, 0.45, 0.6, 0.75 and 0.9. The system was then allowed to equilibrate for an elapsed time of
$1\times 10^{5}\,{\it\delta}t$
, obtaining a configuration at mechanical equilibrium. Reasonable agreement among the three values is observed, with deviations expected due to the smeared interfacial regions for the LBM. The measured interfacial curvature and pressure can deviate from the value implied by the contact angle as a consequence. The three independently measured quantities are plotted in figure 2.
Figure 2. Agreement between the geometrically measured contact angle value and the equilibrium interfacial curvature and phase pressures plotted as a function of the solid boundary value
${\it\phi}_{s}$
as determined using a trapped bubble test in a three-dimensional capillary tube with width
$R=18{\it\delta}x$
.
An advantage of in situ analysis is the capability to efficiently extract information about the system dynamics at high temporal resolution. Piston displacement in a cylindrical capillary tube is considered as a demonstrative example and as a validation case. Reservoirs of wetting and non-wetting fluids were established at opposite ends of a capillary tube of length
$L=400{\it\delta}x$
. The total domain size was
$40\times 40\times 400$
, with a capillary tube radius of
$R=18{\it\delta}x$
to match the prior case of a static trapped bubble. Initially the capillary tube was saturated with the wetting fluid, and a pressure boundary condition was applied to induce flow of non-wetting phase into the tube. Pressure values for the inlet and outlet were set such that the pressure difference exceeds the equilibrium values measured for the static case. The geometry is shown in figure 3. Equilibrium measurements are known from figure 2, since the capillary tube diameter was identical for both static and dynamic simulations.
Figure 3. Piston displacement driven by pressure boundary conditions in a cylindrical tube. The radius of the tube was
$R=18{\it\delta}x$
and the length of the tube was
$L=400{\it\delta}x$
.
Simulations of displacement were performed using three different equilibrium contact angles corresponding to
${\it\phi}_{s}=0.0$
, 0.3 and 0.6. Different capillary numbers were simulated by varying the boundary pressure difference and the interfacial tension. The capillary number was defined as
$$\begin{eqnarray}\mathit{Ca}=\frac{{\it\mu}U}{{\it\gamma}_{wn}},\end{eqnarray}$$
where
$U=|\boldsymbol{w}^{\overline{\overline{wns}}}|$
. This is a non-standard definition for the capillary number based on the velocity of the common curve. More typically, one of the phase velocities is used in the definition of the capillary number. However, for a steady displacement in a capillary tube,
$v_{z}^{\overline{w}}=v_{z}^{\overline{n}}=w_{z}^{\overline{\overline{wns}}}=U$
. This equivalence will not hold when the displacement is not steady, such as when the configuration of the
$wn$
interface is changing.
Each of the averages defined in § 2.4 was computed every 1000 time steps throughout the course of the simulation. The rate of change in saturation was determined by fitting a smoothing spline to the saturation values measured from the simulation. At steady state, the change in
$s^{\overline{\overline{w}}}$
is linear in time due to the geometry, with the slope of the line given by
$$\begin{eqnarray}\frac{\partial s^{\overline{\overline{w}}}}{\partial t}=\frac{U}{L}.\end{eqnarray}$$
Partial derivatives with respect to time were computed by fitting a third-order polynomial spline function to the time history of the variable of interest. An analytical relationship between the rate of change in saturation and the rate of the deformation of
${\it\Omega}_{wn}$
can be obtained by applying the transport theorem for a phase
$$\begin{eqnarray}{\it\epsilon}\frac{\partial s^{\overline{\overline{w}}}}{\partial t}=\frac{1}{V}\int _{{\it\Omega}_{wn}}\boldsymbol{w}_{wn}\boldsymbol{\cdot }\boldsymbol{n}_{w}\,\text{d}\mathfrak{r},\end{eqnarray}$$
where
$V$
is the total volume of the domain. This expression is exact provided that the solid is immobile and non-deformable (Gray et al.
Reference Gray, Dye, McClure, Pyrak-Nolte and Miller2015). In figure 4, the rate of change in saturation is plotted against scaled measurements of the phase velocities
$v_{z}^{\overline{w}}$
and
$v_{z}^{\overline{n}}$
, the kinematic velocity of the common curve
$w_{z}^{\overline{\overline{wns}}}$
and (3.7). Good agreement is obtained among the four independently measured values shown. Of the four measures, errors in the common curve velocity
$\boldsymbol{w}^{\overline{\overline{wns}}}$
are largest. These results demonstrate that our approach to compute interface and common curve velocities delivers results that are consistent with other measures in the steady displacement.
Figure 4. Agreement between measured rate of saturation change, average velocity of the phases and kinematic velocity of the common curve.
Dynamic measurements of the contact angle were determined for each steady-state displacement. The dynamic behaviour of the contact angle in a two-fluid system has a long history of study, with the contact angle changing from its equilibrium value for a moving common curve (e.g. Jiang et al. Reference Jiang, Oh and Slattery1979; Li & Slattery Reference Li and Slattery1991; Brochardwyart & DeGennes Reference Brochardwyart and DeGennes1992; Sagis & Slattery Reference Sagis and Slattery1995; Seppecher Reference Seppecher1996; Dhori & Slattery Reference Dhori and Slattery1997; Shikhmurzaev Reference Shikhmurzaev1997; Jacqmin Reference Jacqmin2000; Siebold et al. Reference Siebold, Nardin, Schultz, Walliser and Oppliger2000; Pomeau Reference Pomeau2002; Slattery et al. Reference Slattery, Oh and Fu2004; Blake Reference Blake2006). The multiphase LBM has previously been shown to recover established models for the dynamic contact angle without a need to explicitly set the slip length or other parameters. Instead, the dynamic behaviour arises naturally due to uncompensated stresses in the vicinity of the common curve (Latva-Kokko & Rothman Reference Latva-Kokko and Rothman2007). A linear dependence on the capillary number was proposed by Sheng & Zhou (Reference Sheng and Zhou1992):
$$\begin{eqnarray}\cos {\it\varphi}^{\overline{\overline{ws,wn}}}-\cos {\it\varphi}_{eq}^{\overline{\overline{wn,ws}}}=\mathit{Ca}\log \left(KD/l_{s}\right),\end{eqnarray}$$
where
$K$
is a constant and
$l_{s}/D$
is the slip length. Equation (3.8) can be shown to approximate the well-known slip model derived by Cox (Reference Cox1986). In figure 5, we show that our implementation recovers (3.8); the dynamic contact angle is a linear function of the capillary number, with the same slope obtained for different equilibrium contact angles based on averaged measurements obtained from simulations of piston displacement.
Figure 5. Dynamic contact angle as a function of the capillary number plotted for three different equilibrium contact angles.
The approach used to compute the microscale kinematic velocity of the interface and common curve is particularly useful for dealing with issues associated with parasitic currents in the LBM. Owing to parasitic currents, the phase velocity tends to be inaccurate in the vicinity of the interfaces and common curve. The normal velocity of the interface can be computed accurately according to (2.6). Furthermore, it can be used to determine the velocity of the common curve even though a no-slip boundary condition is applied in the momentum transport solution.
3.2 Porous medium system
To evaluate the accuracy of the simulations in a realistic porous medium geometry, a resolution study was performed using a periodic pack of 229 equally sized spheres in a cubic domain. Simulations were performed at four different resolutions, which corresponded to cubic domain sizes of
$64^{3},~128^{3},~256^{3}$
and
$512^{3}$
. For these simulations, the sphere diameters were
$D=13.6{\it\delta}x$
,
$27.2{\it\delta}x$
,
$54.5{\it\delta}x$
and
$109{\it\delta}x$
, respectively. Fast drainage and imbibition sequences were carried out following Ahrenholz et al. (Reference Ahrenholz, Toelke, Lehmann, Peters, Kaestner, Krafczyk and Durner2008). The system was initially saturated with a wetting fluid phase, and reservoirs of wetting and non-wetting fluids were imposed at opposing ends of the domain. Pressure boundary conditions were then chosen to drive fluid displacement. The simulation time was set as
$$\begin{eqnarray}t_{max}=\frac{{\it\mu}L}{{\it\gamma}^{wn}}t_{ref},\end{eqnarray}$$
where the non-dimensional quantity
$t_{ref}=4075.312$
was chosen to allow sufficient time for the fluid displacement to occur. The boundary pressures were then chosen to be linear functions of time,
$$\begin{eqnarray}p_{{\it\alpha}}^{{\it\Gamma}}(t)=p_{{\it\alpha}}^{{\it\Gamma}}(0)+\frac{p_{{\it\alpha}}^{{\it\Gamma}}(t_{max})-p_{{\it\alpha}}^{{\it\Gamma}}(0)}{t_{max}}\quad \text{for }{\it\alpha}\in \{w,n\}.\end{eqnarray}$$
The initial and final values for the pressure were specified for each of the two reservoirs. Pressure boundary conditions were chosen to ensure identical non-dimensional behaviour for each the four resolutions. For drainage, these boundary pressure values were
$$\begin{eqnarray}\displaystyle & \displaystyle \left[p_{n}^{{\it\Gamma}}(0)-p_{w}^{{\it\Gamma}}(0)\right]D/{\it\gamma}^{wn}=4.7, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left[p_{n}^{{\it\Gamma}}(t_{max})-p_{w}^{{\it\Gamma}}(t_{max})\right]D/{\it\gamma}^{wn}=28.12, & \displaystyle\end{eqnarray}$$
and for imbibition
$$\begin{eqnarray}\displaystyle & \displaystyle \left[p_{n}^{{\it\Gamma}}(0)-p_{w}^{{\it\Gamma}}(0)\right]D/{\it\gamma}^{wn}=12.53, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left[p_{n}^{{\it\Gamma}}(t_{max})-p_{w}^{{\it\Gamma}}(t_{max})\right]D/{\it\gamma}^{wn}=1.6. & \displaystyle\end{eqnarray}$$
The number of time steps required for each process (drainage and imbibition) was then determined for each resolution as
$3\times 10^{5}{\it\delta}t$
,
$6\times 10^{5}{\it\delta}t$
,
$1.2\times 10^{6}{\it\delta}t$
and
$2.4\times 10^{6}{\it\delta}t$
, respectively.
Figure 6. Behaviour of three average quantities as a function of fluid saturation measured during simulated drainage and imbibition performed using four different resolutions: (a) boundary pressure difference; (b) phase pressure difference.
In a traditional experimental set-up, fluid pressures are only measured at the boundary reservoirs. In our simulations, this value was set explicitly; the resulting saturation values are shown in figure 6(a). While the curves indicate a well-behaved trend towards convergence as the resolution increases, the result cannot be considered fully grid-independent even at a resolution
$D=109{\it\delta}x$
. The largest differences are observed during drainage for irreducible wetting phase saturation. This is due to the fact that wetting phase tends to be trapped in very small spaces in the vicinity of grain contacts. Resolving these features requires high resolution, and, as more pendular rings are resolved, the observed irreducible wetting phase saturation increases. Resolving pendular rings does not appear to have a large effect on the subsequent dynamics of imbibition, however. Trapped non-wetting phase, which is obtained after imbibition, tends to occupy large pores that are more easily resolved. Predicted values for the entrapped non-wetting phase agree quite well for
$D=54.5{\it\delta}x$
and
$D=109{\it\delta}x$
, with final residual non-wetting phase occupying 10.1 % and 10.3 % of the total pore volume at these two resolutions.
Macroscopic theory is expressed in terms of the averaged phase pressures
$p_{}^{w}$
and
$p_{}^{n}$
. These values differ from the boundary measurements because internal pressure values vary throughout the domain, and the pressures of features that are disconnected from the boundary need not have the same pressure as the boundary values. The measured difference between the average phase pressures is shown in figure 6(b). Close agreement is observed between results obtained for
$D=27.2{\it\delta}x$
,
$54.5{\it\delta}x$
and
$109{\it\delta}x$
. Since the phase-averaged pressures reflect the dynamics of the system, fluctuations are observed as a result of Haine’s jumps during drainage and snap-off during imbibition. These dynamics are similar for the three most well-resolved cases.
Figure 7. Consistency of the rate of change in saturation and the rate of boundary deformation as a function of resolution for multiphase displacement in a sphere pack: (a) agreement with (3.7); (b) box plot showing the error distribution.
The rate of change in saturation is plotted against the rate of boundary deformation for all four resolutions in figure 7. For a sharp interface, in the absence of numerical error, the two terms will be equal. However, since the interfacial region simulated by the LBM is diffuse, the location of the interface must be determined approximately within the interfacial region. The errors associated with (3.7) therefore decrease as the width of the interfacial region decreases relative to the overall resolution. Figure 7(a) shows consistency between the right- and left-hand sides of (3.7). Agreement improves as the resolution increases. A box plot for the associated error distribution is shown in figure 7(b). The variances for error obtained at the four resolutions were
$8.16\times 10^{-6}$
,
$4.96\times 10^{-6}$
,
$1.56\times 10^{-6}$
and
$7.78\times 10^{-7}$
, respectively, for
$D=13.6{\it\delta}x$
,
$26.2{\it\delta}x$
,
$54.5{\it\delta}x$
and
$109{\it\delta}x$
. In each case, we find that doubling the resolution decreased the variance in the error by slightly more than a factor of two.
4 Conclusions
An approach was implemented to evaluate the microscale kinematic velocity of the fluid–fluid interface and the common curve, which was shown to accurately determine these quantities for simulations where parasitic currents were present. The analysis approach used to extract kinematic information is formulated in terms of the level-set equation, and may also be applied to analyse the kinematics of experimentally generated data. In this work, a multiphase lattice-Boltzmann simulator was used to demonstrate the analysis capabilities for dynamic systems and validate the approach. A range of equilibrium contact angles were considered, and steady displacement in a cylindrical capillary tube was used to verify that the dynamic behaviour of the contact curve could be determined accurately. Under dynamic conditions, the contact angle was shown to be consistent with widely used slip models without the need for additional parameters. These results were obtained by using the common curve velocity directly.
We note that, while higher-order methods are able to reduce the effects of spurious currents, numerics that rely on a diffuse interface are still inherently limited by the errors associated with the interface representation. As with many other numerical methods, LBMs represent the interfacial region by smoothing it across several grid points so that derivatives can be computed. This is not an ideal approximation for real immiscible fluids; a jump condition would provide the best approximation. The accuracy of the coupled approach is limited by interfacial approximations made in the LBM rather than the quadrature scheme used to determine the averages. The thickness of the interfacial region in the LBM results in boundary effects at the solid surface. The effect is most significant for small contact angles.
A resolution study was performed for fast drainage and imbibition sequences in a pack of 229 equally sized spheres. Four resolutions were considered in which each sphere diameter was resolved with 13.6, 26.2, 54.5 and 109 lattice lengths, respectively. The ability to predict the distribution of the entrapped non-wetting phase that forms during imbibition was achieved at a resolution of
$54.5$
lattice lengths. For displacement in the complex geometry, transport theorem results were used to determine the accuracy of the interfacial kinematics. We find that high resolution is essential to recover this behaviour accurately. Convergence was not obtained for the interfacial kinematics until a resolution of
$109$
lattice lengths per sphere diameter was used; this would be considered a very well-resolved simulation based on previously published work.
Acknowledgement
This work was supported by National Science Foundation Grant 0941235, Department of Energy Grant DE-SC0002163, and Army Research Office Grant W911NF-14-1-0287. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research also used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725.
