4.1.1 Conditions at the boundaries of the macroscopic porous domain

Darcy’s equation can be written in terms of the pore pressure at leading order by combining (4.22) and (4.25), to yield

(4.27) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}X_{k}}\left[{\mathcal{K}}_{ki}\left(\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X_{i}}-f_{i}^{(0)}\right)\right]=0.\end{eqnarray}$$

Close to boundaries the assumption of spatial homogeneity drops and a boundary layer appears. For example, at a solid wall (say, placed on a plane of constant $Y$ ) which bounds a porous, isotropic region, the non-penetration condition suggests imposing $\unicode[STIX]{x2202}P^{(0)}/\unicode[STIX]{x2202}Y=0$ . This ‘intuitive’ condition is confirmed by a boundary layer analysis, as indicated by Carraro, Marǔsić-Paloka & Mikelić (Reference Carraro, Marǔsić-Paloka and Mikelić2018).

Particular care must be used in treating the interface between a free-fluid region and a porous domain and the commonly employed condition of pressure continuity at the dividing surface (Ene & Sanchez-Palencia Reference Ene and Sanchez-Palencia1975; Lācis & Bagheri Reference Lācis and Bagheri2017) is an approximation which, according to Carraro et al. (Reference Carraro, Marǔsić-Paloka and Mikelić2018), is acceptable only for the case of isotropic porous media. Recently, Lācis et al. (Reference Lācis, Sudhakar, Pasche and Bagheri2019) have proposed a pressure jump condition obtained from a normal-momentum balance through the dividing surface; this condition, briefly described later in § 7.1.2, extends and complements that proposed by previous researchers.

4.2.1 Properties of the adjoint function of a linear, parabolic problem

Let us consider the function space of all real continuous functions $f(t)$ defined over $0\leqslant t\leqslant T$ . The inner product between functions in this space can be defined by

(4.28) $$\begin{eqnarray}f_{1}\cdot f_{2}:=\int _{0}^{T}f_{1}(t)f_{2}(t)\,\text{d}t.\end{eqnarray}$$

Given a parabolic, linear equation, to be integrated from $t=0$ to $t=T$ and defined in terms of the operator ${\mathcal{L}}$ as

(4.29) $$\begin{eqnarray}\frac{\text{d}f}{\text{d}t}={\mathcal{L}}f+S,\end{eqnarray}$$

with $S$ a source term, it is possible to introduce the adjoint equation, to be integrated backwards, i.e. proceeding from $t=T$ to $t=0$ , as

(4.30) $$\begin{eqnarray}-\frac{\text{d}f^{\dagger }}{\text{d}t}={\mathcal{L}}^{\dagger }f^{\dagger },\end{eqnarray}$$

with the property that

(4.31) $$\begin{eqnarray}f(T)f^{\dagger }(T)=f(0)f^{\dagger }(0)+S\cdot f^{\dagger }.\end{eqnarray}$$

The above stems immediately from the defining property of adjoint operators, i.e.

(4.32) $$\begin{eqnarray}{\mathcal{L}}f\cdot f^{\dagger }=f\cdot {\mathcal{L}}^{\dagger }f^{\dagger }.\end{eqnarray}$$

The adjoint function $f^{\dagger }$ plays the role of a Green’s function for equation (4.29) for, if we integrate (4.30) using $f^{\dagger }(T)=1$ as the terminal condition, the direct solution at the generic final time $t=T$ is simply

(4.33) $$\begin{eqnarray}f(T)=f(0)f^{\dagger }(0)+S\cdot f^{\dagger }.\end{eqnarray}$$

In other words, once the adjoint solution, viz. the Green’s function, $f^{\dagger }(t)$ , is available, the final direct solution $f(T)$ is found by simply multiplying the Green’s function at $t=0$ with any initial condition $f(0)$ , and summing this to the inner product of $f^{\dagger }(t)$ with any source term, $S$ , of the direct equation. The extension to multidimensional problems is straightforward and is addressed below on the actual problem of interest.

4.2.2 Adjoint of the nonlinear time-dependent, microscopic equations

The adjoint system of equations (4.6) and (4.9) can be derived in a similar manner as in the linear case, carrying out integrations in both space and time. The time integral must run from the initial time, set arbitrarily to 0, to a generic final time, $T$ . After a few integrations by parts it is simple to find

(4.34) $$\begin{eqnarray}\displaystyle & & \displaystyle \int _{0}^{T}\left(\frac{\unicode[STIX]{x2202}u_{i}^{\dagger }}{\unicode[STIX]{x2202}x_{i}},p^{(1)}\right)+\left(Re\left[\frac{\unicode[STIX]{x2202}u_{i}^{\dagger }}{\unicode[STIX]{x2202}t}+u_{j}^{(0)}\frac{\unicode[STIX]{x2202}u_{i}^{\dagger }}{\unicode[STIX]{x2202}x_{j}}\right]-\frac{\unicode[STIX]{x2202}p^{\dagger }}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{\dagger }}{\unicode[STIX]{x2202}x_{j}^{2}},u_{i}^{(0)}\right)\,\text{d}t\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{0}^{T}\left(u_{i}^{\dagger },\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X_{i}}-f_{i}^{(0)}\right)\,\text{d}t+(u_{i}^{\dagger },u_{i}^{(0)})|_{t=0}^{t=T}.\end{eqnarray}$$

The auxiliary problem to be solved with periodicity at the boundaries of the RVE and no-slip at the fluid/grain contact surfaces reads

(4.35a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{i}}=0;\quad -Re\left[\frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}t}+u_{j}^{(0)}\frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{j}}\right]=-\frac{\unicode[STIX]{x2202}p^{\dagger (k)}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{j}^{2}}+S_{ki}^{\dagger },\end{eqnarray}$$

with the source term $S_{ki}^{\dagger }$ equal to the Kronecker delta, $\unicode[STIX]{x1D6FF}_{ki}$ , if the problem is steady, and $S_{ki}^{\dagger }=0$ in a time-dependent problem. In this latter case, the terminal conditions at $t=T$ for the adjoint problem will need to be defined. As before, the superscript $(k)$ in the adjoint variables in (4.35) stems from the need to define and solve three different auxiliary problems, $k=1,2$ and 3, for each of the velocity components.

The steady problem

The adjoint momentum equation is simply

(4.36) $$\begin{eqnarray}-Re\,u_{j}^{(0)}\frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{j}}=-\frac{\unicode[STIX]{x2202}p^{\dagger (k)}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{j}^{2}}+\unicode[STIX]{x1D6FF}_{ki},\end{eqnarray}$$

and the resulting macroscopic problem reads

(4.37) $$\begin{eqnarray}\langle u_{k}^{(0)}\rangle =-{\mathcal{K}}_{ki}^{eff}\left[\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X_{i}}-f_{i}^{(0)}\right],\end{eqnarray}$$

with

(4.38) $$\begin{eqnarray}{\mathcal{K}}_{ki}^{eff}=\langle u_{i}^{\dagger (k)}\rangle .\end{eqnarray}$$

The resulting macroscale equation (4.37) is the same as in the creeping flow case (i.e. the Darcy equation (4.22)), but now the adjoint variables (and consequently the permeability) do not depend only on the geometric small-scale properties of the porous medium, but also on the direct flow state; the permeability is thus called effective (or apparent) and noted with the superscript $^{eff}$ . As opposed to the Stokes permeability of equation (4.23), the apparent permeability tensor is in general not symmetric; Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2017) have shown that it can be decomposed into the sum of a symmetric and a skew-symmetric part, the latter component originating from inertial transport. Resolution of (4.36) is indispensable to obtain all components of the apparent permeability tensor, for any given macroscopic forcing $\unicode[STIX]{x2202}p^{(0)}/\unicode[STIX]{x2202}X_{i}-f_{i}^{(0)}$ . For simplicity, from now on, the $f_{i}^{(0)}$ volumetric force will be taken conservative and absorbed into the pressure gradient term.

The unsteady problem

Because of the minus sign in front of the first term of the adjoint momentum equation in (4.35), and the presence of a positive diffusion term on the right-hand side, the only stable direction of evolution for the equation is for time $t$ running from $T$ to 0. It is then convenient to introduce the time variable $\unicode[STIX]{x1D70F}=T-t$ , which runs from $\unicode[STIX]{x1D70F}=0$ to $\unicode[STIX]{x1D70F}=T$ , for the equation to read

(4.39) $$\begin{eqnarray}Re\left[\frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}-u_{j}^{(0)}\frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{j}}\right]=-\frac{\unicode[STIX]{x2202}p^{\dagger (k)}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{j}^{2}}.\end{eqnarray}$$

From (4.34) we are left with

(4.40) $$\begin{eqnarray}0=\int _{0}^{T}\left(u_{i}^{\dagger (k)},\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X_{i}}\right)\,\text{d}t+(u_{i}^{\dagger (k)},u_{i}^{(0)})|_{t=0}^{t=T};\end{eqnarray}$$

as initial condition for the adjoint variable at $\unicode[STIX]{x1D70F}=0$ we can choose either

(4.41)

With the first choice we recover the phase-averaged velocity at any final time $t=T$ , whereas the second choice yields the instantaneous, pointwise velocity distribution in the RVE. In other words, choice 1 yields

(4.42) $$\begin{eqnarray}\langle u_{k}^{(0)}\rangle |_{t=T}=(u_{i}^{\dagger (k)},u_{i}^{(0)})|_{t=0}-\int _{0}^{T}{\mathcal{K}}_{ki}^{eff}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X_{i}}\,\text{d}t,\end{eqnarray}$$

with ${\mathcal{K}}_{ki}^{eff}=\langle u_{i}^{\dagger (k)}\rangle$ the effective dynamic permeability, function of $\unicode[STIX]{x1D70F}=T-t$ . The last integral above is a convolution when the imposed pressure gradient is time dependent. Choice 2 gives

(4.43) $$\begin{eqnarray}u_{k}^{(0)}(T,\boldsymbol{x}^{\prime },\boldsymbol{X})={\mathcal{V}}(u_{i}^{\dagger (k)},u_{i}^{(0)})|_{t=0}-{\mathcal{V}}\int _{0}^{T}\left(u_{i}^{\dagger (k)},\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X_{i}}\right)\,\text{d}t.\end{eqnarray}$$

Both choices illustrate nicely the fact that adjoint fields operate as Green’s functions, expressing the sensitivity to either the initial condition, $u_{i}^{(0)}|_{t=0}$ , or the source term, $\unicode[STIX]{x2202}p^{(0)}/\unicode[STIX]{x2202}X_{i}$ . In the study of the flow through porous or poroelastic media the average velocity is usually sought for, immediately available from choice 1. Recently, the problem of the unsteady seepage of fluid in a porous medium has been posed and solved by Lasseux, Valdés-Parada & Bellet (Reference Lasseux, Valdés-Parada and Bellet2019), within Whitaker’s (Reference Whitaker1996) volume-averaging framework.

4.2.3 Further considerations on the steady case

A microscopic problem similar to that given in equation (4.36) has been given previously by Whitaker (Reference Whitaker1996) in terms of two closure variables, noted below as $\boldsymbol{e}=\{e_{k}\}$ and $\boldsymbol{E}=\{E_{ik}\}$ . The most notable difference with Whitaker’s closure problem is that his stationary closure equations are

(4.44a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}E_{ik}}{\unicode[STIX]{x2202}x_{i}}=0,\quad Re\,u_{j}^{(0)}\frac{\unicode[STIX]{x2202}E_{ik}}{\unicode[STIX]{x2202}x_{j}}=-\frac{\unicode[STIX]{x2202}e_{k}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}E_{ik}}{\unicode[STIX]{x2202}x_{j}^{2}}+\unicode[STIX]{x1D6FF}_{ik},\end{eqnarray}$$

with periodicity for $E_{ik}$ and $e_{k}$ on the cell boundary and zero value for $E_{ik}$ at $\unicode[STIX]{x1D6E4}$ ; the effective permeability is

(4.45) $$\begin{eqnarray}{\mathcal{K}}_{ik}^{eff}=\langle E_{ik}\rangle .\end{eqnarray}$$

In other words, the advective term on the left-hand side of (4.44) contains $u_{j}^{(0)}$ , whereas it contains $-u_{j}^{(0)}$ in model (4.36). This apparent difference is easily understood upon closer scrutiny of the two models, which turn out to be simply the transpose of one another, as shown by Lasseux & Valdés-Parada (private communication, 2018). Upon solving with the adjoint formulation, for $k=p$ , with $p=1,2$ or $3$ , one recovers $K_{p1}^{eff}$ , $K_{p2}^{eff}$ and $K_{p3}^{eff}$ . Conversely, the solution of (4.44) yields for any $k=p$ : $K_{1p}^{eff}$ , $K_{2p}^{eff}$ and $K_{3p}^{eff}$ . The final entries of the apparent permeability tensor in the two cases are the same.

