Darcy's law survival from no-slip to perfect-slip flow in porous media

Abstract

A macroscopic model for perfect-slip flow in porous media is derived in this work, starting from the pore-scale flow problem and making use of an upscaling technique based on the adjoint method and Green's formula. It is shown that the averaged momentum equation has a Darcy form in which the permeability tensor, $\boldsymbol{\mathsf{K}}_{ps}$, is obtained from an associated adjoint (closure) problem that is to be solved on a (periodic) unit cell representative of the structure. Similarly to the classical permeability tensor, $\boldsymbol{\mathsf{K}}$, in the no-slip regime, $\boldsymbol{\mathsf{K}}_{ps}$ is intrinsic to the porous medium under consideration and is shown to be symmetric and positive. Integral relationships between $\boldsymbol{\mathsf{K}}_{ps}$, the partial-slip flow permeability tensor, $\boldsymbol{\mathsf{K}}_{s}$, and $\boldsymbol{\mathsf{K}}$ are derived. Numerical simulations carried out on two-dimensional model porous structures, together with an approximate analytical solution and an empirical correlation for a particular configuration, confirm the validity of the macroscopic model and the relationship between $\boldsymbol{\mathsf{K}}_{ps}$ and $\boldsymbol{\mathsf{K}}$.

1. Introduction

Slip flow induced by nano- or micro-rough patterned surfaces is of interest for many applications in microfluidics, flow of polymers as encountered in composite manufacturing (Trochu et al. 2006) or measurements of polymer phase properties (Yang 1998). Typically, slip can be conceived as the result of the interaction of the flowing fluid and a heterogeneous surface composed of patches where perfect slip (i.e. when the fluid experiences no shear at the interface) can be reasonably assumed whereas no slip occurs in the complementary part of the surface (Lauga & Stone 2003; Lauga, Brenner & Stone 2007; Asmolov & Vinogradova 2012). The perfect-slip zones result, for instance, from gas trapping due to capillary effects when a liquid flows over the surface, whereas no-slip zones correspond to the region where solid–liquid contact occurs. This is typically the mechanism at play on the so-called superhydrophobic surfaces (Vinogradova 1999) that are of major interest in interfacially driven flows in engineering applications, in particular for electro-osmosis (Bocquet & Barrat 2007). Evidently, the overall shear exerted on the flowing fluid is expected to be a decreasing function of the perfect-slip surface fraction down to the limit where this fraction tends to unity.

Implication of these mechanisms in the process of polymer injection in a fibrous preform during liquid composite moulding, for instance, is of major concern as close to perfect slip conditions would ease polymer infusion in the porous material (Trochu et al. 2006). In another context, the analysis of He superfluid flow in porous media would also imply considering a perfect-slip flow condition for the superfluid component (Allain et al. 2010). Whereas partial- and no-slip flows in porous media have been thoroughly analysed, sometimes in a different context involving Knudsen effects (see e.g. Lasseux, Valdés-Parada & Porter 2016), the perfect-slip limit has not been addressed from a formal upscaling point of view. Recently, Geoffre et al. (2021) performed a series of numerical flow experiments in the transverse direction of two-dimensional random arrays of circular obstacles and proposed an empirical relationship between the permeability corresponding to partial slip and those for perfect- and no-slip conditions. This relationship generalises the one that can be formally obtained in the case of mono-disperse parallel cylinders of circular cross-section in the limit of large enough porosity. However, at the moment, it is unclear whether a Darcy-like model is applicable to the perfect-slip flow situation and, if this is the case, there are no direct means for predicting the perfect-slip permeability for an arbitrary given structure.

The purpose of this work is to address the above fundamental questions by deriving the macroscopic momentum balance equation starting from the pore-scale problem for Newtonian creeping flow with a zero-shear condition at the solid–fluid interface. To this end, the work is organised as follows. The perfect-slip flow problem statement at the pore scale as well as the starting assumptions are presented in § 2. Section 3 is dedicated to the upscaling process that is performed by proposing an adjoint problem, followed by the use of Green's formula. The result of upscaling is a Darcy-like model for perfect slip in which the corresponding permeability tensor is determined by the solution of the adjoint problem in a periodic unit cell representative of the structure. This tensor is shown to be intrinsic, i.e. to depend only on the porous structure, and to be symmetric positive (see Appendix A). Relationships between the perfect-slip, partial-slip and no-slip permeability tensors are reported in § 4, showing that the perfect-slip permeability corresponds to the limit of the partial one in the limit of infinite slip length. This is in agreement with the analysis reported in Mácha & Tichý (2014) at the pore scale. In § 5, the macroscopic model is validated through direct numerical simulations and an approximate analytical solution is proposed for sufficiently large porosity values in a simple geometrical configuration, which appropriately reproduces the numerical results. In addition, the role played by the geometry is examined by comparing the results of the transition from no slip to perfect slip in three unit cell configurations. Finally, the corresponding conclusions are provided in § 6.

2. Pore-scale flow formulation

Consider a homogeneous porous medium whose skeleton is composed of a non-deformable solid phase $\sigma$, the pores, of characteristic size $\ell _\beta$, being saturated by a fluid phase $\beta$ having a constant dynamic viscosity $\mu$. In the absence of body forces, the pore-scale incompressible, creeping and steady flow problem for perfect slip in a periodic unit cell, $\mathscr {V}$, representative of the medium (see figure 1) can be stated as follows:

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

(2.1b)$$\begin{gather}\boldsymbol{0} = \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\tilde{p}}-\frac{1}{\mu} \boldsymbol{\nabla} \langle \, p \rangle^\beta\quad \text{in } \mathscr V_\beta, \end{gather}$$

(2.1c)$$\begin{gather}{\rm B.C.} 1\quad \boldsymbol{\mathsf{P}}\boldsymbol{\cdot} (\boldsymbol{n}\boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{v}+\boldsymbol{\nabla}\boldsymbol{v}^{\rm T}))= \boldsymbol{\mathsf{P}}\boldsymbol{\cdot}(\boldsymbol{n} \boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{\tilde{p}})= \textbf{0}\quad \text{at }\mathscr{A}_{\beta \sigma}, \end{gather}$$

(2.1d)$$\begin{gather}{\rm B.C.} 2\quad \boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{v}=0\quad \text{at }\mathscr{A}_{\beta\sigma}, \end{gather}$$

(2.1e)$$\begin{gather}\psi(\boldsymbol{r})= \psi(\boldsymbol{r}+\boldsymbol{l}_i) \quad {\rm for}\ i=1,2,3, \quad \psi = \boldsymbol{v}, \boldsymbol{\mathsf{T}}_{\tilde{p}}, \end{gather}$$

