3.1 Closure problem

A convenient closure is now derived by noticing that the boundary value problem for the deviations (2.14) at the pore scale is made non-homogeneous due to the macroscopic source term $\langle \boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}$  on the right-hand side of the boundary condition (2.14c ). Since this problem is incompressible in nature, $\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}$ can be regarded as a parameter and does not represent a source for $\widetilde{p}_{\unicode[STIX]{x1D6FD}}$ and $\tilde{\boldsymbol{v}}_{\unicode[STIX]{x1D6FD}}$ despite its dependence on $\langle \unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}$ . Due to the linearity of the problem, the solution can be sought as a linear combination of the source, namely

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\boldsymbol{v}}_{\unicode[STIX]{x1D6FD}}=\unicode[STIX]{x1D63E}\boldsymbol{\cdot }\langle \boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}} & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{p}_{\unicode[STIX]{x1D6FD}}=\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FD}}\boldsymbol{c}\boldsymbol{\cdot }\langle \boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}. & \displaystyle\end{eqnarray}$$

Any additive constant in the representations of $\tilde{\boldsymbol{v}}_{\unicode[STIX]{x1D6FD}}$ and $\widetilde{p}_{\unicode[STIX]{x1D6FD}}$ can be shown to be unimportant in the final macroscopic model, just as for the upscaling of the classical incompressible Stokes flow with no slip (Whitaker Reference Whitaker1999). Using these representations in the closure equations (2.14), while treating $\langle \boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}$ as a constant due to the contrast of scale between $\ell _{\unicode[STIX]{x1D6FD}}$ and $L$ , yields the following closure problem in terms of the closure variables $\unicode[STIX]{x1D63E}$ and $\boldsymbol{c}$

(3.3a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D63E}=0\quad \text{in }\mathscr{V}_{\unicode[STIX]{x1D6FD}} & \displaystyle\end{eqnarray}$$

(3.3b ) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\unicode[STIX]{x1D735}\boldsymbol{c}+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D63E}-\frac{1}{V_{\unicode[STIX]{x1D6FD}}}\int _{\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}}\boldsymbol{n}\boldsymbol{\cdot }(-\unicode[STIX]{x1D644}\boldsymbol{c}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D63E})\,\text{d}A\quad \text{in }\mathscr{V}_{\unicode[STIX]{x1D6FD}} & \displaystyle\end{eqnarray}$$

(3.3c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63E}+\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }(\boldsymbol{n}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D63E}+(\unicode[STIX]{x1D735}\unicode[STIX]{x1D63E})^{T1}))=-\unicode[STIX]{x1D644}\quad \text{at }\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}} & \displaystyle\end{eqnarray}$$

(3.3d ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \boldsymbol{c}\rangle ^{\unicode[STIX]{x1D6FD}}=0 & \displaystyle\end{eqnarray}$$

(3.3e ) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\begin{array}{@{}ll@{}}\displaystyle \text{Periodicity} & \displaystyle \unicode[STIX]{x1D63E}(\boldsymbol{r}+\boldsymbol{l}_{i})=\unicode[STIX]{x1D63E}(\boldsymbol{r}),\\ & \displaystyle \boldsymbol{c}(\boldsymbol{r}+\boldsymbol{l}_{i})=\boldsymbol{c}(\boldsymbol{r}),\quad i=1,2,3.\end{array}\right\} & \displaystyle\end{eqnarray}$$

$T1$

$T$

At this point, a simple transformation, similar to that employed while upscaling the incompressible creeping flow with no slip leading to Darcy’s law, can be used to obtain a pure differential form of this closure (Barrère, Gipouloux & Whitaker Reference Barrère, Gipouloux and Whitaker1992; Whitaker Reference Whitaker1999). Letting

(3.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D646}_{s}^{-1}=-\unicode[STIX]{x1D700}^{-1}\frac{1}{V_{\unicode[STIX]{x1D6FD}}}\int _{\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}}\boldsymbol{n}\boldsymbol{\cdot }(-\unicode[STIX]{x1D644}\boldsymbol{c}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D63E})\,\text{d}A & & \displaystyle\end{eqnarray}$$

defining $\unicode[STIX]{x1D63F}$ and $\boldsymbol{d}$ as

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}=\unicode[STIX]{x1D700}^{-1}(\unicode[STIX]{x1D63E}+\unicode[STIX]{x1D644})\boldsymbol{\cdot }\unicode[STIX]{x1D646}_{s} & \displaystyle\end{eqnarray}$$

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{d}=\unicode[STIX]{x1D700}^{-1}\boldsymbol{c}\boldsymbol{\cdot }\unicode[STIX]{x1D646}_{s} & \displaystyle\end{eqnarray}$$

(3.7a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}=0\quad \text{in }\mathscr{V}_{\unicode[STIX]{x1D6FD}} & \displaystyle\end{eqnarray}$$

(3.7b ) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\unicode[STIX]{x1D735}\boldsymbol{d}+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D63F}+\unicode[STIX]{x1D644}\quad \text{in }\mathscr{V}_{\unicode[STIX]{x1D6FD}} & \displaystyle\end{eqnarray}$$

(3.7c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}=-\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }(\boldsymbol{n}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}+(\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F})^{T1}))\quad \text{at }\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}} & \displaystyle\end{eqnarray}$$

(3.7d ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \boldsymbol{d}\rangle ^{\unicode[STIX]{x1D6FD}}=0 & \displaystyle\end{eqnarray}$$

(3.7e ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \unicode[STIX]{x1D63F}\rangle =\unicode[STIX]{x1D646}_{s} & \displaystyle\end{eqnarray}$$

(3.7f ) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\begin{array}{@{}ll@{}}\displaystyle \text{Periodicity} & \displaystyle \unicode[STIX]{x1D63F}(\boldsymbol{r}+\boldsymbol{l}_{i})=\unicode[STIX]{x1D63F}(\boldsymbol{r}),\\ & \displaystyle \boldsymbol{d}(\boldsymbol{r}+\boldsymbol{l}_{i})=\boldsymbol{d}(\boldsymbol{r}),\quad i=1,2,3.\end{array}\right\} & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D646}_{s}$

$\unicode[STIX]{x1D63F}$

$\unicode[STIX]{x1D63F}$

$\langle \unicode[STIX]{x1D63E}\rangle ^{\unicode[STIX]{x1D6FD}}=0$

$\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}$

3.2 Closed macroscopic model

When the representation of the deviations given in (3.1) and (3.2) are reported in the unclosed from of the averaged momentum equation (2.11), one readily obtains a Darcy-like macroscopic law

(3.8) $$\begin{eqnarray}\displaystyle 0=-\unicode[STIX]{x1D735}\langle p_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}-\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D646}_{s}^{-1}\boldsymbol{\cdot }\langle \boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\rangle +\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D700}^{-1}\unicode[STIX]{x1D6FB}^{2}\langle \boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\rangle , & & \displaystyle\end{eqnarray}$$

which, along with (2.10) and (2.12), forms the closed averaged model for the slightly compressible slip flow considered in this work. In (3.8), $\unicode[STIX]{x1D646}_{s}$ , defined in (3.7e ), is identified as the apparent slip-corrected permeability.

The last term on the right-hand side of (3.8) is referred to as the Brinkman correction. An order of magnitude analysis can be used to show that this macroscale viscous diffusion term has a negligible contribution in the bulk of the porous medium. Indeed, as can be inferred from the boundary condition in (2.14c )