Figure 12.

Streamlines of the steady, adjoint velocity field $u_{i}^{\dagger (1)}$ in the unit cell, built from the solution of (4.36) for $k=1$ , isotropic (a) and anisotropic (b) grains. The direct flow is driven diagonally (at $45^{\circ }$ ) by a macroscopic pressure gradient in both cases.

To illustrate the above we consider two different types of inclusion within a two-dimensional unit cell. The first inclusion is cylindrical while the second is highly irregular, they are shown in figure 12. The direct fields in the two cases are computed by solving equations (4.6) and (4.9) (the latter in its steady form) in the periodic cell, subject to a pressure gradient inclined at $45^{\circ }$ with respect to the horizontal axis of the unit cell. In particular, in the first case (cylindrical grains with porosity of the medium $\unicode[STIX]{x1D703}=0.80$ ) we take $\unicode[STIX]{x2202}p^{(0)}/\unicode[STIX]{x2202}X=\unicode[STIX]{x2202}p^{(0)}/\unicode[STIX]{x2202}Y=-74.32$ , and in the second (asymmetric inclusions with porosity $\unicode[STIX]{x1D703}=0.734$ ) we take $\unicode[STIX]{x2202}p^{(0)}/\unicode[STIX]{x2202}X=\unicode[STIX]{x2202}p^{(0)}/\unicode[STIX]{x2202}Y=-639.4$ . In the symmetric case the microscopic Reynolds number is $Re=80.7$ while in the second case it is $782$ ; the filtration velocity components in the two cases are

(4.46) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\mathit{isotropic\;case}:~\langle u^{(0)}\rangle =\langle v^{(0)}\rangle =0.6666,\\ \mathit{anisotropic\;case}:~\langle u^{(0)}\rangle =0.8467,\quad \langle v^{(0)}\rangle =9.720\times 10^{-2}.\end{array}\right\}\end{eqnarray}$$

Once the direct fields are computed, they are used in (4.36) to obtain the adjoint velocity and pressure. When $k=1$ the adjoint velocity components are denoted by $u_{i}^{\dagger (1)}$ and from them we can draw adjoint streamlines, represented in figure 12 for the two configurations at hand. For the case of isotropic grains it is found through averaging: ${\mathcal{K}}_{11}^{eff}={\mathcal{K}}_{22}^{eff}=6.331\times 10^{-3}$ ; ${\mathcal{K}}_{12}^{eff}={\mathcal{K}}_{21}^{eff}=2.635\times 10^{-3}$ . Using Darcy’s law with the apparent permeability tensor we thus have

(4.47) $$\begin{eqnarray}\langle u^{(0)}\rangle =\langle v^{(0)}\rangle =-{\mathcal{K}}_{11}^{eff}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X}-{\mathcal{K}}_{12}^{eff}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}Y}=0.6665,\end{eqnarray}$$

in excellent agreement with the direct simulation result.

The irregular grain geometry produces the following result for the apparent permeability tensor:

(4.48) $$\begin{eqnarray}\left(\begin{array}{@{}cc@{}}8.65\times 10^{-4} & 4.59\times 10^{-4}\\ -1.71\times 10^{-4} & 3.23\times 10^{-4}\end{array}\right)\end{eqnarray}$$

while, for comparison purposes, the intrinsic (Stokes) permeability components ${\mathcal{K}}_{ki}$ are

(4.49) $$\begin{eqnarray}\left(\begin{array}{@{}cc@{}}5.53\times 10^{-3} & -2.25\times 10^{-4}\\ -2.25\times 10^{-4} & 6.71\times 10^{-3}\end{array}\right).\end{eqnarray}$$

Using Darcy equation with the effective permeability components we recover

(4.50) $$\begin{eqnarray}\displaystyle & \displaystyle \langle u^{(0)}\rangle =-{\mathcal{K}}_{11}^{eff}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X}-{\mathcal{K}}_{12}^{eff}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}Y}=0.8466, & \displaystyle\end{eqnarray}$$

(4.51) $$\begin{eqnarray}\displaystyle & \displaystyle \langle v^{(0)}\rangle =-{\mathcal{K}}_{21}^{eff}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X}-{\mathcal{K}}_{22}^{eff}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}Y}=9.720\times 10^{-2}, & \displaystyle\end{eqnarray}$$

which, again, agree very well with the results from the direct simulation.

The presence of an inertial term in the macroscopic Darcy equation can be made explicit by writing ${\mathcal{K}}_{ki}^{eff}$ as

(4.52) $$\begin{eqnarray}{\mathcal{K}}_{ki}^{eff}=({\mathcal{F}}_{jk}+\unicode[STIX]{x1D6FF}_{jk})^{-1}{\mathcal{K}}_{ji},\end{eqnarray}$$

with ${\mathcal{F}}_{jk}$ the non-symmetric, $u_{j}^{(0)}$ -dependent Forchheimer tensor; equation (4.37) then becomes the Darcy–Forchheimer equation (Whitaker Reference Whitaker1996)

(4.53) $$\begin{eqnarray}\langle u_{k}^{(0)}\rangle =-{\mathcal{K}}_{ki}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X_{i}}-{\mathcal{F}}_{ki}\langle u_{i}^{(0)}\rangle .\end{eqnarray}$$

Clearly, as $Re\rightarrow 0$ Darcy equation is recovered; when $Re$ increases, the so-called regime of weak inertia is entered and ${\mathcal{F}}_{ki}\langle u_{i}^{(0)}\rangle \sim \langle u_{k}^{(0)}\rangle \langle u_{i}^{(0)}\rangle \langle u_{i}^{(0)}\rangle$ (Mei & Auriault Reference Mei and Auriault1991; Firdaouss, Guermond & Le Quéré Reference Firdaouss, Guermond and Le Quéré1997). Upon further increasing $Re$ , the strong inertia regime is attained with ${\mathcal{F}}_{ki}\langle u_{i}^{(0)}\rangle$ proportional to the square of the seepage velocity (Whitaker Reference Whitaker1996; Lasseux, Abbasian Arani & Ahmadi Reference Lasseux, Abbasian Arani and Ahmadi2011). For yet larger values of the Reynolds numbers, the flow in the porous medium undergoes successive bifurcations, before transition to a turbulent state. The onset of unsteady motion within the pores of the medium depends on the geometry, the dimensions and the arrangement of the grains within the porous domain.

4.2.4 Possible solution strategies

In practical calculations, to account for the dependence of ${\mathcal{K}}_{ki}^{eff}$ on the microstate when convection terms are non-negligible, both direct and adjoint states are needed. For any given geometry of the inclusions, the parameters to be considered are the orientation of the driving pressure force and the local Reynolds number.

Extensive parametric calculations under steady conditions have been reported by Lasseux et al. (Reference Lasseux, Abbasian Arani and Ahmadi2011). The two-dimensional case was considered, with solid inclusions of either circular or square shape; interestingly, Lasseux et al. (Reference Lasseux, Abbasian Arani and Ahmadi2011) considered both ordered and disordered arrangements of grains, in the latter case also including grains of varying sizes (and this implied the use of periodic RVEs of rather large dimensions). For the case of ordered structures the Forchheimer tensor ${\mathcal{F}}_{jk}$ was found to be non-symmetric, except when the applied pressure gradient was directed along a symmetry axis of the representative elementary cell. This stems from the fact that, for an arbitrary orientation of $\unicode[STIX]{x2202}p^{(0)}/\unicode[STIX]{x2202}X_{i}$ , the force on the grain due to inertia is not necessarily aligned with the mean flow. Another interesting result was that the range of Reynolds numbers in which the weak inertia form of the Forchheimer equation occurred was more narrow when disordered structures were considered, possibly explaining why experiments with natural media often did not exhibit such a regime.

An alternative strategy, requiring consideration of only adjoint variables, has been proposed by Valdés-Parada, Lasseux & Bellet (Reference Valdés-Parada, Lasseux and Bellet2016). It consists in solving the system

(4.54a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{i}}=0;\quad Re\,u_{l}^{\dagger (j)}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X_{l}}\frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{j}}=-\frac{\unicode[STIX]{x2202}p^{\dagger (k)}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{j}^{2}}+\unicode[STIX]{x1D6FF}_{ik},\end{eqnarray}$$

which holds since the microscopic equation

(4.55) $$\begin{eqnarray}u_{j}^{(0)}=-u_{l}^{\dagger (j)}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X_{l}},\end{eqnarray}$$

applies in the unit cell (Mei & Vernescu Reference Mei and Vernescu2010).

A less involved approximation consists in expressing the direct velocity field in the RVE as the sum of a mean part plus a periodic deviation of zero mean, i.e.

(4.56) $$\begin{eqnarray}u_{j}^{(0)}(x_{i},X_{i})=\langle u_{j}^{(0)}\rangle ^{f}(X_{i})+u_{j}^{\prime }(x_{i},X_{i})=\unicode[STIX]{x1D703}^{-1}\langle u_{j}^{(0)}\rangle (X_{i})+u_{j}^{\prime }(x_{i},X_{i}),\end{eqnarray}$$

and then neglecting the fluctuations, i.e. assuming $u_{j}^{(0)}=\unicode[STIX]{x1D703}^{-1}\langle u_{j}^{(0)}\rangle$ to replace the actual direct field in the microscopic problem (4.36) with its intrinsic average. The results of this Oseen-like procedure deteriorate as the Reynolds number increases, as shown in the representative results of figure 13 for the isotropic, regular medium sketched in the image, when compared to results obtained from the direct solver (labelled as ‘DNS’ in the figure). For $\unicode[STIX]{x1D703}=0.8$ and $Re=0$ the dimensionless permeability of the medium is $1.941\times 10^{-2}$ ; as inertia becomes prominent the permeability decreases as initially shown by Edwards et al. (Reference Edwards, Shapiro, Bar-Yoseph and Shapira1990). The results by these authors are however numerically under-resolved; fine grid direct simulations, displayed with asterisks in the figure, yield slightly larger values of the apparent permeability component, particularly for $Re$ exceeding 100. Values identical to those of the direct solver, to within graphical accuracy, are found when solving (4.54). Conversely, the Oseen-like approach produces approximate results which almost coincide with those by Edwards et al. (Reference Edwards, Shapiro, Bar-Yoseph and Shapira1990) over the whole range of $Re$ examined.

Figure 13.

Variation of the apparent permeability with the microscopic Reynolds number for a two-dimensional medium composed by regularly arranged circular grains.

Another possible approach to treat flows with inertial effects through porous media consists in creating a database of direct and adjoint states, possibly using a metamodel based on kriging interpolation, or any other suitable interpolation technique, to fill the data gaps which inevitably arise when discretizing a very large space of parameters, thus generating a look-up table of permeabilities to be employed, for example, when the pressure gradient changes intensity and/or direction, or when the porosity is not homogeneous. This approach has been followed by Luminari et al. (Reference Luminari, Airiau and Bottaro2018) for a medium formed by a three-dimensional staggered array of cylindrical fibres. Interestingly, for such a case the apparent permeability tensor remains strongly diagonally dominant, regardless of the orientation and intensity of the driving force, a fact which must be correlated to the transversely isotropic nature of the medium considered.

6.1.1 The small roughness limit

When the wall roughness is of sufficiently small amplitude for the flow in its proximity to be dominated by viscosity, e.g. when ${\mathcal{R}}$ is of order $\unicode[STIX]{x1D716}$ or smaller, (6.4) and (6.5) reduce at leading order to the Stokes problem

(6.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{(0)}}{\unicode[STIX]{x2202}x_{i}}=0, & \displaystyle\end{eqnarray}$$

(6.7) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{(0)}}{\unicode[STIX]{x2202}x_{j}^{2}}, & \displaystyle\end{eqnarray}$$

and the RVE is the unit cell. The equations must be supplemented by no-slip conditions at the wall, $y_{wall}$ , and periodicity at the cell boundaries in $x$ and $z$ . The conditions at the upper boundary of the domain will arise as part of the solution in terms of outer, macroscopic variables.

We multiply equations (6.6) and (6.7) by $p^{\dagger }$ and $u_{i}^{\dagger }$ , then sum and integrate in space, over ${\mathcal{V}}_{f}$ . Using Lagrange–Green identity it is easy to find

(6.8) $$\begin{eqnarray}0=\left(p^{(0)},\frac{\unicode[STIX]{x2202}u_{i}^{\dagger }}{\unicode[STIX]{x2202}x_{i}}\right)+\left(u_{i}^{(0)},\left[-\frac{\unicode[STIX]{x2202}p^{\dagger }}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{\dagger }}{\unicode[STIX]{x2202}x_{j}^{2}}\right]\right)+\text{b.t.}\end{eqnarray}$$