(2.1f)$$\begin{gather}\langle \tilde{p}\rangle^\beta =0. \end{gather}$$

In the above equations, $\boldsymbol {v}$ denotes the fluid velocity and $\tilde {p}$ the pressure deviation resulting from the decomposition $\tilde p= p -\langle \, p\rangle ^\beta$ (Gray 1975), in which $\langle \, p\rangle ^\beta$ is the intrinsic average of the pressure. The intrinsic average of a quantity $\psi$ taking values in the $\beta$-phase is given by $\langle \psi \rangle ^\beta =({1}/{V_\beta })\int _{\mathscr {V}_\beta } \psi \,{\rm d} V$, where $\mathscr {V}_\beta$ (of measure $V_\beta$) denotes the portion of the periodic unit cell $\mathscr {V}$ (of measure $V$ and size $\ell _c$) occupied by the $\beta$-phase. Moreover, the stress-like tensor, $\mu \boldsymbol{\mathsf{T}}_{\tilde {p}}$, is defined as $\mu \boldsymbol{\mathsf{T}}_{\tilde {p}} = -\tilde {p}\boldsymbol{\mathsf{I}} +\mu (\boldsymbol {\nabla } \boldsymbol {v}+\boldsymbol {\nabla } \boldsymbol {v}^{\rm T})$. In the perfect-slip condition expressed in (2.1c), $\boldsymbol{\mathsf{P}}=\boldsymbol{\mathsf{I}}-\boldsymbol {n}\boldsymbol {n}$ is the projection tensor on the tangential plane at the fluid–solid interfaces inside the medium denoted by $\mathscr {A}_{\beta \sigma }$, $\boldsymbol{\mathsf{I}}$ being the identity tensor and $\boldsymbol {n}$ the unit normal vector at $\mathscr {A}_{\beta \sigma }$, pointing out of the $\beta$-phase. Note that (2.1c), which represents zero shear at $\mathscr {A}_{\beta \sigma }$, does not make reference to any partial-slip model at this interface. Moreover, in (2.1e), $\boldsymbol {l}_i$ is the periodic lattice vector of the unit cell in the $i$th direction and $\boldsymbol {r}$ is a vector locating any point inside the fluid phase with respect to a fixed origin of the coordinate system. In the course of the development, the superficial average, $\langle \psi \rangle$, of $\psi$ is also used, which is defined as $\langle \psi \rangle =({1}/{V})\int _{\mathscr {V}_\beta } \psi \,{\rm d} V$, i.e. $\langle \psi \rangle =\varepsilon \langle \psi \rangle ^\beta$, $\varepsilon =V_\beta /V$ being the porosity of the medium.

Figure 1. Sketch of a porous medium of characteristic length $L$ that is represented as an array of periodic unit cells $\mathscr {V}$ containing a fluid phase (the $\beta$-phase with characteristic length $\ell _\beta$) and a solid phase (the $\sigma$-phase with characteristic length $\ell _\sigma$).

As a usual prerequisite to carrying out upscaling, it is assumed that characteristic length scales are well separated, i.e. $\ell _c \ll L$, $L$ being the size of the macroscopic system (see figure 1). As a consequence of this scale hierarchy, average quantities can be regarded as constants within the unit cell, so that $\langle \tilde {p}\rangle ^\beta = 0$, as indicated in (2.1f), which is required in order for the above problem to be well posed.

3. Upscaling

The purpose of upscaling is to provide the average mass and momentum balance equations that filter the non-redundant information contained in the flow problem within the unit cell and affect the result to a single point of the effective medium. For mass conservation, the superficial averaging operator can be applied to (2.1a). By making use of the spatial averaging theorem (Whitaker 1999), which reads

(3.1)\begin{equation} \langle\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\psi}\rangle= \boldsymbol{\nabla}\boldsymbol{\cdot}\langle\boldsymbol{\psi}\rangle+\frac{1}{V} \int_{\mathscr{A}_{\beta\sigma}}\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\psi}\,{\rm d} A, \end{equation}

with $\boldsymbol {\psi }=\boldsymbol {v}$, and taking into account the impervious boundary condition given in (2.1d), this leads to the classical macroscopic mass balance equation,

(3.2)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\langle \boldsymbol{v}\rangle=0. \end{equation}

3.1. Adjoint problem and Green's formula

In order to derive the macroscopic momentum balance equation, it is convenient to introduce the closure variables, vector $\boldsymbol {d}$ and second-order tensor $\boldsymbol{\mathsf{D}}$, which solve the following adjoint problem associated with (2.1) (see Bottaro (2019) for applications of the adjoint method):

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

(3.3b)$$\begin{gather}\boldsymbol{0} = \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\boldsymbol{d}} +\boldsymbol{\mathsf{I}}\quad \text{in } \mathscr{V}_\beta, \end{gather}$$

(3.3c)$$\begin{gather}{\rm B.C.} 1\quad \boldsymbol{\mathsf{P}}\boldsymbol{\cdot}(\boldsymbol{n}\boldsymbol{\cdot} (\boldsymbol{\nabla}\boldsymbol{\mathsf{D}}+\boldsymbol{\nabla}\boldsymbol{\mathsf{D}}^{T1}))= \boldsymbol{\mathsf{P}}\boldsymbol{\cdot}(\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\boldsymbol{d}})= \textbf{0}\quad \text{at }\mathscr{A}_{\beta \sigma}, \end{gather}$$

(3.3d)$$\begin{gather}{\rm B.C.} 2\quad \boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\mathsf{D}}=\boldsymbol{0}\quad \text{at}\ \mathscr{A}_{\beta \sigma}, \end{gather}$$

(3.3e)$$\begin{gather}\psi(\boldsymbol{r})= \psi(\boldsymbol{r}+\boldsymbol{l}_i)\quad {\rm for}\ i=1,2,3, \quad \psi = \boldsymbol{\mathsf{D}}, \boldsymbol{\mathsf{T}}_{\boldsymbol{d}}, \end{gather}$$

(3.3f)$$\begin{gather}\boldsymbol{d} =\textbf{0}\quad \text{at }\boldsymbol{r}_{0}, \end{gather}$$

(3.3g)$$\begin{gather}\boldsymbol{\mathsf{T}}_{\boldsymbol{d}}={-}\boldsymbol{\mathsf{I}} \boldsymbol{d}+\boldsymbol{\nabla}\boldsymbol{\mathsf{D}}+ \boldsymbol{\nabla}\boldsymbol{\mathsf{D}}^{T1}. \end{gather}$$