(3.9) $$\begin{eqnarray}\displaystyle \tilde{\boldsymbol{v}}_{\unicode[STIX]{x1D6FD}}=\boldsymbol{O}\left(\frac{\langle v_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}}{1+\unicode[STIX]{x1D709}\overline{Kn}}\right), & & \displaystyle\end{eqnarray}$$

where $\langle v_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}$ represents the leading order of $\langle \boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}$ . This shows that

(3.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63E}=\boldsymbol{O}((1+\unicode[STIX]{x1D709}\overline{Kn})^{-1}) & & \displaystyle\end{eqnarray}$$

and consequently, from the definition of $\unicode[STIX]{x1D646}_{s}$ in (3.4), the order of magnitude of the Darcy term in (3.8) can be estimated as

(3.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D646}_{s}^{-1}\boldsymbol{\cdot }\langle \boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\rangle =\boldsymbol{O}(\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D700}^{-2}\ell _{\unicode[STIX]{x1D6FD}}^{-1}(1+\unicode[STIX]{x1D709}\overline{Kn})^{-1}a_{v}\langle v_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}). & & \displaystyle\end{eqnarray}$$

In this estimate, $a_{v}=A_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}/V$ represents the interfacial area per unit volume of the medium that is expected to scale as $\ell _{\unicode[STIX]{x1D6FD}}^{-1}$ . Similarly, the order of magnitude of the Brinkman term can be estimated to be

(3.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D700}^{-1}\unicode[STIX]{x1D6FB}^{2}\langle \boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\rangle =\boldsymbol{O}(\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FD}}L^{-2}\langle v_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}), & & \displaystyle\end{eqnarray}$$

which, when compared to (3.11), clearly indicates that the Brinkman term is completely negligible in the context of slip flow, i.e. when $\unicode[STIX]{x1D709}\overline{Kn}\lesssim 0.1$ . The macroscopic model can hence be written, to within an approximation $\boldsymbol{O}(\ell _{\unicode[STIX]{x1D6FD}}/L)$ , as

(3.13a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\langle \unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\langle \unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}\langle \boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}})=0 & \displaystyle\end{eqnarray}$$

(3.13b ) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\unicode[STIX]{x1D735}\langle p_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}-\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D646}_{s}^{-1}\boldsymbol{\cdot }\langle \boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\rangle & \displaystyle\end{eqnarray}$$

(3.13c ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}=F(\langle p_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}). & \displaystyle\end{eqnarray}$$

It must be noted that the apparent slip-corrected permeability, $\unicode[STIX]{x1D646}_{s}$ , in the Darcy-like form of the macroscopic momentum equation (3.13b ) is non-intrinsic as it lumps together topological properties, intrinsic to the microstructure of the medium, and slip effects that depend on the slip length $\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}$ . Moreover, a detailed analysis, provided in appendix A, indicates that, unlike the intrinsic permeability that is the fundamental coefficient in the classical Darcy’s law when no slip occurs, $\unicode[STIX]{x1D646}_{s}$ is not symmetric in the general case. Quasi-symmetry is observed under the constraint

(3.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D709}\overline{Kn}\ll \boldsymbol{O}(\unicode[STIX]{x1D700}) & & \displaystyle\end{eqnarray}$$

that is also detailed and illustrated in appendix A, completing a previous analysis of this property (Skjetne & Auriault Reference Skjetne and Auriault1999). An important consequence is that the full characterization of the tensor $\unicode[STIX]{x1D646}_{s}$ from a measurement point of view is significantly simplified when the constraint in (3.14) is satisfied.

In § 3.3 below, the validity of the macroscopic model and of the associated closure problem (3.7) yielding $\unicode[STIX]{x1D646}_{s}$ is verified by comparison with analytic solutions in a simple configuration and with direct numerical simulations on a more complex medium.

3.3 Validation

The finite element solver Comsol Multiphysics 4.4 was used to carry out the numerical solutions of both the pore-scale equations (2.1) and the closure problem (3.7). Meshing tests were carried out to reach convergence and led to the use of nearly 13 000 domain mesh elements with more than 500 boundary elements. This mesh was built by adding a boundary layer at the solid–fluid interface with a stretching factor of 1.2. As a first validation, a simple structure made of a periodic square pattern of parallel circular cylinders was considered, as represented in figure 2.

Figure 2. Unit cell for the square lattice of cylinders of circular cross-section.

The slip flow orthogonal to the cylinder axes is such that $\unicode[STIX]{x1D646}_{s}=k_{s}\unicode[STIX]{x1D644}$ , where $k_{s}$ can be computed from the solution of the projection of the closure problem either on $\boldsymbol{e}_{x}$ or $\boldsymbol{e}_{y}$ . Due to the unit cell symmetry, the projection leads to a boundary value problem equivalent to the initial incompressible pore-scale problem when symmetry on the velocity and Dirichlet boundary conditions on the pressure are applied on faces respectively parallel and orthogonal to the mean flow direction. Taking the projection of the closure problem (3.7) on $\boldsymbol{e}_{x}$ , the equivalence is obvious when $\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\boldsymbol{e}_{x}$ and $\boldsymbol{d}\boldsymbol{\cdot }\boldsymbol{e}_{x}$ are respectively identified as $\unicode[STIX]{x1D707}\boldsymbol{v}_{\unicode[STIX]{x1D6FD}}/\Vert \unicode[STIX]{x1D735}\langle p_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}\Vert$ and $\widetilde{p}_{\unicode[STIX]{x1D6FD}}/\Vert \unicode[STIX]{x1D735}\langle p_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}\Vert$ in the incompressible version of the pore-scale flow problem (2.1). The validation of the macroscopic model through a comparison of $k_{s}$ obtained either from the closure problem solution or from a DNS of the initial boundary value problem is hence trivial in that case. Nevertheless, further validation can be made by comparing numerical results on $k_{s}$ with analytic predictions obtained on a Chang’s unit cell (Chai et al. Reference Chai, Lu, Shi and Guo2011; Lasseux et al. Reference Lasseux, Valdes Parada, Ochoa Tapia and Goyeau2014) which has been shown to be a reliable estimate for the periodic structure under consideration.

The closure problem was solved on the unit cell of figure 2 for $\unicode[STIX]{x1D700}=0.8$ . Examples of the dimensionless $xx$ component of $\unicode[STIX]{x1D63F}$ and $x$ component of $\boldsymbol{d}$ are reported in figure 3 for two values of the cell Knudsen number $\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }=\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}/\ell$ . Note that $\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }=\unicode[STIX]{x1D709}\overline{Kn}\ell _{\unicode[STIX]{x1D6FD}}^{\ast }=\unicode[STIX]{x1D709}\overline{Kn}\ell _{\unicode[STIX]{x1D6FD}}/\ell$ . The closure fields represented in figures 3(a) and 3(b) are, in fact, the $x$ components of the velocity field made dimensionless by $\ell ^{2}\Vert \unicode[STIX]{x1D735}\langle p_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}\Vert /\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FD}}$ whereas the closure fields in figures 3(c) and 3(d) are the pressure deviation fields made dimensionless by $\ell \Vert \unicode[STIX]{x1D735}\langle p_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}\Vert$ .