with b.t. denoting boundary terms coming from integration by parts. The adjoint variables are taken to solve the system

(6.9) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\dagger }}{\unicode[STIX]{x2202}x_{i}}=0, & \displaystyle\end{eqnarray}$$

(6.10) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\unicode[STIX]{x2202}p^{\dagger }}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{\dagger }}{\unicode[STIX]{x2202}x_{j}^{2}}, & \displaystyle\end{eqnarray}$$

together with ‘no-slip’ conditions for $u_{i}^{\dagger }$ at $y=y_{wall}$ and periodicity on the lateral boundaries of the unit cell for both adjoint velocity components and adjoint pressure. The equation

(6.11) $$\begin{eqnarray}\text{b.t.}=\int _{\unicode[STIX]{x1D6FA}}-p^{(0)}v^{\dagger }+v^{(0)}p^{\dagger }+u_{i}^{\dagger }\frac{\unicode[STIX]{x2202}u_{i}^{(0)}}{\unicode[STIX]{x2202}y}-u_{i}^{(0)}\frac{\unicode[STIX]{x2202}u_{i}^{\dagger }}{\unicode[STIX]{x2202}y}\,\text{d}x\,\text{d}z=0,\end{eqnarray}$$

must be satisfied, with all variables evaluated at $y\rightarrow +\infty$ . We still have some freedom in choosing the adjoint conditions at the upper boundary of the microscopic domain and can thus set up two distinct auxiliary problems, which will be denoted below by the superscript $(k)$ on the adjoint variables, with $k$ which can be equal to either 1 or 3

(6.12a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}y}=\unicode[STIX]{x1D6FF}_{ik},\quad v^{\dagger (k)}:=u_{2}^{\dagger (k)}=0,\quad \text{at }y\rightarrow +\infty ;\end{eqnarray}$$

also the index $i$ can take the values 1 and 3 to denote, respectively, the components of the adjoint vector along $x$ and $z$ , so that from (6.11)–(6.12) the surface-averaged velocities become

(6.13) $$\begin{eqnarray}\displaystyle & \displaystyle \lim _{y\rightarrow +\infty }\lceil u^{(0)}\rceil :=\lim _{y\rightarrow +\infty }\lceil u_{1}^{(0)}\rceil =\lim _{y\rightarrow +\infty }\left\lceil u^{\dagger (1)}\frac{\unicode[STIX]{x2202}u^{(0)}}{\unicode[STIX]{x2202}y}+w^{\dagger (1)}\frac{\unicode[STIX]{x2202}w^{(0)}}{\unicode[STIX]{x2202}y}+p^{\dagger (1)}v^{(0)}\right\rceil ; & \displaystyle\end{eqnarray}$$

(6.14) $$\begin{eqnarray}\displaystyle & \displaystyle \lim _{y\rightarrow +\infty }\lceil w^{(0)}\rceil :=\lim _{y\rightarrow +\infty }\lceil u_{3}^{(0)}\rceil =\lim _{y\rightarrow +\infty }\left\lceil u^{\dagger (3)}\frac{\unicode[STIX]{x2202}u^{(0)}}{\unicode[STIX]{x2202}y}+w^{\dagger (3)}\frac{\unicode[STIX]{x2202}w^{(0)}}{\unicode[STIX]{x2202}y}+p^{\dagger (3)}v^{(0)}\right\rceil . & \displaystyle\end{eqnarray}$$

As already stated after (4.17), the superscript $(k)$ next to the name of an adjoint variable stands for the row of a matrix, and not for some asymptotic order. Thus, in the present case, $u_{i}^{\dagger (k)}$ is the component of a $2\times 2$ matrix.

In practice, the upper boundary conditions are enforced at $y=y_{\infty }$ sufficiently far away from the wall. The operator $\lceil \cdot \rceil$ introduced above denotes averaging over a planar surface $\unicode[STIX]{x1D6FA}$ of the RVE at any given position $y$ , i.e.

(6.15) $$\begin{eqnarray}\lceil a\rceil =\frac{1}{\unicode[STIX]{x1D6FA}}\int _{\unicode[STIX]{x1D6FA}_{f}}a\,\text{d}x\,\text{d}z,\end{eqnarray}$$

with $\unicode[STIX]{x1D6FA}_{f}$ the area occupied by the fluid; $\unicode[STIX]{x1D6FA}$ and $\unicode[STIX]{x1D6FA}_{f}$ might differ, for example when averaging over a plane which cuts through roughness elements. This definition is the planar counterpart of volumetric phase averaging (4.11).

The third constraint is simply

(6.16) $$\begin{eqnarray}\lceil v^{(0)}\rceil =0,\end{eqnarray}$$

at any $y$ , because of mass conservation and periodicity.

Provided $y_{\infty }$ is taken large enough, it is found that the adjoint variables there are constant over $\unicode[STIX]{x1D6FA}$ ; such values will be denoted by a $\infty$ subscript. Equation (6.13) expressed in outer variables, defined so that $U_{i}=\unicode[STIX]{x1D716}\lceil u_{i}^{(0)}\rceil$ , reduces to

(6.17) $$\begin{eqnarray}U|_{Y_{\infty }}=\unicode[STIX]{x1D716}\left(u_{\infty }^{\dagger (1)}\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Y}+w_{\infty }^{\dagger (1)}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Y}\right)\biggm\vert_{Y_{\infty }},\end{eqnarray}$$

and similarly for $W|_{Y_{\infty }}$ . The simple Taylor expansion

(6.18) $$\begin{eqnarray}U|_{Y_{\infty }}=U|_{Y=0}+Y_{\infty }\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Y}\biggm\vert_{Y=0}+\cdots \,,\end{eqnarray}$$

permits us to transfer the boundary conditions to the reference surface in $Y=0$ , to recover the classical Navier condition at $Y=0$ , correct to order $\unicode[STIX]{x1D716}$ , i.e.

(6.19) $$\begin{eqnarray}U_{s}=U|_{Y=0}=\unicode[STIX]{x1D716}\left(\unicode[STIX]{x1D706}_{x}\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Y}+\unicode[STIX]{x1D706}_{xz}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Y}\right)\biggm\vert_{Y=0},\end{eqnarray}$$

with $\unicode[STIX]{x1D706}_{x}$ and $\unicode[STIX]{x1D706}_{xz}$ dimensionless slip lengths defined by

(6.20a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{x}=u_{\infty }^{\dagger (1)}-y_{\infty },\quad \unicode[STIX]{x1D706}_{xz}=w_{\infty }^{\dagger (1)}.\end{eqnarray}$$

Analogous reasoning leads easily to the transverse slip velocity at $Y=0$ , given by

(6.21) $$\begin{eqnarray}W_{s}=W|_{Y=0}=\unicode[STIX]{x1D716}\left(\unicode[STIX]{x1D706}_{zx}\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}Y}+\unicode[STIX]{x1D706}_{z}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Y}\right)\biggm\vert_{Y=0},\end{eqnarray}$$

with the slip lengths defined by

(6.22a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{zx}=u_{\infty }^{\dagger (3)},\quad \unicode[STIX]{x1D706}_{z}=w_{\infty }^{\dagger (3)}-y_{\infty },\end{eqnarray}$$

whereas $V$ vanishes at order $\unicode[STIX]{x1D716}$ in $Y=0$ (cf. (6.16)). In compact form we can write the in-plane components of the velocity at the fictitious wall in $Y=0$ as

(6.23) $$\begin{eqnarray}\left(\begin{array}{@{}c@{}}U_{s}\\ W_{s}\end{array}\right)=\unicode[STIX]{x1D716}\unicode[STIX]{x1D726}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}Y}\left(\begin{array}{@{}c@{}}U\\ W\end{array}\right)\biggm\vert_{Y=0},\end{eqnarray}$$

with

(6.24) $$\begin{eqnarray}\unicode[STIX]{x1D726}=\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D706}_{x} & \unicode[STIX]{x1D706}_{xz}\\ \unicode[STIX]{x1D706}_{zx} & \unicode[STIX]{x1D706}_{z}\end{array}\right),\end{eqnarray}$$

the slip tensor; $\unicode[STIX]{x1D726}$ is diagonal for symmetric wall patterns, in particular for isotropic or orthotropic surface roughnesses, in which case $\unicode[STIX]{x1D706}_{xz}=\unicode[STIX]{x1D706}_{zx}=0$ . This is the case of riblets, for example, and in this simpler limit our result coincides with that derived by other methods by Bechert & Bartenwerfer (Reference Bechert and Bartenwerfer1989) and Luchini, Manzo & Pozzi (Reference Luchini, Manzo and Pozzi1991).

Another interesting contribution on this problem is due to Achdou, Pironneau & Valentin (Reference Achdou, Pironneau and Valentin1998); these authors went up to second order in $\unicode[STIX]{x1D716}$ , still in the limit of negligible ${\mathcal{R}}$ , and focussed on two-dimensional problems. Thus, they identified the nonlinear correction to the first-order Navier condition, $U_{s}=\unicode[STIX]{x1D716}\unicode[STIX]{x1D706}_{x}\unicode[STIX]{x2202}U/\unicode[STIX]{x2202}Y|_{Y=0}$ , always maintaining the wall-normal velocity equal to zero, and analysed numerically a few steady, laminar cases, comparing model results to direct simulations which fully accounted for the surface topography. The second-order approximation, which required the solution of two additional microscopic closure problems in the unit cell, consistently provided an excellent match to the direct simulations. However, in all the cases tested, the first-order term gave already very satisfactory comparisons to the feature-resolving simulations, and this raises the question of whether there is any substantial advantage in deriving the second-order correction for (only) the longitudinal velocity.

6.1.2 An alternative derivation of the slip velocity

Since we are interested in expressing the slip velocity at some fictitious smooth surface as a function of an external, macroscopic shear, we can assume that the microscopic equations are driven by two shear stress components, along $x$ and $z$ , which will be called $\unicode[STIX]{x1D70F}_{x}$ and $\unicode[STIX]{x1D70F}_{z}$ , functions of $X$ and $Z$ ; $\unicode[STIX]{x1D70F}_{x}$ and $\unicode[STIX]{x1D70F}_{z}$ are positioned at some generic coordinate ${\mathcal{Y}}$ (above, or at the most coinciding with, the upper edge of the wall protrusions) via the use of the delta function. The direct equations (6.6) and (6.7) now become

(6.25) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{(0)}}{\unicode[STIX]{x2202}x_{i}}=0, & \displaystyle\end{eqnarray}$$

(6.26) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{(0)}}{\unicode[STIX]{x2202}x_{j}^{2}}+\unicode[STIX]{x1D6FF}(y-{\mathcal{Y}})S_{i}, & \displaystyle\end{eqnarray}$$

with $S_{i}$ the components of the volume source term $\boldsymbol{S}=(\unicode[STIX]{x1D70F}_{x},0,\unicode[STIX]{x1D70F}_{z})$ . The boundary conditions at $y\rightarrow \infty$ are $\unicode[STIX]{x2202}u^{(0)}/\unicode[STIX]{x2202}y=v^{(0)}=\unicode[STIX]{x2202}w^{(0)}/\unicode[STIX]{x2202}y=0$ ; the other conditions are, as before, no slip at $y=y_{wall}$ and periodicity on opposing vertical planes. The adjoint variables, i.e. the Green’s functions, are taken to satisfy the system

(6.27) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{i}}=0, & \displaystyle\end{eqnarray}$$

(6.28) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\unicode[STIX]{x2202}p^{\dagger (k)}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{j}^{2}}+\unicode[STIX]{x1D6FF}(y-{\mathcal{Y}})\unicode[STIX]{x1D6FF}_{ki}, & \displaystyle\end{eqnarray}$$

with the same boundary conditions as in the direct problem for variables with corresponding names. As before, the index $k$ can only take values 1 and 3. For example, when $k=1$ the last term on the right-hand side of (6.28) is the vector $[\unicode[STIX]{x1D6FF}(y-{\mathcal{Y}}),0,0]$ . The Lagrange–Green identity yields

(6.29) $$\begin{eqnarray}\lceil u_{k}^{(0)}(x,{\mathcal{Y}},z,X_{j})\rceil =\lceil u_{i}^{\dagger (k)}(x,{\mathcal{Y}},z)\rceil S_{i},\end{eqnarray}$$

from which the macroscopic slip velocity components when ${\mathcal{Y}}=0$ follow