The constraint given in (3.3f) is necessary in order to make the problem well posed, and it replaces (2.1f). Although it could have also been formulated as $\langle \boldsymbol {d}\rangle ^\beta = \textbf {0}$, (3.3f) may be preferred as it is less computationally demanding. In addition, in (3.3c) and (3.3g), the superscript $T1$ denotes the transpose of a third-order tensor that permutes its first and second indices, i.e. $(\boldsymbol {\nabla } \boldsymbol{\mathsf{D}})_{ijk}^{T1}=(\boldsymbol {\nabla } \boldsymbol{\mathsf{D}})_{jik}$ in Gibbs’ notation. In (3.3g), $\boldsymbol{\mathsf{I}}\boldsymbol {d}$ represents the dyadic product between the identity tensor $\boldsymbol{\mathsf{I}}$ and vector $\boldsymbol {d}$, i.e. $\boldsymbol{\mathsf{I}}\boldsymbol {d}\equiv \boldsymbol{\mathsf{I}}\otimes \boldsymbol {d}$. However, for the sake of simplicity in notation, the dyadic product sign is omitted in the rest of this work since this does not cause any ambiguity. It is worth adding that the closure problem given in (3.3) can be conceived as a simplified version of the one reported by Penta & Merodio (2017) within a different physical context using the homogenisation approach (Hornung 2012).

To derive the upscaled form of the momentum balance equation, it is of interest to consider the following Green's formula that is valid for any arbitrary scalar field $a$, vector fields $\boldsymbol {a}$ and $\boldsymbol {b}$, and second-order tensor field $\boldsymbol{\mathsf{B}}$, both $\boldsymbol {a}$ and $\boldsymbol{\mathsf{B}}$ being solenoidal and all the variables spatially periodic. This formula is given by (see the derivation in the appendix of Lasseux & Valdés-Parada 2022; Sánchez-Vargas, Valdés-Parada & Lasseux 2022)

(3.4)\begin{equation} \langle \boldsymbol{a}\boldsymbol{\cdot} (\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\boldsymbol{b}}) - (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_a) \boldsymbol{\cdot} \boldsymbol{\mathsf{B}}\rangle= \frac{1}{V} \int_{\mathscr{A}_{\beta\sigma}}[\boldsymbol{a} \boldsymbol{\cdot} (\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{\boldsymbol{b}})- \boldsymbol{n} \boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_a \boldsymbol{\cdot} \boldsymbol{\mathsf{B}}]\,{\rm d} A. \end{equation}

Here, $\boldsymbol{\mathsf{T}}_a=-a\boldsymbol{\mathsf{I}} + \boldsymbol {\nabla } \boldsymbol {a}+\boldsymbol {\nabla } \boldsymbol {a}^{{\rm T}}$ and $\boldsymbol{\mathsf{T}}_{\boldsymbol {b}}=-\boldsymbol{\mathsf{I}}\textbf {b} + \boldsymbol {\nabla } \boldsymbol{\mathsf{B}}+\boldsymbol {\nabla } \boldsymbol{\mathsf{B}}^{T1}$. The use of the adjoint method with Green's formula is convenient as it allows the formal derivation of the upscaled model. Moreover, the adjoint problem is formally obtained as the volume integral of the associated Green's functions, as shown by Lasseux, Valdés-Parada & Bottaro (2021). In this reference it is also shown that if more than one source term is present in the flow problem, it is only necessary to solve one closure problem. Owing to these features, the adjoint method with Green's formula is used in this work. Identical results would be obtained while employing the traditional asymptotic homogenisation or volume averaging technique, with, however, the present approach having the advantage of being much more compact in the development.

3.2. Macroscopic momentum balance equation

To relate the physical and adjoint problems, Green's formula given in (3.4) can be used taking $a=\tilde {p}$, $\boldsymbol {a}=\boldsymbol {v}$, $\boldsymbol {b}=\boldsymbol {d}$ and $\boldsymbol{\mathsf{B}}=\boldsymbol{\mathsf{D}}$. Considering (2.1b) and (3.3b), this yields

(3.5)\begin{equation} {-}\langle \boldsymbol{v}\rangle-\frac{1}{\mu}\langle\boldsymbol{\mathsf{D}}\rangle^{\rm T} \boldsymbol{\cdot}\boldsymbol{\nabla}\langle \, p\rangle^\beta=\frac{1}{V}\int_{\mathscr{A}_{\beta\sigma}} [\boldsymbol{v}\boldsymbol{\cdot}(\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\boldsymbol{d}})-\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\tilde{p}}\boldsymbol{\cdot} \boldsymbol{\mathsf{D}}]\,{\rm d} A. \end{equation}

Since $\boldsymbol {v}$ and $\boldsymbol{\mathsf{D}}$ are tangential at $\mathscr {A}_{\beta \sigma }$, $\boldsymbol {v}=\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol{\mathsf{P}}$ and $\boldsymbol{\mathsf{D}}=\boldsymbol{\mathsf{P}}\boldsymbol {\cdot }\boldsymbol{\mathsf{D}}$. When these relationships are substituted back into the area integral on the right-hand side of (3.5), taking into account that $\boldsymbol{\mathsf{P}}$ is a symmetric tensor, one can write this term as $({1}/{V})\int _{\mathscr {A}_{\beta \sigma }} [\boldsymbol {v}\boldsymbol {\cdot } \boldsymbol{\mathsf{P}}\boldsymbol {\cdot }(\boldsymbol {n}\boldsymbol {\cdot } \boldsymbol{\mathsf{T}}_{\boldsymbol {d}})- \boldsymbol{\mathsf{P}}\boldsymbol {\cdot }(\boldsymbol {n} \boldsymbol {\cdot }\boldsymbol{\mathsf{T}}_{\tilde {p}})\boldsymbol {\cdot } \boldsymbol{\mathsf{D}}]\,{\rm d} A$, which is zero due to boundary conditions (2.1c) and (3.3c). The macroscopic momentum equation for perfect-slip flow finally follows from (3.5), and takes the form

(3.6)\begin{equation} \langle \boldsymbol{v}\rangle={-}\frac{1}{\mu}\boldsymbol{\mathsf{K}}_{ps} \boldsymbol{\cdot}\boldsymbol{\nabla}\langle \, p\rangle^\beta, \end{equation}

with the perfect-slip permeability tensor given by $\boldsymbol{\mathsf{K}}_{ps}=\langle \boldsymbol{\mathsf{D}}\rangle$. This definition takes into account the fact that $\boldsymbol{\mathsf{K}}_{ps}$ is a symmetric tensor, the proof of which is provided in Appendix A, together with proof of its positivity. It must be noted that $\boldsymbol{\mathsf{K}}_{ps}$ is intrinsic, i.e. is a characteristic of the porous structure only. It is expected that this intrinsic effective coefficient is larger than the no-slip intrinsic permeability and that it corresponds to the limit as the slip length tends to infinity in partial slip as clearly shown below.