Figure 3. Dimensionless closure variable fields computed on the unit cell of figure 2, $\unicode[STIX]{x1D700}=0.8$ : (a) $D_{xx}/\ell ^{2}$ , $\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }=10^{-4}$ ; (b) $D_{xx}/\ell ^{2}$ , $\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }=0.1$ ; (c) $d_{x}/\ell$ , $\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }=10^{-4}$ ; (d) $d_{x}/\ell$ , $\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }=0.1$ .

Numerical results on $k_{s}^{\ast }=k_{s}/\ell ^{2}$ are reported in figure 4 $a$ versus $\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }=\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}/\ell$ , for $\unicode[STIX]{x1D700}=0.8$ . As discussed elsewhere (Lasseux et al. Reference Lasseux, Valdes Parada, Ochoa Tapia and Goyeau2014), a Chang’s unit cell, composed of the solid cylinder immersed in a circular fluid shell, may be used, along with a zero vorticity boundary condition at the outer edge, to derive an approximate analytic solution. For this relatively large value of the porosity, it is given by

(3.15) $$\begin{eqnarray}\displaystyle k_{s}^{\ast } & = & \displaystyle \frac{1}{8\unicode[STIX]{x03C0}\left(1+2\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }\sqrt{\displaystyle \frac{\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D719}}}\right)}\left(-\ln \unicode[STIX]{x1D719}-\frac{3}{2}+2\unicode[STIX]{x1D719}-\frac{\unicode[STIX]{x1D719}^{2}}{2}\right.\nonumber\\ \displaystyle \displaystyle & & \displaystyle +\,\left.2\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }\sqrt{\frac{\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D719}}}\left(-\ln \unicode[STIX]{x1D719}-\frac{1}{2}+\frac{\unicode[STIX]{x1D719}^{2}}{2}\right)\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D719}=1-\unicode[STIX]{x1D700}$ .

Figure 4. Dimensionless apparent slip-corrected permeability $k_{s}^{\ast }=k_{s}/\ell ^{2}$ for the unit cell of figure 2, $\unicode[STIX]{x1D700}=0.8$ , $k^{\ast }=k/\ell ^{2}\simeq 0.01938$ , $k$ being the intrinsic permeability: (a) $k_{s}^{\ast }$ versus $\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }=\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}/\ell$ , comparison between the computed values obtained from the solution of the closure problem (3.7) and the analytic prediction from (3.15); (b) $k_{s}^{\ast }/k^{\ast }$ versus $\unicode[STIX]{x1D709}\overline{Kn}=\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }/\ell _{\unicode[STIX]{x1D6FD}}^{\ast }$ , $\ell _{\unicode[STIX]{x1D6FD}}^{\ast }$ being the dimensionless distance between parallel plates that would exhibit the same intrinsic permeability, i.e.  $\ell _{\unicode[STIX]{x1D6FD}}^{\ast }=2\sqrt{3k^{\ast }/\unicode[STIX]{x1D700}}$ .

The agreement of the prediction from (3.15) with our numerical results is excellent, as shown in figure 4. The relative error on $k_{s}^{\ast }$ , taking the computed value as the reference, is less than $0.4\,\%$ over the interval $10^{-4}\leqslant \unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }\leqslant 0.13$ . Within this range of $\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }$ , $\unicode[STIX]{x1D709}\overline{Kn}$ remains smaller than ${\sim}0.23$ when $\ell _{\unicode[STIX]{x1D6FD}}^{\ast }=\ell _{\unicode[STIX]{x1D6FD}}/\ell$ is estimated from the slit aperture of the equivalent bundle of regularly spaced parallel plates given by $\ell _{\unicode[STIX]{x1D6FD}}^{\ast }=2\sqrt{3k^{\ast }/\unicode[STIX]{x1D700}}$ , $k^{\ast }=k/\ell ^{2}$ being the dimensionless intrinsic permeability, i.e. the permeability without slip. As clearly indicated in figure 4 $b$ , representing $k_{s}^{\ast }/k^{\ast }$ versus $\unicode[STIX]{x1D709}\overline{Kn}$ , $k_{s}^{\ast }$ exhibits a nonlinear dependence on the Knudsen number even though the flow is expected to remain in the slip regime at the microscale. This important remark will be analysed in more detail in § 4.2.

An additional validation is now investigated in the case of a more complex 2-D structure featuring no special symmetry properties and for which no analytic solution is available. In that case, validation of both the upscaling procedure and macroscopic model can be achieved through the comparison of DNS of the pore-scale flow problem (2.1) with the numerical solution of the closure problem (3.7).

The structure at scale $L$ for DNS is obtained by duplicating $n$ times, in the direction of the applied macroscopic pressure gradient, a periodic unit cell of size $\ell$ made of randomly placed cylinders of circular cross-sections with randomly chosen radii within a prescribed interval $[r_{min},r_{max}]$ according to a log-normal distribution. Placement of cylinder centres is constrained to ensure periodicity in the $x$ and $y$ directions and is repeated until a target porosity is achieved, while controlling cylinders overlapping.

Two random unit cells, of porosity $\unicode[STIX]{x1D700}=0.375$ and $\unicode[STIX]{x1D700}=0.804$ , were generated using respectively $r_{min}=0.023\ell$ , $r_{max}=0.062\ell$ and $r_{min}=0.016\ell$ , $r_{max}=0.031\ell$ . They were used to compute the $xx$ dimensionless component, $k_{sxx}^{\ast }=k_{sxx}/\ell ^{2}$ , of the apparent slip-corrected permeability tensor from the closure problem solution. In addition, DNS were performed on structures at scale $L$ obtained from these generic unit cells with $n$ , the number of unit cell replicas, ranging from $3$ to $9$ . Constant pressures were prescribed at the entrance ( $x=0$ ) and exit ( $x=L$ ) of the structure while periodic boundary conditions were applied in the $y$ direction. Superficial and intrinsic averages of $\boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\boldsymbol{\cdot }\boldsymbol{e}_{x}$ and $p_{\unicode[STIX]{x1D6FD}}$ were performed on the central unit cell $(n+1)/2$ to compute $k_{sxx}^{\ast }$ . Computations were carried out for $10^{-4}\leqslant \unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }\leqslant 10^{-2}$ .

Figure 5. Slip-flow DNS on 2-D packs of cylinders. (a) Large-scale structures obtained from periodic replicas along $\boldsymbol{e}_{x}$ of random unit cells. Velocity maps in the $x$ direction are superimposed. (b) Random unit cells of respective porosity $\unicode[STIX]{x1D700}=0.375$ and $\unicode[STIX]{x1D700}=0.804$ used to generate large-scale structures along with $x$ -velocity maps. (c) Pressure profiles at the left edge of the central unit cells.

The DNS procedure is illustrated in figure 5 where we have represented an example of the two $L$ -scale structures, along with the two generic unit cells on which maps of $\boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\boldsymbol{\cdot }\boldsymbol{e}_{x}$ were superimposed. In this figure, we have also reported the pressure profiles computed on the left vertical edge of the central unit cell. The evident non-uniformity of $p_{\unicode[STIX]{x1D6FD}}$ explains why DNS cannot be carried out on a single unit cell with predefined entrance/exit constant Dirichlet boundary conditions on the pressure.

As shown in table 1 and in figure 6 $a$ , results on $k_{sxx}^{\ast }$ obtained from DNS and the closure problem solution are in excellent agreement since the relative error is less than $1.3\,\%$ and decreases when $n$ increases, confirming the validity of the macroscopic model and associated closure.

