Macroscopic model for generalised Newtonian inertial two-phase flow in porous media

Abstract

A closed macroscopic model for quasi-steady, inertial, incompressible, two-phase generalised Newtonian flow in rigid and homogeneous porous media is formally derived. The model consists of macroscopic equations for mass and momentum balance as well as an expression for the macroscopic pressure difference between the two fluid phases. The model is obtained by upscaling the pore-scale equations, employing a methodology based on volume averaging, the adjoint method and Green's formulation, only assuming the existence of a representative elementary volume and the separation of scales between the microscale and the macroscale. The average mass equations coincide with those for Newtonian flow. The macroscopic momentum balance equation in each phase expresses the seepage velocity in terms of a dominant and a coupling Darcy-like term, a contribution from interfacial tension effects and another one from interfacial inertia. Finally, the expression of the macroscopic pressure difference is obtained in terms of the macroscopic pressure gradient and body force in each phase, and interfacial terms that account for capillary effects and inertia, if present when the interface is not stationary. All terms involved in the macroscale equations are predicted from the solution of adjoint closure problems in periodic representative domains. Numerical predictions from the upscaled models are compared with direct numerical simulations for two-dimensional configurations, considering flow of a Newtonian non-wetting fluid and a Carreau wetting fluid. Excellent agreement between the two approaches confirms the pertinence of the macroscopic models derived here.

1. Introduction

Two-phase non-Newtonian flow in porous media is involved in numerous applications such as the production of scaffolds for biomedical applications (Elbert 2011), remediation of polluted soils (Philippe et al. 2020), in trickle bed reactors (Giri & Majumder 2014) and packed bed bioreactors (Sen, Nath & Bhattacharjee 2017), gas and liquid flow through porous shale reservoirs (Singh & Cai 2019), polymer injection for composite materials manufacturing (Wittemann et al. 2019), food engineering (Liu et al. 2019), gas–liquid and liquid–liquid flow in packed columns for extraction processes and distillation (Chhabra & Basavaraj 2019) and soil flooding for enhanced oil recovery (Sorbie 1991; Nilsson et al. 2013), to cite only a few. Despite this extensive interest, flow modelling at the macroscopic scale relying on a physically sound basis remains challenging and, most of the time, empirical approaches are employed. This undoubtedly stems from the nonlinear character of the fluid rheology, combined with the complexity inherent to fluid-phase dynamics and the porous medium description.

Triggered by the presence of non-Newtonian flows in enhanced oil recovery, the earliest reported analyses were made essentially at the macroscopic scale, employing postulated extensions of Darcy's law. In this regard, two main types of heuristic models have been proposed in the literature. The first one relies on an extension of that employed for one-phase flow of a power-law fluid (Pascal 1983). It consists of a modified version of Darcy's law relating the seepage velocity, elevated at a power dependent on the fluid rheology, to the macroscopic pressure gradient (eventually shifted by a threshold value for a yield-stress fluid) making use of an effective dynamic viscosity. The second type of model uses the same generalised Darcy's law as for Newtonian fluids, albeit the nonlinearity is lumped into the viscosity coefficients.

During the 1980s, Pascal reported a series of works concerning non-Newtonian two-phase flow in porous media, making use of the first type of macroscopic model. In particular, the frontal displacement of a Newtonian fluid by a power-law yield-stress fluid, both being incompressible, was studied in the work of Pascal (1984) and further extended to compressible and non-Newtonian conditions for both phases by Pascal & Pascal (1988). In these works, where the displaced fluid was assumed motionless behind the front, the question of the non-Newtonian effect on the potential instability of the front was raised (further addressed by Pascal (1986)), and the existence of a pressure front, ahead of the saturation front in the compressible case, was identified. Its modification when both phases are mobile behind the saturation front was also analysed shortly after (Pascal 1990). In all these studies, capillary effects were neglected, in accordance with the frontal displacement approach.

The same type of model (i.e. generalised Darcy's law with a power-law relationship for the (weakly) non-Newtonian phase) was considered in a work dedicated to the analytical and numerical solution of simultaneous incompressible two-phase flow of a power-law wetting fluid and a Newtonian non-wetting fluid, including capillary pressure effects (Xi & Shangping 1987). Using power-series expansions and a perturbation method, an analytical solution was obtained. It was shown that, when the non-Newtonian wetting phase is dilatant, an increment in the fractional fluid flow is observed, whereas the opposite occurs when this phase is pseudoplastic.

The second type of heuristic model was used by Wu, Pruess & Witherspoon (1991) to develop analytical solutions for one-dimensional displacement of a Newtonian fluid by a non-Newtonian one in porous media. Using a Buckley–Leverett approach (neglecting the capillary pressure gradient), they showed how the saturation profile and the front displacement were affected by the rheological parameters. Later on, Wu & Pruess (1998) extended their analysis to three dimensions, using a finite difference scheme. The same type of model was recently considered by Katiyar et al. (2020).

A somewhat general tendency in the study of two-phase non-Newtonian flow since the beginning of this century, has been the use of pore-scale approaches, favoured by the development of microfluidics, computational means and imaging techniques. In this regard, Tsakiroglou (2004) performed experiments of a shear-thinning non-wetting phase displacing a Newtonian wetting fluid in glass-etched quasi-two-dimensional porous media. By performing inverse modelling, the relative permeabilities and capillary pressure curves were identified, on the basis of the second heuristic model mentioned above. With this approach, a nonlinear dependency of the non-wetting phase relative permeability upon the capillary number was reported.

The analysis of viscous fingering resulting from Saffman–Taylor instability, mainly motivated by the study of sweep efficiency in oil recovery processes, has been the focus of several works, also following a pore-scale approach. The displacement of a non-Newtonian (Bingham) fluid by a Newtonian one, using two-dimensional pore-network modelling in a square array of tubes having uniformly distributed radii, was investigated by Tian & Yao (1999). They highlighted the transition from the diffusion-limited aggregation regime to a compact regime. More recently, Shi & Tang (2016) reported lattice Boltzmann pore-scale simulations incorporating a phase-field method to study viscous fingering during displacement of a non-Newtonian fluid by a Newtonian one in a staggered two-dimensional pattern of square inclusions. Their results showed that displacement becomes unstable as the capillary number, ratio of displaced to displacing fluids viscosity and contact angle increase. In a later work, Wang et al. (2019) considered a shear-thinning (power-law) fluid displaced by a Newtonian one, using lattice Boltzmann simulations in a two-dimensional pattern of staggered circular solid inclusions. They showed that, as the power-law index is reduced (i.e. increasing the shear-thinning character of the fluid), viscous fingering is stabilised resulting in thicker fingers and better displacement efficiency. In addition, the porous medium microstructure has been reported to play a decisive role in the fingering phenomenon, as evidenced by Shende, Niasar & Babaei (2021). Using also the lattice Boltzmann method, Xie, Lv & Wang (2018) investigated the cases of shear-thinning and shear-thickening displacing fluids at the microscale in both homogeneous and heterogeneous porous structures. In the homogeneous case, those authors reported that the non-Newtonian character of the flow did not result in a noticeable difference in the favourable displacement regime. However, in the unfavourable displacement regime, shear-thinning fluids did not enhance diversion effects (i.e. the capability of favouring flow of the displacing fluid in lower-permeability regions with respect to the high-permeability zones). This counter-intuitive result was attributed to the viscosity ratio and porous medium properties.

Although pore-scale approaches are of great interest in identifying and deciphering complex conjugated mechanisms, their use is limited in terms of the size of the medium that can be investigated. In addition, they are not dedicated to the derivation of general models at the macroscopic (Darcy) scale. Moreover, surface tension effects may play a relevant role in the prediction of the seepage velocities as suggested by the numerical simulations reported by Druetta & Picchioni (2020) and Wang et al. (2020), and recently formally identified for two-phase Newtonian flow (Lasseux & Valdés-Parada 2022). However, these effects are not present in these commonly employed macroscopic models, and despite the fact that their use is an appealing approach that may ease flow description, this is a serious limitation. Indeed, even if inverse modelling is used to predict the coefficients present in these models from experimental data, their origin remains hypothesised. Although these issues were suggested in Pearson & Tardy (2002), formal models for macroscopic two-phase non-Newtonian flow in porous media are still lacking.

While many of the works cited in the previous paragraphs deal with displacement fronts and fingering effects, where there may not be separation of length scales between the microscale and the macroscale, it is nevertheless relevant to derive macroscopic models in regions sufficiently far away from the front displacement zone (if present). Furthermore, even if upscaled models were available in this zone, they should reproduce the results from models applicable sufficiently far away from their influence. In many other circumstances, two-phase flow in porous media may occur while the two fluids coexist everywhere in the pore space, as for instance in two fluid-phase reactors (Safinski & Adesina 2012). Hence, the objective of this work is to derive upscaled models in porous media regions where it is acceptable to assume separation of length scales between the microscale and the macroscale, and this is detailed in the following sections. In § 2, the governing equations at the pore scale are described, along with the associated starting assumptions. Section 3 is dedicated to the derivation of the macroscopic mass and momentum balance equations. In addition, in § 4, the expression of the macroscopic pressure difference is reported. These macroscopic models are obtained following the approach recently employed by Lasseux & Valdés-Parada (2022, 2023), which makes use of a simplified version of the volume-averaging technique along with an adjoint method and an integral (Green's) formulation. The average models are summarised in § 5 and they are evaluated and validated by means of direct numerical simulations (DNS) in § 6. In this section, an analysis of the effects of the fluid rheology, inertia and phase distribution over the macroscopic variables and effective-medium coefficients is provided. Conclusions are finally reported in § 7.

2. Pore-scale model

2.1. Starting assumptions and system configuration

The incompressible, inertial and quasi-steady flow under consideration is that of two non-Newtonian fluid phases, $\beta$ and $\gamma$, through a rigid and homogeneous porous medium (for which the solid skeleton is denoted as the $\sigma$ phase). By quasi-steady it is meant that the shapes and positions of the fluid–fluid interfaces within the pores are not stationary, although momentum transfer in the bulk of each phase is considered steady. Furthermore, the derivations that follow are applicable in situations in which both fluid phases circulate in either a continuous or discontinuous manner. Under these conditions, the mass and momentum balance equations at the pore scale are written as (in the following, $\alpha =\beta,\gamma$ is used to represent either one of the fluid phases)

(2.1a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}_\alpha=0, \quad\text{in }\mathscr{V}_\alpha, \end{gather}$$

(2.1b)$$\begin{gather}\rho_\alpha\boldsymbol{v}_\alpha\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{v}_\alpha=\rho_\alpha\boldsymbol{b}_\alpha + \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{p_\alpha} \quad\text{in }\mathscr{V}_\alpha, \end{gather}$$

where ${\boldsymbol{\mathsf{T}}}_{p_\alpha }$ is the total stress tensor in the $\alpha$ phase, defined by

(2.1c)\begin{equation} {\boldsymbol{\mathsf{T}}}_{p_\alpha}={-}{\boldsymbol{\mathsf{I}}}p_\alpha+\mu(\varGamma_\alpha) ( \boldsymbol{\nabla} \boldsymbol{v}_\alpha+\boldsymbol{\nabla} \boldsymbol{v}^T_\alpha ). \end{equation}

In these equations, $\boldsymbol {v}_\alpha$, $\rho _\alpha$ and $p_\alpha$ are, respectively, the velocity, density and pressure of the $\alpha$ phase, whereas ${\boldsymbol{\mathsf{I}}}$ is the identity tensor and $\boldsymbol {b}_\alpha$ represents the body force per unit mass, which is assumed to be constant within $\mathscr {V}_\alpha$. Moreover, in the viscous part of the total stress tensor, a generalised Newton's law of viscosity is used, which means that the apparent viscosity, $\mu (\varGamma _\alpha )$, is a function of the shear-rate modulus, $\varGamma _\alpha$, defined as

(2.1d)\begin{equation} \varGamma_\alpha =\sqrt{ \tfrac{1}{2} [ ( \boldsymbol{\nabla} \boldsymbol{v}_\alpha+ \boldsymbol{\nabla} \boldsymbol{v}^T_\alpha ): ( \boldsymbol{\nabla} \boldsymbol{v}_\alpha+ \boldsymbol{\nabla} \boldsymbol{v}^T_\alpha )]}. \end{equation}

While writing (2.1b), it is assumed that, albeit convective acceleration may be significant in accordance with a Reynolds number not exceedingly small compared with unity, the frequency parameter in the $\alpha$ phase, defined as $\rho _\alpha \ell _{{ref} \alpha }^2/(\mu _{{ref} \alpha } t_{ref})$, remains small compared with 1 ($\ell _{{ref} \alpha }$, $\mu _{{ref} \alpha }$ and $t_{ref}$ being, respectively, a reference length and viscosity in the $\alpha$ phase and a reference time; these definitions are further discussed in § 6 once a particular rheological model is chosen).

2.2. Flow problem in a periodic unit cell

Following a classical assumption in upscaling, the analysis is restricted to systems in which there is a separation of length scales. This means that the following developments correspond to two-phase flow far away from a displacement front, if present, or within the entire system if the two phases coexist mixed as exemplified in figure 1. In the former case, it is assumed that the porous medium (of characteristic length $\mathscr {L}$) can be subdivided into regions (of characteristic length $L<\mathscr {L}$) that have a hierarchical nature. In these circumstances, an averaging domain, $\mathscr {V}$ (of measure $V$ and characteristic size $r_0\ll L$), representative of the pore-scale configuration and transport processes (at least for momentum as discussed below) can be assumed to exist. For this reason, $r_0$ must also be much larger than the characteristic pore-scale length, $\max (\ell _\beta,\ell _\gamma )$. This concept, which is introduced in the form of a periodic unit cell at the beginning in the homogenisation technique (Auriault, Boutin & Geindreau 2009) or at the closure level in the volume-averaging method (Whitaker 1999), is necessary for a local description. At this point, it is pertinent to mention that the periodicity assumption is a commodity more than a necessity in the upscaling approach used here. It allows the closure problems to be solved only in a unit cell, thus leading to local effective-medium coefficients. Furthermore, the upscaled models can be used in practice even in situations when the real porous medium is not periodic. In other words, periodicity is a convenient assumption for the closure problem definition and it has a minor impact at the macroscopic level. In the following, the averaging domain is denoted as $\mathscr {V}$ and the space occupied by each fluid phase is $\mathscr {V}_\alpha$ (of measure $V_\alpha$).

Figure 1. (a) Sketch of a porous medium saturated with two fluid phases, $\beta$ (in blue) and $\gamma$ (in white), flowing in the presence of a displacement front or as coexisting mixed phases including the scales, averaging domain and characteristic lengths. (b) Example of a periodic unit cell, of side length $\ell$, illustrating the position vectors locating the barycentres of each fluid phase ($\boldsymbol {x}_\beta$ and $\boldsymbol {x}_\gamma$), as well as an example of the position vectors that locate the fluid–fluid interface with respect to a fixed system of coordinates ($\boldsymbol {r}_{\beta \gamma }$) and with respect to the phase barycentres ($\boldsymbol {z}_{\beta \gamma }$ and $\boldsymbol {z}_{\gamma \beta }$).

In terms of the averaging domain, the superficial and intrinsic averages of a pore-scale quantity, $\psi _\alpha$, defined at any point within the $\alpha$ phase, located by $\boldsymbol {r}_{\alpha }$ with respect to a fixed system of coordinates, are respectively defined by the following integral operators:

(2.2a)$$\begin{gather} \langle\psi_\alpha\rangle_\alpha|_{\boldsymbol{x}_\alpha}= \frac{1}{V}\int_{\mathscr{V}_\alpha} \psi_\alpha|_{\boldsymbol{r}_\alpha}\, {\rm d} V, \end{gather}$$

(2.2b)$$\begin{gather}\langle \psi_\alpha \rangle^\alpha |_{\boldsymbol{x}_\alpha} = \varepsilon_\alpha^{{-}1} \langle \psi_\alpha \rangle_\alpha |_{\boldsymbol{x}_\alpha} = \frac{1}{V_\alpha}\int_{\mathscr{V}_\alpha} \psi_\alpha|_{\boldsymbol{r}_\alpha}\, {\rm d} V. \end{gather}$$

The two resulting averages are assigned at the barycentre of each phase, $\boldsymbol {x}_\alpha =\langle \boldsymbol {r}_\alpha \rangle ^\alpha$, of $\mathscr {V}_\alpha$. The volume fraction of the $\alpha$ phase is denoted by $\varepsilon _\alpha =V_\alpha /V=\varepsilon S_\alpha$, where $\varepsilon =(V_\beta +V_\gamma )/V$ is the porosity of the porous medium and $S_\alpha$ the $\alpha$ phase saturation. From the above definition, a pore-scale quantity can be expressed in terms of its intrinsic average and deviation by means of the following decomposition (Gray 1975):

(2.3)\begin{equation} \psi_\alpha |_{\boldsymbol{r}_\alpha} = \langle \psi_\alpha\rangle^\alpha |_{\boldsymbol{r}_\alpha} + \tilde{\psi}_\alpha |_{\boldsymbol{r}_\alpha}.\end{equation}

Furthermore, the following Taylor series expansion can be used in order to evaluate $\langle \psi _\alpha \rangle ^\alpha$ at the phase barycentre:

(2.4)\begin{equation} \langle \psi_\alpha\rangle^\alpha |_{\boldsymbol{r}_\alpha} =\langle \psi_\alpha\rangle^\alpha |_{\boldsymbol{x}_\alpha} +\boldsymbol{z}_\alpha\boldsymbol{\cdot} \boldsymbol{\nabla} \langle \psi_\alpha\rangle^\alpha |_{\boldsymbol{x}_\alpha} +\tfrac{1}{2} \boldsymbol{z}_\alpha\boldsymbol{z}_\alpha: \boldsymbol{\nabla}\boldsymbol{\nabla} \langle \psi_\alpha\rangle^\alpha |_{\boldsymbol{x}_\alpha}+\cdots ,\end{equation}

where $\boldsymbol {z}_\alpha =\boldsymbol {r}_\alpha -\boldsymbol {x}_\alpha$ locates the same point as $\boldsymbol {r}_\alpha$ but with respect to $\boldsymbol {x}_\alpha$ as illustrated in figure 1(b). Note that, by definition, $\langle \boldsymbol {z}_\alpha \rangle ^\alpha = \boldsymbol {0}$ and $\boldsymbol {z}_\alpha = \boldsymbol {O}(r_0)$. If this expansion is used for $\boldsymbol{\nabla} \langle p_\alpha \rangle ^\alpha$, since this quantity can be assumed to be slowly varying (i.e. its spatial variations are characterised by $L$), it is reasonable to assume, on the basis of the length-scale constraint $r_0\ll L$, that

(2.5)\begin{equation} \boldsymbol{\nabla}\langle p_\alpha\rangle^\alpha|_{\boldsymbol{r}_\alpha}\simeq \boldsymbol{\nabla}\langle p_\alpha\rangle^\alpha|_{\boldsymbol{x}_\alpha}\equiv \boldsymbol{\nabla}\langle p_\alpha\rangle^\alpha, \end{equation}

which means that $\boldsymbol{\nabla} \langle p_\alpha \rangle ^\alpha$ can be considered as constant at the scale $r_0$.

The pore-scale momentum balance equation is hence written in a periodic unit cell under the following form:

(2.6)\begin{equation} \rho_\alpha\boldsymbol{v}_\alpha\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{v}_\alpha={-}\boldsymbol{\nabla} \langle p_\alpha\rangle^\alpha +\rho_\alpha\boldsymbol{b}_\alpha + \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\alpha} \quad\text{in }\mathscr{V}_\alpha,\end{equation}

with

(2.7)\begin{equation} {\boldsymbol{\mathsf{T}}}_{\tilde{p}_\alpha}={-}{\boldsymbol{\mathsf{I}}} \tilde{p}_\alpha+\mu(\varGamma_\alpha) ( \boldsymbol{\nabla} \boldsymbol{v}_\alpha+\boldsymbol{\nabla} \boldsymbol{v}^T_\alpha ). \end{equation}

To complete the statement of the flow problem, it is assumed that no mass transport and no slip take place at the solid–fluid interfaces, $\mathscr {A}_{\alpha \sigma }$, which leads to the following boundary condition:

(2.8a)\begin{equation} \boldsymbol{v}_\alpha=\boldsymbol{0}, \quad \text{at }\mathscr{A}_{\alpha \sigma}. \end{equation}

Moreover, the same assumption is imposed at the fluid–fluid interface, $\mathscr {A}_{\beta \gamma }$. Hence, the mass balance boundary condition at this interface can be written as

(2.8b)\begin{equation} \boldsymbol{v}_\beta= \boldsymbol{v}_\gamma = \boldsymbol{w} \quad\text{at }\mathscr{A}_{\beta\gamma}, \end{equation}

where $\boldsymbol {w}$ represents the speed of displacement of $\mathscr {A}_{\beta \gamma }$. In this work, no triple-phase contact is considered. Nevertheless, the pore-scale flow model statement considered here has been used by many authors (e.g. Auriault & Sanchez-Palencia 1986; Whitaker 1986, 1994; Lasseux, Quintard & Whitaker 1996; Lasseux, Ahmadi & Arani 2008; Picchi & Battiato 2018) and serves as a useful archetype for future studies that would address this issue. Indeed, there is no definite model that allows one to take into account the physics in the contact line neighbourhood, and this deserves a specific analysis that is beyond the scope of the present work.

Assuming that no specific viscosity effects occur at the fluid–fluid interface, the stress jump at $\mathscr {A}_{\beta \gamma }$ is compensated by capillary forces, leading to the following boundary condition:

(2.8c)\begin{equation} \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot} ({\boldsymbol{\mathsf{T}}}_{p_\beta}-{\boldsymbol{\mathsf{T}}}_{p_\gamma})= \nabla_s\gamma+2H\gamma\boldsymbol{n}_{\beta\gamma} \quad\text{at }\mathscr{A}_{\beta\gamma}.\end{equation}