(6.30a,b ) $$\begin{eqnarray}\displaystyle U_{s}=\unicode[STIX]{x1D716}(\unicode[STIX]{x1D706}_{x}\unicode[STIX]{x1D70F}_{x}+\unicode[STIX]{x1D706}_{xz}\unicode[STIX]{x1D70F}_{z}),\quad W_{s}=\unicode[STIX]{x1D716}\,(\unicode[STIX]{x1D706}_{zx}\unicode[STIX]{x1D70F}_{x}+\unicode[STIX]{x1D706}_{z}\unicode[STIX]{x1D70F}_{z}), & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D706}_{x}=\lceil u^{\dagger (1)}(x,0,z)\rceil$ , $\unicode[STIX]{x1D706}_{xz}=\lceil w^{\dagger (1)}(x,0,z)\rceil$ , $\unicode[STIX]{x1D706}_{zx}=\lceil u^{\dagger (3)}(x,0,z)\rceil$ and $\unicode[STIX]{x1D706}_{z}=\lceil w^{\dagger (3)}(x,0,z)\rceil$ . In a numerical application, provided the delta function is well approximated and discretized, the results for the components of the slip tensor $\unicode[STIX]{x1D726}$ when computed as described here coincide with those reported in (6.24); equation (6.30) coincides with (6.23), aside from formally higher-order terms related to $\unicode[STIX]{x2202}V/\unicode[STIX]{x2202}X|_{Y=0}$ and $\unicode[STIX]{x2202}V/\unicode[STIX]{x2202}Z|_{Y=0}$ .

6.1.3 Accounting for near-wall inertia

If ${\mathcal{R}}$ is order one, using the same approach just described the leading-order equations become

(6.31) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{(0)}}{\unicode[STIX]{x2202}x_{i}}=0, & \displaystyle\end{eqnarray}$$

(6.32) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{R}}\left[\frac{\unicode[STIX]{x2202}u_{i}^{(0)}}{\unicode[STIX]{x2202}t}+u_{j}^{(0)}\frac{\unicode[STIX]{x2202}u_{i}^{(0)}}{\unicode[STIX]{x2202}x_{j}}\right]=-\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{(0)}}{\unicode[STIX]{x2202}x_{j}^{2}}+\unicode[STIX]{x1D6FF}(y-{\mathcal{Y}})S_{i}, & \displaystyle\end{eqnarray}$$

subject to the same boundary conditions as in § 6.1.2. In this case, the volumetric source term $S_{i}$ can, in principle, depend also on time. We now multiply equations (6.31) and (6.32) by $p^{\dagger }$ and $u_{i}^{\dagger }$ , sum and integrate in time, from $t=0$ to some final time $t=T$ , and in space, over ${\mathcal{V}}_{f}$ . We impose for the backward-in-time auxiliary system to satisfy

(6.33) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{i}}=0, & \displaystyle\end{eqnarray}$$

(6.34) $$\begin{eqnarray}\displaystyle & \displaystyle -{\mathcal{R}}\left(\frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}t}+u_{j}^{(0)}\frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{j}}\right)=-\frac{\unicode[STIX]{x2202}p^{\dagger (k)}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{j}^{2}}, & \displaystyle\end{eqnarray}$$

with the same adjoint boundary conditions used in § 6.1.2, i.e.  $\unicode[STIX]{x2202}u^{\dagger (k)}/\unicode[STIX]{x2202}y=v^{\dagger (k)}=\unicode[STIX]{x2202}w^{\dagger (k)}/\unicode[STIX]{x2202}y=0$ when $y\rightarrow \infty$ and $u_{i}^{\dagger }=0$ at $y=y_{wall}$ . We are left with

(6.35) $$\begin{eqnarray}{\mathcal{R}}\int _{{\mathcal{V}}_{f}}u_{i}^{\dagger (k)}u_{i}\biggm\vert_{t=T}\,\text{d}{\mathcal{V}}={\mathcal{R}}\int _{{\mathcal{V}}_{f}}u_{i}^{\dagger (k)}u_{i}\biggm\vert_{t=0}\,\text{d}{\mathcal{V}}+\int _{0}^{T}\int _{{\mathcal{V}}_{f}}u_{i}^{\dagger (k)}\unicode[STIX]{x1D6FF}(y-{\mathcal{Y}})S_{i}\,\text{d}{\mathcal{V}}\,\text{d}t.\end{eqnarray}$$

Setting the terminal condition for the adjoint problem as

(6.36) $$\begin{eqnarray}u_{i}^{\dagger (k)}|_{t=T}=\unicode[STIX]{x1D6FF}_{ki}\unicode[STIX]{x1D6FF}(y-{\mathcal{Y}}),\end{eqnarray}$$

(6.35) gives

(6.37) $$\begin{eqnarray}\lceil u_{k}^{(0)}\rceil |_{t=T,y={\mathcal{Y}}}=\frac{1}{{\mathcal{R}}}\int _{0}^{T}\lceil u_{i}^{\dagger (k)}\rceil |_{y={\mathcal{Y}}}S_{i}\,\text{d}t+\frac{{\mathcal{V}}}{\unicode[STIX]{x1D6FA}}(u_{i}^{\dagger (k)},u_{i}^{(0)})|_{t=0}.\end{eqnarray}$$

In terms of outer variables, at the upper edge of the rough pattern (i.e. when ${\mathcal{Y}}=0$ ) we thus have

(6.38) $$\begin{eqnarray}\displaystyle & \displaystyle U_{s}(T,X,Z)=\frac{\unicode[STIX]{x1D716}}{{\mathcal{R}}}\int _{0}^{T}\unicode[STIX]{x1D706}_{1i}S_{i}\,\text{d}t+\unicode[STIX]{x1D716}\frac{{\mathcal{V}}}{\unicode[STIX]{x1D6FA}}(u_{i}^{\dagger (1)},u_{i}^{(0)})|_{t=0}, & \displaystyle\end{eqnarray}$$

(6.39) $$\begin{eqnarray}\displaystyle & \displaystyle W_{s}(T,X,Z)=\frac{\unicode[STIX]{x1D716}}{{\mathcal{R}}}\int _{0}^{T}\unicode[STIX]{x1D706}_{3i}S_{i}\,\text{d}t+\unicode[STIX]{x1D716}\frac{{\mathcal{V}}}{\unicode[STIX]{x1D6FA}}(u_{i}^{\dagger (3)},u_{i}^{(0)})|_{t=0}, & \displaystyle\end{eqnarray}$$

with the dynamic slip lengths given by

(6.40) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D706}_{x}:=\unicode[STIX]{x1D706}_{11}=\lceil u_{1}^{\dagger (1)}\rceil |_{y=0},\quad \unicode[STIX]{x1D706}_{xz}:=\unicode[STIX]{x1D706}_{13}=\lceil u_{3}^{\dagger (1)}\rceil |_{y=0},\\ \unicode[STIX]{x1D706}_{zx}:=\unicode[STIX]{x1D706}_{31}=\lceil u_{1}^{\dagger (3)}\rceil |_{y=0},\quad \text{and}\quad \unicode[STIX]{x1D706}_{z}:=\unicode[STIX]{x1D706}_{33}=\lceil u_{3}^{\dagger (3)}\rceil |_{y=0}.\end{array}\right\}\end{eqnarray}$$

We observe, once more, that nothing has been said yet about the vertical velocity condition at the surface in $y=0$ , aside from stating that, at leading order, $\lceil v^{(0)}\rceil =0$ .

6.1.4 Sample results for the slip lengths

An example of isosurfaces of $u^{\dagger (1)}$ and $w^{\dagger (1)}$ in a unit cell, for the case of a quasi-isotropic rough wall, is displayed in figure 15, for the case of a near-wall flow with negligible inertial effects. The results have been computed by the procedure described in § 6.1.1. Two microscopic vertical axes are defined; one is $y$ and, as usual, has its origin at the upper rim of the roughness elements. The other axis is ${\tilde{y}}$ and it runs through the rough surface, in such a way that $\int _{\unicode[STIX]{x1D6FA}}{\tilde{y}}_{wall}\,\text{d}x\,\text{d}z=0$ . The highest roughness peak is found at ${\tilde{y}}=0.1685$ and the lowest trough is at ${\tilde{y}}=-0.1144$ . For the specific roughness pattern considered here, we have $u_{\infty }^{\dagger (1)}=3.9482$ and $w_{\infty }^{\dagger (1)}=\unicode[STIX]{x1D706}_{xz}=5.60\times 10^{-4}$ at ${\tilde{y}}_{\infty }=4$ . The solution of the adjoint system, when setting $k=3$ in (6.12), is $u_{\infty }^{\dagger (3)}=\unicode[STIX]{x1D706}_{zx}=5.60\times 10^{-4}$ and $w_{\infty }^{\dagger (3)}=3.9459$ .

Figure 15.

(a) Rough wall used to illustrate the procedure. Isosurfaces of $u^{\dagger \,(1)}$ (b) and $w^{\dagger \,(1)}$ (c) in the unit cell, arising from the solution of the adjoint equations (6.9)–(6.10). The vertical axis ${\tilde{y}}$ in the figure has its origin on a mid-plane through the roughness.

By (6.20) and (6.22) it is immediately obtained that the diagonal components of $\unicode[STIX]{x1D726}$ are negative if measured with respect to ${\tilde{y}}$ , i.e. the two virtual origins, for the longitudinal and the transverse motion, sit above the ${\tilde{y}}=0$ axis, with reverse flow in the region underneath the virtual origin (cf. figure 16). This causes numerical problems because of loss of ellipticity (Achdou et al. Reference Achdou, Pironneau and Valentin1998). In an actual simulation the origin of the vertical coordinate should be shifted, to make it coincide, e.g. with the position of the roughness crests (in ${\tilde{y}}=0.1685$ ) as shown in the sketch of figure 16, thus using $y$ as wall-normal coordinate in place of ${\tilde{y}}$ . For the particular case of figure 15 the slip tensor becomes

(6.41) $$\begin{eqnarray}\unicode[STIX]{x1D726}=\left(\begin{array}{@{}cc@{}}0.1167 & 5.6\times 10^{-4}\\ 5.6\times 10^{-4} & 0.1144\end{array}\right),\end{eqnarray}$$

and the structure of the matrix indicates clearly the weakly anisotropic nature of the surface considered. When computing the entries of $\unicode[STIX]{x1D726}$ with the alternative approach described in § 6.1.2 the same results are obtained.

Figure 16.

Sketch of a generic near-wall region in a streamwise–wall-normal plane. In a macroscopic description the velocity decreases linearly towards a virtual origin, indicated in the image with a black dot. The control volume drawn within dashed lines is meant to show that, locally, a wall-normal velocity (in this case at $Y=0$ ) arises when there is a gradient (in $X$ or $Z$ ) of the velocity components parallel to the $Y=0$ plane.

6.2.1 Evaluating the transfer matrix $\unicode[STIX]{x1D63C}$

We must now give the coefficients of the matrix $\unicode[STIX]{x1D63C}$ . In the case of a steady, linear system (cf. (6.25)–(6.26)) the solution depends on the volume forcing $\boldsymbol{S}$ as

(6.53) $$\begin{eqnarray}u_{i}^{(0)}=u_{j}^{\dagger (i)}S_{j},\end{eqnarray}$$

(cf. (6.29)), so that the solution of the adjoint system (6.27)–(6.28) forced by a delta function positioned in $y=0$ yields directly the required transfer matrix. In other words, the components of $\unicode[STIX]{x1D63C}$ needed to evaluate the transpiration length, $m_{ik}$ , are

(6.54) $$\begin{eqnarray}a_{ij}=u_{j}^{\dagger (i)}.\end{eqnarray}$$

This procedure was applied by Lācis et al. (Reference Lācis, Sudhakar, Pasche and Bagheri2019) for the case of a turbulent channel flow with a regularly patterned wall. A texture of aligned rectangular cuboids was placed on the lower surface of the channel with $\unicode[STIX]{x1D716}=0.2$ and height $\hat{k}$ of the cuboids equal to $0.04L$ . For this configuration our own computations give $\unicode[STIX]{x1D706}_{x}=\unicode[STIX]{x1D706}_{z}=5.724\times 10^{-2}$ (and also $\unicode[STIX]{x1D706}_{xz}=\unicode[STIX]{x1D706}_{zx}=0$ , as expected because of the isotropy of the wall pattern), so that $b_{11}=b_{33}=17.47$ and $b_{13}=b_{31}=0$ . The value of the slip length coincides with that computed by Lācis et al. (Reference Lācis, Sudhakar, Pasche and Bagheri2019). Our calculations further yield $\int _{y_{wall}}^{0}\lceil a_{11}\rceil \,\text{d}y=\int _{y_{wall}}^{0}\lceil a_{33}\rceil \,\text{d}y=4.596\times 10^{-3}$ . Using definition (6.52), we finally obtain $m_{11}=m_{33}=8.029\times 10^{-2}$ , very close to the value quoted by Lācis et al. (Reference Lācis, Sudhakar, Pasche and Bagheri2019) for the same geometry.