Table 1. Comparison of the predictions of $k_{sxx}^{\ast }$ from DNS with those from volume averaging. $n$ is the number of unit cell replicas. $k_{sxx}^{\ast }$ values are those obtained from DNS taking the average at the $(n+1)/2$ unit cell. The error per cent is computed relative to the DNS predictions.

Figure 6. (a) Comparison of computed values of $k_{sxx}^{\ast }$ obtained from DNS and volume averaging for the two unit cells of figure 5(b). (b) Dependence of $k_{sxx}^{\ast }/k_{xx}^{\ast }$ on $\unicode[STIX]{x1D709}\overline{Kn}=\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }\ell /\ell _{\unicode[STIX]{x1D6FD}}=\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }/\ell _{\unicode[STIX]{x1D6FD}}^{\ast }$ for $\unicode[STIX]{x1D709}\overline{Kn}\lesssim 0.2$ . Here $\unicode[STIX]{x1D709}\overline{Kn}$ is estimated from the dimensionless distance between plane parallel plates, $\ell _{\unicode[STIX]{x1D6FD}}^{\ast }=2\sqrt{3k^{\ast }/\unicode[STIX]{x1D700}}$ , that would lead to the same intrinsic permeability $k^{\ast }=k_{xx}^{\ast }$ . $k_{xx}^{\ast }\simeq 2.44\times 10^{-6}$ ( $\unicode[STIX]{x1D700}=0.375$ );  $k_{xx}^{\ast }\simeq 1.60\times 10^{-4}$ ( $\unicode[STIX]{x1D700}=0.804$ ).

In figure 6 $b$ , a clear nonlinear dependence of $k_{sxx}^{\ast }$ on $\unicode[STIX]{x1D709}\overline{Kn}$ , most pronounced for the largest value of the porosity, is again highlighted for the two unit cells of figure 5 $b$ , although $\unicode[STIX]{x1D709}\overline{Kn}$ remains smaller than approximately $0.23$ , a range typical of the slip-flow regime as justified by some measurements for even larger Knudsen numbers in straight channels (Harley et al. Reference Harley, Huang, Bau and Zemel1995). This confirms the behaviour reported above for the regular square pattern of cylinders.

At this point it is worth mentioning that, in the seminal work of Klinkenberg (Reference Klinkenberg1941), the linear dependence of the apparent permeability with the Knudsen number was regarded as an approximation resulting from a representation of the porous medium as a system of straight capillaries. As a matter of fact, Klinkenberg (Reference Klinkenberg1941) noticed some nonlinearities in his experimental results with the same tendency as the one observed here. However, the porous media that he considered were relatively tight with intrinsic permeabilities that could go as low as 2.4 mDarcy. As a consequence, the porosity is expected to be quite small and this may explain why the observed nonlinearity remained not too drastic in his experimental results. As shown in figure 6(b), our predictions indicate that the nonlinearity decreases as porosity decreases and becomes almost insignificant for a porosity of 0.25, which is in a qualitative agreement with Klinkenberg’s results.

Before considering in more detail the above mentioned nonlinearity, the impact of an incomplete version of the slip boundary condition, in which the $\unicode[STIX]{x1D735}\boldsymbol{v}_{\unicode[STIX]{x1D6FD}}^{T}$ term is omitted, shall be illustrated. When such an incomplete boundary condition is employed, the analytic form of the apparent slip-corrected permeability of the periodic square pattern of parallel circular cylinders, for a sufficiently large porosity ( $\unicode[STIX]{x1D700}\gtrsim 0.55$ ), is given by

(3.16) $$\begin{eqnarray}\displaystyle k_{s}^{\ast } & = & \displaystyle \frac{1}{8\unicode[STIX]{x03C0}\left(1+2\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }\sqrt{\frac{\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D719}}}\right)}\left(-\ln \unicode[STIX]{x1D719}-\frac{3}{2}+2\unicode[STIX]{x1D719}-\frac{\unicode[STIX]{x1D719}^{2}}{2}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast }\sqrt{\frac{\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D719}}}\left(-\ln \unicode[STIX]{x1D719}+\frac{1}{2}-2\unicode[STIX]{x1D719}+\frac{3\unicode[STIX]{x1D719}^{2}}{2}\right)\right)\end{eqnarray}$$

instead of (3.15) above, with again $\unicode[STIX]{x1D719}=1-\unicode[STIX]{x1D700}$ . This last expression underestimates $k_{s}^{\ast }$ obtained from (3.15). The two slip-flow corrections can be extracted from (3.15) and (3.16) as $k_{s}^{\ast }/k^{\ast }-1$ , where $k^{\ast }=1/8\unicode[STIX]{x03C0}(-\ln \unicode[STIX]{x1D719}-(3/2)+2\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}^{2}/2)$ is again the dimensionless intrinsic permeability and the relative error between the two, taking the former as the reference, is given by $(-2\ln \unicode[STIX]{x1D719}-3+4\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}^{2})/(4(1-2\unicode[STIX]{x1D719}+\unicode[STIX]{x1D719}^{2}))$ . This relative % error, independent of $\unicode[STIX]{x1D709}\overline{Kn}$ , is represented versus the porosity in figure 7 $a$ showing an increasingly significant error with $\unicode[STIX]{x1D700}$ , as it can reach approximately $62\,\%$ for $\unicode[STIX]{x1D700}=0.9$ .

Figure 7. Relative error on the slip correction while using an incomplete slip boundary condition (i.e. omitting the $\unicode[STIX]{x1D735}\boldsymbol{v}_{\unicode[STIX]{x1D6FD}}^{T}$ term). (a) Unit cell of figure 2. The error is computed from the analytic expressions (3.15) and (3.16). (b) Random pack of parallel cylinders (see figure 5 b), $\unicode[STIX]{x1D700}=0.804$ . The error is computed from the closure problem solutions using complete and incomplete boundary conditions at $\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}$ .

When the $\unicode[STIX]{x1D735}\boldsymbol{v}_{\unicode[STIX]{x1D6FD}}^{T}$ term is omitted in the pore-scale slip boundary condition, the closure problem (3.7) is unchanged, except the term $(\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F})^{T1}$ that is no longer present in (3.7c ). This form of the closure was solved on the unit cell of the random pack of parallel cylinders (see figure 5 $b$ ) for $\unicode[STIX]{x1D700}=0.804$ , yielding $k_{sxx}^{\ast }$ that turns out to be larger than that obtained with the complete boundary condition. The slip correction can again be estimated from $k_{sxx}^{\ast }/k^{\ast }-1$ computed with and without the complete shear rate in the boundary condition and the relative error can be formed taking the former as the reference. This % error, represented versus $\unicode[STIX]{x1D709}\overline{Kn}$ in figure 7 $b$ , remains small at exceedingly small values of the Knudsen number (approximately $7\,\%$ for $\unicode[STIX]{x1D709}\overline{Kn}\simeq 0.04$ ), but strongly increases with $\unicode[STIX]{x1D709}\overline{Kn}$ as it reaches approximately $26\,\%$ for $\unicode[STIX]{x1D709}\overline{Kn}\simeq 0.2$ . This clearly evidences that the complete strain rate must be kept in the expression of the slip at the solid boundary. More insight on the impact of the form of the boundary condition and on the nonlinear behaviour mentioned above will be provided in § 4, following a reformulation of the closure problem.