After application of Gray's decomposition to the pressure in each phase, the above equation can be equivalently written as

(2.8d)\begin{equation} \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot} ({\boldsymbol{\mathsf{T}}}_{\tilde{p}_\beta}-{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma}) = \boldsymbol{n}_{\beta\gamma} (\langle p_\beta\rangle^\beta -\langle p_\gamma\rangle^\gamma )_{\boldsymbol{r}_{\beta\gamma}}+ \nabla_s\gamma+2H\gamma\boldsymbol{n}_{\beta\gamma} \quad\text{at }\mathscr{A}_{\beta\gamma}.\end{equation}

In the above equations, $\boldsymbol {n}_{\beta \gamma }$ is the unit normal vector at $\mathscr {A}_{\beta \gamma }$, directed from the $\beta$ phase towards the $\gamma$ phase, $\boldsymbol {r}_{\beta \gamma }$ locates points at $\mathscr {A}_{\beta \gamma }$ with respect to a fixed system of coordinates (see figure 1b), $H$ is the mean curvature of the interface and $\gamma$ is the interfacial tension between the two fluids. In this boundary condition, possible spatial variations of $\gamma$ along $\mathscr {A}_{\beta \gamma }$ are considered, expressed by $\nabla _s \gamma$, where $\nabla _s$ is the surface gradient operator defined as $\nabla _s=({\boldsymbol{\mathsf{I}}}-\boldsymbol {n}_{\beta \gamma } \boldsymbol {n}_{\beta \gamma })\boldsymbol{\cdot} \boldsymbol{\nabla}$. Note that $2H=-\nabla _s\boldsymbol{\cdot} \boldsymbol {n}_{\beta \gamma }$.

Using Taylor series expansions for the average pressure terms (see (2.4)) about the phase barycentres, $\boldsymbol {x}_\alpha$, recalling that $\boldsymbol{\nabla} \langle p_\alpha \rangle ^\alpha$ can be assumed constant within the unit cell, the boundary condition given in (2.8d) can be rewritten as

(2.8e)\begin{align} \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot} ({\boldsymbol{\mathsf{T}}}_{\tilde{p}_\beta}-{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma}) &= \boldsymbol{n}_{\beta\gamma} ( \langle p_\beta\rangle^\beta|_{\boldsymbol{x}_\beta} -\langle p_\gamma\rangle^\gamma|_{\boldsymbol{x}_\gamma} )-\boldsymbol{n}_{\beta\gamma} (\boldsymbol{x}_{\beta}-\boldsymbol{x}_{\gamma})\boldsymbol{\cdot}\boldsymbol{\nabla} \langle p_\beta \rangle^\beta\nonumber\\ &\quad +\nabla_s\gamma+2H\gamma\boldsymbol{n}_{\beta\gamma} \quad\text{at }\mathscr{A}_{\beta\gamma}. \end{align}

While writing (2.8e), use was made of the fact that, if a periodic unit cell representative of the local process is employed for the averaging domain, the macroscopic pressure gradients are equal (see the proof in Appendix B in Lasseux & Valdés-Parada (2022)).

The problem statement in the unit cell is completed with the following periodicity conditions:

(2.8f)\begin{equation} \psi (\boldsymbol{r}+\boldsymbol{l}_i)=\psi (\boldsymbol{r}), \quad i=1,2,3; \ \psi = \boldsymbol{v}_{\alpha}, \tilde{p}_{\alpha}, \end{equation}

and with the average constraint

(2.8g)\begin{equation} \langle \tilde{p}_\alpha\rangle^\alpha=0,\end{equation}

which results from applying the intrinsic averaging operator to (2.3) (with $\psi _\alpha = p_\alpha$) and the use of the Taylor series expansion (2.4), while recalling that $\langle \boldsymbol {z}_\alpha \rangle ^\alpha =\boldsymbol {0}$ and that $\boldsymbol{\nabla} \langle p_\alpha \rangle ^\beta$ is taken as a constant in the unit cell.

To finalise this section, it is worth mentioning that the assumptions imposed in this work do not constitute a severe restriction in terms of application of the model. In fact, some systems where it is reasonable to adopt these assumptions (and hence where the macroscopic models presented below are applicable) are the following: incompressible gas–liquid and liquid–liquid packed columns for extraction processes and distillation (Chhabra & Basavaraj 2019), incompressible gas and liquid flow through porous shale reservoirs (Singh & Cai 2019), water and non-Newtonian soil flooding for enhanced oil recovery (Sorbie 1991; Nilsson et al. 2013), as well as industrial applications in packed bed reactors (Giri & Majumder 2014) and bioreactors (Sen et al. 2017).

3. Mass and momentum macroscopic equations

The method to derive the macroscopic equations takes elements of the volume-averaging method, the adjoint technique and Green's integral formulation. The essential steps to carry out the developments are reported in the following paragraphs.

3.1. Averaged mass equations

Derivation of the macroscopic mass equation in each phase was reported in previous works dedicated to two-phase flow in homogeneous porous media (Whitaker 1986, 1994; Lasseux et al. 1996; Lasseux & Valdés-Parada 2022) and remains unchanged in the present configuration. The upscaled equation can be written as

(3.1)\begin{equation} \frac{\partial\varepsilon_\alpha}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot}\langle\boldsymbol{v}_\alpha\rangle_\alpha=0.\end{equation}

This expression, which is obtained without assuming that $\mathscr {V}$ is periodic, takes into account the temporal changes of the $\alpha$ phase volume fraction. Time dependency at the macroscale arises from the possible non-stationary character of $\mathscr {A}_{\beta \gamma }$, implicitly making the pore-scale problem time-dependent, even if temporal acceleration is neglected with respect to viscous effects, as already discussed in § 2. Indeed, if $\mathscr {V}$ is assumed periodic, $({1}/{V})\int _{\mathscr {A}_{\beta \gamma }}\boldsymbol {n}_{\beta \gamma } \boldsymbol{\cdot} \boldsymbol {v}_\alpha \, {\rm d} A=\langle \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol {v}_\alpha \rangle _\alpha$, which, after application of the superficial averaging operator to (2.1a), allows concluding that $\boldsymbol{\nabla} \boldsymbol{\cdot} \langle \boldsymbol {v}_\alpha \rangle _\alpha =0$. An equivalent result is obtained if the fluid–fluid interface is stationary, which means