4. Relationship between perfect-, partial- and no-slip flow permeability tensors

4.1. Recap of the macroscopic partial-slip and no-slip flow models

When flow is in the (partial-)slip regime, obeying a first-order slip condition at $\mathscr {A}_{\beta \sigma }$, with slip length $\ell _s$, the pore-scale flow problem in a periodic unit cell is still defined by (2.1), with the difference that the two interfacial boundary conditions given in (2.1c) and (2.1d) are replaced by (see Lauga et al. 2007; Lasseux et al. 2014, 2016)

(4.1)\begin{equation} {\rm B.C.}\quad \boldsymbol{v}={-}\ell_s\boldsymbol{\mathsf{P}}\boldsymbol{\cdot}(\boldsymbol{n} \boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{v}+\boldsymbol{\nabla}\boldsymbol{v}^{\rm T}))={-}\ell_s\boldsymbol{\mathsf{P}}\boldsymbol{\cdot}(\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\tilde{p}})\quad \text{at }\mathscr{A}_{\beta \sigma}.\end{equation}

Here, the slip flow condition is formulated assuming the slip length is a scalar. However, extension to the case where this quantity is tensorial is straightforward.

As in § 3.1, the corresponding adjoint (closure) problem for the velocity is now defined in terms of the closure variables $\boldsymbol {d}_s$ and $\boldsymbol{\mathsf{D}}_s$ instead of $\boldsymbol {d}$ and $\boldsymbol{\mathsf{D}}$ that solve the problem given by (3.3), with the difference that B.C.1 and B.C.2 are replaced by

(4.2)\begin{equation} {\rm B.C.}\quad \boldsymbol{\mathsf{D}}_s={-}\ell_s \boldsymbol{\mathsf{P}}\boldsymbol{\cdot} (\boldsymbol{n}\boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{\mathsf{D}}_s+ \boldsymbol{\nabla}\boldsymbol{\mathsf{D}}^{T1}_s))={-}\ell_s\boldsymbol{\mathsf{P}}\boldsymbol{\cdot}(\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\boldsymbol{d}_s})\quad \text{at }\mathscr{A}_{\beta \sigma}. \end{equation}

Using the same procedure as the one described in § 3.2, i.e. Green's formula provided in (3.4) with $a=\tilde {p}$, $\boldsymbol {a}=\boldsymbol {v}$, $\boldsymbol {b}=\boldsymbol {d}_s$ and $\boldsymbol{\mathsf{B}}=\boldsymbol{\mathsf{D}}_s$, together with (2.1b) and (3.3b) (with $\boldsymbol{\mathsf{D}}_s$ and $\boldsymbol {d}_s$ respectively replacing $\boldsymbol{\mathsf{D}}$ and $\boldsymbol {d}$), leads to

(4.3)\begin{equation} {-}\langle \boldsymbol{v}\rangle-\frac{1}{\mu}\langle\boldsymbol{\mathsf{D}}_s\rangle^{\rm T} \boldsymbol{\cdot}\boldsymbol{\nabla}\langle \, p\rangle^\beta= \frac{1}{V}\int_{\mathscr{A}_{\beta\sigma}} [\boldsymbol{v}\boldsymbol{\cdot}(\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\boldsymbol{d}_s})- \boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{\tilde{p}}\boldsymbol{\cdot} \boldsymbol{\mathsf{D}}_s]\,{\rm d} A. \end{equation}

When the boundary conditions expressed in (4.1) and (4.2) are taken into account, the area integral term on the right-hand side of the above expression can be rewritten as $-({\ell _s}/{V})\int _{\mathscr {A}_{\beta \sigma }} [\boldsymbol{\mathsf{P}}\boldsymbol {\cdot }(\boldsymbol {n}\boldsymbol {\cdot } \boldsymbol{\mathsf{T}}_{\tilde {p}})\boldsymbol {\cdot }(\boldsymbol {n}\boldsymbol {\cdot } \boldsymbol{\mathsf{T}}_{\boldsymbol {d}_s})- \boldsymbol {n}\boldsymbol {\cdot }\boldsymbol{\mathsf{T}}_{\tilde {p}}\boldsymbol {\cdot } \boldsymbol{\mathsf{P}}\boldsymbol {\cdot }(\boldsymbol {n}\boldsymbol {\cdot } \boldsymbol{\mathsf{T}}_{\boldsymbol {d}_s})]\,{\rm d} A$. Since $\boldsymbol{\mathsf{P}}$ is a symmetric tensor, $\boldsymbol{\mathsf{P}}\boldsymbol {\cdot }(\boldsymbol {n} \boldsymbol {\cdot }\boldsymbol{\mathsf{T}}_{\tilde {p}})= \boldsymbol {n}\boldsymbol {\cdot }\boldsymbol{\mathsf{T}}_{\tilde {p}}\boldsymbol {\cdot } \boldsymbol{\mathsf{P}}$, which indicates that this area integral is zero. Hence, (4.3) yields the macroscopic (partial-)slip flow momentum equation that is given by

(4.4)\begin{equation} \langle \boldsymbol{v}\rangle={-}\frac{1}{\mu}\boldsymbol{\mathsf{K}}_s\boldsymbol{\cdot} \boldsymbol{\nabla}\langle \, p\rangle^\beta, \end{equation}

where the apparent slip permeability tensor is defined as $\boldsymbol{\mathsf{K}}_s=\langle \boldsymbol{\mathsf{D}}_s\rangle$, taking into account that this is a symmetric (and positive) tensor (see the proof in Appendix A). The tensor $\boldsymbol{\mathsf{K}}_s$ is apparent in the sense that it depends on $\ell _s$.

Statements of both the flow (2.1) and adjoint (3.3) problems are still valid when no slip occurs, in which circumstances one only needs to consider $\ell _s=0$. In that case, the adjoint problem given in (3.3) is written in terms of closure variables $\boldsymbol {d}_0$ and $\boldsymbol{\mathsf{D}}_0$ instead of $\boldsymbol {d}$ and $\boldsymbol{\mathsf{D}}$, and the interfacial boundary condition is now

(4.5)\begin{equation} {\rm B.C.}\quad \boldsymbol{\mathsf{D}}_0= \boldsymbol{0}\quad \text{at }\mathscr{A}_{\beta \sigma}. \end{equation}