4.1 Expansion in $\unicode[STIX]{x1D709}\overline{Kn}$

The development is carried out on the dimensionless form of the closure given by the starred dimensionless quantities $\boldsymbol{d}^{\ast }=\boldsymbol{d}/\ell _{\unicode[STIX]{x1D6FD}}$ , $\unicode[STIX]{x1D63F}^{\ast }=\unicode[STIX]{x1D63F}/\ell _{\unicode[STIX]{x1D6FD}}^{2}$ , $\unicode[STIX]{x1D6FB}^{\ast }=\ell _{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D735}$ as

(4.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D63F}^{\ast }=0\quad \text{in }\mathscr{V}_{\unicode[STIX]{x1D6FD}} & \displaystyle\end{eqnarray}$$

(4.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\unicode[STIX]{x1D6FB}^{\ast }\boldsymbol{d}^{\ast }+\unicode[STIX]{x1D6FB}^{\ast 2}\unicode[STIX]{x1D63F}^{\ast }+\unicode[STIX]{x1D644}\quad \text{in }\mathscr{V}_{\unicode[STIX]{x1D6FD}} & \displaystyle\end{eqnarray}$$

(4.1c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}^{\ast }=-\unicode[STIX]{x1D709}\overline{Kn}(\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }(\boldsymbol{n}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FB}^{\ast }\unicode[STIX]{x1D63F}^{\ast }+(\unicode[STIX]{x1D6FB}^{\ast }\unicode[STIX]{x1D63F}^{\ast })^{T1}))\quad \text{at }\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}} & \displaystyle\end{eqnarray}$$

(4.1d ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \boldsymbol{d}^{\ast }\rangle ^{\unicode[STIX]{x1D6FD}}=0 & \displaystyle\end{eqnarray}$$

(4.1e ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \unicode[STIX]{x1D63F}^{\ast }\rangle =\unicode[STIX]{x1D646}_{s}\ell _{\unicode[STIX]{x1D6FD}}^{-2} & \displaystyle\end{eqnarray}$$

(4.1f ) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\begin{array}{@{}ll@{}}\displaystyle \text{Periodicity} & \displaystyle \unicode[STIX]{x1D63F}^{\ast }(\boldsymbol{r}^{\ast }+\boldsymbol{l}_{i}^{\ast })=\unicode[STIX]{x1D63F}^{\ast }(\boldsymbol{r}^{\ast }),\\ & \displaystyle \boldsymbol{d}^{\ast }(\boldsymbol{r}^{\ast }+\boldsymbol{l}_{i}^{\ast })=\boldsymbol{d}^{\ast }(\boldsymbol{r}^{\ast }),\quad i=1,2,3.\end{array}\right\} & \displaystyle\end{eqnarray}$$

Since $\unicode[STIX]{x1D709}\overline{Kn}$ remains smaller than unity in the context of slip flow, $\boldsymbol{d}^{\ast }$ and $\unicode[STIX]{x1D63F}^{\ast }$ can be developed up to the $m\text{th}$ order in Maclaurin series expansions that we shall write as

(4.2) $$\begin{eqnarray}\displaystyle \boldsymbol{d}^{\ast }=\boldsymbol{d}_{0}^{\ast }+\mathop{\sum }_{j=1}^{m}(\unicode[STIX]{x1D709}\overline{Kn})^{j}\boldsymbol{e}_{j}^{\ast }+R_{dm} & & \displaystyle\end{eqnarray}$$

and

(4.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63F}^{\ast }=\unicode[STIX]{x1D63F}_{0}^{\ast }+\mathop{\sum }_{j=1}^{m}(\unicode[STIX]{x1D709}\overline{Kn})^{j}\unicode[STIX]{x1D640}_{j}^{\ast }+R_{Dm}, & & \displaystyle\end{eqnarray}$$

where $R_{dm}$ and $R_{Dm}$ are the residuals at the $m\text{th}$ order. When these expansions are introduced back into (4.1), the closure problem can be split in a series of closures that ensue from the identification at the successive orders in $\unicode[STIX]{x1D709}\overline{Kn}$ . Returning to the dimensional form, the closure problem at the $0\text{th}$ order takes the form

0th order:

(4.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}_{0}=0\quad \text{in }\mathscr{V}_{\unicode[STIX]{x1D6FD}} & \displaystyle\end{eqnarray}$$

(4.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\unicode[STIX]{x1D735}\boldsymbol{d}_{0}+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D63F}_{0}+\unicode[STIX]{x1D644}\quad \text{in }\mathscr{V}_{\unicode[STIX]{x1D6FD}} & \displaystyle\end{eqnarray}$$

(4.4c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}_{0}=0\quad \text{at }\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}} & \displaystyle\end{eqnarray}$$

(4.4d ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \boldsymbol{d}_{0}\rangle ^{\unicode[STIX]{x1D6FD}}=0 & \displaystyle\end{eqnarray}$$

(4.4e ) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\begin{array}{@{}ll@{}}\displaystyle \text{Periodicity} & \displaystyle \unicode[STIX]{x1D63F}_{0}(\boldsymbol{r}+\boldsymbol{l}_{i})=\unicode[STIX]{x1D63F}_{0}(\boldsymbol{r}),\\ & \displaystyle \boldsymbol{d}_{0}(\boldsymbol{r}+\boldsymbol{l}_{i})=\boldsymbol{d}_{0}(\boldsymbol{r}),\quad i=1,2,3.\end{array}\right\} & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D646}$

(4.5) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D63F}_{0}\rangle =\unicode[STIX]{x1D646}. & & \displaystyle\end{eqnarray}$$

Taking this relationship into account, along with the expansion of (4.3) and the definition of $\unicode[STIX]{x1D646}_{s}$ in (3.7e ), allows writing, at the $m\text{th}$ order

(4.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D646}_{s}\simeq \unicode[STIX]{x1D646}\boldsymbol{\cdot }\left(\unicode[STIX]{x1D644}+\mathop{\sum }_{j=1}^{m}(\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}})^{j}\unicode[STIX]{x1D64E}_{j}\right). & & \displaystyle\end{eqnarray}$$

Here $\unicode[STIX]{x1D64E}_{j}$ ( $j=1,\ldots ,m$ ) is the macroscopic coefficient obtained from the closure problem at the $j\text{th}$ order given, in its dimensional form, by

$j\text{th}$ order ( $j=1,\ldots ,m$ ):

(4.7a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}_{j}=0\quad \text{in }\mathscr{V}_{\unicode[STIX]{x1D6FD}} & \displaystyle\end{eqnarray}$$

(4.7b ) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\unicode[STIX]{x1D735}\boldsymbol{d}_{j}+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D63F}_{j}\quad \text{in }\mathscr{V}_{\unicode[STIX]{x1D6FD}} & \displaystyle\end{eqnarray}$$

(4.7c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}_{j}=-(\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }(\boldsymbol{n}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}_{j-1}+(\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}_{j-1})^{T1}))\quad \text{at }\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}} & \displaystyle\end{eqnarray}$$