(3.2)\begin{equation} \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}\boldsymbol{w}=0, \quad \text{at } \mathscr{A}_{\beta\gamma}. \end{equation}

In that case, it follows from the general transport theorem that ${\partial \varepsilon _\alpha }/{\partial t}=0$. This allows concluding that

(3.3)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\langle\boldsymbol{v}_\alpha\rangle_\alpha=0, \quad \mathscr{V}\ \text{periodic or } \mathscr{A}_{\beta\gamma} \ \text{stationary}. \end{equation}

3.2. Averaged momentum equations

To obtain the macroscopic momentum equation in each phase, the upscaling approach detailed in Lasseux & Valdés-Parada (2022) is used. In brief, it consists of the following steps. First, propose pertinent adjoint (or closure) problems in a periodic unit cell, which can be subsequently related to the flow problem by means of Green's formula. Second, substitute the corresponding differential equations and boundary conditions to carry out the necessary algebraic manipulations to obtain the upscaled models. In this way, the derivations commence by introducing the four second-order tensors, ${\boldsymbol{\mathsf{D}}}_{\alpha \kappa }$, and vectors, $\boldsymbol {d}_{\alpha \kappa }(\alpha,\kappa=\beta,\gamma)$, that solve the following adjoint closure problems in a periodic unit cell (cf. Bottaro 2019):

Problem I

(3.4a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\alpha \beta}=\boldsymbol{0}, \quad\text{in }\mathscr{V}_\alpha, \end{gather}$$

(3.4b)$$\begin{gather}-\frac{\rho_\alpha}{\mu_{{ref}\alpha}}\boldsymbol{v}_\alpha \boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\mathsf{D}}}_{\alpha \beta}=\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{{d}}_{\alpha\beta}} +\delta^K_{\alpha\beta} {\boldsymbol{\mathsf{I}}}, \quad\text{in }\mathscr{V}_\alpha, \end{gather}$$

(3.4c)$$\begin{gather}\boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot} \mu_{{ref}\beta} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{{d}}_{\beta\beta}} =\boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}\mu_{{ref}\gamma} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{{d}}_{\gamma\beta}} , \quad\text{at }\mathscr{A}_{\beta\gamma}, \end{gather}$$

(3.4d)$$\begin{gather}{\boldsymbol{\mathsf{D}}}_{\beta\beta}={\boldsymbol{\mathsf{D}}}_{\gamma\beta}, \quad\text{at }\mathscr{A}_{\beta\gamma}, \end{gather}$$

(3.4e)$$\begin{gather}{\boldsymbol{\mathsf{D}}}_{ \alpha\beta}=\boldsymbol{0}, \quad\text{at }\mathscr{A}_{ \alpha\sigma}, \end{gather}$$

(3.4f)$$\begin{gather}\psi (\boldsymbol{r}+\boldsymbol{l}_i)=\psi (\boldsymbol{r}), \quad i=1,2,3; \ \psi = {\boldsymbol{\mathsf{D}}}_{\alpha\beta}, \boldsymbol{d}_{\alpha\beta}, \end{gather}$$

(3.4g)$$\begin{gather}\boldsymbol{d}_{\alpha\beta} =\boldsymbol{0}, \quad \text{at }\boldsymbol{r}_{\alpha}^0. \end{gather}$$

Problem II

(3.5a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\alpha\gamma}=\boldsymbol{0}, \quad\text{in }\mathscr{V}_\alpha, \end{gather}$$

(3.5b)$$\begin{gather}-\frac{\rho_\alpha}{\mu_{{ref}\alpha}}\boldsymbol{v}_\alpha \boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\mathsf{D}}}_{\alpha\gamma}=\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{\boldsymbol{d}}_{\alpha\gamma}}+ \delta^K_{\alpha\gamma}{\boldsymbol{\mathsf{I}}} , \quad\text{in }\mathscr{V}_\alpha, \end{gather}$$

(3.5c)$$\begin{gather}\boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}\mu_{{ref}\beta} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{{d}}_{\beta\gamma}} =\boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}\mu_{{ref}\gamma} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{{d}}_{\gamma\gamma}}, \quad\text{at }\mathscr{A}_{\beta\gamma}, \end{gather}$$

(3.5d)$$\begin{gather}{\boldsymbol{\mathsf{D}}}_{\beta\gamma}={\boldsymbol{\mathsf{D}}}_{\gamma\gamma}, \quad\text{at }\mathscr{A}_{\beta\gamma}, \end{gather}$$

(3.5e)$$\begin{gather}{\boldsymbol{\mathsf{D}}}_{ \alpha\gamma}=\boldsymbol{0}, \quad\text{at }\mathscr{A}_{ \alpha\sigma}, \end{gather}$$

(3.5f)$$\begin{gather}\psi (\boldsymbol{r}+\boldsymbol{l}_i)=\psi (\boldsymbol{r}), \quad i=1,2,3; \ \psi = {\boldsymbol{\mathsf{D}}}_{\alpha\gamma},\boldsymbol{d}_{\alpha\gamma}, \end{gather}$$

(3.5g)$$\begin{gather}\boldsymbol{d}_{\alpha\gamma} =\boldsymbol{0}, \quad \text{at }\boldsymbol{r}_{\alpha}^0. \end{gather}$$

In the above equations, $\boldsymbol {r}_{\alpha }^0$ locates an arbitrary point in the $\alpha$ phase and $\delta _{\alpha \kappa }^K$ is the Kronecker delta. In addition,

(3.6)\begin{equation} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{{d}}_{\alpha\kappa}}={-}{\boldsymbol{\mathsf{I}}}\boldsymbol{d}_{\alpha\kappa}+ \frac{\mu(\varGamma_\alpha)}{\mu_{{ref}\alpha}} (\boldsymbol{\nabla}{\boldsymbol{\mathsf{D}}}_{\alpha\kappa}+\boldsymbol{\nabla}{\boldsymbol{\mathsf{D}}}_{\alpha \kappa}^{T1}), \end{equation}

where the superscript $T1$ denotes the transpose of a third-order tensor that permutes its first two indices, i.e. $(\boldsymbol{\nabla} {\boldsymbol{\mathsf{D}}}_{\alpha \kappa })_{ijk}^{T1} =(\boldsymbol{\nabla} {\boldsymbol{\mathsf{D}}}_{\alpha \kappa })_{jik}$. Note that the boundary conditions given in (3.4g) and (3.5g) are necessary in order to well pose the closure problems. A practical solution of closure problems I and II, written under the forms given in (3.4) and (3.5), entails a pre-determination of the pore-scale velocity and viscosity fields over the corresponding periodic domain. These adjoint problems can be derived from the associated velocity Green's function pair problems, as detailed in Lasseux & Valdés-Parada (2022) for Newtonian creeping two-phase flow. It is worth noting that, in that case, the same closure problems as those previously reported in Auriault & Sanchez-Palencia (1986), Whitaker (1994), Lasseux et al. (1996) and Picchi & Battiato (2018) were retrieved, albeit the macroscale model is different. Analogously, the closure problems reported in (3.4) and (3.5) are the adjoint versions of those derived in Lasseux et al. (2008) (see (114) and (115) therein) in the case of Newtonian and inertial two-phase flow in homogeneous porous media.

The next step in the derivations consists of relating the associated closure problems with the flow problem with the help of a Green's formula. Considering a constant scalar $d$ and arbitrary and sufficiently regular scalar fields $a$ and $c$, vector fields $\boldsymbol {a}$, $\boldsymbol {b}$ and $\boldsymbol {c}$ and a second-order tensor field ${\boldsymbol{\mathsf{B}}}$, together with the following definitions:

(3.7a)$$\begin{gather} {\boldsymbol{\mathsf{T}}}_{a} ={-}{\boldsymbol{\mathsf{I}}}a + c (\boldsymbol{\nabla} \boldsymbol{a}+ \boldsymbol{\nabla} \boldsymbol{a}^T), \end{gather}$$

(3.7b)$$\begin{gather}{\boldsymbol{\mathsf{T}}}_{{\boldsymbol{b}}} ={-}{\boldsymbol{\mathsf{I}}}\boldsymbol{b}+ c(\boldsymbol{\nabla} {\boldsymbol{\mathsf{B}}}+\boldsymbol{\nabla} {\boldsymbol{\mathsf{B}}}^{T1}), \end{gather}$$

this formula is written (see the proof in Appendix A in Sánchez-Vargas, Valdés-Parada & Lasseux (2022))

(3.8)\begin{align} &\int_{\mathscr{V}_\alpha } [ \boldsymbol{a} \boldsymbol{\cdot} ( d \boldsymbol{c} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{\mathsf{B}}} +\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{b}} ) - ({-}d \boldsymbol{c}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{a} + \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{a}) \boldsymbol{\cdot} {\boldsymbol{\mathsf{B}}} -\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{a} \boldsymbol{b} \nonumber\\ &\quad +d (\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{c} \boldsymbol{a})\boldsymbol{\cdot} {\boldsymbol{\mathsf{B}}} +a\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\mathsf{B}}} ]\, {\rm d} V\nonumber\\ &\qquad = \int_{\mathscr{A}_\alpha } [\boldsymbol{a} \boldsymbol{\cdot} ( \boldsymbol{n}_\alpha\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{\boldsymbol{b}} )- \boldsymbol{n}_\alpha \boldsymbol{\cdot} ({-}d \boldsymbol{c}\boldsymbol{a} +{\boldsymbol{\mathsf{T}}}_{{a}} ) \boldsymbol{\cdot} {\boldsymbol{\mathsf{B}}}]\, {\rm d} A. \end{align}

Here, $\mathscr {A}_\alpha$ denotes the surfaces bounding $\mathscr {V}_\alpha$.

Attention is now focused on the derivation of the average momentum balance in the $\beta$ phase. To this end, the above Green's formula is employed, taking $\boldsymbol {a}=\boldsymbol {v}_\beta$, $a= \tilde {p}_\beta$, ${\boldsymbol{\mathsf{B}}}={\boldsymbol{\mathsf{D}}}_{\beta \beta }$, $\boldsymbol {b}=\mu _{{ref}\beta }\boldsymbol {d}_{\beta \beta }$, $c=\mu (\varGamma _\beta )$, $\boldsymbol {c}=\boldsymbol {v}_\beta$ and $d=\rho _\beta$. Once the corresponding differential equations and boundary conditions from the pore-scale problem and closure problem I are substituted in the ensuing relationship, keeping in mind the definition of the superficial averaging operator and the fact that $\boldsymbol{\nabla} \langle p_\alpha \rangle ^\alpha$ and $\rho _\alpha \boldsymbol {b}_\alpha$ are constant in $\mathscr {V}_\alpha$, the following expression for $\langle \boldsymbol {v}_\beta \rangle _\beta$ is obtained:

(3.9a)\begin{align} \langle\boldsymbol{v}_\beta\rangle_\beta &={-}\frac{1}{\mu_{{ref}\beta}} (\boldsymbol{\nabla}\langle p_\beta\rangle^\beta -\rho_\beta\boldsymbol{b}_\beta ) \boldsymbol{\cdot} \langle{\boldsymbol{\mathsf{D}}}_{\beta\beta}\rangle_{\beta} + \frac{1}{\mu_{{ref}\beta}V}\int_{\mathscr{A}_{\beta\gamma}} ( 2\gamma H \boldsymbol{n}_{\beta \gamma}+\nabla_s \gamma )\boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\beta\beta}\, {\rm d} A \nonumber\\ &\quad -\frac{1}{\mu_{{ref}\beta}V} \int_{\mathscr{A}_{\beta\gamma} } [\mu_{{ref}\gamma}\boldsymbol{w} \boldsymbol{\cdot} ( \boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{\boldsymbol{d}}_{\gamma\beta}} )-\boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot} (-\rho_\beta\boldsymbol{w}\boldsymbol{w}+{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma} )\boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\gamma\beta}]\, {\rm d} A\nonumber\\ &\quad+\frac{1}{\mu_{{ref}\beta}V} \int_{\mathscr{A}_{\beta\gamma}}\boldsymbol{n}_{\beta \gamma} \boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\beta\beta }\, {\rm d} A (\langle p_\beta \rangle^\beta |_{\boldsymbol{x}_\beta} -\langle p_\gamma \rangle^\gamma |_{\boldsymbol{x}_\gamma} -(\boldsymbol{x}_\beta -\boldsymbol{x}_\gamma) \boldsymbol{\cdot} \boldsymbol{\nabla} \langle p_\beta \rangle^\beta). \end{align}

Since ${\boldsymbol{\mathsf{D}}}_{\beta \beta }$ is a solenoidal field that is periodic and zero at $\mathscr {A}_{\beta \sigma }$, one can make use of the divergence theorem to conclude that $\int _{\mathscr {V}_\beta }\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\beta \beta }\, {\rm d} V=\int _{\mathscr {A}_{\beta \gamma }}\boldsymbol {n}_{\beta \gamma } \boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\beta \beta }\, {\rm d} A=\boldsymbol {0}$. This leads to simplifying the expression of $\langle \boldsymbol {v}_\beta \rangle _\beta$ into the following form:

(3.9b)\begin{align} \langle\boldsymbol{v}_\beta\rangle_\beta &={-}\frac{1}{\mu_{{ref}\beta}} (\boldsymbol{\nabla}\langle p_\beta\rangle^\beta -\rho_\beta\boldsymbol{b}_\beta ) \boldsymbol{\cdot} \langle{\boldsymbol{\mathsf{D}}}_{\beta\beta}\rangle_{\beta} +\frac{1}{\mu_{{ref}\beta}V}\int_{\mathscr{A}_{\beta\gamma}} ( 2\gamma H \boldsymbol{n}_{\beta \gamma}+\nabla_s \gamma )\boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\beta\beta}\, {\rm d} A \nonumber\\ &\quad -\frac{1}{\mu_{{ref}\beta}V} \int_{\mathscr{A}_{\beta\gamma} } [\mu_{{ref}\gamma}\boldsymbol{w} \boldsymbol{\cdot} ( \boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{\boldsymbol{d}}_{\gamma\beta}} ) -\boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot} (-\rho_\beta\boldsymbol{w}\boldsymbol{w}+{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma} )\boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\gamma\beta}]\, {\rm d} A. \end{align}

This last relationship is still non-local due to the presence of $\tilde {p}_\gamma$ in the last interfacial integral term. To progress towards a local form, Green's formula (3.8) can be employed once more taking $\boldsymbol {a}=\boldsymbol {v}_\gamma$, $a= \tilde {p}_\gamma$, ${\boldsymbol{\mathsf{B}}}={\boldsymbol{\mathsf{D}}}_{\gamma \beta }$, $\boldsymbol {b}=\mu _{{ref}\gamma }\boldsymbol {d}_{\gamma \beta }$, $c=\mu (\varGamma _\gamma )$, $\boldsymbol {c}=\boldsymbol {v}_\gamma$ and $d=\rho _\gamma$. Using the definition of the superficial averaging operator, the resulting expression can be written as follows:

(3.10)\begin{align} &\frac{1}{V}\int_{\mathscr{A}_{\beta\gamma} } [\mu_{{ref}\gamma}\boldsymbol{w} \boldsymbol{\cdot}( \boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{\boldsymbol{d}}_{\gamma\beta}} )- \boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot} (-\rho_\gamma\boldsymbol{w}\boldsymbol{w}+{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma} )\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\gamma\beta}]\, {\rm d} A \nonumber\\ &\quad =(\boldsymbol{\nabla} \langle p_\gamma\rangle^\gamma-\rho_\gamma\boldsymbol{b}_\gamma )\boldsymbol{\cdot} \langle{\boldsymbol{\mathsf{D}}}_{\gamma\beta}\rangle_\gamma. \end{align}

Substituting this expression in (3.9b), while recalling the interfacial boundary condition given in (3.4d), yields

(3.11)\begin{align} \langle \boldsymbol{v}_\beta\rangle_\beta &={-}\frac{1}{\mu_{{ref}\beta}} (\boldsymbol{\nabla}\langle p_\beta\rangle^\beta -\rho_\beta\boldsymbol{b}_\beta ) \boldsymbol{\cdot} \langle {\boldsymbol{\mathsf{D}}}_{\beta\beta}\rangle_\beta -\frac{1}{\mu_{{ref}\beta}}(\boldsymbol{\nabla} \langle p_\gamma\rangle^\gamma-\rho_\gamma\boldsymbol{b}_\gamma )\boldsymbol{\cdot} \langle {\boldsymbol{\mathsf{D}}}_{\gamma\beta} \rangle_\gamma\nonumber\\ &\quad +\frac{1}{\mu_{{ref}\beta}V}\int_{\mathscr{A}_{\beta\gamma}} ( 2\gamma H \boldsymbol{n}_{\beta \gamma}+\nabla_s \gamma )\boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}_{\beta\beta} \, {\rm d} A +\frac{\rho_\gamma -\rho_\beta}{\mu_{{ref}\beta}V} \int_{\mathscr{A}_{\beta\gamma} } \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}\boldsymbol{w} \boldsymbol{w}\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\beta\beta}\,{\rm d} A. \end{align}

The last area integral term on the right-hand side of the above equation is related to the displacement of $\mathscr {A}_{\beta \gamma }$ when inertia is present.

The above procedure can be repeated in a similar way for the $\gamma$ phase. The resulting equation can be written as follows:

(3.12)\begin{align} \langle\boldsymbol{v}_\gamma\rangle_\gamma &={-}\frac{1}{\mu_{{ref}\gamma}}(\boldsymbol{\nabla}\langle p_\gamma\rangle^\gamma -\rho_\gamma\boldsymbol{b}_\gamma ) \boldsymbol{\cdot} \langle {\boldsymbol{\mathsf{D}}}_{\gamma\gamma}\rangle_\gamma - \frac{1}{\mu_{{ref}\gamma}}(\boldsymbol{\nabla}\langle p_\beta\rangle^\beta -\rho_\beta\boldsymbol{b}_\beta ) \boldsymbol{\cdot} \langle {\boldsymbol{\mathsf{D}}}_{\beta\gamma}\rangle_\beta\nonumber\\ &\quad+\frac{1}{\mu_{{ref}\gamma}V}\int_{\mathscr{A}_{\beta\gamma} } ( 2\gamma H \boldsymbol{n}_{\beta\gamma}+\nabla_s \gamma)\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\gamma\gamma}\, {\rm d} A +\frac{\rho_\gamma -\rho_\beta}{\mu_{{ref}\gamma}V} \int_{\mathscr{A}_{\beta\gamma} }\boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot}\boldsymbol{w}\boldsymbol{w}\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\gamma\gamma}\, {\rm d} A. \end{align}

At this point, it is convenient to adopt the following nomenclature ($\alpha,\kappa =\beta, \gamma$):

(3.13)\begin{equation} {\boldsymbol{\mathsf{H}}}_{\alpha\kappa}=\frac{\mu_{{ref} \kappa}}{\mu_{{ref}\alpha}}\langle{\boldsymbol{\mathsf{D}}}_{\kappa \alpha}\rangle_\kappa^T,\end{equation}

${\boldsymbol{\mathsf{H}}}_{\alpha \alpha }$ and ${\boldsymbol{\mathsf{H}}}_{\alpha \kappa }$ being respectively identified as the classical dominant and coupling apparent permeability tensors in the $\alpha$ phase. This allows writing the macroscopic momentum balance equations in the following compact form:

(3.14a)\begin{align} \langle\boldsymbol{v}_\alpha\rangle_\alpha &={-} \frac{{\boldsymbol{\mathsf{H}}}_{\alpha\alpha}}{\mu_{{ref}\alpha}} \boldsymbol{\cdot}(\boldsymbol{\nabla}\langle p_\alpha\rangle^\alpha -\rho_\alpha \boldsymbol{b}_\alpha ) -\frac{{\boldsymbol{\mathsf{H}}}_{\alpha\kappa}}{\mu_{{ref}\kappa}}\boldsymbol{\cdot} (\boldsymbol{\nabla}\langle p_\kappa\rangle^\kappa -\rho_\kappa\boldsymbol{b}_\kappa )\nonumber\\ &\quad+\frac{1}{\mu_{{ref}\alpha}V}\int_{\mathscr{A}_{\beta\gamma} } ( 2\gamma H \boldsymbol{n}_{\beta\gamma}+\nabla_s \gamma)\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\alpha\alpha}\, {\rm d} A +\frac{\rho_\gamma -\rho_\beta}{\mu_{{ref}\alpha}V} \int_{\mathscr{A}_{\beta\gamma} } \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}\boldsymbol{w} \boldsymbol{w}\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\alpha\alpha}\, {\rm d} A. \end{align}

This is the first important result of the present work. The first two terms on the right-hand side of (3.14a) are the Darcy-like terms. The first one accounts for the viscous resistance in the phase under concern, while the second one results from the viscous coupling between the two fluid phases. These terms include bulk inertial effects, justifying the apparent and non-intrinsic character of the permeability tensors ${\boldsymbol{\mathsf{H}}}_{\alpha \kappa }$. In addition to the Darcy-like terms, a contribution of the interfacial effects is also present, which is accounted for by the last two terms in (3.14a). The first one, dealing with interfacial tension effects, is formally the same as that identified in a previous work dedicated to creeping Newtonian two-phase flow in porous media (Lasseux & Valdés-Parada 2022). The second one results from inertial effects. Note that, both inertial and non-Newtonian effects are reflected in the values of the closure variable ${\boldsymbol{\mathsf{D}}}_{\alpha \alpha }$ at $\mathscr {A}_{\beta \gamma }$ (and also in $H$). The first interfacial term was not present in the development of the average model proposed by Lasseux et al. (2008) in the case of two-phase inertial Newtonian flow. This is because a simplification of the contribution of the capillary forces was assumed. Such an assumption lacks careful justification, and this was extensively discussed in the case of creeping Newtonian two-phase flow (see Lasseux & Valdés-Parada (2022), pp. 14–15). When the fluid–fluid interface is stationary, or flow remains in the creeping regime in the $\alpha$ phase, this equation reduces to

(3.14b)\begin{align} \langle\boldsymbol{v}_\alpha\rangle_\alpha &={-}\frac{{\boldsymbol{\mathsf{H}}}_{\alpha\alpha}}{\mu_{{ref}\alpha}} \boldsymbol{\cdot}(\boldsymbol{\nabla}\langle p_\alpha\rangle^\alpha -\rho_\alpha\boldsymbol{b}_\alpha ) - \frac{{\boldsymbol{\mathsf{H}}}_{\alpha\kappa}}{\mu_{{ref}\kappa}} \boldsymbol{\cdot}(\boldsymbol{\nabla}\langle p_\kappa\rangle^\kappa -\rho_\kappa\boldsymbol{b}_\kappa )\nonumber\\ &\quad +\frac{1}{\mu_{{ref}\alpha}V}\int_{\mathscr{A}_{\beta\gamma} } ( 2\gamma H \boldsymbol{n}_{\beta\gamma}+\nabla_s \gamma)\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\alpha\alpha}\, {\rm d} A,\nonumber\\ &\quad\quad \quad \quad \quad \text{creeping flow in the }\alpha \text{ phase or}\ \mathscr{A}_{\beta\gamma} \ \text{stationary}. \end{align}

Note that (3.14b) is valid in the creeping flow regime, even if $\boldsymbol {n}_{\beta \gamma }\boldsymbol{\cdot} \boldsymbol {w}\ne 0$. To that extent, the model reported in Lasseux & Valdés-Parada (2022) is valid even if the fluid–fluid interface is not in its stationary configuration.

All the terms in the macroscopic momentum equations are determined from the solution of the two (adjoint) closure problems I and II, which requires knowledge of the interface position, $\mathscr {A}_{\beta \gamma }$, and of its speed of displacement, $\boldsymbol {w}$, if inertia is significant and the interface has not reached a stationary configuration. Note that these closure problems also require accounting for the velocity field in each phase even under creeping flow conditions for the computation of the viscosity coefficient $\mu (\varGamma _\alpha )$.

3.3. Symmetry and reciprocity relationships for ${\boldsymbol{\mathsf{H}}}_{\alpha \alpha }$ and ${\boldsymbol{\mathsf{H}}}_{\alpha \kappa }$ ($\alpha, \gamma =\beta, \gamma$, $\alpha \ne \kappa$)

In Appendix A, a symmetry analysis of ${\boldsymbol{\mathsf{H}}}_{\alpha \alpha }$ is provided, along with a reciprocity relationship between ${\boldsymbol{\mathsf{H}}}_{\beta \gamma }$ and ${\boldsymbol{\mathsf{H}}}_{\gamma \beta }$. The results are respectively given in (A7) and (A10). The relationship (A7) shows, first, that ${\boldsymbol{\mathsf{H}}}_{\alpha \alpha }$ is not symmetric and, second, that its skew-symmetric part only results from inertial effects. In addition, the reciprocity relationship given in (A10) indicates that the two coupling apparent permeability tensors are not independent and that the right-hand side also only results from inertia. Indeed, in the absence of inertial effects, the right-hand sides of both (A7) and (A10) are zero. This allows concluding that, under creeping flow conditions, the dominant permeability tensors are symmetric, as for Newtonian flow, and the coupling permeability tensors satisfy the same relationship as for $\mu (\varGamma _\alpha )={Cst}=\mu _{{ref}\alpha }$. This extends the results reported in Lasseux et al. (1996) and Lasseux & Valdés-Parada (2017) (section II.E) for Newtonian two-phase flow, showing that the non-Newtonian character of the flow does not alter the properties of the permeability tensors in the absence of inertia. It also extends the result obtained for one-phase incompressible Newtonian and inertial flow in porous media for which the apparent permeability tensor was shown to be non-symmetric, with its skew-symmetric part only due to inertia (Lasseux & Valdés-Parada 2017, § II.A).

To complete the macroscopic description of the two-phase flow process, it is of interest to derive an expression for the macroscpic pressure difference. This is the objective of the development that follows.

4. Macroscopic pressure difference

The derivation steps are the same as those followed in Lasseux & Valdés-Parada (2023) in which a model for the macroscopic pressure difference was proposed for two-phase creeping Newtonian flow in homogeneous porous media. In order to determine the expression for the macroscopic pressure difference, i.e. $\langle p_\beta \rangle ^\beta |_{\boldsymbol {x}_\beta }-\langle p_\gamma \rangle ^\gamma |_{\boldsymbol {x}_\gamma }$, it is convenient to introduce the scalar and vector fields, $f_\alpha$ and $\boldsymbol {f}_\alpha$. These closure variables are associated with the pore-scale fluid pressures and solve a new adjoint (closure) problem that is defined as follows:

(4.1a)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{f}_\beta = \frac{1}{V_\beta}, \quad\text{in }\mathscr{V}_\beta, \end{gather}$$

(4.1b)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{f}_\gamma ={-}\frac{1}{V_\gamma}, \quad\text{in }\mathscr{V}_\gamma, \end{gather}$$

(4.1c)$$\begin{gather}-\rho_\alpha\boldsymbol{v}_\alpha\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{f}_\alpha=\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{{f}_\alpha}, \quad\text{in }\mathscr{V}_\alpha, \end{gather}$$

(4.1d)$$\begin{gather}\boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{{f}_\beta}= \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{{f}_\gamma}, \quad\text{at }\mathscr{A}_{\beta\gamma}, \end{gather}$$