6.2.2 Testing the rough-wall boundary conditions

For the simple geometry examined here, (6.51) becomes

(6.55) $$\begin{eqnarray}V|_{Y=0}=-\unicode[STIX]{x1D716}m_{11}\left[\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}X}+\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Z}\right]\biggm\vert_{Y=0},\end{eqnarray}$$

i.e. from continuity

(6.56) $$\begin{eqnarray}V|_{Y=0}=\unicode[STIX]{x1D716}m_{11}\frac{\unicode[STIX]{x2202}V}{\unicode[STIX]{x2202}Y}\biggm\vert_{Y=0}.\end{eqnarray}$$

The physical significance of the transpiration coefficient in this case is thus immediately evident: $m_{11}$ plays the same role as $\unicode[STIX]{x1D706}_{x}$ and $\unicode[STIX]{x1D706}_{z}$ , i.e. the vertical velocity component can be extrapolated to zero in $Y=-\unicode[STIX]{x1D716}m_{11}$ . The existence of a vertical Navier coefficient was previously hypothesized by Gómez-de-Segura et al. (Reference Gómez-de-Segura, Fairhall, MacDonald, Chung and García-Mayoral2018), but no strategy was put forward to assign its value.

The average velocity distribution computed by employing the equivalent wall conditions (6.23) and (6.56) agrees quite well with a geometry-resolving simulation, as shown first by Lācis et al. (Reference Lācis, Sudhakar, Pasche and Bagheri2019), whereas when fluid transpiration at the wall is neglected $|\unicode[STIX]{x0394}U^{+}|$ is underestimated. We have repeated the simulations by Lācis et al. (Reference Lācis, Sudhakar, Pasche and Bagheri2019) with the coefficients given in § 6.2.1, and the results, scaled in wall units using the friction velocity based on the total stress at $Y=0$ , are shown in figure 17. The mean velocity plotted on the vertical axis is relative to the slip velocity on the roughness’ crests, thus eliminating the need to define a zero velocity plane (Orlandi & Leonardi Reference Orlandi and Leonardi2006). The result embodied by this figure demonstrates the importance of adding a nominally smaller-order wall-normal velocity term to reproduce mean profiles over a rough surface. It has been shown by Lācis et al. (Reference Lācis, Sudhakar, Pasche and Bagheri2019) that also the r.m.s. values of the fluctuating velocity components are well captured by the approach which includes transpiration through the fictitious wall.

Figure 17.

Mean velocity profiles for a channel flow at $Re_{\unicode[STIX]{x1D70F}}=180$ . The top curve (empty circles) corresponds to a smooth wall. The dashed line represents the result of a direct simulation over a wall patterned with cuboids, modelled using Navier’s slip and $V|_{Y=0}=0$ . The solid line corresponds to results obtained using both slip and wall transpiration, for the same rough wall. This latter curve is rather close to the results of the geometry-resolving direct simulation by Lācis et al. (Reference Lācis, Sudhakar, Pasche and Bagheri2019), represented with filled square symbols.

The configuration tested corresponds to the lower transitionally rough regime (the roughness Reynolds number is $k^{+}=7.1$ and $\unicode[STIX]{x0394}U^{+}\approx -2.4$ ); as such, we are justified in employing slip and transpiration conditions with coefficients evaluated on the basis of microscopic Stokes problems. When advection through the roughness is non-negligible and ${\mathcal{R}}$ cannot be set to zero, an Oseen-like strategy, such as that described in § 4.2.4, might be attempted. This is presented in the next section, while in § 6.4 we will describe a few interesting cases in which microscopic wall features produce sizeable effects.

6.4.1 Riblets

The prototypical case of small-amplitude, anisotropic roughness is constituted of riblets (figure 3 b), minute wall grooves elongated along the mean flow direction. The state of the art in riblets research is comprehensively reviewed by Bannier (Reference Bannier2016), including a discussion on the effects of riblets size, the influence of pressure gradients, the laminar–turbulent transition, Reynolds and Mach number effects, three-dimensional riblets, etc. Here we limit ourselves to longitudinal riblets for which the slip tensor is simply

(6.64) $$\begin{eqnarray}\unicode[STIX]{x1D726}=\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D706}_{x} & 0\\ 0 & \unicode[STIX]{x1D706}_{z}\end{array}\right).\end{eqnarray}$$

Whereas riblets are ineffective in laminar flow as a drag reducing agent (Luchini Reference Luchini1995), they represent an interesting passive control approach in the turbulent case, provided they remain immersed within the viscous sublayer, for the flow to be in the lower transitionally rough regime. Riblets are one of the few technologies which have been tested with some degree of success in both laboratory experiments and in aero/hydrodynamic macroscale settings, without having transitioned yet to operational fleets of aircraft or ships for reasons of maintenance and effectiveness over time. Several types of protrusion shapes have been analysed, including triangular, trapezoidal, blade, U-shaped, sinusoidal and notched-peak configurations, reaching a maximum reduction in skin friction drag of $10\,\%$ (cf. figure 21). The best performances are achieved for spanwise periodicity $s^{+}$ close to fifteen wall units, which translates to a spacing of $30{-}70~\unicode[STIX]{x03BC}\text{m}$ for typical aeronautical configurations.

Figure 21.

Experimental points for skin friction drag reduction for triangular and blade riblets as a function of spanwise spacing, $s^{+}$ , scaled in viscous wall units; the dashed lines represent qualitatively the predictions of a linear Stokes model. Adapted from Bechert et al. (Reference Bechert, Bruse, Hage, Van Der Hoeven and Hoppe1997).

The agreed-upon argument is that riblets create a gap between the virtual origin of the mean, longitudinal flow and that of the transverse flow, as by the sketch in figure 22. The turbulent mean flow is primarily affected by the longitudinal slip length, $\unicode[STIX]{x1D706}_{x}$ , whereas turbulent fluctuations also see the transverse slip length, $\unicode[STIX]{x1D706}_{z}$ . Assuming that turbulent fluctuations see their effective wall only at $\unicode[STIX]{x1D706}_{z}$ , the difference in slip lengths becomes the only physically significant length scale parameter. Luchini et al. (Reference Luchini, Manzo and Pozzi1991) were the first to argue that the properties of ribbed walls must be a function of $\unicode[STIX]{x0394}\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}_{x}-\unicode[STIX]{x1D706}_{z}$ , and that variations in $\unicode[STIX]{x0394}\unicode[STIX]{x1D706}$ are responsible for a displacement of the logarithmic layer, yielding a reduction of the skin friction coefficient, $C_{f}=2\hat{\unicode[STIX]{x1D70F}}/\unicode[STIX]{x1D70C}\hat{U} _{avg}^{2}=2(u_{\ast }/\hat{U} _{avg})^{2}$ . Such a variation is approximated by

(6.65) $$\begin{eqnarray}\frac{\unicode[STIX]{x0394}C_{f}}{C_{f_{0}}}=-\frac{\unicode[STIX]{x0394}U^{+}}{(2C_{f_{0}})^{-1/2}+(2\unicode[STIX]{x1D705})^{-1}},\end{eqnarray}$$

with $\unicode[STIX]{x0394}C_{f}=C_{f}-C_{f_{0}}$ , and $C_{f_{0}}$ the reference coefficient for the flow over a smooth wall under the same external conditions (Luchini Reference Luchini and Désideri1996; García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2011). Equation (6.65) applies for small changes of $C_{f}$ , i.e. for small roughness functions $\unicode[STIX]{x0394}U^{+}$ , in the linear, dashed tracts of figure 21. Because of linearity it is

(6.66) $$\begin{eqnarray}\unicode[STIX]{x0394}U^{+}\simeq \unicode[STIX]{x1D707}_{0}\unicode[STIX]{x0394}\unicode[STIX]{x1D706}^{+},\end{eqnarray}$$

with $\unicode[STIX]{x0394}\unicode[STIX]{x1D706}^{+}$ scaled in wall units, and Luchini (Reference Luchini and Désideri1996) has shown that, to a good approximation, $\unicode[STIX]{x1D707}_{0}=1.$ In turn, $\unicode[STIX]{x0394}\unicode[STIX]{x1D706}^{+}$ is a function of $s^{+}$ , dimensionless periodicity of the riblets, as shown for several shapes by Luchini et al. (Reference Luchini, Manzo and Pozzi1991). The linear dependence of $\unicode[STIX]{x0394}C_{f}$ on $\unicode[STIX]{x0394}\unicode[STIX]{x1D706}^{+}$ (respectively, $s^{+}$ ) persists up to $\unicode[STIX]{x0394}\unicode[STIX]{x1D706}^{+}$ equal to approximately 2 (respectively, $s^{+}\simeq 10$ , cf. dashed lines in figure 21), i.e. as long as the texture length scales are small compared to the characteristic scales of near-wall turbulence. Beyond some thresholds on $\unicode[STIX]{x0394}\unicode[STIX]{x1D706}^{+}$ and $s^{+}$ the linear behaviour deteriorates. When the spacing of the riblets is above a value of the order of $20\div 30$ wall units, quasi-two dimensional spanwise rollers appear, related to a Kelvin–Helmholtz instability of the mean velocity profile (García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2011). This occurrence, as we will see in the following, is common to the case of the flow over porous and poroelastic coatings, and is associated with a drag increase.

Figure 22.

Sketch of streamwise (a) and spanwise (b) velocity profiles at the crest and trough of a ribbed wall. The virtual origins are indicated with open circles.

Numerical modelling of the flow over walls with riblets by the use of equivalent, homogenized conditions has consistently used $V=0$ at $Y=0$ . The possibility of a transpiration condition has been suggested, but never much elaborated upon. On the basis of the findings reported in § 6.2, we can now state that to second order in $\unicode[STIX]{x1D716}$ a suitable boundary condition is

(6.67) $$\begin{eqnarray}V|_{Y=0}=-\unicode[STIX]{x1D716}\left[m_{11}\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}X}+m_{33}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}Z}\right]\biggm\vert_{Y=0}.\end{eqnarray}$$

The coefficients are

(6.68a,b ) $$\begin{eqnarray}m_{11}=\left[\int _{y_{wall}}^{0}\lceil a_{11}\rceil \,\text{d}y\right]b_{11}\quad \text{and}\quad m_{33}=\left[\int _{y_{wall}}^{0}\lceil a_{33}\rceil \,\text{d}y\right]b_{33},\end{eqnarray}$$

with

(6.69) $$\begin{eqnarray}\unicode[STIX]{x1D63D}=\unicode[STIX]{x1D726}^{-1}=\left(\begin{array}{@{}cc@{}}1/\unicode[STIX]{x1D706}_{x} & 0\\ 0 & 1/\unicode[STIX]{x1D706}_{z}\end{array}\right).\end{eqnarray}$$

The two transpiration coefficients $m_{11}$ and $m_{33}$ are in general different, which means that the transpiration velocity cannot take the simple form (6.56); the presence of two additional length-related parameters modifies at higher order the picture which has emerged over the years on the role of the length scale $\unicode[STIX]{x0394}\unicode[STIX]{x1D706}$ .

Riblets can be optimized for given flow conditions (e.g. for cruise flight of an aircraft in still air), but their use can be detrimental in off-design conditions (for example, during take-off or landing, or in gusty winds). Three-dimensional variations of the conventional flow-parallel riblet shape have been proposed in the recent past, including streamwise-sinusoidal and zigzag riblets (Kramer et al. Reference Kramer, Grüneberger, Thiele, Wassen, Hage and Meyer2010; Grüneberger et al. Reference Grüneberger, Kramer, Wassen, Hage, Meyer, Thiele, Tropea and Bleckmann2012), and grooves of variable height (McClure, Smith & Baker Reference McClure, Smith and Baker2010), with limited success. In general, three-dimensional riblets do not outperform their two-dimensional, streamwise-invariant, counterpart. A possible, as yet little-explored alternative to potentially extend the range of flow parameters for which riblets, or any other kind of protrusions at the wall, exert a beneficial effect consists in rendering them capable to deform elastically, adapting to the flow. Trapezoidal riblets made of high-elongation elastomeric materials have been patented by Rawlings & Burg (Reference Rawlings and Burg2016), with the claim that their optimized structural design provides the capability for riblets to be ‘thinner, lower weight and more aerodynamically efficient’. To date, the only fluid–solid coupling analysis of riblets is due to Zampogna et al. (Reference Zampogna, Naqvi, Magnaudet and Bottaro2019c ); a similar upscaling strategy, within the adjoint framework, is described in the following section.

6.4.2 Deformable surface protrusions

The study of the fluid flow interacting with microscopic surface features made of a linearly elastic material and attached to a rigid substrate starts from the equations of motion for the fluid and Cauchy’s equations for the solid, i.e. in dimensional form