(4.7d ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \boldsymbol{d}_{j}\rangle ^{\unicode[STIX]{x1D6FD}}=0 & \displaystyle\end{eqnarray}$$

(4.7e ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \unicode[STIX]{x1D63F}_{j}\rangle =\unicode[STIX]{x1D646}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}_{j} & \displaystyle\end{eqnarray}$$

(4.7f ) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\begin{array}{@{}ll@{}}\displaystyle \text{Periodicity} & \displaystyle \unicode[STIX]{x1D63F}_{j}(\boldsymbol{r}+\boldsymbol{l}_{i})=\unicode[STIX]{x1D63F}_{j}(\boldsymbol{r}),\\ & \displaystyle \boldsymbol{d}_{j}(\boldsymbol{r}+\boldsymbol{l}_{i})=\boldsymbol{d}_{j}(\boldsymbol{r}),\quad i=1,2,3.\end{array}\right\} & \displaystyle\end{eqnarray}$$

$j\text{th}$

$\boldsymbol{d}_{j}$

$\unicode[STIX]{x1D63F}_{j}$

$\boldsymbol{e}_{j}=\ell _{\unicode[STIX]{x1D6FD}}\boldsymbol{e}_{j}^{\ast }$

$\unicode[STIX]{x1D640}_{j}=\ell _{\unicode[STIX]{x1D6FD}}^{2}\unicode[STIX]{x1D640}_{j}^{\ast }$

$j=1,\ldots ,m$

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{d}_{j}=\boldsymbol{e}_{j}/(\ell _{\unicode[STIX]{x1D6FD}})^{j} & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}_{j}=\unicode[STIX]{x1D640}_{j}/(\ell _{\unicode[STIX]{x1D6FD}})^{j}. & \displaystyle\end{eqnarray}$$

$m\text{th}$

$\unicode[STIX]{x1D709}\overline{Kn}$

(4.10a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D700}\frac{\unicode[STIX]{x2202}\langle \unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\langle \unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}\langle \boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\rangle )=0 & \displaystyle\end{eqnarray}$$

(4.10b ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \boldsymbol{v}_{\unicode[STIX]{x1D6FD}}\rangle \simeq -\frac{1}{\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FD}}}\unicode[STIX]{x1D646}\boldsymbol{\cdot }\left(\unicode[STIX]{x1D644}+\mathop{\sum }_{j=1}^{m}(\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}})^{j}\unicode[STIX]{x1D64E}_{j}\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}\langle p_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}} & \displaystyle\end{eqnarray}$$

(4.10c ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}=F(\langle p_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}) & \displaystyle\end{eqnarray}$$

$(\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}})^{j}\unicode[STIX]{x1D64E}_{j}$

$j\text{th}$

$\unicode[STIX]{x1D646}$

$\unicode[STIX]{x1D64E}_{j}$

$1/\langle p_{\unicode[STIX]{x1D6FD}}\rangle ^{\unicode[STIX]{x1D6FD}}$

One must be clear that the expansion carried out in $\unicode[STIX]{x1D709}\overline{Kn}$ at the closure level is an alternative representation of the original closure problem given by (3.7) that is basically local, and this is quite a different approach from the one followed by Skjetne & Auriault (Reference Skjetne and Auriault1999) that requires an additional constraint on the Knudsen number. The approach, in this last reference, potentially allows us to take into account non-local effects that could play a role while studying gas slip momentum transport near macroscopic boundaries.

Under the form of (4.10b ), the macroscopic momentum equation clearly indicates that the apparent slip-corrected permeability, up to the first order, remains linear in $\unicode[STIX]{x1D709}\overline{Kn}$ . As a consequence, the nonlinear behaviour observed in figures 4(b) and 6(b) can only be captured by higher-order terms in the macroscopic model, and this is a major result of the present development. In this perspective, a more thorough analysis of the macroscopic slip-flow correction, beyond the first order is proposed in § 4.2.

4.2 Effective coefficients

All the closure problems providing the tensors $\unicode[STIX]{x1D646}$ and $\unicode[STIX]{x1D64E}_{j}$ have an incompressible Stokes structure, the latter being coupled to each other and to the $0\text{th}$ -order problem through the slip-like boundary condition, which makes them non-homogeneous. Therefore, a unique numerical procedure can be used to determine all the macroscopic coefficients in a simple manner.

The same computational tool as the one employed for the determination of the apparent slip-corrected permeability, $\unicode[STIX]{x1D646}_{s}$ , (see §§ 3.1 and 3.3), was used to solve the closure problems (4.4) and (4.7) up to $m=3$ on the model structure represented by the unit cell of figure 2 and $\unicode[STIX]{x1D700}$ ranging from $0.25$ to $0.8$ . Due to the symmetries in this particular case, all the macroscopic coefficients are spherical tensors ( $\unicode[STIX]{x1D646}=k\unicode[STIX]{x1D644}$ , $\unicode[STIX]{x1D64E}_{j}=s_{j}\unicode[STIX]{x1D644}$ ), thus requiring the solution of the projection of the corresponding closure problems on $\boldsymbol{e}_{x}$ (or $\boldsymbol{e}_{y}$ ) only. Fields of the $xx$ component of the tensors $\unicode[STIX]{x1D63F}_{0}$ , $\unicode[STIX]{x1D63F}_{1}$ , $\unicode[STIX]{x1D63F}_{2}$ and of the $x$ components of the vectors $\boldsymbol{d}_{0}$ , $\boldsymbol{d}_{1}$ , $\boldsymbol{d}_{2}$ are reported in their dimensionless forms in figure 8.

Figure 8. Dimensionless closure variable fields up to the second order computed on the unit cell of figure 2, $\unicode[STIX]{x1D700}=0.8$ ; (a) $D_{0xx}/\ell ^{2}$ ; (b) $D_{1xx}/\ell ^{2}$ ; (c) $D_{2xx}/\ell ^{2}$ ; (d) $d_{0x}/\ell$ ; (e) $d_{1x}/\ell$ ; (f) $d_{2x}/\ell$ .

In figure 9(a), we have represented $k^{\ast }=k/\ell ^{2}$ versus $\unicode[STIX]{x1D700}$ , which shows an excellent agreement with our previous results (Lasseux et al. Reference Lasseux, Valdes Parada, Ochoa Tapia and Goyeau2014) as well as with predictions available in the literature in the range of small (Bruschke & Advani Reference Bruschke and Advani1993) and large (Kuwabara Reference Kuwabara1959) porosities.

Figure 9. Dimensionless effective coefficients for slip flow versus $\unicode[STIX]{x1D700}$ . Unit cell of figure 2. (a) Intrinsic permeability. (b) Slip-flow corrective terms $k^{\ast }s_{j}^{\ast }$ for $j=1$ – $3$ . Note that $s_{2}^{\ast }$ is negative.

Dimensionless slip-flow corrective terms, $k^{\ast }s_{j}^{\ast }=k^{\ast }(\ell )^{j}s_{j}$ are represented versus $\unicode[STIX]{x1D700}$ in figure 9 $b$ for $j=1$ to $3$ . In this figure, we have also reported the $O(\unicode[STIX]{x1D709}\overline{Kn})$ slip-flow corrective term $k^{\ast }s^{\ast }$ appearing in the classical macroscopic model as the one reported earlier (Lasseux et al. Reference Lasseux, Valdes Parada, Ochoa Tapia and Goyeau2014). A perfect agreement with the first-order term of the present model can be observed, showing the equivalence of the two approaches restricted to the first order in $\unicode[STIX]{x1D709}\overline{Kn}$ . Higher-order terms in the present model are such that $s_{j}$ is positive when $j$ is odd and negative otherwise.