The macroscopic model corresponds to the classical Darcy law given by (4.4) in which $\boldsymbol{\mathsf{K}}_s$ is replaced by the intrinsic no-slip permeability tensor $\boldsymbol{\mathsf{K}}$, defined as $\boldsymbol{\mathsf{K}}=\langle \boldsymbol{\mathsf{D}}_0\rangle$.

4.2. The special case of rectilinear flow

In the case of slip flow in a straight tube of radius $R$ and length L, the component of (4.4) in the axial direction (i.e. the $z$-direction) is given by

(4.6a,b)\begin{equation} \langle v_z \rangle ={-}\frac{K_{szz}}{\mu} \frac{{\rm d}\langle \, p\rangle^\beta }{{\rm d} z},\quad K_{szz} =\frac{R^2}{8\mu L} \left(1+\frac{4\ell_s}{R}\right). \end{equation}

Clearly, in this case, the apparent permeability diverges when $\ell_s \rightarrow \infty$, which may be understood as perfect-slip conditions. This is consistent with the fact that the pore-scale flow solution indicates that $v_z$ is constant, whereas ${\rm d} p/{\rm d} z=0$, thus leading to $\langle v_z\rangle$ being a finite constant. A similar result is obtained in any configuration allowing rectilinear flow. The physical explanation of this particular feature is the following and stems from the fact that, in the case of a straight channel, the velocity exhibits a piston-like solution. Even if the fluid is viscous in nature, the flow is such that there is no viscous dissipation, neither at the solid–fluid interfaces (due to perfect slip) nor in the bulk, because fluid layers are flowing parallel to each other. As a result, the fluid flows as a rigid body with a constant and uniform velocity given by the imposed flow rate (not by a pressure gradient which is evanescent). Consequently, the permeability in that case is undefined. In any other non-parallel flow situation, although there is no shear at $\mathscr {A}_{\beta \sigma }$, shear is induced in the bulk since streamlines are not straight. This yields a finite permeability.

4.3. Relationship between $\boldsymbol{\mathsf{K}}_{ps}$, $\boldsymbol{\mathsf{K}}_s$ and $\boldsymbol{\mathsf{K}}$

In order to relate $\boldsymbol{\mathsf{K}}_{ps}$ and $\boldsymbol{\mathsf{K}}_s$ (and $\boldsymbol{\mathsf{K}}$), it is convenient to consider the following Green's formula that is valid for any two arbitrary vector fields, $\boldsymbol {a}$ and $\boldsymbol {b}$, and second-order solenoidal tensor fields, $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{B}}$. When all the variables are spatially periodic, this formula is expressed as (Lasseux, Zaouter & Valdés-Parada 2023)

(4.7)\begin{equation} \langle \boldsymbol{\mathsf{A}}^{\rm T} \boldsymbol{\cdot} (\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\boldsymbol{b}}) - (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{\boldsymbol{a}})^{\rm T} \boldsymbol{\cdot}\boldsymbol{\mathsf{B}}\rangle= \frac{1}{V}\int_{\mathscr{A}_{\beta\sigma}} [\boldsymbol{\mathsf{A}}^{\rm T} \boldsymbol{\cdot}(\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\boldsymbol{b}})-(\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\boldsymbol{a}})^{\rm T} \boldsymbol{\cdot}\boldsymbol{\mathsf{B}}]\,{\rm d} A, \end{equation}

where, as before, $\boldsymbol{\mathsf{T}}_{\boldsymbol {a}}=-\boldsymbol{\mathsf{I}}\boldsymbol {a} + \boldsymbol {\nabla }\boldsymbol{\mathsf{A}}+\boldsymbol {\nabla } \boldsymbol{\mathsf{A}}^{T1}$ and $\boldsymbol{\mathsf{T}}_{\boldsymbol {b}}=-\boldsymbol{\mathsf{I}}\boldsymbol {b} + \boldsymbol {\nabla } \boldsymbol{\mathsf{B}}+\boldsymbol {\nabla } \boldsymbol{\mathsf{B}}^{T1}$. When the above formula is used with $\boldsymbol {a}=\boldsymbol {d}_s$, $\boldsymbol{\mathsf{A}}=\boldsymbol{\mathsf{D}}_s$, $\boldsymbol {b}=\boldsymbol {d}$ and $\boldsymbol{\mathsf{B}}=\boldsymbol{\mathsf{D}}$, and once the momentum-like equations for both ($\boldsymbol {d}$,$\boldsymbol{\mathsf{D}}$) and ($\boldsymbol {d}_s$,$\boldsymbol{\mathsf{D}}_s$) are taken into account, one obtains

(4.8)\begin{equation} {-}\boldsymbol{\mathsf{K}}_s+\boldsymbol{\mathsf{K}}_{ps}= \frac{1}{V} \int_{\mathscr{A}_{\beta\sigma}} [\boldsymbol{\mathsf{D}}_s^{\rm T} \boldsymbol{\cdot} (\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{\boldsymbol{d}})- (\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{\boldsymbol{d}_s})^{\rm T} \boldsymbol{\cdot}\boldsymbol{\mathsf{D}}]\,{\rm d} A. \end{equation}

Since $\boldsymbol{\mathsf{D}}_s=\boldsymbol{\mathsf{P}}\boldsymbol {\cdot }\boldsymbol{\mathsf{D}}_s$ at $\mathscr {A}_{\beta \sigma }$ and $\boldsymbol{\mathsf{P}}$ is a symmetric tensor, the first term in the area integral of the above equation can be equivalently rewritten as $\boldsymbol{\mathsf{D}}_s^{\rm T} \boldsymbol {\cdot } (\boldsymbol {n}\boldsymbol {\cdot } \boldsymbol{\mathsf{T}}_{\boldsymbol {d}})=\boldsymbol{\mathsf{D}}_s^{\rm T} \boldsymbol {\cdot } \boldsymbol{\mathsf{P}}\boldsymbol {\cdot }(\boldsymbol {n} \boldsymbol {\cdot }\boldsymbol{\mathsf{T}}_{\boldsymbol {d}})$, which is obviously zero due to the boundary condition given in (3.3c). Using the fact that $\boldsymbol{\mathsf{D}}=\boldsymbol{\mathsf{P}}\boldsymbol {\cdot }\boldsymbol{\mathsf{D}}$ at $\mathscr {A}_{\beta \sigma }$ and the boundary condition (4.2) in the remaining area integral term of (4.8), together with the notation $\boldsymbol{\mathsf{W}}_s=\boldsymbol {\nabla }\boldsymbol{\mathsf{D}}_s+ \boldsymbol {\nabla }\boldsymbol{\mathsf{D}}^{T1}_s$, the relationship between $\boldsymbol{\mathsf{K}}_{ps}$ and $\boldsymbol{\mathsf{K}}_{s}$ follows from (4.8) as