(6.70a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\hat{U} _{i}}{\unicode[STIX]{x2202}\hat{X}_{i}}=0,\quad \unicode[STIX]{x1D70C}_{f}\left(\frac{\unicode[STIX]{x2202}\hat{U} _{i}}{\unicode[STIX]{x2202}\hat{t}}+\hat{U} _{j}\frac{\unicode[STIX]{x2202}\hat{U} _{i}}{\unicode[STIX]{x2202}\hat{X}_{j}}\right)=-\frac{\unicode[STIX]{x2202}\hat{P}}{\unicode[STIX]{x2202}\hat{X}_{i}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\hat{X}_{j}}[2\unicode[STIX]{x1D707}\hat{\unicode[STIX]{x1D700}}_{ij}(\hat{\boldsymbol{U}})],\end{eqnarray}$$

(6.71) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{s}\frac{\unicode[STIX]{x2202}^{2}\hat{V}_{i}}{\unicode[STIX]{x2202}\hat{t}^{2}}=\frac{\unicode[STIX]{x2202}\hat{\unicode[STIX]{x1D70E}}_{ij}}{\unicode[STIX]{x2202}\hat{X}_{j}}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\hat{X}_{j}}[\hat{{\mathcal{C}}}_{ijkl}\hat{\unicode[STIX]{x1D700}}_{kl}(\hat{\boldsymbol{V}})],\end{eqnarray}$$

with the operator $\hat{\unicode[STIX]{x1D700}}_{kl}$ defined as in (5.5) and applicable to both the fluid velocity, $\hat{\boldsymbol{U}}$ , and the solid deformation, $\hat{\boldsymbol{V}}$ . The fluid and solid equations are coupled through the matching of velocities and tractions across the microscopic fluid–solid boundary in ${\hat{Y}}={\hat{Y}}_{wall}(\hat{X},\hat{Z},\hat{t})$ , viz.

(6.72) $$\begin{eqnarray}\hat{U} _{i}=\frac{\unicode[STIX]{x2202}\hat{V}_{i}}{\unicode[STIX]{x2202}\hat{t}},\end{eqnarray}$$

and

(6.73) $$\begin{eqnarray}-\hat{P}n_{i}+2\unicode[STIX]{x1D707}\hat{\unicode[STIX]{x1D700}}_{ij}(\hat{\boldsymbol{U}})n_{j}=\hat{\unicode[STIX]{x1D70E}}_{ij}n_{j}.\end{eqnarray}$$

At the lower boundary of the unit cell ( ${\hat{Y}}\rightarrow -\infty$ ) we simply set $\hat{V}_{i}=0$ , under the assumption that the thick elastic coating is anchored onto a rigid substrate; at the upper boundary ( ${\hat{Y}}\rightarrow \infty$ ) we match with the outer shear flow, which requires setting the velocity and length scales in the region just above the upper boundary. Before doing this, it is however important to establish scaling relations within the unit cell, to normalize the microscale equations. Upon assuming that the continuum coating, made up by fluid and solid, is characterized by a frequency, $\unicode[STIX]{x1D627}$ , sufficiently large for dynamic effects to be felt at leading order, we can write that in the fluid domain

(6.74) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{f}{\mathcal{U}}\unicode[STIX]{x1D627}\sim \frac{\unicode[STIX]{x0394}P}{l}\sim \unicode[STIX]{x1D707}\frac{{\mathcal{U}}}{l^{2}},\end{eqnarray}$$

with ${\mathcal{U}}$ the velocity scale, $\unicode[STIX]{x0394}P$ the pressure scale and $l$ the microscopic length scale. From the above, we can choose the velocity scale to be

(6.75) $$\begin{eqnarray}{\mathcal{U}}=\frac{\unicode[STIX]{x0394}Pl}{\unicode[STIX]{x1D707}}.\end{eqnarray}$$

We also have a relation between the microscale $l$ and the frequency $\unicode[STIX]{x1D627}$ , which states that, for viscous effects to balance inertia, $l$ must be of the order of the Stokes layer thickness, i.e.

(6.76) $$\begin{eqnarray}l=\sqrt{\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D627}}}.\end{eqnarray}$$

The small displacement of the surface micropattern is assumed to occur coherently over a macroscopic length $L$ , as in the case of honami waves of canopies under the effect of wind. By balancing inertia and diffusion in Cauchy’s equation for the solid, we have

(6.77) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D627}^{2}\sim \frac{E}{L^{2}},\end{eqnarray}$$

so that the macroscale $L$ can be taken to coincide with the elastic wavelength (cf. (5.9)), i.e.

(6.78) $$\begin{eqnarray}L=\frac{1}{\unicode[STIX]{x1D627}}\sqrt{\frac{E}{\unicode[STIX]{x1D70C}_{s}}}.\end{eqnarray}$$

The kinematic condition (6.72) at the interface is useful since it permits us to relate the solid displacement to the fluid velocity, so that the length scale of the deformation, $\hat{\boldsymbol{V}}$ , is chosen equal to

(6.79) $$\begin{eqnarray}\frac{{\mathcal{U}}}{\unicode[STIX]{x1D627}}=\frac{\unicode[STIX]{x0394}Pl}{\unicode[STIX]{x1D707}\unicode[STIX]{x1D627}}.\end{eqnarray}$$

The dimensionless variables are thus written (and expanded) as follows:

(6.80) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle t=\unicode[STIX]{x1D627}\hat{t},\quad \boldsymbol{x}=\frac{\hat{\boldsymbol{X}}}{l},\quad p=p^{(0)}+\unicode[STIX]{x1D716}p^{(1)}+\cdots =\frac{\hat{P}}{\unicode[STIX]{x0394}P},\\ \displaystyle u_{i}=u_{i}^{(0)}+\unicode[STIX]{x1D716}u_{i}^{(1)}+\cdots =\frac{\unicode[STIX]{x1D707}\hat{U} _{i}}{\unicode[STIX]{x0394}Pl},\quad v_{i}=v_{i}^{(0)}+\unicode[STIX]{x1D716}v_{i}^{(1)}+\cdots =\frac{\unicode[STIX]{x1D707}\hat{V}_{i}\unicode[STIX]{x1D627}}{\unicode[STIX]{x0394}Pl}.\end{array}\right\}\end{eqnarray}$$

Substituting (6.80) into the equations for fluid and solid phases, using $X_{i}=\unicode[STIX]{x1D716}x_{i}$ and (4.4), assuming that the microscale Reynolds number is small, ${\mathcal{R}}=(\unicode[STIX]{x1D70C}_{f}{\mathcal{U}}l/\unicode[STIX]{x1D707})=O(\unicode[STIX]{x1D716})$ , the governing equations at different orders become

$O(\unicode[STIX]{x1D716}^{-2})$ :

(6.81) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}^{(-2)}}{\unicode[STIX]{x2202}x_{j}}=0,\end{eqnarray}$$

$O(\unicode[STIX]{x1D716}^{-1})$ :

(6.82) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}^{(-1)}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}^{(-2)}}{\unicode[STIX]{x2202}X_{j}}=0,\end{eqnarray}$$

$O(\unicode[STIX]{x1D716}^{0})$ :

(6.83) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}^{2}v_{i}^{(0)}}{\unicode[STIX]{x2202}t^{2}}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}^{(0)}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}^{(-1)}}{\unicode[STIX]{x2202}X_{j}}, & \displaystyle\end{eqnarray}$$

(6.84) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{(0)}}{\unicode[STIX]{x2202}x_{i}}=0, & \displaystyle\end{eqnarray}$$

(6.85) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{(0)}}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{(0)}}{\unicode[STIX]{x2202}x_{j}^{2}}, & \displaystyle\end{eqnarray}$$

with the following kinematic and dynamic conditions at the interface in $y=y_{wall}(x,z,t)$ :

(6.86) $$\begin{eqnarray}\displaystyle & \displaystyle u_{i}^{(0)}=\frac{\unicode[STIX]{x2202}v_{i}^{(0)}}{\unicode[STIX]{x2202}t}, & \displaystyle\end{eqnarray}$$

(6.87) $$\begin{eqnarray}\displaystyle & \displaystyle u_{i}^{(1)}=\frac{\unicode[STIX]{x2202}v_{i}^{(1)}}{\unicode[STIX]{x2202}t}, & \displaystyle\end{eqnarray}$$

(6.88) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{ij}^{(-2)}n_{j}=0, & \displaystyle\end{eqnarray}$$

(6.89) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{ij}^{(-1)}n_{j}=0, & \displaystyle\end{eqnarray}$$

(6.90) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D70C}_{s}}{\unicode[STIX]{x1D70C}_{f}}\unicode[STIX]{x1D70E}_{ij}^{(0)}n_{j}=-p^{(0)}n_{i}+2\unicode[STIX]{x1D700}_{ij}(\boldsymbol{u}^{(0)})n_{j}. & \displaystyle\end{eqnarray}$$

As before, the unit normal vector $\boldsymbol{n}$ points into the fluid; also, in the equations above we have used the following notations

(6.91) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{ij}^{(-2)}=C_{ijkl}\unicode[STIX]{x1D700}_{kl}(\boldsymbol{v}^{(0)}), & \displaystyle\end{eqnarray}$$

(6.92) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{ij}^{(-1)}=C_{ijkl}[\unicode[STIX]{x1D700}_{kl}(\boldsymbol{v}^{(1)})+{\mathcal{E}}_{kl}(\boldsymbol{v}^{(0)})], & \displaystyle\end{eqnarray}$$

(6.93) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{ij}^{(0)}=C_{ijkl}[\unicode[STIX]{x1D700}_{kl}(\boldsymbol{v}^{(2)})+{\mathcal{E}}_{kl}(\boldsymbol{v}^{(1)})], & \displaystyle\end{eqnarray}$$

with the microscopic and macroscopic dimensionless strain operators, $\unicode[STIX]{x1D700}_{kl}(\cdot )$ and ${\mathcal{E}}_{kl}(\cdot )$ , as in (5.12) and (5.20). We immediately notice that $v_{i}^{(0)}=v_{i}^{(0)}(t,X_{j})$ and $\unicode[STIX]{x1D70E}_{ij}^{(-2)}=0$ on account of (6.81) and (6.88). This allows us to express the fluid mass conservation and momentum equations at leading order in the unit cell as a function of the fluid velocity relative to the solid skeleton, as

(6.94a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{i}}(u_{i}^{(0)}-\dot{v}_{i}^{(0)})=0,\quad \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(u_{i}^{(0)}-\dot{v}_{i}^{(0)})=-\ddot{v}_{i}^{(0)}-\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x_{j}^{2}}(u_{i}^{(0)}-\dot{v}_{i}^{(0)}).\end{eqnarray}$$

The leading-order equation for the solid is $\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}^{(-1)}/\unicode[STIX]{x2202}x_{j}=0$ .

Figure 23.

Sketch of representative volume element of volume ${\mathcal{V}}$ ( ${\mathcal{V}}={\mathcal{V}}_{s}+{\mathcal{V}}_{f}$ ) drawn within dashed lines and projected onto the $(x,y)$ plane, for the flow past a rough, linearly elastic and impermeable surface. Periodicity of all variables holds on opposing planes parallel to the $y$ -axis.

At this point we can apply the procedure described in § 6.1.3, assuming that the fluid momentum equation is driven by a stress term, $\boldsymbol{S}=(\unicode[STIX]{x1D70F}_{x},0,\unicode[STIX]{x1D70F}_{z})$ , which only depends on macroscopic independent variables and time. This forcing is applied at a generic position $y={\mathcal{Y}}$ in ${\mathcal{V}}_{f}$ by the use of the delta function. The boundary conditions for $y\rightarrow \pm \infty$ are shown in figure 23, together with a sketch of the RVE. Three test functions, $p^{\dagger }$ , $u_{i}^{\dagger }$ and $v_{i}^{\dagger }$ , periodic along the horizontal directions, are introduced, and adjoint problems are built by employing the inner product defined in (4.12). However, since the unit cell includes portions filled with either fluid or solid, the spatial domain of integration in the inner product is either ${\mathcal{V}}_{f}$ or ${\mathcal{V}}_{s}$ , depending on whether we consider fluid-based or solid-based variables (cf. the analysis in § 5). Integrating also in time we have