As can be inferred from the nonlinear dependence of $k_{s}$ on the Knudsen number reported above, the first-order macroscopic slip correction might become relatively inaccurate in some circumstances. This is highlighted in figure 10 where we have reported the apparent slip-corrected permeability, $k_{s}^{\ast }$ , versus $\unicode[STIX]{x1D709}\overline{Kn}$ for $\unicode[STIX]{x1D700}=0.25$ , $0.4$ , $0.6$ and $0.8.$ , keeping $\unicode[STIX]{x1D709}\overline{Kn}$ smaller than ${\sim}0.19$ so that slip flow is expected to remain physically relevant. The apparent slip-corrected permeability was computed from the closure problem (3.7), on the one hand, and from (4.6) (i.e. from the solutions of problems (4.4) and (4.7)) for $m$ up to $3$ , on the other hand. Figure 10 shows that the inaccuracy of the classical first-order approximation increases with $\unicode[STIX]{x1D709}\overline{Kn}$ and becomes much more significant while increasing $\unicode[STIX]{x1D700}$ .

Figure 10. Comparison of the dimensionless apparent slip-corrected permeability, $k_{s}^{\ast }=k_{s}/\ell ^{2}$ , computed from the solution of the closure problem (3.7), on the one hand, and estimated at the first, second and third orders with (4.6) and the solutions of expanded closure problems (4.4) and (4.7), on the other hand. Periodic model structure of figure 2: (a) $\unicode[STIX]{x1D700}=0.25$ ; (b) $\unicode[STIX]{x1D700}=0.4$ ; (c) $\unicode[STIX]{x1D700}=0.6$ ; (d) $\unicode[STIX]{x1D700}=0.8$ . Here $\overline{Kn}$ is the macroscale Knudsen number based on the fluid-phase characteristic length scale: $\unicode[STIX]{x1D709}\overline{Kn}=\unicode[STIX]{x1D709}(\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}/\ell _{\unicode[STIX]{x1D6FD}})$ .

For a more quantitative picture of the above, the relative error on the slip correction in $k_{s}$ at the $j\text{th}$ order, given by $(|(k_{s}/k)-1-\sum _{j=1}^{m}(\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}})^{j}s_{j}|)/((k_{s}/k)-1)=(|k_{s}^{\ast }-k^{\ast }(1+\sum _{j=1}^{m}(\unicode[STIX]{x1D709}\overline{\unicode[STIX]{x1D706}}_{\unicode[STIX]{x1D6FD}}^{\ast })^{j}s_{j}^{\ast })|)/(k_{s}^{\ast }-k^{\ast })$ , is reported, up to the third order, versus $\unicode[STIX]{x1D709}\overline{Kn}$ in figure 11, for the four values of $\unicode[STIX]{x1D700}$ . When $k_{s}$ is estimated at the first order, the relative error on the slip correction is ${\sim}3.7\,\%$ for $\unicode[STIX]{x1D700}=0.25$ and $\unicode[STIX]{x1D709}\overline{Kn}\sim 0.19$ but reaches ${\sim}20.3\,\%$ for $\unicode[STIX]{x1D700}=0.6$ and $\unicode[STIX]{x1D709}\overline{Kn}\sim 0.1$ . The estimation of the slip correction is significantly improved while taking into account the second-order term, or further, the third-order term. In fact, for $m=3$ and the structure under consideration, the relative error on the slip correction falls to approximately $0.6\,\%$ for $\unicode[STIX]{x1D700}=0.25$ and $\unicode[STIX]{x1D709}\overline{Kn}\sim 0.19$ or $1.7\,\%$ for $\unicode[STIX]{x1D700}=0.6$ and $\unicode[STIX]{x1D709}\overline{Kn}\sim 0.1$ . This is clearly illustrated in figure 11. If a different Knudsen number, say $\overline{Kn}_{\unicode[STIX]{x1D70E}}$ , based on the characteristic size, $\ell _{\unicode[STIX]{x1D70E}}$ , of the solid phase is used instead of $\overline{Kn}$ , one would observe that the nonlinearity appears for $\unicode[STIX]{x1D709}\overline{Kn}_{\unicode[STIX]{x1D70E}}=\ell _{\unicode[STIX]{x1D6FD}}/\ell _{\unicode[STIX]{x1D70E}}\unicode[STIX]{x1D709}\overline{Kn}\simeq 0.02$ , a value smaller than, or, when $\unicode[STIX]{x1D700}$ is large, almost equal to $\unicode[STIX]{x1D709}\overline{Kn}$ , as $\ell _{\unicode[STIX]{x1D6FD}}/\ell _{\unicode[STIX]{x1D70E}}$ takes the values $1.98\times 10^{-2}$ , $1.49\times 10^{-1}$ , $4.21\times 10^{-1}$ and $1.07$ for $\unicode[STIX]{x1D700}=0.25$ , $0.4$ , $0.6$ and $0.8$ respectively.

Figure 11. Relative error on the slip-flow correction in $k_{s}$ estimated at the first, second and third orders taking the slip correction extracted from $k_{s}$ computed from the solution of the closure problem (3.7) as the reference (see complete expression in the text). Periodic model structure of figure 2: (a) $\unicode[STIX]{x1D700}=0.25$ ; (b) $\unicode[STIX]{x1D700}=0.4$ ; (c) $\unicode[STIX]{x1D700}=0.6$ ; (d) $\unicode[STIX]{x1D700}=0.8$ .

One must however be clear about the physical meaning of the macroscopic slip-correction terms for $j\geqslant 2$ obtained from the development leading to (4.10b ). In fact, it must be kept in mind that the pore-scale physical model from which the macroscale balance equations derive involves a slip boundary condition and a momentum equation that are both first-order accurate in $\unicode[STIX]{x1D709}Kn$ . As a consequence, the macroscopic slip-correction terms beyond the first order in $\unicode[STIX]{x1D709}\overline{Kn}$ account for the contribution of microstructural effects, i.e. the signature of pore topology on the averaged slip effect that can significantly depend upon local constrictions, enlargements and curvature of pore walls through the strain rate at $\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}$ . One can note, in particular, that the macroscopic model and associated closures provide an exact macroscopic solution in the case of slip flow in a bundle of straight capillary tubes or slits made of plane parallel plates. Indeed, this solution exactly corresponds to that obtained by averaging the pore-scale flow solution. In this situation with no curvature at the pore scale in the direction of the flow, which is 1-D, the macroscopic slip correction remains at the first order that is exactly obtained from the closure problem (4.7) for $j=1$ , all the higher-order closure problems yielding $\unicode[STIX]{x1D63F}_{j}=0$ ( $j\geqslant 2$ ). To be more precise, specific general properties of the slip boundary condition and the consequences on the closure problems are highlighted in the following section.