(4.1e)$$\begin{gather}\boldsymbol{f}_{\beta}=\boldsymbol{f}_\gamma, \quad\text{at }\mathscr{A}_{\beta\gamma}, \end{gather}$$

(4.1f)$$\begin{gather}\boldsymbol{f}_{ \alpha}=\boldsymbol{0}, \quad\text{at }\mathscr{A}_{ \alpha\sigma}, \end{gather}$$

(4.1g)$$\begin{gather}\psi (\boldsymbol{r}+\boldsymbol{l}_i)=\psi (\boldsymbol{r}), \quad i=1,2,3; \ \psi = \boldsymbol{f}_\alpha,{f}_\alpha, \end{gather}$$

(4.1h)$$\begin{gather}{f}_\alpha={0}, \quad \text{at }\boldsymbol{r}_{\alpha}^0. \end{gather}$$

In this problem, in which ${\boldsymbol{\mathsf{T}}}_{f_\alpha }=-f_\alpha {\boldsymbol{\mathsf{I}}}+\mu (\varGamma _\alpha )(\boldsymbol{\nabla} \boldsymbol {f}_\alpha +\boldsymbol{\nabla} \boldsymbol {f}_\alpha ^T)$ is the associated stress-like second-order tensor, the closure variables $\boldsymbol {f}_\alpha$ are non-solenoidal. In fact, this new closure problem involves the source terms in the mass-like equations and can be inferred from the pressure Green's functions pair problem, as reported by Choi & Dong (2021) for one-phase creeping flow. It clearly differs from the two adjoint problems associated with the velocities reported in § 3.2, where the source terms are located in the momentum-like equations (3.4b) and (3.5b).

In order to derive the relationship between the new adjoint closure variables and the pore-scale velocities and pressures, the following Green's formula is employed (see the derivation in Sánchez-Vargas et al. 2022):

(4.2)\begin{align} &\int_{\mathscr{V}_\alpha} [ \boldsymbol{a}\boldsymbol{\cdot} ( d \boldsymbol{c}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{b} +\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{b}} ) -({-}d\boldsymbol{c}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{a} +\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{a}} ) \boldsymbol{\cdot} \boldsymbol{b}-b\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{a}\nonumber\\ &\quad +d (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{c} \boldsymbol{a})\boldsymbol{\cdot} \boldsymbol{b} + a \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{b}]\, {\rm d} V \nonumber\\ &\qquad = \int_{\mathscr{A}_\alpha} [ \boldsymbol{a}\boldsymbol{\cdot} ( \boldsymbol{n}_\alpha\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{b}}) - ( \boldsymbol{n}_\alpha \boldsymbol{\cdot} ({-}d\boldsymbol{ca}+{\boldsymbol{\mathsf{T}}}_{{a}} )) \boldsymbol{\cdot} \boldsymbol{b} ]\, {\rm d} A, \end{align}

with $\boldsymbol{\mathsf{T}}_a$ defined in (3.7a) and

(4.3)\begin{equation} {\boldsymbol{\mathsf{T}}}_{b} ={-}{\boldsymbol{\mathsf{I}}}b + c (\boldsymbol{\nabla} \boldsymbol{b}+ \boldsymbol{\nabla} \boldsymbol{b}^T). \end{equation}

In these two equations, $d$ is a constant scalar, whereas $a$, $b$, $c$ and $\boldsymbol {a}$, $\boldsymbol {b}$, $\boldsymbol {c}$ are respectively three arbitrary scalar and three arbitrary vector fields having sufficient regularity.

When this formula is employed with $\boldsymbol {a}=\boldsymbol {c}=\boldsymbol {v}_\beta$, $c=\mu (\varGamma _\beta )$, $a=\tilde {p}_\beta$, $\boldsymbol {b}=\boldsymbol {f}_\beta$, $b=f_\beta$ and $d=\rho _\beta$, and by making use of the differential equations and boundary conditions for the pore-scale flow and the adjoint closure problem, the resulting expression can be written as follows:

(4.4)\begin{align} \langle p_\gamma\rangle^\gamma|_{\boldsymbol{x}_\gamma}- \langle p_\beta\rangle^\beta|_{\boldsymbol{x}_\beta}&={-}\int_{\mathscr{A}_{\beta\gamma}} [ \boldsymbol{w}\boldsymbol{\cdot} ( \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{f}_\gamma} ) - ( \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}(-\rho_\beta \boldsymbol{w}\boldsymbol{w}+{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma} )) \boldsymbol{\cdot} \boldsymbol{f}_\beta]\, {\rm d} A \nonumber\\ & \quad -(\boldsymbol{\varphi}_\beta+\boldsymbol{x}_\beta- \boldsymbol{x}_\gamma)\boldsymbol{\cdot}\boldsymbol{\nabla}\langle p_\beta\rangle^\beta +\rho_\beta\boldsymbol{b}_\beta\boldsymbol{\cdot}\boldsymbol{\varphi}_\beta +s_{\beta\gamma}. \end{align}

To arrive at this result, the fact that $\int _{\mathscr {A}_{\beta \gamma }} \boldsymbol {n}_{\beta \gamma }\boldsymbol{\cdot} \boldsymbol {f}_\beta \, {\rm d} A=1$ was employed. Moreover, the following nomenclature was adopted:

(4.5a)$$\begin{gather} \boldsymbol{\varphi}_\alpha = \int_{\mathscr{V}_\alpha} \boldsymbol{f}_\alpha\, {\rm d} V, \end{gather}$$

(4.5b)$$\begin{gather}s_{\beta\gamma } = \int_{\mathscr{A}_{\beta\gamma}} (\nabla_s\gamma+2H\gamma\boldsymbol{n}_{\beta\gamma} )\boldsymbol{\cdot} \boldsymbol{f}_\beta\, {\rm d} A. \end{gather}$$

Let now attention be directed to the $\gamma$ phase. The Green's formula given in (4.2) is used again, setting $\boldsymbol {a}=\boldsymbol {c}=\boldsymbol {v}_\gamma$, $c=\mu (\varGamma _\gamma )$, $a=\tilde {p}_\gamma$, $\boldsymbol {b}=\boldsymbol {f}_\gamma$, $b=f_\gamma$ and $d=\rho _\gamma$, and, by making use of the corresponding differential equations and boundary conditions, this leads to

(4.6)\begin{align} ( \boldsymbol{\nabla}\langle p_\beta\rangle^\beta-\rho_\gamma\boldsymbol{b}_\gamma)\boldsymbol{\cdot} \boldsymbol{\varphi}_\gamma =\int_{\mathscr{A}_{\beta\gamma}} [ \boldsymbol{w}\boldsymbol{\cdot} ( \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{f}_\gamma} ) - ( \boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot} (-\rho_\gamma\boldsymbol{w}\boldsymbol{w}+{\boldsymbol{\mathsf{T}}}_{\tilde{p}_\gamma} ) ) \boldsymbol{\cdot} \boldsymbol{f}_\beta]\, {\rm d} A. \end{align}

When this result is substituted back into (4.4), the following expression for the macroscopic pressure difference, ${\rm \Delta} \mathcal {P} \equiv \langle p_\gamma \rangle ^\gamma |_{\boldsymbol {x}_\gamma }-\langle p_\beta \rangle ^\beta |_{\boldsymbol {x}_\beta }$, is obtained:

(4.7a)\begin{align} {\rm \Delta}\mathcal{P} &={-}(\boldsymbol{\varphi}_\beta+ \boldsymbol{\varphi}_\gamma+\boldsymbol{x}_\beta- \boldsymbol{x}_\gamma)\boldsymbol{\cdot}\boldsymbol{\nabla}\langle p_\beta\rangle^\beta +\rho_\beta \boldsymbol{b}_\beta\boldsymbol{\cdot}\boldsymbol{\varphi}_\beta +\rho_\gamma \boldsymbol{b}_\gamma\boldsymbol{\cdot}\boldsymbol{\varphi}_\gamma +s_{\beta\gamma} \nonumber\\ &\quad +(\rho_\gamma-\rho_\beta) \int_{\mathscr{A}_{\beta\gamma}}\boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot}\boldsymbol{w}\boldsymbol{w}\boldsymbol{\cdot}\boldsymbol{f}_\beta\, {\rm d} A. \end{align}

This is the second salient result of the developments reported in this work, which extends the analysis detailed in Lasseux & Valdés-Parada (2023) to inertial and non-Newtonian flow. When there is no inertia and $\mu (\varGamma _\alpha )$ is constant, the expression for ${\rm \Delta} \mathcal {P}$ exactly coincides with the one derived in that reference.

From (4.7a), it is clear that the average pressure difference depends on the macroscopic pressure gradient and bulk force in each phase as well as on capillary effects, implicitly involving the shape, position and area of $\mathscr {A}_{\beta \gamma }$. It also depends on inertia, when present, if the fluid–fluid interface is not stationary, as reflected by the last term in (4.7a). All these dependencies can be determined from the solution of the adjoint closure problem given in (4.1), providing the fields of $\boldsymbol {f}_\alpha$ that are affected by inertia and the non-Newtonian character of the fluids. The first three terms and the last one on the right-hand side of (4.7a) only result from dynamic effects contributing to the average pressure difference. Nevertheless, dynamic effects are also implicitly present in $s_{\beta \gamma }$ through $\boldsymbol {f}_\beta$, which operates as a weighting factor. When inertia is insignificant, or when $\mathscr {A}_{\beta \gamma }$ is stationary, (4.7a) simplifies to

(4.7b)\begin{align} {\rm \Delta}\mathcal{P} ={-}(\boldsymbol{\varphi}_\beta+ \boldsymbol{\varphi}_\gamma+\boldsymbol{x}_\beta- \boldsymbol{x}_\gamma)\boldsymbol{\cdot}\boldsymbol{\nabla}\langle p_\beta\rangle^\beta &+\rho_\beta \boldsymbol{b}_\beta\boldsymbol{\cdot}\boldsymbol{\varphi}_\beta +\rho_\gamma \boldsymbol{b}_\gamma\boldsymbol{\cdot}\boldsymbol{\varphi}_\gamma +s_{\beta\gamma}, \nonumber\\ & \quad \text{creeping flow or} \mathscr{A}_{\beta\gamma} \ \text{stationary}. \end{align}

Finally, when there is no macroscopic forcing, no interfacial tension gradient and the fluid–fluid interface is in its stationary position, $H$ is constant and ${\rm \Delta} \mathcal {P}$ reduces to the classical Laplace relationship, ${\rm \Delta} \mathcal {P} =\mathcal {P}_{cL}=2H\gamma$. When these conditions are not met, ${\rm \Delta} \mathcal {P}$ may significantly differ from this equilibrium relationship as previously suggested by Marle (1981), Hassanizadeh & Gray (1993) and Gray & Miller (2011) and evidenced in Lasseux & Valdés-Parada (2023).

5. Macroscopic model summary

The set of averaged equations that constitutes the macroscopic models derived here is now summarised. On the basis of the existence and pertinence of a periodic unit cell, justified by the length-scale separation, the macroscale mass and momentum balance equations for incompressible, quasi-steady, inertial, generalised Newtonian flow in rigid and homogeneous porous media can be written in the following manner:

(5.1a)\begin{align} &\frac{\partial\varepsilon_\alpha}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\boldsymbol{v}_\alpha\rangle_\alpha=0, \end{align}

(5.1b)\begin{align} & \langle\boldsymbol{v}_\alpha\rangle_\alpha ={-} \frac{{\boldsymbol{\mathsf{H}}}_{\alpha\alpha}}{\mu_{{ref}\alpha}} \boldsymbol{\cdot}(\boldsymbol{\nabla}\langle p_\alpha\rangle^\alpha -\rho_\alpha\boldsymbol{b}_\alpha ) - \frac{{\boldsymbol{\mathsf{H}}}_{\alpha\kappa}}{\mu_{{ref}\kappa}} \boldsymbol{\cdot}(\boldsymbol{\nabla}\langle p_\kappa\rangle^\kappa -\rho_\kappa \boldsymbol{b}_\kappa )\nonumber\\ &\qquad\quad +\frac{1}{\mu_{{ref}\alpha}V}\int_{\mathscr{A}_{\beta\gamma} } ( 2\gamma H \boldsymbol{n}_{\beta\gamma}+\nabla_s \gamma) \boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\alpha\alpha}\, {\rm d} A +\frac{\rho_\gamma -\rho_\beta}{\mu_{{ref}\alpha}V} \int_{\mathscr{A}_{\beta\gamma} } \boldsymbol{n}_{\beta\gamma} \boldsymbol{\cdot}\boldsymbol{w}\boldsymbol{w}\boldsymbol{\cdot}{\boldsymbol{\mathsf{D}}}_{\alpha \alpha}\, {\rm d} A. \end{align}

In addition, the expression that predicts the macroscopic pressure difference is given by

(5.1c)\begin{align} {\rm \Delta}\mathcal{P} &=\langle p_\gamma\rangle^\gamma|_{\boldsymbol{x}_\gamma}-\langle p_\beta\rangle^\beta|_{\boldsymbol{x}_\beta}={-}(\boldsymbol{\varphi}_\beta+\boldsymbol{\varphi}_\gamma+ \boldsymbol{x}_\beta-\boldsymbol{x}_\gamma) \boldsymbol{\cdot}\boldsymbol{\nabla}\langle p_\beta\rangle^\beta \nonumber\\ &\quad +\rho_\beta \boldsymbol{b}_\beta\boldsymbol{\cdot}\boldsymbol{\varphi}_\beta +\rho_\gamma \boldsymbol{b}_\gamma\boldsymbol{\cdot}\boldsymbol{\varphi}_\gamma \nonumber\\ & \quad +s_{\beta\gamma} +(\rho_\gamma-\rho_\beta)\int_{\mathscr{A}_{\beta\gamma}} \boldsymbol{n}_{\beta\gamma}\boldsymbol{\cdot}\boldsymbol{w} \boldsymbol{w}\boldsymbol{\cdot}\boldsymbol{f}_\beta\, {\rm d} A. \end{align}