(6.95) $$\begin{eqnarray}\displaystyle 0 & = & \displaystyle \int _{0}^{T}\left(p^{\dagger },\frac{\unicode[STIX]{x2202}(u_{i}^{(0)}-\dot{v}_{i}^{(0)})}{\unicode[STIX]{x2202}x_{i}}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\left(u_{i}^{\dagger },-\frac{\unicode[STIX]{x2202}(u_{i}^{(0)}-\dot{v}_{i}^{(0)})}{\unicode[STIX]{x2202}t}-\ddot{v}_{i}^{(0)}-\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}(u_{i}^{(0)}-\dot{v}_{i}^{(0)})}{\unicode[STIX]{x2202}x_{j}^{2}}+S_{i}\unicode[STIX]{x1D6FF}(y-{\mathcal{Y}})\right)\nonumber\\ \displaystyle & & \displaystyle +\,\left(v_{i}^{\dagger },\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\{{\mathcal{C}}_{ijkl}[\unicode[STIX]{x1D700}_{kl}(\boldsymbol{v}^{(1)})+{\mathcal{E}}_{kl}(\boldsymbol{v}^{(0)})]\}\right)\,\text{d}t.\end{eqnarray}$$

The auxiliary system in the RVE is

(6.96a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{i}}=0,\quad -\frac{\unicode[STIX]{x2202}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x2202}p^{\dagger (k)}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}u_{i}^{\dagger (k)}}{\unicode[STIX]{x2202}x_{j}^{2}},\end{eqnarray}$$

(6.97) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}[{\mathcal{C}}_{ijkl}\unicode[STIX]{x1D700}_{kl}(\boldsymbol{v}^{\dagger (m)})]=\unicode[STIX]{x1D6FF}_{im}.\end{eqnarray}$$

The superscript $(k)$ in system (6.96) reflects the fact that two problems must be solved in ${\mathcal{V}}_{f}$ , one for $k=1$ and the second for $k=3$ . The superscript $(m)$ in (6.97) indicates that three problems must be solved in ${\mathcal{V}}_{s}$ , for $m=1,2$ and 3.

Integrating by parts the right-hand side of (6.95), imposing that

(6.98a,b ) $$\begin{eqnarray}u_{i}^{\dagger (k)}=0;\quad {\mathcal{C}}_{ijkl}\unicode[STIX]{x1D700}_{kl}(\boldsymbol{v}^{\dagger (m)})n_{j}=0,\quad \text{on }y=y_{wall},\end{eqnarray}$$

(6.99) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}u_{1}^{\dagger (k)}}{\unicode[STIX]{x2202}y}=u_{2}^{\dagger (k)}=\frac{\unicode[STIX]{x2202}u_{3}^{\dagger (k)}}{\unicode[STIX]{x2202}y}=0,\quad \text{on }y=y_{\infty },\end{eqnarray}$$

and

(6.100) $$\begin{eqnarray}\boldsymbol{v}^{\dagger (m)}=0,\quad \text{on }y=-y_{\infty },\end{eqnarray}$$

we obtain, using the surface average notation $\lceil \cdot \rceil$ defined in (6.15) and imposing (6.36) at $t=T$ ,

(6.101) $$\begin{eqnarray}\displaystyle & & \displaystyle \lceil u_{k}^{(0)}\rceil |_{t=T,y={\mathcal{Y}}}-{\dot{v}_{k}}^{(0)}|_{t=T}=\int _{0}^{T}\lceil u_{i}^{\dagger (k)}\rceil |_{y={\mathcal{Y}}}S_{i}\,\text{d}t\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{{\mathcal{V}}}{\unicode[STIX]{x1D6FA}}\left[(u_{i}^{\dagger (k)},u_{i}^{(0)}-\dot{v}_{i}^{(0)})\biggm\vert_{t=0}\;-\;\int _{0}^{T}(u_{i}^{\dagger (k)},\ddot{v}_{i}^{(0)})\,\text{d}t\right],\qquad\end{eqnarray}$$

together with

(6.102) $$\begin{eqnarray}v_{m}^{(1)}={\mathcal{C}}_{ijkl}{\mathcal{E}}_{kl}(\boldsymbol{v}^{(0)})\frac{\unicode[STIX]{x2202}v_{i}^{\dagger (m)}}{\unicode[STIX]{x2202}x_{j}}.\end{eqnarray}$$

The two time integrals on the right-hand side of (6.101) are convolutions because of the temporal behaviour of the auxiliary, microscopic problem in ${\mathcal{V}}_{f}$ . Equation (6.102) furnishes the order-one solid displacement as a function of the gradient of a microscopic tensorial coefficient of rank 2 (which we can solve for) and the macroscopic solid strain at order zero. Knowledge of $\boldsymbol{v}^{(1)}$ is necessary to address the solution of Cauchy’s equations at leading order, which we rewrite here for clarity

(6.103) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\{{\mathcal{C}}_{ijkl}[\unicode[STIX]{x1D700}_{kl}(\boldsymbol{v}^{(1)})+{\mathcal{E}}_{kl}(\boldsymbol{v}^{(0)})]\}=0.\end{eqnarray}$$

We can now move to macroscopic variables; the two instantaneous (at $t=T$ ) slip velocity components become

(6.104) $$\begin{eqnarray}\displaystyle U_{s} & := & \displaystyle (U_{1})_{s}=\unicode[STIX]{x1D716}\left\{{\dot{v}_{1}}^{(0)}+\int _{0}^{T}\unicode[STIX]{x1D706}_{x}\unicode[STIX]{x1D70F}_{x}+\unicode[STIX]{x1D706}_{xz}\unicode[STIX]{x1D70F}_{z}\,\text{d}t\right\}\nonumber\\ \displaystyle & & \displaystyle +\,2\unicode[STIX]{x1D716}y_{\infty }\left[(u_{i}^{\dagger (1)},u_{i}^{(0)}-\dot{v}_{i}^{(0)})|_{t=0}-\int _{0}^{T}(u_{i}^{\dagger (1)},\ddot{v}_{i}^{(0)})\,\text{d}t\right],\end{eqnarray}$$

and

(6.105) $$\begin{eqnarray}\displaystyle W_{s} & := & \displaystyle (U_{3})_{s}=\unicode[STIX]{x1D716}\left\{{\dot{v}_{3}}^{(0)}+\int _{0}^{T}\unicode[STIX]{x1D706}_{zx}\unicode[STIX]{x1D70F}_{x}+\unicode[STIX]{x1D706}_{z}\unicode[STIX]{x1D70F}_{z}\,\text{d}t\right\}\nonumber\\ \displaystyle & & \displaystyle +\,2\unicode[STIX]{x1D716}y_{\infty }\left[(u_{i}^{\dagger (3)},u_{i}^{(0)}-\dot{v}_{i}^{(0)})|_{t=0}-\int _{0}^{T}(u_{i}^{\dagger (3)},\ddot{v}_{i}^{(0)})\,\text{d}t\right],\end{eqnarray}$$

with the dynamic slip coefficients formally identical to those given in (6.40).

In the absence of transpiration at the wall, when horizontal gradients of the slip velocity components can be neglected, the remaining kinematic boundary condition for the direct problem is simply $\lceil u_{2}^{(0)}\rceil =\dot{v}_{2}^{(0)}$ , which in macroscopic variables reads

(6.106) $$\begin{eqnarray}V|_{Y=0}:=U_{2}|_{Y=0}=\unicode[STIX]{x1D716}\dot{v}_{2}^{(0)}.\end{eqnarray}$$

Cell-averaged balance equations for the fluid/solid composite are still necessary to close the problem for the unknowns $\lceil u_{i}^{(0)}\rceil$ and $v_{i}^{(0)}$ , in a manner similar to the procedure in §§ 5.1, 5.2 and 5.3 for the case of a poroelastic medium away from boundaries with a free-fluid domain. We refer to Zampogna et al. (Reference Zampogna, Naqvi, Magnaudet and Bottaro2019c ) for the derivation of macroscopic equilibrium laws for the mixture in the unit cell.

6.4.3 Superhydrophobic and lubricant-infused coatings

Another interesting problem which requires a close microscopic look is that of a surface able to capture and maintain a layer of lubricant fluid, interposed between the solid matrix and the working fluid. If the working fluid is a liquid, the most efficient such lubricant is a gas or the vapour; in this case we speak of hydrophobic or superhydrophobic materials. If the lubricant fluid is a liquid, we speak of SLIPS (slippery liquid-infused porous surfaces), according to the acronym coined by Joanna Aizenberg, or of LIS (lubricant-impregnated surfaces) as suggested by Kripa K. Varanasi.

Figure 24.

Drag reduction in a channel with superhydrophobic walls for $Re_{\unicode[STIX]{x1D70F}}=180$ . The top right image illustrates the microscopic surface pattern, and the light-coloured regions refer to the gas pockets. The solid line in the $\unicode[STIX]{x0394}C_{f}/C_{f_{0}}$ versus $s^{+}$ plot comes from numerical results obtained using the Navier slip condition, while filled dots are used to indicate direct simulation results found when modelling the walls with alternating patches of no shear and no slip. Finally, the red dashed line corresponds to the analytical approximation (6.65)–(6.66). Adapted from Luchini (Reference Luchini, Choi and García-Mayoral2015).

Much work has been conducted on the flow over hydrophobic and superhydrophobic surfaces and several reviews exist (Quéré Reference Quéré2008; Rothstein Reference Rothstein2010). Such surfaces achieve their water repelling properties by combining the wall micro-pattern with a low surface energy coating. The large majority of the numerical models dedicated to studying the flow of a liquid over a superhydrophobic wall have considered the gas plastron as flat and undeformable (i.e. the surface tension is assumed to be very large), with the liquid able to freely slip over it (i.e. the dynamic viscosity of the gas is negligible compared to that of the liquid). The first computation of slip lengths under these assumptions, for the case of streamwise elongated microgrooves, was conducted by Philip (Reference Philip1972). Since then many other surface textures have been considered, including circular and square posts, cavities and grid-like patterns. A result common to all the surface textures examined is that the slip length is a decreasing function of the solid area fraction of the surface. This is an immediate consequence of the no-shear condition imposed at the water/gas interface, and means that larger gas pockets can produce larger drag reduction. Significant numerical results for flat, shear-free interfaces have been reported by Luchini (Reference Luchini, Choi and García-Mayoral2015). For the case of grooves aligned along the mean flow, as sketched in figure 24, Luchini modelled superhydrophobic walls both as a flat wall with alternating patches of no shear and no slip, and by employing the Navier condition (6.23) with the analytically derived slip lengths. In both cases the transpiration velocity was set to zero. Figure 24 shows that the direct numerical simulation results with the two approaches almost coincide up to $s^{+}$ equal to 30, demonstrating the accuracy which can be attained by the simple homogenized condition when the periodicity of the wall pattern is not too large. In any case, a large pattern periodicity exposes a real superhydrophobic surface to deformation of the interface and to depletion of the air layer, so that the results obtained for large $s^{+}$ are mostly of academic interest. Another significant result obtained by Luchini (Reference Luchini, Choi and García-Mayoral2015) is the accuracy of the analytical approximation embodied by equations (6.65) and (6.66), reported in the figure as a dashed line. Such an agreement extends beyond that found for riblets, cf. figure 21. The effect of a deforming air–water interface has been addressed in the steady, Stokes limit by Alinovi & Bottaro (Reference Alinovi and Bottaro2018) (see also references therein), the main conclusion being that $\unicode[STIX]{x0394}\unicode[STIX]{x1D706}$ increases when the interface protrudes outside of the representative microcavity. The interaction between turbulence and the gas pockets has been studied with direct simulations by Seo, García-Mayoral & Mani (Reference Seo, García-Mayoral and Mani2015); these authors looked at posts and grooves at the wall and treated the deforming interface by solving the Young–Laplace equation, after the pressure distribution had been determined. Later, Seo, García-Mayoral & Mani (Reference Seo, García-Mayoral and Mani2018) extended their numerical study and proposed a threshold criterion for the failure of superhydrophobic surfaces.

Figure 25.

Artist’s impression of the fabrication process of SLIPS/LIS. When a porous/textured layer (a) is impregnated by a low surface energy, chemically inert lubricant (b), a slippery surface is created. When the SLIPS/LIS is mildly inclined (c), a drop of water rolls easily down the surface (Smith et al. Reference Smith, Dhiman, Anand, Reza-Garduno, Cohen, McKinley and Varanasi2013). Image courtesy of Joanna Aizenberg, Harvard University.

The fragility of the gas plastron is a primary concern; it can collapse because of hydrostatic pressure (when the superhydrophobic surface is immersed in water) or can be carried away by the action of shear forces (in this sense, longitudinal grooves are particularly prone to destabilization). Some techniques have been proposed to maintain the gas layer, including entraining outside air or using electrolysis to replenish the plastron, but none has made it into actual underwater applications yet. Probably, the most effective alternative consists in using some liquid lubricant instead of air or vapour, employing an oil with some affinity to the porous, or textured, solid. It is the case of SLIPS/LIS (cf. the images in figure 25). Fabricating textures for creating SLIPS/LIS requires the same steps as fabricating superhydrophobic textures: the solid must be characterized by a low surface energy and, for the lubricant impregnated surface to be stable, also the surface energy of the lubricant liquid needs to be low.