Before moving on, it is important to mention that some results obtained from a numerical solution of an approximate version of the Boltzmann equations (i.e. Bhatnagar Gross and Krook (BGK) model) under a linearized approach, predict a linear dependence of the apparent permeability with the Knudsen number over the whole range of $\unicode[STIX]{x1D709}\overline{Kn}$ investigated here for $\unicode[STIX]{x1D700}=0.8$ (these results were provided by an anonymous reviewer and are available from the authors). The contrast with the present results would suggest that the range of application of the slip-flow model should be drastically reduced as the porosity increases. This result is however difficult to justify physically. Moreover, the linear dependence of $k_{s}^{\ast }$ with the Knudsen number resulting from this linearized BGK approach seems to hold for Knudsen number values up to approximately 0.5, fitting the predictions of Chang’s unit cell solution that is however expanded at the first order in the limit of $\unicode[STIX]{x1D709}\overline{Kn}\ll 1$ ; the reason for such a surprising behaviour remains unclear. Some possible sources for explaining this linear behaviour are the use of a linearized BGK model, along with the numerical method used to solve it.

The above remarks should not be considered as a definitive argument to accept a nonlinear dependence of the apparent permeability with the Knudsen number. In this regard, the finding of this behaviour should motivate the performance of slightly compressible rarefied gas flow experiments in highly permeable porous media.

4.3 Some important features of the slip boundary condition

An incomplete slip boundary condition, in which the $\unicode[STIX]{x1D735}\boldsymbol{v}_{\unicode[STIX]{x1D6FD}}^{T}$ term has been dropped, has often been used for porous media flow (Skjetne & Auriault Reference Skjetne and Auriault1999; Pavan & Chastanet Reference Pavan and Chastanet2011), and we have shown, in accordance with some previous references (Barber et al. Reference Barber, Sun, Gu and Emerson2004; Lockerby et al. Reference Lockerby, Reese, Emerson and Barber2004; Panzer et al. Reference Panzer, Liu and Einzel1992), that this physically unjustified form can lead to a significant error in some circumstances (see § 3.3). Together with the nonlinear dependence of the slip correction on $\unicode[STIX]{x1D709}\overline{Kn}$ , this suggests a close attention to the slip-like boundary condition of (4.7c ) in the $j\text{th}$ -order closure problem.

The analysis starts with an alternative expression of the second term on the right-hand side of (4.7c ) given by

(4.11) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }(\boldsymbol{n}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}_{j-1})^{T1}) & = & \displaystyle (\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}_{j-1})\nonumber\\ \displaystyle & & \displaystyle -\,(\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}_{j-1})\quad \text{at }\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}.\end{eqnarray}$$

For all $j$ , $\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}_{j-1}=0$ at $\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}$ since $\unicode[STIX]{x1D63F}_{j-1}$ is either $0$ ( $j=1$ ) or purely tangential to $\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}$ at this interface ( $j\geqslant 2$ ). Consequently, $\unicode[STIX]{x1D735}(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}_{j-1})$ is purely normal to $\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}$ and this allows rewriting (4.11) as

(4.12) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }(\boldsymbol{n}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}_{j-1})^{T1})=-\unicode[STIX]{x1D735}_{s}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}_{j-1}\quad \text{at }\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D735}_{s}$ denotes the surface differentiation operator defined by $\unicode[STIX]{x1D735}_{s}\equiv (\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ .

For $j=1$ , this result can be further simplified due to the no-slip boundary condition on $\unicode[STIX]{x1D63F}_{0}$ (see (4.4c )) leading to

(4.13) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }(\boldsymbol{n}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}_{0})^{T1})=0\quad \text{at }\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}} & & \displaystyle\end{eqnarray}$$

and this proves that the $\unicode[STIX]{x1D735}\boldsymbol{v}_{\unicode[STIX]{x1D6FD}}^{T}$ term in the slip boundary condition does not play any role up to the first order of the expansion in $\unicode[STIX]{x1D709}\overline{Kn}$ at the closure level. At higher orders, (4.12) cannot be further simplified and indicates that the contribution of $(\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }(\boldsymbol{n}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}_{j-1})^{T1})$ depends on $\unicode[STIX]{x1D735}_{s}\boldsymbol{n}$ that is related to the curvature of $\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}$ . In 2-D, for instance, it can be easily demonstrated (see appendix B) that $\unicode[STIX]{x1D735}_{s}\boldsymbol{n}=-\unicode[STIX]{x1D705}(\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})$ where $\unicode[STIX]{x1D705}=-\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }\boldsymbol{n}$ is exactly the curvature of $\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}$ . This clearly shows that the $(\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}_{j-1})^{T1}$ term does play a role in the slip-like boundary condition (4.7c ) of the closure problem at the $j\text{th}$ order, $j\geqslant 2$ , in the general case. Its contribution vanishes only in some special geometrical situations where the flow is 1-D and the curvature of $\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}$ is zero in the flow direction. This last remark applies to the boundary condition of (2.1d ) at the pore scale and justifies why the incomplete form of the slip boundary condition can be used in the case of straight tubes, slits or channels that have been extensively used for the study of Knudsen effects in the slip regime (Fishman & Hetsroni Reference Fishman and Hetsroni2005; Lauga & Cossu Reference Lauga and Cossu2005; Shen et al. Reference Shen, Chen, Crone and Anaya-Dufresne2007).

The analysis is pursued with an alternative expression of the first term on the right-hand side of (4.7c ), which can be written as

(4.14) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}_{j-1})=\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(\boldsymbol{n}\unicode[STIX]{x1D63F}_{j-1})+\unicode[STIX]{x1D705}\unicode[STIX]{x1D63F}_{j-1}\quad \text{at }\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}. & & \displaystyle\end{eqnarray}$$

This finally allows us to express the slip-like boundary condition in (4.7c ) for the $j\text{th}$ -order closure problem as

(4.15) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63F}_{j}=-\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(\boldsymbol{n}\unicode[STIX]{x1D63F}_{j-1})+(\unicode[STIX]{x1D735}_{s}\boldsymbol{n}-\unicode[STIX]{x1D705}\unicode[STIX]{x1D644})\boldsymbol{\cdot }\unicode[STIX]{x1D63F}_{j-1}\quad \text{at }\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}. & & \displaystyle\end{eqnarray}$$

In 2-D, this last expression takes the form $\unicode[STIX]{x1D63F}_{j}=-\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(\boldsymbol{n}\unicode[STIX]{x1D63F}_{j-1})-2\unicode[STIX]{x1D705}\unicode[STIX]{x1D63F}_{j-1}$ . For $j=1$ , it simplifies to

(4.16) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63F}_{1}=-(\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}_{0})=-\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(\boldsymbol{n}\unicode[STIX]{x1D63F}_{0})\quad \text{at }\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}. & & \displaystyle\end{eqnarray}$$

Equations (4.15) and (4.16) show that the slip-like boundary condition in the $j\text{th}$ -order closure problem contains an explicit dependence on the curvature of $\mathscr{A}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}$ that is however filtered out at the first order in $\unicode[STIX]{x1D709}\overline{Kn}$ . This explains why the classical first-order macroscopic slip-flow model remains inaccurate for systems presenting a significant curvature of the solid–fluid interface. For the periodic square pattern of circular cylinders, it clearly justifies why the error on the slip correction, when restricted to the first order, increases with the porosity, as curvature scales as $(1-\unicode[STIX]{x1D700})^{-1/2}$ in that case. This also justifies the necessity of taking into account higher-order terms in the macroscopic model as in (4.10b ).

Finally, it must be noted that the analysis carried out in this section may be straightforwardly applied to an expansion in the Knudsen number performed on the original initial boundary value problem given in (2.1). In that way, the slip boundary condition properties at the successive orders in $Kn$ may be generalized far beyond the strict context of porous media flow, i.e. for any situation allowing such an expansion in $Kn$ .