(4.9)\begin{align} \boldsymbol{\mathsf{K}}_{ps}-\boldsymbol{\mathsf{K}}_s &={-}\frac{1}{V}\int_{\mathscr{A}_{\beta\sigma}} (\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\boldsymbol{d}_s})^{\rm T} \boldsymbol{\cdot} \boldsymbol{\mathsf{D}}\,{\rm d} A={-}\frac{1}{V}\int_{\mathscr{A}_{\beta\sigma}} (\boldsymbol{\mathsf{P}}\boldsymbol{\cdot}(\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\boldsymbol{d}_s}))^{\rm T} \boldsymbol{\cdot} \boldsymbol{\mathsf{D}}\,{\rm d} A \nonumber\\ &={-}\frac{1}{V}\int_{\mathscr{A}_{\beta\sigma}} (\boldsymbol{\mathsf{P}} \boldsymbol{\cdot}(\boldsymbol{n} \boldsymbol{\cdot}\boldsymbol{\mathsf{W}}_s))^{\rm T}\boldsymbol{\cdot} \boldsymbol{\mathsf{D}}\,{\rm d} A={-}\frac{1}{V}\int_{\mathscr{A}_{\beta\sigma}} (\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\mathsf{W}}_s)^{\rm T} \boldsymbol{\cdot} \boldsymbol{\mathsf{D}}\,{\rm d} A. \end{align}

In the limit of perfect slip, i.e. $\ell _s\rightarrow \infty$, $\boldsymbol{\mathsf{D}}_s\rightarrow \boldsymbol{\mathsf{D}}$ and $\boldsymbol{\mathsf{K}}_s\rightarrow \boldsymbol{\mathsf{K}}_{ps}$, the area integral term in (4.9) is zero due to the boundary condition given in (3.3c) and the above result is an identity, as expected. More importantly, in the limit of no slip, $\ell _s\rightarrow 0$, $\boldsymbol{\mathsf{D}}_s\equiv \boldsymbol{\mathsf{D}}_0$ and $\boldsymbol{\mathsf{K}}_s\rightarrow \boldsymbol{\mathsf{K}}$ so that the relationship between the two intrinsic permeability tensors, $\boldsymbol{\mathsf{K}}_{ps}$ and $\boldsymbol{\mathsf{K}}$, becomes

(4.10)\begin{gather} \boldsymbol{\mathsf{K}}_{ps}-\boldsymbol{\mathsf{K}} ={-}\frac{1}{V}\int_{\mathscr{A}_{\beta\sigma} } (\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{\boldsymbol{d}_0})^{\rm T} \boldsymbol{\cdot} \boldsymbol{\mathsf{D}}\,{\rm d} A ={-}\frac{1}{V}\int_{\mathscr{A}_{\beta\sigma}} [\boldsymbol{n}\boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{\mathsf{D}}_0+ \boldsymbol{\nabla}\boldsymbol{\mathsf{D}}_0^{T1})]^{\rm T} \boldsymbol{\cdot}\boldsymbol{\mathsf{D}}\,{\rm d} A. \end{gather}

A physical interpretation of this result can be provided by noticing that the solutions for the pressure deviation and the velocity in the periodic unit cell are $\tilde {p}=-\boldsymbol {d}\boldsymbol {\cdot }\boldsymbol {\nabla }\langle \, p\rangle ^\beta$ and $\boldsymbol {v}=-({1}/{\mu })\boldsymbol{\mathsf{D}}\boldsymbol {\cdot }\boldsymbol {\nabla }\langle \, p\rangle ^\beta$ in the perfect-slip case (respectively, $\tilde {p}_0=-\boldsymbol {d}_0\boldsymbol {\cdot }\boldsymbol {\nabla }\langle \, p\rangle ^\beta$ and $\boldsymbol {v}_0=-({1}/{\mu })\boldsymbol{\mathsf{D}}_0\boldsymbol {\cdot }\boldsymbol {\nabla }\langle \, p\rangle ^\beta$ in the no-slip case); see chapter 4 in Whitaker 1999. Pre- and post-multiplying the first of (4.10) by $-\boldsymbol {\nabla }\langle \, p\rangle ^\beta /\mu$ leads to

(4.11)\begin{equation} (\langle \boldsymbol{v}\rangle -\langle \boldsymbol{v}_0 \rangle) \boldsymbol{\cdot} \boldsymbol{\nabla} \langle \, p\rangle^\beta= a_v \langle (\boldsymbol{n}\boldsymbol{\cdot} \mu \boldsymbol{\mathsf{T}}_{p_0}) \boldsymbol{\cdot} \boldsymbol{v}\rangle_{\beta\sigma},\end{equation}

with $\mu \boldsymbol{\mathsf{T}}_{p_0} = -p_0\boldsymbol{\mathsf{I}} +\mu (\boldsymbol {\nabla } \boldsymbol {v}_0+\boldsymbol {\nabla } \boldsymbol {v}_0^{\rm T})$. Here, $a_v=A_{\beta \sigma }/V$ is the interfacial area per unit volume of the medium and $\langle \psi \rangle _{\beta \sigma }= ({1}/{A_{\beta \sigma }}) \int _{\mathscr {A}_{\beta \sigma }}\psi \,{\rm d} A$ denotes the interfacial average of $\psi$. The above equation indicates that the excess macroscopic velocity in the direction of the macroscopic pressure gradient induced by perfect slip (with respect to no slip) is exactly the interfacial average of the projection of the (perfect-)slip velocity (at $\mathscr {A}_{\beta \sigma }$) onto the no-slip stress exerted by the solid skeleton on the fluid. The same conclusion generalises if partial slip is meant instead of perfect slip, the latter representing the limiting case $\ell _s\rightarrow \infty$. A similar interpretation also holds between perfect slip and partial slip. Finally, it should be noted that all the above relationships remain valid in the case of rectilinear flow.

5. Illustrative results on two-dimensional porous structures

To illustrate the developments reported above, the apparent permeabilities for no slip, partial slip and perfect slip are predicted by solving the adjoint problems firstly on a periodic unit cell for a simple two-dimensional structure, consisting of a square pattern of parallel cylinders of circular cross-section like the one depicted in figure 2(a). In addition, direct numerical simulations (DNS) are performed in the perfect-slip regime over a finite-size domain composed of three juxtaposed unit cells of this pattern in the $x$-direction, considering periodicity in the $y$-direction and specifying the normal stress along $\boldsymbol {e}_x$ in such a way that each unit cell is subject to a unitary macroscopic pressure gradient. Secondly, two other classical structures of parallel cylinders are considered, namely, the staggered and hexagonal lattices (see figure 4).