In fact, the above equations can be simplified under creeping flow conditions or when the fluid–fluid interface is stationary as, under such circumstances, the last interfacial integral term in both (5.1b) and (5.1c) vanishes. Otherwise, knowledge of the displacement velocity of the fluid–fluid interfaces, $\boldsymbol {w}$, is required. The macroscopic equations (5.1) are written in terms of effective-medium quantities (see (3.13) and (4.5)) that can be predicted from the solution of the adjoint closure problems derived in the previous sections, given by (3.4) and (3.5) for the macroscopic momentum balance equation and (4.1) for the average pressure difference. These problems require the location of the fluid–fluid interfaces, as well as the velocity field in each phase in the periodic unit cell, even under creeping flow conditions, due to the apparent viscosity coefficient $\mu (\varGamma _\alpha )$. The validity of these models, as well as the dependence with the rheological parameters, are investigated in § 6 for a particular non-Newtonian viscosity model.

To conclude this section, it is worth recalling that the macroscopic models derived here involve the use of no-slip conditions at the solid–fluid interfaces. If such condition is replaced by a Navier-slip condition and triple-phase contact line effects are disregarded, it is reasonable to propose, from previous works on one-phase flow with slip (cf. Lasseux, Valdés-Parada & Porter 2016), that the structure of macroscopic models remains unchanged as the slip effects play a role only at the closure problems level. In addition, many works have reported on hysteretic effects when studying two-phase flow in porous media. For example, this is the case when the capillary pressure is assumed to only depend on the saturation. In this regard, the macroscopic models reported here clearly show that both macroscale fluid velocities and pressure difference depend on the macroscopic forcings, saturation, shapes and positions of the interfaces, interfacial tension effects as well as inertial effects, when present. This more exhaustive dependence is prone to alleviate the conjectured hysteretic behaviour.

6. Numerical simulations

The objective of this section is twofold. Firstly, average velocities and pressure difference predicted from the upscaled models are compared with the corresponding average quantities resulting from DNS of the pore-scale equations, for validation purposes. Secondly, the influence of the non-Newtonian and inertial character of the flow on the predictions of the effective-medium coefficients and macroscale variables is analysed.

The macroscale models derived in the previous two sections were obtained without assuming any particular geometrical configuration in the periodic unit cell nor any constitutive relationship for the apparent viscosity $\mu (\varGamma _\alpha )$, hence keeping generality. Here, results are reported considering that the $\gamma$ phase is Newtonian and that a Carreau model is applicable for $\mu (\varGamma _\beta )$. This model involves four parameters, namely the infinite and zero shear-rate viscosities $\mu _\infty$ and $\mu _0$, the relaxation time $\lambda$ and the power-law index $n$, and can be expressed as follows:

(6.1)\begin{equation} \mu (\varGamma_\beta) = \mu_\infty + (\mu_0 -\mu_\infty) [1 +(\lambda \varGamma_\beta)^2 ]^{(n-1)/2}. \end{equation}

The particular two-phase flow configurations used for the numerical simulations are sketched in figure 2 and consist of square (figure 2a,b) and random (figure 2c) patterns of solid cylinders with circular cross-section in periodic unit cells of side length $\ell$. In the numerical simulations that follow, the $\beta$-phase saturation is varied in the geometries depicted in figures 2(a) and 2(b) (albeit the porosity is fixed at $\varepsilon =0.8$), whereas, in the case of figure 2(c), $S_\beta = 0.346$ and $\varepsilon =0.915$. Although the models derived in the previous sections are applicable over a wider set of conditions, for the sake of simplicity in the analysis, here the flow is considered to be two-dimensional. In addition, fluid phase distributions are considered such that the (wetting) $\beta$ phase flows under continuous and/or discontinuous forms of layers attached to each solid cylinder (note that the latter case mimics residual fluid trapping), whereas the (non-wetting) $\gamma$ phase flows in the remaining core of the pore space. Moreover, in the case of figure 2(c), the non-wetting fluid is also partially in contact with the solid phase, as would be the case for a mixed-wet medium. The two phases are flowing under the macroscopic forcing given by

(6.2)\begin{equation} \boldsymbol{\nabla}\langle p_\alpha\rangle^\alpha=\frac{\partial\langle p_\alpha\rangle^\alpha}{\partial x}\boldsymbol{e}_x={-}h\boldsymbol{e}_x,\end{equation}

assuming $\boldsymbol {b}_\alpha =\boldsymbol {0}$ and a constant interfacial tension, i.e. $\nabla _s \gamma =\boldsymbol {0}$. Moreover, the analysis comprises both the laminar and creeping flow regimes. In the latter case, the terms involving $\boldsymbol {n}\boldsymbol{\cdot} \boldsymbol {w}$ in the upscaled models can be safely discarded as explained in the previous sections. It is worth mentioning that, in the absence of inertia, the apparent permeability tensors (${\boldsymbol{\mathsf{H}}}_{\alpha \kappa }$) reduce to the permeability tensors ${\boldsymbol{\mathsf{K}}}_{\alpha \kappa }$. In addition, in the case depicted in figure 2(a), for each wetting-phase saturation value, $S_\beta$, the shape and position of the fluid–fluid interface are taken to correspond to the steady pore-scale flow solution when both phases are Newtonian. They were obtained using a boundary integral element method, as reported by Lasseux & Valdés-Parada (2022). Obviously, this interface shape and position cannot be assumed to correspond, in general, to steady-state conditions when the wetting phase is non-Newtonian. Fortunately, as explained above, the upscaled models derived here are applicable under unsteady conditions as discussed in § 2.1. For this reason, the shapes and positions of $\mathscr {A}_{\beta \gamma }$ obtained for Newtonian two-phase flow can also be used in the present analysis in order to test the robustness of the models for situations in which the fluid–fluid interfaces are not stationary. With the same reasoning, the shapes and positions used for the configurations in figure 2(b,c) can be arbitrarily set. Indeed, the unit cell concept allows considering three-dimensional porous media geometries and phase distributions that can be as complex as desired and that would be obtained from laboratory images, albeit without compromising the periodicity condition. Furthermore, some statistical distribution of flow properties can be determined by applying the upscaling approach developed here over a range of specific microscale realisations, and then combining these in different arrangements to extend the model applicability. A thorough investigation of the many geometrical configurations that can be considered together with particular applications, although interesting, lies beyond the scope of this work.

Figure 2. Sketch of the two-dimensional periodic unit cells, of side length $\ell$, corresponding to square (a,b) and random (c) patterns of parallel cylinders of circular cross-section. The wetting phase (blue) can be continuous (a), discontinuous (b) or both (c), whereas the non-wetting phase (white) is flowing continuously in all cases. In (c), both fluids are in contact with the $\sigma$ phase.

For the analysis that follows, it is convenient to reformulate both the microscale and macroscale models, as well as the closure problems, in their dimensionless form as presented in Appendix B. Under this form, the flow and closure problems depend on the porosity, $\varepsilon$, the wetting phase saturation, $S_\beta$, the viscosity ratio, $\mu ^*={\mu _{{ref}\beta }}/{\mu _{{ref}\gamma }}$, the power-law index, $n$, the dimensionless relaxation time, $\lambda ^*=\lambda h\ell /\mu _{{ref}\gamma }$, the Reynolds numbers, $Re_\alpha = \rho _\alpha \ell ^2/(t_{{ref}\alpha } \mu _{{ref}\alpha })$, and the capillary number, $Ca=h \ell ^2/\gamma$.

Following previous works (e.g. Airiau & Bottaro 2020; Sánchez-Vargas et al. 2022), it is assumed that $\mu _\infty =0$. Thus, the remaining degrees of freedom in the Carreau model are $\mu _0$ ($= \mu _{{ref}\beta }$), $n$ and $\lambda$.

Reported values for the power-law index of a variety of shear-thinning fluids range mostly from 0.2 to 0.8, whereas the dimensional relaxation time can vary over several orders of magnitude, from $\boldsymbol {O}(0.001)$ s to $\boldsymbol {O}(100)$ s (Chhabra & Uhlherr 1980; Yasuda, Armstrong & Cohen 1981; Bird, Armstrong & Hassager 1987; Piau, Kissi & Tremblay 1988; Marcovich et al. 2004; Kwon et al. 2008; Liu et al. 2022). In order to reduce the number of degrees of freedom, the analysis is carried out fixing the following parameters in the unit cells depicted in figure 2(a,b): $\varepsilon =0.8$, $\mu ^*=0.1$ and $Ca=10$, while $S_\beta$ is varied from 0.4 up to 0.75, $n$ ranging between 0.4 and 1, and $\lambda ^*$ between 0 and 100. Some examples of fluids with a power-law index of about 0.4 include blood with 65 % haematocrit (Kwon et al. 2008), plant-based biopolymers derived from soy, wheat or pea protein (Hayashi, Hayakawa & Fujio 1993; Klüver & Meyer 2014) and some ultrahigh-viscosity silica nanoparticle suspensions (Liu et al. 2022). Finally, 2 %–5 % alginate solutions in distilled water (Rezende et al. 2009; Belalia & Djelali 2014) are examples of shear-thinning fluids with $n=0.8$. In addition, $\lambda ^*$ depends on the dimensional value of the relaxation time, and on the magnitude of the applied pressure gradient. Therefore, $\lambda ^*=100$ represents a fluid with a strong non-Newtonian character, subject to large pressure gradients as can be encountered, for instance, in polymer and resin injection (Meng 2019; Liu et al. 2022).

To carry out the closure problem solutions, as well as the DNS, the finite-element software Comsol Multiphysics 6.1 was employed, using a direct PARDISO solver. In addition, typical mesh refinement tests were performed in order to ensure that the results were independent of this degree of freedom. The numerical simulations were carried out in two steps. The first one consisted of the computation of the pore-scale velocity and pressure fields in a unit cell. These results allow, on the one hand, computing the dimensionless average velocity and average pressure difference, which constitute the DNS results. On the other hand, the wetting-phase viscosity becomes available, together with the velocity fields, $\boldsymbol {v}_\alpha$ ($\alpha =\beta,\gamma$), within each phase and at the fluid–fluid interface (i.e. $\boldsymbol {w}$) when inertia is present, for the closure problem solution carried out in the second step. In the following sections, an analysis of the prediction of the effective parameters and the validation of the macroscopic models are proposed for both momentum (§ 6.1) and average pressure difference (§ 6.2).

6.1. Analysis of the macroscopic velocity

This part of the work is dedicated first to the study of the influence of rheological parameters under creeping flow conditions for the configuration in which the $\beta$ phase is continuous (see figure 2a). The use of creeping flow conditions in the first set of simulations is justified both to contrast with previously reported results (Lasseux & Valdés-Parada 2022), and also to isolate the effects of rheology. Secondly, inertial effects are studied in the three configurations depicted in figure 2. In all cases, emphasis is laid upon validation with DNS results.

6.1.1. Influence of the rheological parameters in continuous fluid phases

In figure 3, a first set of results corresponding to the pore-scale flow streamlines, superimposed to the velocity magnitude, are represented, taking $n=0.4$, $\lambda ^*=0.5$ and four values of $S_\beta$ ranging from 0.4 to 0.7. For the Newtonian case (figure 3a,c,e,g), the results reproduce those reported in figure 5 of Lasseux & Valdés-Parada (2022) and are recalled here for comparison purposes. As expected, $\mathscr {A}_{\beta \gamma }$ coincides with a streamline in that case. This feature is not observed in figure 3(b,d, f,h), which correspond to a situation where the $\beta$ phase obeys a Carreau model. In addition, the viscosity parameters employed here, corresponding to a shear-thinning wetting phase, lead to a very significant increase of the velocity by more than an order of magnitude in both phases in comparison with the case in which the $\beta$ phase is Newtonian. The change in the flow structure observed when the $\beta$ phase is non-Newtonian results from both the change in rheology and the fact that the fluid–fluid interface does not correspond to its steady configuration. In all cases, the largest wetting-phase velocity is located in the constriction between the solid cylinder and $\mathscr {A}_{\beta \gamma }$, which corresponds to a low-pressure zone. In contrast, the portions of the wetting phase where eddies are present correspond to those where the viscosity is the largest and correlate to large-pressure zones. Consequently, as the saturation increases, the low-pressure zones increase in size, thus favouring the overall fluid flow. The same effect is observed when $\lambda ^*$ increases, or when $n$ decreases, which is consistent with the shear-thinning nature of the Carreau fluid considered here.

Figure 3. Fields of the velocity magnitude and streamlines (in white) under Newtonian (a,c,e,g) and non-Newtonian (b,d, f,h) flow conditions for the wetting $\beta$ phase, taking four values of the $\beta$-phase saturation ($S_\beta$): (a,b) $S_\beta =0.4$; (c,d) $S_\beta =0.5$; (e, f) $S_\beta =0.6$; (g,h) $S_\beta =0.7$. The fluid–fluid interface is represented as a black curve and the central solid cylinder in white. Simulations for non-Newtonian flow (b,d, f,h) correspond to a Carreau fluid considering the power-law index $n=0.4$ and the dimensionless relaxation time $\lambda ^*=5$. Due to symmetry, results are only reported for the top half of the unit cell depicted in figure 2(a). Other parameters: $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1$.

At this point, it is worth pondering the contribution of the different terms present in the dimensionless macroscopic flow model as expressed in (B7). In table 1, data are reported to evaluate both the Darcy-like part and the interfacial part of $\langle v_{\alpha x}^*\rangle _\alpha$ predicted by this model, taking $n=0.4$ and $\lambda ^*=5$, for the same wetting-phase saturation values considered in figure 3. These predictions yield relative errors with respect to DNS of less than 1 % for all the values of $S_\beta$ considered here. Moreover, the Darcy-like terms are found to be, at least, one order of magnitude larger than the interfacial terms for the set of parameters considered here. Furthermore, since, for each given saturation, the macroscopic pressure gradient is along $\boldsymbol {e}_x$ (see (6.2)), $\langle v_{\alpha y}^*\rangle _\alpha$ is numerically verified to be zero, so, from this point forth, focus is laid upon $\langle v_{\alpha x}^*\rangle _\alpha$. Results in table 1 also show that the contributions from $K^*_{\beta \beta xx}$, $K^*_{\beta \gamma xx}$ and $K^*_{\gamma \beta xx}$ increase with $S_\beta$. However, for $K^*_{\gamma \gamma xx}$ the same observation is not applicable and this needs to be further analysed. To address this point, figures 4 and 5 report the predictions of the permeability coefficients and the average velocities for different wetting-phase saturations, while varying $n$ and $\lambda ^*$. In these figures, results are presented in terms of the $xx$ component $k_{r\alpha \kappa }$ of the relative permeability tensor, defined as