Multiple applications can be envisioned for SLIPS/LIS. These include promoting dropwise condensation, reducing ice adhesion, enhancing optical transmission, anti-fouling, friction reduction, etc. (Wong et al. Reference Wong, Kang, Tang, Smythe, Hatton, Grinthal and Aizenberg2011; Solomon et al. Reference Solomon, Subramanyam, Farnham, Khalil, Anand, Varanasi, Ras and Marmur2017). As far as drag reduction is concerned, the drainage of the infused fluid by the shearing action of the working fluid is a potential drawback. Liu et al. (Reference Liu, Wexler, Schönecker and Stone2016) have shown that, for a given external fluid, a lower lubricant viscosity allows more lubricant to be retained. In the case of recirculating motion (such as in a Taylor–Couette cell, for example), with azimuthally periodic grooves, the problem of shear-induced drainage is absent; for this configuration Van Buren & Smits (Reference Van Buren and Smits2017) carried out torque measurements to find drag reductions of up to 35 % in the turbulent regime. The level of drag reduction for varying parameters was roughly comparable to that obtained in the superhydrophobic case, and appeared to peak for a pattern wavelength $s^{+}$ approximately equal to 70. Above such a value, drag was seen to increase again and this appeared to be related to the partial water wetting of the textured surface.

The reason why SLIPS/LIS reduce skin friction drag, by a similar mechanism and by a comparable amount than superhydrophobic surfaces, is surprising when one considers that the viscosity ratio between water and a lubricating liquid is $O(1)$ , while it is approximately 50 when the lubricant is air. This issue has been addressed recently by Arenas et al. (Reference Arenas, García, Fu, Orlandi, Hultmark and Leonardi2019) through direct numerical simulations for a few texture types, in the approximation of flat, undeformable interfaces. The water/lubricant interface was shown to damp wall-normal velocity fluctuations, limiting the flux of momentum inside the micro-cavities and consequently reducing the Reynolds stress and the dissipation within. In hindsight, the role of the vertical velocity at the crests of the roughness, and the fact that its absence has the effect of reducing drag, is also evident from figure 17, upon comparing the two lower curves (dashed line versus solid line).

7.1.1 Testing the slip condition

To verify the appropriateness of (7.13) we consider the simple case of the two-dimensional incompressible Stokes flow in a plane channel bounded from above (in $Y=2$ ) by an impermeable wall and from below ( $Y=0$ ) by the isotropic porous medium characterized in figure 26. The $Y=0$ axis is tangent to the uppermost inclusion (see sketch of the porous medium in figure 28 a), and the medium extends up to $Y=-0.5$ . The small parameter $\unicode[STIX]{x1D716}$ which defines the size of the pores is chosen equal to $0.1$ and a direct simulation which captures all details of the flow through the pores for $Y<0$ is conducted as reference. The simulation considers a minimal domain in $X$ , of length equal to $\unicode[STIX]{x1D716}$ on account of the fully developed nature of the motion. The scales used to normalize the continuity and momentum equations are $L$ for length, the bulk speed $U_{avg}$ for velocity and $\unicode[STIX]{x1D70C}U_{avg}^{2}$ for pressure, so that the Reynolds number is defined by $\text{Re}=\unicode[STIX]{x1D70C}U_{avg}L/\unicode[STIX]{x1D707}$ . The motion is driven by a constant pressure gradient $\text{d}p^{(0)}/\text{d}X$ chosen so that $\text{Re}=100$ . With the given normalization, the slip velocity takes the dimensionless form

(7.14) $$\begin{eqnarray}U_{s}=\unicode[STIX]{x1D716}\unicode[STIX]{x1D706}_{x}\frac{\text{d}U}{\text{d}Y}\biggm\vert_{Y\rightarrow 0^{+}}\;-\;\unicode[STIX]{x1D716}^{2}{\mathcal{K}}_{x}^{itf}\text{Re}\frac{\text{d}p^{(0)}}{\text{d}X},\end{eqnarray}$$

Figure 28.

Continuous lines: streamwise velocity distribution in a channel bounded from below by a porous domain, accounting for the motion through the pores when $Y<0$ . (a,c) Show close-ups around the Darcy–Stokes dividing line. Empty circles in all frames represent either the exact solution of (7.15) or the solution of Darcy equation. The dimensionless horizontal filtration velocity is equal to $4.40\times 10^{-4}$ .

so that the laminar flow solution in the fluid domain reads

(7.15) $$\begin{eqnarray}U(Y)=\text{Re}\left\{\frac{\text{d}p^{(0)}}{\text{d}X}\left[\frac{Y^{2}}{2}-\unicode[STIX]{x1D716}^{2}{\mathcal{K}}_{x}^{itf}\right]+A[Y+\unicode[STIX]{x1D716}\unicode[STIX]{x1D706}_{x}]\right\},\end{eqnarray}$$

with the constant $A$ given by

(7.16) $$\begin{eqnarray}A=-\frac{\text{d}p^{(0)}}{\text{d}X}\left[\frac{2-\unicode[STIX]{x1D716}^{2}{\mathcal{K}}_{x}^{itf}}{2+\unicode[STIX]{x1D716}\unicode[STIX]{x1D706}_{x}}\right].\end{eqnarray}$$

The solution $U$ arising from the direct simulation is displayed in figure 28 as a function of $Y$ , together with the distribution (7.15). The agreement with the exact solution is excellent in the free-fluid region $(0<Y<2)$ and also the slip velocity evaluated at $Y=0$ agrees to within $10^{-5}$ between the two cases, in particular $U_{s}=1.912\times 10^{-2}$ in the simulation and $U_{s}=1.913\times 10^{-2}$ from equation (7.14). As expected the contribution to the slip velocity from the Navier term is dominant, $\unicode[STIX]{x1D716}\unicode[STIX]{x1D706}_{x}\text{d}U/\text{d}Y|_{Y=0}=1.893\times 10^{-2}$ , and this term alone is sufficient to produce a very satisfactory result, as argued by Saffman (Reference Saffman1971). (The error committed in neglecting the term proportional to the mean pressure gradient is equal to $1\,\%$ in the simple case presented here to illustrate the technique, but it is expected to change upon varying the ratio between the pore size and the macroscopic length scale, the porosity and the permeability of the coating and, most of all, the value of the macroscopic Reynolds number.) Darcy’s solution within the porous domain (figure 28 a) matches well the exact numerical solution, which oscillates in $Y$ because of the alternation of pores and inclusions, except in two layers of thickness $O(\unicode[STIX]{x1D716})$ , one right below the dividing line (at $Y=0$ ) and the other right above the lower boundary of the porous medium (at $Y=-0.5$ ).

Whereas the theory just outlined provides satisfactory results, this is also due to the simplicity of the test case considered. More strenuous tests of the theory should include rough, permeable interfaces, possibly including fluid transpiration across, in the manner described in § 6.2, three-dimensional pressure gradients, anisotropic porous media and inertial effects.

7.1.2 The pressure jump condition

In non-idealized situations a pressure condition at the dividing surface is needed whenever Darcy’s equation must be solved in order to describe the flow in the porous medium. An interesting effective condition has been recently proposed by Lācis et al. (Reference Lācis, Sudhakar, Pasche and Bagheri2019), briefly summarized below.

The first step consists in writing a dimensionless, microscopic problem in the unit cell spanning across the interface; such a problem is the same as (7.5)–(7.6), with the addition of a shear term on the right-hand side of the momentum equation equal to $\unicode[STIX]{x1D6FF}(y)S_{i}$ (cf. § 6.1.2). Beyond periodicity along $x$ and $z$ , the other boundary conditions specify zero shear stress at $y=+y_{\infty }$ and $(u_{i},p)|_{y=-y_{\infty }}=(u_{i},p)|_{y=-y_{\infty }+1}$ (the latter on account of the periodicity of the fields on the boundaries of a unit cell, deep within the porous medium). The solution is linearly dependent on the macroscopic pressure gradient and on the volume forcing $\boldsymbol{S}$ , i.e.

(7.17) $$\begin{eqnarray}\displaystyle & \displaystyle u_{i}=-\tilde{a}_{ij}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X_{j}}+\tilde{b}_{ij}S_{j}, & \displaystyle\end{eqnarray}$$

(7.18) $$\begin{eqnarray}\displaystyle & \displaystyle p=-\tilde{c}_{j}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X_{j}}+\tilde{d}_{j}S_{j}, & \displaystyle\end{eqnarray}$$

with (7.17) leading directly to (7.14). Two microscopic systems can be set up and read

pore pressure related problem

(7.19a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\tilde{a}_{ij}}{\unicode[STIX]{x2202}x_{i}}=0,\quad -\frac{\unicode[STIX]{x2202}\tilde{c}_{j}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}\tilde{a}_{ij}}{\unicode[STIX]{x2202}x_{k}^{2}}+H(-y)\unicode[STIX]{x1D6FF}_{ij}=0,\end{eqnarray}$$

outer shear related problem

(7.20a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\tilde{b}_{ij}}{\unicode[STIX]{x2202}x_{i}}=0,\quad -\frac{\unicode[STIX]{x2202}\tilde{d}_{j}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}^{2}\tilde{b}_{ij}}{\unicode[STIX]{x2202}x_{k}^{2}}+\unicode[STIX]{x1D6FF}(y)\unicode[STIX]{x1D6FF}_{ij}=0,\end{eqnarray}$$

with boundary conditions easily available from those of the original equations for $(u_{i},p)$ . Once solutions of the auxiliary problems above are found, and given that a macroscopic pressure condition is eventually needed to couple the Stokes system above the interface to the Darcy system below, the next step consists in the intrinsic volume averaging of the pressure in the unit cells right above and right below the interface, i.e.

(7.21) $$\begin{eqnarray}\displaystyle & \displaystyle P^{+}:=\langle p\rangle _{+}^{f}=\int _{0}^{1}\lceil p\rceil \,\text{d}y=-c_{j}^{+}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X_{j}}+d_{j}^{+}S_{j}, & \displaystyle\end{eqnarray}$$

(7.22) $$\begin{eqnarray}\displaystyle & \displaystyle P^{-}:=\langle p\rangle _{-}^{f}=\frac{1}{{\mathcal{V}}_{f}}\int _{0}^{1}\int _{-1}^{0}\int _{0}^{1}p\,\text{d}x\,\text{d}y\,\text{d}z=-c_{j}^{-}\frac{\unicode[STIX]{x2202}p^{(0)}}{\unicode[STIX]{x2202}X_{j}}+d_{j}^{-}S_{j}. & \displaystyle\end{eqnarray}$$

The microscopic pressure $p$ in (7.22) is set to zero in the points where solid inclusions are present, in the unit cell right below (and tangent to) the Stokes–Darcy dividing plane; the coefficients of (7.21)–(7.22) are

(7.23) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle c_{j}^{+}=\int _{0}^{1}\lceil \tilde{c}_{j}\rceil \,\text{d}y,\quad d_{j}^{+}=\int _{0}^{1}\lceil \tilde{d}_{j}\rceil \,\text{d}y,\\ \displaystyle c_{j}^{-}=\frac{1}{{\mathcal{V}}_{f}}\int _{0}^{1}\int _{-1}^{0}\int _{0}^{1}\tilde{c}_{j}\,\text{d}x\,\text{d}y\,\text{d}z,\quad d_{j}^{-}=\frac{1}{{\mathcal{V}}_{f}}\int _{0}^{1}\int _{-1}^{0}\int _{0}^{1}\tilde{d}_{j}\,\text{d}x\,\text{d}y\,\text{d}z.\end{array}\right\}\end{eqnarray}$$

Eventually, the macroscopic pressure jump is $\unicode[STIX]{x0394}P=P^{+}-P^{-}$ , a linear function of the mean velocity within the porous medium (through $\unicode[STIX]{x2202}p^{(0)}/\unicode[STIX]{x2202}X_{j}$ ) and of the Navier slip speed (through $S_{j}$ ). The macroscopic pressure is thus discontinuous at the interface, as already noted by other researchers in the past (see, e.g. Carraro et al. (Reference Carraro, Marǔsić-Paloka and Mikelić2018) and references therein). Whereas it might seem contrived to enforce a macroscopic pressure jump across a surface through which the microscopic pressure $p$ is continuous, sample calculations by Lācis et al. (Reference Lācis, Sudhakar, Pasche and Bagheri2019) demonstrate that imposing such a discontinuity gives solutions in excellent agreement with feature-resolving direct numerical simulation results, whereas the more common choice (i.e. Ene & Sanchez-Palencia Reference Ene and Sanchez-Palencia1975; Lācis & Bagheri Reference Lācis and Bagheri2017) of pressure continuity at the Darcy–Stokes dividing surface yields rather large differences.