Figure 2. Sketch of a (a) periodic and (b) non-periodic (Chang's) unit cell corresponding to a representation of the porous medium geometry as a square pattern of parallel cylinders of circular cross-section.

5.1. Square pattern of parallel cylinders

For convenience, the flow and closure problems are made dimensionless using $\ell _c$ (i.e. the unit cell side length), $\ell _c\|\boldsymbol {\nabla }\langle \, p\rangle ^\beta\|$ and $\ell _c^2\|\boldsymbol {\nabla }\langle \, p\rangle ^\beta\|/\mu$ for the reference length, pressure and velocity, respectively. Because of the symmetry of the structure, all the permeability tensors under concern are spherical and interest is focused only on their dimensionless $xx$-component, namely $K_{psxx}^*=K_{psxx}/\ell _c^2$, $K_{sxx}^*=K_{sxx}/\ell _c^2$ and $K_{xx}^*=K_{xx}/\ell _c^2$ for the perfect-slip, partial-slip and no-slip permeabilities, respectively. Solutions of the flow and closure problems are computed using the finite element software Comsol Multiphysics 6.1 or a boundary element method (BEM) with constant elements (Pozrikidis 1992), performing a prior mesh convergence check for both methods.

For large enough values of $\varepsilon$, $K_{psxx}^*$ can be approximated analytically by replacing the unit cell of figure 3(a) by Chang's unit cell such as the one sketched in figure 2(b) (see also figure 1.12 in Whitaker 1999), consisting of a fluid domain whose outer boundary is an artificial circle of radius $R_c$ centred on the solid cylinder such that the porosity considered for the periodic unit cell of interest is conserved. At $r=R_c$, the flow is assumed not to be significantly perturbed by the solid obstacle (requiring $\varepsilon$ to be large enough), together with a zero vorticity (Kuwabara 1959). In this way, formulating the flow with the stream function in cylindrical coordinates and solving the resulting bi-harmonic equation leads to

(5.1)\begin{equation} K_{psxx}^*=\frac{1}{16{\rm \pi}}[{-}2 \ln(1-\varepsilon) -1+(1-\varepsilon)^2].\end{equation}

The above expression is expected to reproduce well the numerical solution corresponding to transverse flow across inline arrays of cylinders for sufficiently large porosity values. Using the same approach, the corresponding analytical expression for the permeability under no-slip conditions is (see, for example, (8-4.23) in Happel & Brenner 1981)

(5.2)\begin{equation} K_{xx}^*=\frac{1}{16{\rm \pi}} [{-}2\ln (1-\varepsilon)-3+4(1-\varepsilon)-(1-\varepsilon)^2].\end{equation}

These two limit values of the permeability can be summarised in the equation (Chai et al. 2011)

(5.3)\begin{equation} K_{sxx}^*=\frac{K_{psxx}^*+\dfrac{1}{2}\dfrac{\ell_\sigma/2}{\ell_s} K_{xx}^*}{1+ \dfrac{1}{2}\dfrac{\ell_\sigma/2}{\ell_s}}, \end{equation}

with $\ell_\sigma$ being the cylinder diameter (see figure 2a). This last result was empirically generalised by Geoffre et al. (2021) for transverse flow through a random arrangement of parallel cylinders of circular cross-section with random diameters to

(5.4)\begin{equation} K_{sxx}^*=K_{xx}^*+\frac{K_{psxx}^*-K_{xx}^*}{1+\dfrac{1}{2}\dfrac{\bar{r}}{\ell_s}}, \end{equation}

with $\bar {r}$ being the cylinders' mean radius. The difference between this equation and (5.3) is that in the latter, $K_{psxx}^*$ and $K_{xx}^*$ respectively come from (5.1) and (5.2), whereas in the former, these values are the intrinsic transverse permeabilities for perfect- and no-slip conditions in the geometry under consideration and its use is not restricted to large porosity values.

Figure 3. (a) System under consideration for DNS with periodicity in the $y$-direction and a pressure gradient applied along $\boldsymbol {e}_x$. The colour bar corresponds to the dimensionless pore-scale velocity magnitude in perfect-slip flow condition. (b) Dimensionless perfect-slip, $K_{psxx}^*$, and no-slip, $K_{xx}^*$, permeabilities versus porosity. Here, $K_{psxx}^*$ is computed either from $\langle D_{xx}^*\rangle$, from (4.10) or from (5.1). Numerical results are obtained from DNS, Comsol Multiphysics 6.1 or a BEM. (c) Partial-slip permeability normalised by the no-slip permeability versus $\ell _s^*$. The dashed line represents the normalised perfect-slip permeability. The dotted line represents the empirical correlation (5.4). Here $\varepsilon =0.8$ and $K_{xx}^*=0.01941$.

Results on the perfect-slip, $K_{psxx}^*$, and no-slip, $K_{xx}^*$, permeabilities for porosity values, $\varepsilon$, ranging from $0.3$ to $0.95$, are reported in figure 3(b). In the perfect-slip case, results from DNS (see a representation of the velocity magnitude field in figure 3a) are obtained after averaging the $x$-component of the dimensionless velocity over the central unit cell. As expected, the structure of the medium under consideration is such that the edges $x^*=1$ and $x^*=2$ are isobars so that DNS could be restricted to a single unit cell on which the problem exactly stated as in (2.1) is solved.

As shown in figure 3(b), results on $K_{psxx}^*$ obtained from DNS and closure problems (3.3) are in excellent agreement, with the maximum relative error, considering DNS as the reference, being less than 0.12 %. Moreover, results obtained from Comsol Multiphysics and BEM are also in perfect agreement, with a relative error (taking BEM as the reference) smaller than 0.09 % for perfect slip and 0.07 % for no slip. The results are also perfectly consistent with the prediction from (5.1), which remains less than 1.3 % in error compared with the BEM results, taking the latter as the reference, for $\varepsilon >0.7$. Finally, $K_{psxx}^*$, computed either as $\langle D_{xx}^*\rangle =\langle D_{xx}/\ell _c^2\rangle$ or from the dimensionless version of (4.10), perfectly match, with the relative error, considering the former as the reference, remaining less than 0.05 % over the whole range of $\varepsilon$. All these results validate the macroscopic model for perfect slip and the relationship (4.10) between the perfect-slip and no-slip intrinsic permeabilities.

To complete the results, the partial-slip permeability, $K_{sxx}^*=\langle D_{sxx}^*\rangle$, computed from the solution of the closure problem for partial slip for $\varepsilon =0.8$ and normalised by $K_{xx}^*= \langle D_{0xx}^*\rangle$, are represented in figure 3(c) versus $\ell _s^*$ ranging from $10^{-4}$ to $10^4$. The normalised value of $K_{psxx}^*$ for this porosity is also reported in this figure, showing that the asymptotic value of $K_{sxx}^*$ in the limit $\ell _s^*\rightarrow \infty$ is in excellent agreement with this intrinsic value. Indeed, the relative error between $K_{sxx}^*$ for $\ell _s^*=10^4$ and $K_{psxx}^*$ is 0.04 %.

Finally, the values obtained from the relationship given in (5.4) are reported in figure 3(c) and show agreement with the model reported in the present work.

5.2. Staggered and hexagonal arrays of cylinders

To conclude the numerical analysis, the closure problem solution is extended to periodic unit cells for staggered and hexagonal arrays of cylinders with circular cross-section like those depicted in figure 4, keeping the same dimensionless forms as those employed in § 5.1. The corresponding predictions of the transition of the apparent permeability from no slip to perfect slip are reported in figure 5, both normalised by $\ell _c^2$ (a–d) and by the dimensionless intrinsic permeability under no-slip conditions (i–iv), for several porosity values. In this figure, the results from solving the closure problem in an inline array of cylinders and those resulting from the analytical solution in (5.3) (for $\varepsilon >0.6$) are also included for comparison purposes. Regarding these results, the following comments are in order:

1. (a) When normalised by $\ell _c^2$, the apparent permeability for a staggered array of obstacles is always smaller than the one corresponding to hexagonal or inline arrays of cylinders. This is to be expected since, in the staggered configuration, the streamlines are more deformed than in the two other arrays. Nevertheless, when the apparent permeability is normalised by $K_{xx}^*$, the predictions for the staggered array are quite close to those obtained with an inline array of cylinders for all the porosity values considered here.

2. (b) The value of $K_{sxx}^*$ for the hexagonal array is larger than $K_{sxx}^*$ for the inline array only for $\varepsilon \le 0.5$. Note that, when the apparent permeability is normalised by $K_{xx}^*$, the predicted permeability for the hexagonal array is always smaller than those of both inline and staggered lattices. Moreover, for sufficiently large porosity values (i.e. $\varepsilon \ge 0.9$), all the predictions tend to collapse onto a single curve.

3. (c) As expected, $K_{psxx}^*$ increases with $\varepsilon$ for all the structures. Nevertheless, $K_{xx}^*$ increases more rapidly, and from the ratio between (5.1) and (5.2), it can be concluded that ${K_{psxx}^*}/{K_{xx}^*}\rightarrow 1$ when $\varepsilon \rightarrow 1$.

4. (d) All the results reported in figure 5 were verified to be well reproduced by the empirical equation (5.4). This can be confirmed upon substitution into this equation of the values of $K_{xx}^*$ and $K_{psxx}^*$ reported in table 1 for the three geometrical configurations considered here. Furthermore, the analytical solution given in (5.3) reasonably reproduces the results for $K_{sxx}^*$ for an inline array of cylinders for $\varepsilon =0.7$ and $0.9$.

Figure 4. Sketch of configurations and unit cells consisting of (a) staggered and (b) hexagonal arrays of cylinders of circular cross-section and uniform radius.

Figure 5. Comparison of the predictions of the $xx$ component of the partial-slip permeability tensor with the dimensionless slip length $\ell _s^*$ resulting from solving the adjoint closure problem in inline, staggered and hexagonal arrays of solid obstacles. In (a–d) the results are presented normalised by $\ell _c^2$, i.e. $K_{sxx}^*$, whereas in (i)–(iv) they are normalised by the dimensionless intrinsic permeability under no-slip conditions, $K_{xx}^*$. Porosity values are 0.3 ((a), (i)), 0.5 ((b), (ii)), 0.7 ((c), (iii)) and 0.9 ((d), (iv)). In (c), (d), (iii) and (iv), results from the analytical solution given in (5.3) are included.

Table 1. Predictions of the $xx$ components of the no-slip and perfect-slip permeability tensors normalised by $\ell _c^2$ resulting from solving the corresponding closure problems in periodic unit cells for inline, staggered and hexagonal arrays of cylinders with circular cross-section for several porosity values.

Certainly, the unit cell concept allows consideration of more elaborate geometries in both two and three dimensions. Performing an exhaustive analysis of these geometries, although interesting, surpasses the scope of the present work and shall be considered in future studies.

6. Conclusions

The closed macroscopic model for steady, Newtonian, incompressible, creeping flow in the perfect-slip regime in homogeneous porous media is formally derived in this work, making use of the adjoint method and Green's formula applied to the pore-scale flow problem. It is shown that the averaged momentum equation corresponds to Darcy's law with a permeability tensor, $\boldsymbol{\mathsf{K}}_{ps}$, determined from the solution of a new adjoint (closure) problem on a (periodic) unit cell, representative of the porous material under consideration. This permeability tensor is characteristic of the medium as it is intrinsic to the pore-scale structure. Moreover, $\boldsymbol{\mathsf{K}}_{ps}$ is shown to be symmetric and (semi-definite) positive. Relationships between $\boldsymbol{\mathsf{K}}_{ps}$ and the permeability tensor in the partial-slip flow regime, as well as the classical no-slip intrinsic permeability tensor, $\boldsymbol{\mathsf{K}}$, are developed. Validity of the macroscopic model is assessed from comparison with DNS, and is further supported by an approximated analytical prediction of the permeability considering an inline array of cylinders, for which the relationships between $\boldsymbol{\mathsf{K}}_{ps}$ and $\boldsymbol{\mathsf{K}}$ are also verified. Moreover, the role of the porous medium structure is investigated by comparing the transition permeability curves from no slip to perfect slip in three geometrical configurations, noticing relevant differences for small enough porosity values (i.e. for $\varepsilon \le 0.5$). Nevertheless, when the slip permeability is normalised by the no-slip permeability, it is found that the geometry becomes unimportant in the limit of sufficiently large porosities (i.e. $\varepsilon \ge 0.9$). Moreover, the perfect-slip permeability reaches the same value as the no-slip permeability in the limit of exceedingly large porosities.

Finally, the results reported here confirm that the no-slip and perfect-slip permeabilities respectively represent the lower and upper bounds for Newtonian incompressible creeping flow. They are intrinsic to the porous material and respectively correspond to the maximum and minimum resistances to flow for a given system. This work opens the way for further comparison with experimental data that could be obtained from aqueous one-phase flow in porous materials exhibiting superhydrophobic properties, for instance.