(6.3)\begin{equation} k_{r\alpha\kappa}=\frac{K_{\alpha\kappa xx}}{K}, \quad \alpha,\kappa=\beta,\gamma, \end{equation}

where ${K}_{\alpha \kappa xx}$ is the $xx$ component of the permeability tensor ${\boldsymbol{\mathsf{K}}}_{\alpha \kappa }$ and $K$ is the intrinsic permeability, which is a scalar since, for the case under consideration, the structure is orthotropic.

Table 1. Predictions of the macroscopic flow model for different saturations of the non-Newtonian wetting phase, $S_\beta$. The error is defined as $\text {Error}_\alpha ~(\%)= 100 \times |\langle v_{\alpha x}\rangle _{\alpha DNS}-\langle v_{\alpha x}\rangle _\alpha |/\langle v_{\alpha x}\rangle _{\alpha DNS}$. Here, $K^*_{\alpha \kappa xx}$ stands for the xx component of the corresponding dimensionless permeability tensor and $\alpha _\alpha ^*={2\mu _{{ref}\gamma }}/({\mu _{{ref}\alpha } Ca V^*})\int _{\mathscr {A}_{\beta \gamma } } H^* \boldsymbol {n}_{\beta \gamma }\boldsymbol{\cdot} {\boldsymbol{\mathsf{D}}}^*_{\alpha \alpha }\, {\rm d} A^*\boldsymbol{\cdot} \boldsymbol {e}_x$. Also, $\varepsilon =0.8$, $Ca=10$, $\mu ^*=0.1$, $n=0.4$ and $\lambda ^*=5$. The results correspond to the geometry depicted in figure 2(a).

Figure 4. The $x$ component of the dimensionless average velocity in the $\beta$ phase (a) and in the $\gamma$ phase(b) versus the wetting-phase saturation $S_\beta$. (c–f) The xx component of the dominant ($k_{r\alpha \alpha }$) and coupling ($k_{r\alpha \kappa }$, $\alpha \neq \kappa$) relative permeabilities versus $S_\beta$, taking the intrinsic permeability as the reference. Results correspond to a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase with $\lambda ^*=5$. $\bullet$ (red), Newtonian case ($n=1$); $\blacklozenge$, $n=0.8$; $\bigstar$, $n=0.7$; $\blacktriangledown$, $n=0.6$; $\blacktriangleright$, $n=0.5$; $\blacksquare$, $n=0.4$. In (a,b), results from DNS and the upscaled model are reported with open circles and plain symbols, respectively. Other parameters: $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1$. The results correspond to the geometry depicted in figure 2(a).

Figure 5. The $x$ component of the dimensionless average velocity in the $\beta$ phase (a) and in the $\gamma$ phase(b) versus the wetting-phase saturation $S_\beta$. (c–f) The xx component of the dominant ($k_{r\alpha \alpha }$) and coupling ($k_{r\alpha \kappa }$, $\alpha \neq \kappa$) relative permeabilities versus $S_\beta$, taking the intrinsic permeability as the reference. Results correspond to a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase with $n=0.5$. $\bullet$ (red), Newtonian case ($\lambda ^*=0$); $\blacktriangle$, $\lambda ^*=1$; $\blacksquare$, $\lambda ^*=5$; $\blacklozenge$, $\lambda ^*=25$; $\blacktriangledown$, $\lambda ^*=50$; $\bigstar$, $\lambda ^*=100$. In (a,b), results from DNS and the upscaled model are reported with open circles and plain symbols, respectively. Other parameters: $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1$. The results correspond to the geometry depicted in figure 2(a).

In figure 4, an analysis is presented fixing $\lambda ^*$ and varying the power-law index $n$ in a range of values representative of different non-Newtonian fluids. Regarding these results, the following remarks are in order.

• • From figure 4(a,b), it is clear that the values of the x component of the velocity, resulting from averaging the pore-scale simulation fields and predicted from solving the upscaled model, show excellent agreement for both fluid phases. In fact, the relative error of the upscaled model predictions with respect to those from pore-scale simulations is less than 5 % in all cases. This comment is also applicable to the results reported in figure 5(a,b), for $0\leq \lambda ^*\leq 100$.

• • Clearly, as the non-Newtonian character of the wetting phase is more pronounced, the average flow in both phases is favoured, in agreement with the shear-thinning model under consideration. The increment in both phases is very important as it can reach two orders of magnitude with respect to the Newtonian case. This is attributed to the increase of the relative permeability xx component and this is in agreement with the observations made in figure 3.

• • As $n$ decreases, the relative permeabilities and the average velocities become more sensitive to its variations. This is more evident for $n<0.5$.

• • The graphs in figure 4(e, f), describing the xx component of the coupling relative permeabilities, are related by $k_{r\gamma \beta }= \mu ^* k_{r\beta \gamma }$. This proportionality can be deduced from (A10) under creeping flow conditions.

The second degree of freedom in this analysis is the dimensionless relaxation time, $\lambda ^*$. Results corresponding to the sensitivity of macroscale fluid velocities and relative permeabilities to this parameter are reported in figure 5 for several values of the wetting-phase saturation. Regarding these results, the following observations are pertinent.

• • Consistent with the results shown in figure 4, flow within the unit cell is improved as the shear-thinning character of the wetting phase is increased. Once more, both the relative permeabilities and macroscale velocities increase by roughly two orders of magnitude as $\lambda ^*$ increases from 1 to 100. Note that the flow is more sensitive to variations of $\lambda ^*$ in the range $0\le \lambda ^*\le 5$.

• • Similar to the analysis of the impact of $n$, it must be noted that $k_{r\beta \gamma }/k_{r\gamma \beta }=10$ (see figure 5e, f) and the similarity with figure 4(e, f). This confirms the expected reciprocity relationship reported in § 3.3 (see also (A10)), in agreement with the fact that $\mu ^*=0.1$.

The above analysis shows the important role played by both the dominant and coupling Darcy-like terms. For the specific combination of values of the capillary number and viscosity ratio, together with the porous medium geometry and interface configuration chosen here, there is no relevant contribution from surface tension effects (i.e. $\boldsymbol {\alpha }^*_\alpha$, cf. (B8a)) to the prediction of the macroscopic fluid velocities. However, this should not be considered as a general result as the contribution of this interfacial term was shown to become even dominant in some circumstances, in particular while decreasing $Ca$ or increasing $\mu ^*$ (see table 1 in Lasseux & Valdés-Parada (2022)). The contribution of this term, along with the inertial interfacial term, $\boldsymbol {\omega }^*_\alpha$ (cf. (B8)), is further explored next when inertia is present.

6.1.2. Influence of inertial effects

Let now attention be directed to the influence of inertial effects over the predictions of the macroscale velocity and effective-medium coefficients. The results presented below require specifying the Reynolds number values, $Re_\alpha$, $\alpha =\beta,\gamma$, in both fluid phases. For this reason, $Re_\beta$ was arbitrarily set equal to one-tenth of $Re_\gamma$. Despite the fact that in the case of figure 2(c) the absence of symmetry of the unit cell makes the apparent permeability tensors fully populated, in the following, only the $xx$ components of these tensors are reported, normalised with the $xx$ component of the intrinsic permeability. This is,

(6.4)\begin{equation} h_{r\alpha \kappa} = \frac{H_{\alpha \kappa xx}}{K}, \quad \alpha, \kappa = \beta, \gamma. \end{equation}

Accordingly, only the $x$ components of the average velocities are reported.

In figures 6 and 7, the predictions of the average velocities and relative permeabilities versus $S_\beta$, taking four values of $Re_\gamma$, resulting from solving the closure problems in the geometries depicted in figures 2(a) and 2(b), are respectively reported for $n=0.5$ and $\lambda ^*=5$ in the Carreau model. From these results, it is clear that the influence of the rheology, combined with inertia, yields a complex dependence of the average velocities and effective-medium coefficients on $S_\beta$ for the Reynolds number values considered here. An important remark must be made for the situation in which the $\beta$ phase is trapped (see figure 2b). It is worth noticing in that case, that $\langle v^*_{\beta x}\rangle _\beta$ is not zero, despite the fact that it is trapped. However, this is not surprising as this is a signature of the fact that the interface is not at its steady position (i.e. the interface position can be thought of as being captured at a given time at which this phase is rearranging, before equilibrium is reached). Therefore, there is a net local momentum transport in the $\beta$ phase in the $x$ direction. This should not be confused with the fact that, in this configuration, the $\beta$ phase remains macroscopically trapped in the porous medium, hence yielding no mass flow rate at the macroscopic system boundaries.

Figure 6. The $x$ component of the dimensionless average velocity in the $\beta$ phase (a) and in the $\gamma$ phase(b) versus the wetting-phase saturation $S_\beta$. (c–f) The xx component of the dominant ($h_{r\alpha \alpha }$) and coupling ($h_{r\alpha \kappa }$, $\alpha \neq \kappa$) apparent relative permeabilities versus $S_\beta$, taking the $xx$ component of the intrinsic permeability tensor as the reference. Results correspond to a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase with $n=0.5$ and $\lambda ^*=5$ for the geometry depicted in figure 2(a) in which the $\beta$ phase is continuous. $\bullet$ (red), creeping flow regime ($Re_\gamma =0$); $\blacktriangle$, $Re_\gamma =10$; $\blacktriangledown$, $Re_\gamma =50$; $\blacklozenge$ (grey), $Re_\gamma =100$. In (a,b), results from DNS and the upscaled model are reported with open circles and plain symbols, respectively. Other parameters: $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1$.

Figure 7. The $x$ component of the dimensionless average velocity in the $\beta$ phase (a) and in the $\gamma$ phase(b) versus the wetting-phase saturation $S_\beta$. (c–f) The xx component of the dominant ($h_{r\alpha \alpha }$) and coupling ($h_{r\alpha \kappa }$, $\alpha \neq \kappa$) apparent relative permeabilities versus $S_\beta$, taking the $xx$ component of the intrinsic permeability tensor as the reference. Results correspond to a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase with $n=0.5$ and $\lambda ^*=5$ for the geometry depicted in figure 2(b) in which the $\beta$ phase is discontinuous. $\bullet$ (red), creeping flow regime ($Re_\gamma =0$); $\blacktriangle$, $Re_\gamma =10$; $\blacktriangledown$, $Re_\gamma =50$; $\blacklozenge$ (grey), $Re_\gamma =100$. In (a,b), results from DNS and the upscaled model are reported with open circles and plain symbols, respectively. Other parameters: $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1$.

The effects of $\alpha ^*_{\kappa x}$ are not reported here since they are negligible when the $\beta$ phase is continuous (cf. figure 2a), and they are exactly zero when it is discontinuous (cf. figure 2b), since $H$ is constant in that case. Note that the dependence of the average velocities and $h^*_{r\alpha \alpha }$ ($\alpha =\beta,\gamma$) on the Reynolds number is not monotonic over the reported range of saturations. This can be attributed to pore-scale balance between viscous dissipation and inertial acceleration in both fluid phases and to the fact that the fluid–fluid interface is not at its stationary position. Nevertheless, as expected, when the $\beta$ phase is continuous, the values of both the average velocities and $h^*_{r\alpha \alpha }$ are larger than those corresponding to the case in which the $\beta$ phase is trapped. This contrast is more pronounced in the $\gamma$ phase, as a result of the fact that, when the $\beta$ phase is trapped, the wetting phase exerts an extra drag on the surrounding non-wetting phase flow. Finally, for $h^*_{r\alpha \kappa }$ ($\alpha, \kappa =\beta,\gamma$, $\alpha \ne \kappa$), the simple reciprocity relationship observed in the previous section is no longer valid, especially when the $\beta$ phase is discontinuous (see figure 7e, f). This is because the terms on the right-hand side of (A10) are no longer vanishingly small compared with their counterparts on the left-hand side when inertia is present. Despite this, the predictions for $h_{r\beta \gamma }$ remain one order of magnitude larger than those for $h_{r\gamma \beta }$ as was observed in the previous paragraphs under creeping flow conditions.

The last part of the analysis of inertial effects is focused on the configuration depicted in figure 2(c). In this case, $\varepsilon =0.915$, $S_\beta =0.346$, $\mu ^*=0.1$, $n=0.5$, $\lambda ^*=5$ and $Ca=10$, so the results for the velocity and the relative (dominant and coupling) apparent permeabilities reported in figure 8 are presented as functions of the Reynolds number in the $\gamma$ phase. These results show that larger values of $Re_\gamma$ are required to appreciate the effect of inertia, in contrast to the results for the configurations in figures 2(a) and 2(b). This is attributed to the viscous drag exerted by the solid on both fluid phases, and this was not the case in the other geometries for which the $\gamma$ phase experienced lubrication effects. The predictions for $\langle v_{\alpha x}^*\rangle _\alpha$ are qualitatively alike in both fluid phases, although $\langle v_{\gamma x}^*\rangle _\gamma$ and $h_{r\gamma \gamma }$ are one order of magnitude larger than their counterparts in the $\beta$ phase. The dependence of $\langle v_{\alpha x}^*\rangle _\alpha$ on the Reynolds number is consistent with observations made in inertial one-phase flow in porous media (cf. Lasseux, Arani & Ahmadi 2011). Finally, the coupling permeability terms maintain the same properties observed in the previous geometries, i.e. $h_{r\beta \gamma }$ is one order of magnitude larger than $h_{r\gamma \beta }$, albeit their dependence upon the Reynolds number is not the same. It is worth adding that the results presented in figures 6–8 show an excellent agreement between DNS and the upscaled model predictions for the average fluid velocities, thus validating the upscaled model in the presence of inertia.

