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 (Reference Bottaro2019) 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 (Reference Penta and Merodio2017) within a different physical context using the homogenisation approach (Hornung Reference Hornung2012).

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 Reference Lasseux and Valdés-Parada2022; Sánchez-Vargas, Valdés-Parada & Lasseux Reference Sánchez-Vargas, Valdés-Parada and Lasseux2022)

(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 (Reference Lasseux, Valdés-Parada and Bottaro2021). 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.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. Reference Lauga, Brenner and Stone2007; Lasseux et al. Reference Lasseux, Valdés-Parada, Ochoa-Tapia and Goyeau2014, Reference Lasseux, Valdés-Parada and Porter2016)

(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 Reference Lasseux, Zaouter and Valdés-Parada2023)

(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 Reference Whitaker1999. 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.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 Reference Pozrikidis1992), 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 Reference Whitaker1999), 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 Reference Kuwabara1959). 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 Reference Happel and Brenner1981)

(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. Reference Chai, Lu, Shi and Guo2011)

(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. (Reference Geoffre, Ghestin, Moulin, Bruchon and Drapier2021) 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.