Figure 8. The $x$ component of the dimensionless average velocity in the $\beta$ phase (a) and in the $\gamma$ phase(b) versus $Re_\gamma$. (c–f) The xx component of the dominant ($h_{r\alpha \alpha }$) and coupling ($h_{r\alpha \kappa }$, $\alpha \neq \kappa$) apparent relative permeabilities versus $Re_\gamma$, taking the $xx$ component of the intrinsic permeability as the reference. Results correspond to $Ca=10$ and $\mu ^*=0.1$ for a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase with $n=0.5$ and $\lambda ^*=5$ for the geometry depicted in figure 2(c) in which the $\beta$ phase is both continuous and discontinuous and both phases are in contact with the solid phase $\sigma$. In (a,b), results from DNS and the upscaled model are reported with open and filled circles, respectively. Other parameters are $\varepsilon=0.915$ and $S_\beta=0.346$.

To conclude this part of the analysis, it is pertinent to examine the contributions to $\langle v_{\kappa x}^*\rangle _\kappa$ ($\kappa =\beta,\gamma$) from the $x$ component of the interfacial terms $\boldsymbol {\alpha }_\kappa ^*$ and $\boldsymbol {\omega }_\kappa ^*$ ($\kappa =\beta,\gamma$). In figure 9, these contributions are reported normalised with $\langle v_{\kappa x}^*\rangle _\kappa$ ($\kappa =\beta,\gamma$) in both fluid phases against $Re_\gamma$. These results show that, for the particular configuration under study, the contributions from $\boldsymbol {\omega }_\kappa ^*$ can be neglected for $Re_\gamma \le 10$ with respect to those from $\boldsymbol {\alpha }_\kappa ^*$. Nevertheless, $\boldsymbol {\omega }_\kappa ^*$ becomes especially relevant for $Re_\gamma >10^3$ as it can be roughly equal to the macroscopic velocity values of both fluid phases. In contrast, the contributions from $\boldsymbol {\alpha }_\kappa ^*$, representing between 3 % and 10 % of the average velocities, remain relevant in the entire range of Reynolds number values considered here.

Figure 9. Per cent contribution to $\langle v_{\kappa x}^*\rangle _\kappa$ ($\kappa =\beta,\gamma$) of the interfacial term $x$ component defined in (B8) in (a) the $\beta$ phase and (b) the $\gamma$ phase (i.e. $\psi _\kappa ^{\%} = 100 \times | \psi _{\kappa x}^*|/\langle v_{\kappa x}^*\rangle _\kappa$, $\psi =\alpha,\omega$, $\kappa = \beta,\gamma$) as functions of the Reynolds number in the $\gamma$ phase ($Re_\gamma$). Results correspond to $Ca=10$ and $\mu ^*=0.1$ for a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase with $n=0.5$ and $\lambda ^*=5$ for the geometry depicted in figure 2(c). Other parameters are $\varepsilon=0.915$ and $S_\beta=0.346$.

6.2. Analysis of the average pressure difference

An analysis similar to that reported above for the average velocities is now performed for the macroscopic pressure difference model, as given in (B10). For the sake of brevity, only the case depicted in figure 2(a) in the creeping regime is considered, keeping $\nabla ^*\langle p^*_\alpha \rangle ^\alpha$ as given by (B2). In figure 10(a), the predictions of ${\rm \Delta} \mathcal {P}^*$ from the dimensionless upscaled equation (B10) in which Reynolds numbers are zero are compared with the DNS results, ${\rm \Delta} \mathcal {P}^*_{DNS}$, for different wetting-phase saturations $S_\beta$ and values of the power-law index $n$ ranging from 0.4 to 1.0, fixing $\lambda ^*=5$. Clearly, the derived model predictions perfectly match the pore-scale numerical simulations. In all cases, the error %, calculated as $100 \times |{\rm \Delta} \mathcal {P}^*_{DNS}-{\rm \Delta} \mathcal {P}^*|/ {\rm \Delta} \mathcal {P}^*_{DNS}$, remains smaller than 5 %. This is also true for the results reported in figure 10(b) corresponding to values of the dimensionless relaxation time $\lambda ^*$ ranging from 0 to 50, keeping $n=0.5$. Results for ${\rm \Delta} \mathcal {P}^*$ for larger values of $\lambda ^*$ do not significantly change, and, therefore, are not reported in figure 10(b).

Figure 10. Dimensionless macroscopic pressure difference ${\rm \Delta} \mathcal {P}^*$ versus the wetting-phase saturation $S_\beta$ resulting from the upscaled model (filled symbols) and from DNS (open circles). Results correspond to the geometry depicted in figure 2(a) with a Newtonian $\gamma$ phase and a Carreau fluid for the $\beta$ phase. In (a), $\lambda ^*=5$ and $\bullet$ (red), Newtonian case ($n=1$); $\blacklozenge$, $n=0.8$; $\bigstar$, $n=0.7$; $\blacktriangledown$, $n=0.6$; $\blacktriangleright$, $n=0.5$; $\blacksquare$, $n=0.4$. In (b) $n=0.5$ and $\bullet$ (red), Newtonian case ($\lambda ^*=0$); $\blacktriangle$, $\lambda ^*=1$; $\blacksquare$, $\lambda ^*=5$; $\blacklozenge$, $\lambda ^*=25$; $\blacktriangledown$, $\lambda ^*=50$. Other parameters: $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1.$

Increasing the shear-thinning effect of the $\beta$ phase, by either decreasing $n$ or increasing $\lambda ^*$, yields a decrease of the macroscopic pressure difference as shown in figures 10(a) and 10(b), respectively. This is particularly noticeable for $S_\beta <0.6$, whereas, for larger values of the saturation, the impact of the wetting-phase rheology is almost insignificant. Moreover, in agreement with the velocity fields reported in figure 3, it is verified that the average pressure in each phase decreases while increasing the wetting-phase saturation. This is consistent with the fact that the size of the recirculation zones in the wetting phase decreases with $S_\beta$. Nevertheless, the above does not necessarily explain why the macroscopic pressure difference tends to zero as the wetting-phase saturation increases. To address this issue, it is worth recalling that ${\rm \Delta} \mathcal {P}^*$ results from the contribution of both dynamic (i.e. the macroscopic forcing) and capillary effects, represented by $\psi ^*$ and $s_{\beta \gamma }^*$, respectively. Therefore, the impact of the rheological parameters $n$ and $\lambda ^*$ on $\psi ^*$ and $s_{\beta \gamma }^*$, representing the (dimensionless) dynamic and capillary contributions, respectively given by (B11a) and (B11b), is now analysed.

The effective quantities $\psi ^*$ and $s_{\beta \gamma }^*$, obtained from the closure problem solution with $n$ ranging from 1 to 0.4, while fixing $\lambda ^*=5$, are reported versus $S_\beta$ in figures 11(a) and 11(b), whereas, in figures 11(c) and 11(d), the same representation is provided for results obtained by varying $\lambda ^*$ from 0 to 50, fixing $n=0.5$. From figures 11(a) and 11(c), it can be noted that the contribution of the dynamic effects $\psi ^*$ to ${\rm \Delta} \mathcal {P}^*$ is always negative, contrary to the dimensionless capillary effect contribution $s_{\beta \gamma }^*$ (see figures 11b and 11d). From these figures, it is clear that both $\psi ^*$ and $s_{\beta \gamma }^*$ tend to zero as $S_\beta$ increases. This results in a macroscopic pressure difference that also tends to zero when $S_\beta$ takes large enough values (i.e. $>0.6$). In addition, results in figures 11(a) and 11(c) show that $\psi ^*$ increases in magnitude when the shear-thinning effect in the $\beta$ phase becomes more pronounced, whereas the opposite holds for $s_{\beta \gamma }^*$ (see figures 11b and 11d). It must be noted that both dynamic and capillary effects play a significant role in the prediction of the average pressure difference for the conditions under consideration here. This is in contrast with what was concluded for the average velocities for which the interfacial term is insignificant under the same conditions.

Figure 11. Dynamic (i.e. macroscopic forcing) $\psi ^*=\varphi ^*_{\beta x}+\varphi ^*_{\gamma x}+x^*_\beta -x^*_\gamma$ and capillary effect $s_{\beta \gamma }^*$ contributions versus the wetting-phase saturation $S_\beta$. In (a,b) $\lambda ^*=5$ and $\bullet$ (red), Newtonian case ($n=1$); $\blacklozenge$, $n=0.8$; $\bigstar$, $n=0.7$; $\blacktriangledown$, $n=0.6$; $\blacktriangleright$, $n=0.5$; $\blacksquare$, $n=0.4$. In (c,d) $n=0.5$ and $\bullet$ (red), Newtonian case ($\lambda ^*=0$); $\blacktriangle$, $\lambda ^*=1$; $\blacksquare$, $\lambda ^*=5$; $\blacklozenge$, $\lambda ^*=25$; $\blacktriangledown$, $\lambda ^*=50$. In all cases, $\varepsilon =0.8$, $Ca=10$ and $\mu ^*=0.1$ in the geometrical configuration of figure 2(a).

To conclude this section, it is worth remarking that the main objective of the present work is the derivation of the macroscopic models summarised in § 5. Therefore, an exhaustive analysis of all the degrees of freedom involved in the closure problems solution was not carried out. Indeed, the numerical analysis presented here is only illustrative, as the upscaled models are applicable under broader flow conditions. Moreover, rather than reproducing the conditions of a particular experimental system, the interest lies here in analysing the contributions of each term constituting the upscaled models. Importantly, in all the situations considered in this work, the upscaled models for both the average velocities and average pressure difference are validated by the pore-scale direct numerical simulations. It must be emphasised that this holds even if the fluid–fluid interface is not in its stationary configuration, if it is continuous or not through the unit cell and in potentially mixed-wet conditions. This should be regarded as a first validation step and motivates further comparisons with laboratory experiments.

7. Conclusions

Closed upscaled models for generalised Newtonian two-phase, quasi-steady, incompressible and inertial fluid flow in rigid and homogeneous porous media, including possible interfacial tension gradient (i.e. Marangoni) effects, are formally derived from the governing equations at the pore scale employing a simplified volume-averaging technique coupled to an adjoint method and Green's formulation. No other assumptions than those classically adopted for upscaling, i.e. separation between the microscopic and macroscopic characteristic length scales and existence of a representative elementary volume, are required in the derivation. The models consist of macroscopic mass and momentum balance equations in each phase and an expression of the average pressure difference. The macroscopic mass balance equations are identical to those already reported for the Newtonian case. For momentum, the averaged velocities are given by the classical generalised Darcy's law that takes inertia and viscous coupling into account, but that is complemented by additional terms resulting from the upscaling process. Indeed, the average velocities are shown to depend on the intrinsic average pressure gradient and body force in each phase through Darcy-like terms expressed with dominant and coupling apparent (i.e. inertia-dependent) permeability tensors. The former are symmetric while the latter satisfy a reciprocity relationship. Moreover, two additional terms are involved, including a contribution from interfacial tension effects and another one, accounting for interfacial inertia, that vanishes when fluid–fluid interfaces are in their stationary position. All the effective quantities involved in the expression of the macroscopic velocities are evaluated from the solution of a pair of (adjoint) closure problems on periodic unit cells representative of the microscale process and microstructure of the porous medium. Similarly, the average pressure difference is shown to include terms that involve the macroscopic forcing (i.e. the macroscopic pressure gradient and body force) in each phase, interfacial tension and inertia, the latter vanishing when fluid–fluid interfaces are stationary. The effective quantities for this model are determined from the solution of another closure problem on the periodic unit cell.

Numerical tests were performed on model two-dimensional configurations considering a Newtonian non-wetting phase and a shear-thinning Carreau fluid for the wetting phase. Many different situations were envisaged, for which the wetting phase is either continuous or not, a combination of the two and a mixed-wet type of system where both fluid phases are in contact with the solid phase. The effects of the (wetting-phase) saturation, rheology parameters in the Carreau model and Reynolds numbers were analysed. In all the situations, excellent agreement was achieved between the average velocities and pressure difference obtained from DNS on the one-hand, and predicted by the macroscopic models, after solving the closure problems, on the other-hand, hence validating the models derived here. Indeed, in all the situations, the maximum error is less than 5 %, confirming the performance of the average models. All the effective quantities, and consequently the macroscopic velocity in each phase and the average pressure difference, were shown to be very sensitive to changes in the wetting-phase rheology and the presence or not of inertia. The interfacial terms, resulting from capillary effects and inertia, present in both the macroscopic momentum balance equations and pressure difference relationship, were shown to depend on the rheology and flow conditions and their contributions in the macroscale models were clearly identified and analysed.

To the best of our knowledge, the derivation of these two models in the context of non-Newtonian two-phase inertial flow in porous media represents a novelty. A salient feature lies in the fact that, in addition to the formal derivation, the roles played by the different sources, originating from viscous, interfacial tension and inertial effects, are clearly identified and can be evaluated from the adjoint closure problems. In particular, for a given porous structure, all the effective quantities are expected to depend on the porosity, Reynolds number in each phase, capillary number, (wetting-phase) saturation, ratio of the reference viscosities, rheology in each phase and pressure gradient orientation and, importantly, on the configuration of the fluid–fluid interfaces. The locality of closure problems, which is required to define the effective coefficients involved in the average equations, relies on the identification of a representative periodic unit cell, a concept that is classical during upscaling. Nevertheless, this is not prone to represent a major limitation in the validity of the macroscopic models. However, in the case of a frontal displacement of one phase by another, more attention would have to be paid in the vicinity of the macroscopic front since, in this region, separation of length scales may fail. This situation and other physical conditions remain to be further studied in future works.