2.1 Considered axisymmetric flow fields

Consider, as sketched in figure 1, a Newtonian fluid (liquid or gas), with uniform density $\unicode[STIX]{x1D70C}$ and viscosity $\unicode[STIX]{x1D707}$ , embedding a porous slab domain ${\mathcal{S}}$ with thickness $e>0$ in $-e<z<0$ and two open fluid domains in $z<-e$ and $z>0$ . In the absence of the sphere, the fluid velocity is uniform, i.e. $-V_{0}\boldsymbol{e}_{z}$ in each domain. The fluid pressure is a constant $p_{0}>0$ for $z>0$ , zero for $z<-e$ and obeys the Darcy equation in the slab with uniform Darcy permeability $k^{\ast }>0$ . The pressure $Q$ in the slab therefore reads $Q=\unicode[STIX]{x1D707}V_{0}(z+e)/k^{\ast }$ . Finally, the fluid pressure is continuous across the porous slab boundaries $z=-e$ and $z=0$ , called here the lower and upper surfaces $\unicode[STIX]{x1D6F4}_{l}$ and $\unicode[STIX]{x1D6F4}_{u}$ , respectively. Thus, the so-called background permeation velocity $V_{0}$ in the whole space, imposed by the pressure $p_{0}$ , is given by $V_{0}=p_{0}k^{\ast }/(\unicode[STIX]{x1D707}e)$ .

Figure 1. A solid sphere with radius $R$ settling towards a permeable porous slab ${\mathcal{S}}$ with permeability $k^{\ast }>0$ and thickness $e$ .

A solid sphere, with centre $O^{\prime }$ and radius $R$ , is immersed above the slab in the ambient permeation flow $-V_{0}\boldsymbol{e}_{z}$ . It translates towards the slab at the steady velocity $-V_{p}\boldsymbol{e}_{z}$ . The resulting perturbed fluid flow has a quasi-steady velocity $\boldsymbol{u}-V_{0}\boldsymbol{e}_{z}$ and a pressure $p_{0}+p$ in the fluid domain ${\mathcal{D}}$ above the slab and a quasi-steady velocity $\boldsymbol{u}^{\prime }-V_{0}\boldsymbol{e}_{z}$ and a pressure $Q+p^{\prime }$ in the slab ${\mathcal{S}}$ . As will be shown below, it is irrelevant to consider in the present work the velocity and pressure disturbances occurring below the slab (i.e. in the $z<-e$ domain) and such quantities are thus henceforth discarded. The flow is governed by the quasi-steady Navier–Stokes and Darcy equations in ${\mathcal{D}}$ and ${\mathcal{S}}$ , respectively. Accordingly,

(2.1a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}(\boldsymbol{u}-V_{0}\boldsymbol{e}_{z})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}^{2}\boldsymbol{u}\quad \text{and}\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0\quad \text{in }{\mathcal{D}}, & \displaystyle\end{eqnarray}$$

(2.2a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}^{\prime }=-\frac{k^{\ast }}{\unicode[STIX]{x1D707}}\unicode[STIX]{x1D735}p^{\prime }\quad \text{and}\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\prime }=0\quad \text{in }{\mathcal{S}}. & \displaystyle\end{eqnarray}$$

The problem (2.1)–(2.2) must be supplemented with proper boundary conditions. Denoting by $S$ the sphere boundary, the following far-field condition and no-slip boundary condition on the sphere apply:

(2.3a,b ) $$\begin{eqnarray}(\boldsymbol{u},p)\rightarrow (\mathbf{0},0)\quad \text{at }\infty ,\quad \boldsymbol{u}=(V_{0}-V_{p})\boldsymbol{e}_{z}\quad \text{on }S.\end{eqnarray}$$

Additional conditions on the permeable $z=0$ upper surface $\unicode[STIX]{x1D6F4}_{u}$ are also needed. The continuity of both the velocity component normal to the slab and the pressure holds there

(2.4a,b ) $$\begin{eqnarray}\boldsymbol{u}\cdot \boldsymbol{e}_{z}=\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{e}_{z},\quad p=p^{\prime }\quad \text{on }\unicode[STIX]{x1D6F4}_{u}(z=0).\end{eqnarray}$$

In the present work we also impose on $\unicode[STIX]{x1D6F4}_{u}$ the Saffman (Reference Saffman1971) (SA) slip boundary condition for the velocities $\boldsymbol{u}_{t}=\boldsymbol{u}-(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{e}_{z})\boldsymbol{e}_{z}$ tangent to the slab. This condition reads

(2.5) $$\begin{eqnarray}\text{SA:}\frac{\unicode[STIX]{x2202}\boldsymbol{u}_{t}}{\unicode[STIX]{x2202}z}=\left(\frac{\unicode[STIX]{x1D70E}}{\sqrt{k^{\ast }}}\right)\boldsymbol{u}_{t}\quad \text{on }\unicode[STIX]{x1D6F4}_{u}(z=0),\end{eqnarray}$$

where the dimensionless slip parameter $\unicode[STIX]{x1D70E}$ is positive. In practice, the Darcy law is obtained when the porous slab thickness $e$ is much larger than the typical length of the pores $r_{p}=\sqrt{k^{\ast }}$ . Accordingly, we assume that $e\gg \sqrt{k^{\ast }}$ . In addition, the Saffman condition suggests introducing the so-called slab slip length $l=\sqrt{k^{\ast }}/\unicode[STIX]{x1D70E}$ which means that $\unicode[STIX]{x1D70E}=r_{p}/l$ . The widely employed Navier (Reference Navier1823) slip condition on a solid and impermeable wall with given slip length $l>0$ is retrieved by taking $\unicode[STIX]{x1D70E}=\sqrt{k^{\ast }}/l$ and letting $k^{\ast }$ vanish. The case of a permeable and no-slip slab, considered in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013), is retrieved from the Saffman boundary condition with $k^{\ast }>0$ and $\unicode[STIX]{x1D70E}\rightarrow \infty$ . In summary, the retained assumptions are

(2.6a,b ) $$\begin{eqnarray}\sqrt{k^{\ast }}\ll e;\quad \unicode[STIX]{x1D70E}>0.\end{eqnarray}$$

For symmetry reasons both flows $(\boldsymbol{u},p)$ and $(\boldsymbol{u}^{\prime },p^{\prime })$ are axisymmetric, with no swirl, about the $(O,\boldsymbol{e}_{z})$ axis. We use Cartesian coordinates $(O,x,y,z)$ and cylindrical coordinates $(r,z)$ with the unit local vector $\boldsymbol{e}_{r}$ such that $r=\{x^{2}+{y^{2}\}}^{1/2}$ and $\boldsymbol{x}=z\boldsymbol{e}_{z}+r\boldsymbol{e}_{r}$ . Accordingly, $p(\boldsymbol{x})=p(r,z)$ and $\boldsymbol{u}(\boldsymbol{x})=u_{r}(r,z)\boldsymbol{e}_{r}+w(r,z)\boldsymbol{e}_{z}$ with similar properties for $p^{\prime }$ and $\boldsymbol{u}^{\prime }$ , the cylindrical components of which are $u_{r}^{\prime }$ and $w^{\prime }$ . The stress tensor of the flow field $(\boldsymbol{u},p)$ is denoted as $\unicode[STIX]{x1D748}$ . For symmetry reasons, the force $\boldsymbol{F}_{s}$ experienced by the translating sphere then reads $\boldsymbol{F}_{s}=F_{s}\boldsymbol{e}_{z}$ with

(2.7) $$\begin{eqnarray}F_{s}=\int _{S}\boldsymbol{e}_{z}\boldsymbol{\cdot }[\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}-p_{0}\boldsymbol{n}]\,\text{d}S=\int _{S}\boldsymbol{e}_{z}\boldsymbol{\cdot }\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S,\end{eqnarray}$$

where $\boldsymbol{n}$ denotes the unit normal vector on $S$ pointing into the liquid (see figure 1). The sphere centre is $O^{\prime }$ and each point $M$ on its boundary $S$ is specified by an angle $\unicode[STIX]{x1D703}\in [0,\unicode[STIX]{x03C0}]$ with $\boldsymbol{O}^{\prime }\boldsymbol{M}=R\boldsymbol{n}$ and $\boldsymbol{n}=\sin \unicode[STIX]{x1D703}\boldsymbol{e}_{r}+\cos \unicode[STIX]{x1D703}\boldsymbol{e}_{z}$ . Consequently,

(2.8) $$\begin{eqnarray}F_{s}=2\unicode[STIX]{x03C0}R^{2}\int _{0}^{\unicode[STIX]{x03C0}}\left\{\cos \unicode[STIX]{x1D703}\left[2\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}z}-p\right]+\unicode[STIX]{x1D707}\sin \unicode[STIX]{x1D703}\left[\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}r}+\frac{\unicode[STIX]{x2202}u_{r}}{\unicode[STIX]{x2202}z}\right]\right\}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$

2.2 Lubrication approximation and resulting problem for the pressure

In our notation, the continuity equation in (2.1) becomes the linear relation

(2.9) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}z}=-\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}[ru_{r}],\end{eqnarray}$$

while the momentum equation in (2.1) is nonlinear and thus much less tractable. Now, for a near-contact sphere the well-known lubrication approximation is relevant and makes it possible to solely retain the linear part of this equation. From the assumption that the sphere–slab gap $\unicode[STIX]{x1D6FF}$ is small compared with the sphere radius $R$ , the classical (see Taylor, G. I., private communication cited in Handy & Bircumshaw Reference Handy and Bircumshaw1925) lubrication analysis consists of asymptotically expanding in terms of the small parameter $\unicode[STIX]{x1D716}=\unicode[STIX]{x1D6FF}/R$ the flow velocity components $(u_{r},w)$ and pressure $p$ in the vicinity of the origin $O$ (recall figure 1), i.e. for $z=O(\unicode[STIX]{x1D6FF})$ and $r=O(R)$ . In the present paper the analysis is restricted to the first-order approximation. Let $U$ and $W$ denote the typical magnitudes of the velocity components $u_{r}$ and $u_{z}$ in the lubrication domain. The continuity equation (2.9) first gives $W=\unicode[STIX]{x1D716}U$ . In approximating the momentum equation in (2.1) with boundary condition (2.3), the following Reynolds number $Re_{0}$ arises: $Re_{0}=\unicode[STIX]{x1D70C}|V_{0}-V_{p}|\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D707}$ . We assume that $\unicode[STIX]{x1D716}=\unicode[STIX]{x1D6FF}/R$ is small enough so that $Re_{0}\ll 1$ . Then, the nonlinear term on the left-hand side of the momentum equation (2.1) vanishes in the first approximation and, at the leading order in $\unicode[STIX]{x1D716}$ , this equation becomes

(2.10a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}r}=\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}^{2}u_{r}}{\unicode[STIX]{x2202}z^{2}},\quad \frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}z}=0.\end{eqnarray}$$

In the lubrication domain the following scalings, also obtained using (2.9)–(2.10), hold:

(2.11a-d ) $$\begin{eqnarray}p\sim \unicode[STIX]{x1D707}\frac{U}{R}\left(\frac{R}{\unicode[STIX]{x1D6FF}}\right)^{2},\quad \unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}z}\sim \unicode[STIX]{x1D707}\frac{W}{\unicode[STIX]{x1D6FF}}=\unicode[STIX]{x1D707}\frac{U}{R},\quad \unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}r}\sim \unicode[STIX]{x1D707}\frac{\unicode[STIX]{x1D6FF}U}{R^{2}},\quad \unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}u_{r}}{\unicode[STIX]{x2202}z}\sim \unicode[STIX]{x1D707}\frac{U}{\unicode[STIX]{x1D6FF}}.\end{eqnarray}$$

In the lubrication approximation the leading contribution to the force $F_{s}$ given by (2.8) is due to the part of the sphere located near the slab, i.e. to points of $S$ for which $r=R\sin \unicode[STIX]{x1D703}$ with $\unicode[STIX]{x1D703}$ close to $\unicode[STIX]{x03C0}$ . Furthermore, using (2.11) on the sphere boundary $S$ shows that it is possible to retain only the pressure term on the right-hand side of (2.8). Upon performing the change of variable $r=r(\unicode[STIX]{x1D703})$ , and recalling from (2.10) that the pressure $p$ is independent of $z$ in the lubrication domain gives

(2.12) $$\begin{eqnarray}F_{s}\sim 2\unicode[STIX]{x03C0}\int _{0}^{{\mathcal{R}}}p(r)r\,\text{d}r.\end{eqnarray}$$

In (2.12) the upper limit ${\mathcal{R}}$ is large compared with $\unicode[STIX]{x1D6FF}$ . It will be shown below, after (2.20), that the present (inner) expansion at scale $\unicode[STIX]{x1D6FF}$ is valid for ${\mathcal{R}}\rightarrow \infty$ .

In summary, the task in this work reduces to the determination of the approximated lubrication pressure $p$ and velocities $(u_{r},w)$ governed by the linear equations (2.9) and (2.7) and the far-field behaviour and boundary conditions (2.3)–(2.5).

The conditions (2.5) may be recast into the following form:

(2.13a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}u_{r}}{\unicode[STIX]{x2202}z}=\unicode[STIX]{x1D706}u_{r}\quad \text{for }z=0,\quad \unicode[STIX]{x1D706}=\frac{1}{l}=\frac{\unicode[STIX]{x1D70E}}{\sqrt{k^{\ast }}}.\end{eqnarray}$$

Now, the second equation (2.10) yields $p=p(r)$ . Integrating twice the first equation (2.10) and taking into account (2.13) give the velocity profile

(2.14) $$\begin{eqnarray}u_{r}(r,z)=[1+\unicode[STIX]{x1D706}z]u_{r}(r,0)+\frac{z^{2}}{2\unicode[STIX]{x1D707}}\left[\frac{\text{d}p}{\text{d}r}\right](r).\end{eqnarray}$$

At this stage the velocity $u_{r}(r,0)$ on the interface is still unknown. It is obtained here by enforcing the boundary condition (2.3) on the sphere surface $S$ in the radial direction. As shown in figure 1, this surface is described near the porous slab side by the equation $z=h(r)$ . Therefore, equation (2.3) requires that $u_{r}(r,h)=0$ . With this condition the result for the radial velocity profile is

(2.15) $$\begin{eqnarray}u_{r}(r,z)=\frac{1}{2\unicode[STIX]{x1D707}}\left[\frac{\text{d}p}{\text{d}r}\right](r)\left\{z^{2}-h^{2}\left(\frac{\unicode[STIX]{x1D706}z+1}{\unicode[STIX]{x1D706}h+1}\right)\right\}.\end{eqnarray}$$

Next, integrating the continuity equation (2.9) over $z$ in $[0,h(r)]$ yields

(2.16) $$\begin{eqnarray}w(r,h)=w(r,0)-\frac{1}{r}\int _{0}^{h}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}[ru_{r}]\,\text{d}z.\end{eqnarray}$$

The boundary condition (2.3), when projected on the $z$ direction, also provides $w(r,h)=V_{0}-V_{p}$ while the first condition (2.4) gives $w(r,0)=w^{\prime }(r,0)$ . This latter quantity is given by the Darcy law (see first equation (2.2)), i.e. by $w^{\prime }(r,0)=-(k^{\ast }/\unicode[STIX]{x1D707})(\unicode[STIX]{x2202}p^{\prime }/\unicode[STIX]{x2202}z)_{z=0}$ . As in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013), we assume a linear variation of the pressure $p^{\prime }$ across the slab: $(\unicode[STIX]{x2202}p^{\prime }/\unicode[STIX]{x2202}z)_{z=0}\sim [p^{\prime }(r,0)-p^{\prime }(r,-e)]/e$ . As shown in appendix B, this assumption actually requires that $e\ll R$ ; that is, the porous slab thickness is small compared with the sphere radius. We also consider that the fluid below the slab imposes the pressure on the slab lower surface $\unicode[STIX]{x1D6F4}_{l}(z=-e)$ , i.e. that $p^{\prime }(r,-e)=0$ . Then we obtain the additional link $w^{\prime }(r,0)=-k^{\ast }p^{\prime }(r,0)/(\unicode[STIX]{x1D707}e)$ . Noting that $p^{\prime }(r,0)=p(r,0)=p(r)$ , the combination of (2.16) and (2.15) provides the governing equation for the pressure $p$ in the lubrication domain

(2.17) $$\begin{eqnarray}\frac{1}{12\unicode[STIX]{x1D707}r}\frac{\text{d}}{\text{d}r}\left[\left\{\frac{(4+\unicode[STIX]{x1D706}h)h^{3}}{1+\unicode[STIX]{x1D706}h}\right\}r\frac{\text{d}p}{\text{d}r}\right]=V_{0}-V_{p}+\frac{k^{\ast }p(r)}{\unicode[STIX]{x1D707}e}.\end{eqnarray}$$

Equation (2.17) may be seen as an extension of the classical Reynolds equation (Reynolds Reference Reynolds1886; Nguyen Reference Nguyen2000). There is no mass injection on the $(O,\boldsymbol{e}_{z})$ axis and thus $u_{r}^{\prime }(r,0)=0$ at $r=0$ . By symmetry, $(\text{d}p/\text{d}r)(r=0)$ vanishes. In addition $p$ should vanish for large $r$ . From the small normalized sphere–wall gap $\unicode[STIX]{x1D6FF}/R$ assumption, $h(r)\sim \unicode[STIX]{x1D6FF}+r^{2}/(2R)$ in the vicinity of the origin $O$ , namely for $r\leqslant O(\sqrt{\unicode[STIX]{x1D6FF}R})$ as classically found. As in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013), we henceforth use the lengths $(l_{2},k)$ and the dimensionless parameter $\unicode[STIX]{x1D6FD}$ , length $s$ and pressure $P$ such that

(2.18a-e ) $$\begin{eqnarray}l_{2}=\sqrt{2\unicode[STIX]{x1D6FF}R},\quad k=\frac{k^{\ast }}{e},\quad \unicode[STIX]{x1D6FD}=\frac{24kR}{\unicode[STIX]{x1D6FF}^{2}},\quad s=\frac{r}{l_{2}},\quad p=\frac{24\unicode[STIX]{x1D707}R}{\unicode[STIX]{x1D6FF}^{2}}(V_{p}-V_{0})P.\end{eqnarray}$$

Henceforth, the reduced permeability $\unicode[STIX]{x1D6FD}$ is simply called the permeability while $k$ is the permeance. In addition, we also introduce the dimensionless sphere-to-slab gap $\tilde{\unicode[STIX]{x1D706}}$ :

(2.19) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D706}}=\unicode[STIX]{x1D706}\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D6FF}/l=\frac{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FF}}{\sqrt{k^{\ast }}},\end{eqnarray}$$

and then obtain for $P$ the following dimensionless boundary value problem:

(2.20a,b ) $$\begin{eqnarray}\frac{1}{s}\frac{\text{d}}{\text{d}s}\left\{\left[\frac{(4+\tilde{\unicode[STIX]{x1D706}}H)H^{3}}{1+\tilde{\unicode[STIX]{x1D706}}H}\right]s\frac{\text{d}P}{\text{d}s}\right\}-\unicode[STIX]{x1D6FD}P=-1,\quad \left(\frac{\text{d}P}{\text{d}s}\right)_{s=0}=P(\infty )=0,\end{eqnarray}$$

with $H=1+s^{2}$ the normalized sphere-to-wall distance at position $s$ . The governing problem (2.20) extends that considered in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013) by applying the Saffman slip boundary condition. Not surprisingly, taking $\tilde{\unicode[STIX]{x1D706}}\rightarrow \infty$ in (2.20) recovers the Ramon formulation.

The problem (2.20) is well-posed for large $s$ as may be proven by considering its asymptotic expansion for $s\gg 1$ . Then $H\sim s^{2}$ and the differential equation may be recast in terms of the variable $H$

(2.21) $$\begin{eqnarray}4\frac{\text{d}}{\text{d}H}\left(H^{4}\frac{\text{d}P}{\text{d}H}\right)-\unicode[STIX]{x1D6FD}P=-1.\end{eqnarray}$$

For $H\gg 1$ the first term takes over and the formal solution behaves like $P\sim H^{-2}\sim s^{-4}$ . Thus, the condition at infinity in (2.20) can indeed be applied and the lubrication problem is regular as expected. Moreover, in the integral (2.12) for the force $F_{s}$ the upper limit ${\mathcal{R}}$ is rejected to infinity.

As in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013), the force $F_{s}$ is normalized by that obtained for an impermeable no-slip wall. Using our previous scalings (2.11) and (2.18) for $s$ and $P$ , the normalized force $F$ is

(2.22) $$\begin{eqnarray}F=\frac{\unicode[STIX]{x1D6FF}F_{s}}{6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}(V_{p}-V_{0})R^{2}}=16\int _{0}^{\infty }P(s)s\,\text{d}s.\end{eqnarray}$$

In this paper we focus on the normalized force $F$ which, by virtue of (2.19)–(2.22), depends on the sphere–slab gap $\unicode[STIX]{x1D6FF}$ , the sphere radius $R$ and the porous slab properties $(k^{\ast },e,\unicode[STIX]{x1D70E})$ through the dimensionless parameters $\unicode[STIX]{x1D6FD}$ and $\tilde{\unicode[STIX]{x1D706}}$ . The validity of the present formulation however brings some additional restrictions to the previously mentioned ones in (2.6). First, we must recall the condition, discussed after (2.16), of a sufficiently thin slab with $e\ll R$ . Second, the lubrication approximation is valid for $r\leqslant O(l_{2})$ and the associated domain on the upper slab boundary $z=0$ must contain a large number of pores, i.e.  $r_{p}=\sqrt{k^{\ast }}\ll l_{2}$ . In summary, we require in the present paper that for a slip porous slab

(2.23a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D70E}\geqslant 0,\quad \sqrt{k^{\ast }}\ll e\ll R,\quad \sqrt{k^{\ast }}\ll \sqrt{\unicode[STIX]{x1D6FF}R}.\end{eqnarray}$$

The particular value $\unicode[STIX]{x1D70E}=\infty$ is here employed with $k^{\ast }>0$ to retrieve the no-slip permeable slab (case $l=0$ ) handled in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013).

2.3 Numerical implementation and tests

For the no-slip permeable slab case (Saffman case with $\tilde{\unicode[STIX]{x1D706}}=\infty$ ) the solution of (2.20) for the lubrication pressure $P$ is obtained in closed form. This solution, surprisingly not reported in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013), is given in appendix A. In all other cases, equation (2.20) has to be solved numerically. This is achieved using a truncated domain $0\leqslant s\leqslant L$ and a $O(h^{2})$ second-order finite-difference scheme with $N+1$ nodal points $s_{n}=nh$ (taking $n=0,1,\ldots ,N$ ) and $h=L/N$ . We then calculate the normalized force $F$ taking as $L$ the upper value of the domain of integration in the integral occurring on the right-hand side of (2.22). This procedure is benchmarked for $\unicode[STIX]{x1D6FD}=1,100$ against the force obtained for the no-slip permeable case. That force is obtained here by integrating, using Maple computer algebra software and taking $L=\infty$ , the product $P(s)s$ with $P$ the analytical lubrication pressure displayed in appendix A.

As shown in table 1, the numerical results converge as $L$ increases and become very close to the analytical results. In view of these results we henceforth take $L=80$ which is sufficient for a three-digit accuracy.

3.1 Asymptotic expansion for $\unicode[STIX]{x1D6FD}\gg 1$

As achieved in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013) for the no-slip case, it is worth estimating for the Saffman slip case the resulting normalized lubrication pressure $P$ and force $F$ for either large or small values of the permeability $\unicode[STIX]{x1D6FD}$ . In this subsection we consider the case of a sphere translating towards a near-contact permeable slab $(k=k^{\ast }/e>0)$ in the limit $\unicode[STIX]{x1D6FD}\gg 1$ of large permeability. In such circumstances asymptotic expansions for $P$ and $F$ can be derived by extending the treatment proposed in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013). Curtailing the details, the pressure $P$ takes, for $\unicode[STIX]{x1D6FD}$ large, the following form:

(3.2a-c ) $$\begin{eqnarray}P(s)=\frac{P^{\ast }(\unicode[STIX]{x1D702})}{\unicode[STIX]{x1D6FD}},\quad \unicode[STIX]{x1D702}=\frac{s}{\unicode[STIX]{x1D6FD}^{1/4}},\quad P^{\ast }(\unicode[STIX]{x1D702})\sim P_{0}^{\ast }(\unicode[STIX]{x1D702})+\frac{P_{1}^{\ast }(\unicode[STIX]{x1D702})}{\unicode[STIX]{x1D6FD}^{1/2}}+\frac{P_{2}^{\ast }(\unicode[STIX]{x1D702})}{\unicode[STIX]{x1D6FD}},\end{eqnarray}$$

with, using Maple software, the functions

(3.3a,b ) $$\begin{eqnarray}\displaystyle & & \displaystyle P_{0}^{\ast }(\unicode[STIX]{x1D702})=1-\left[\frac{1+2\unicode[STIX]{x1D702}^{2}}{2\unicode[STIX]{x1D702}^{2}}\right]\text{e}^{-1/(2\unicode[STIX]{x1D702}^{2})},\quad P_{1}^{\ast }(\unicode[STIX]{x1D702})=-\frac{3}{16}\left[\frac{1+\tilde{\unicode[STIX]{x1D706}}}{\tilde{\unicode[STIX]{x1D706}}\unicode[STIX]{x1D702}^{6}}\right]\text{e}^{-1/(2\unicode[STIX]{x1D702}^{2})},\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle P_{2}^{\ast }(\unicode[STIX]{x1D702}) & = & \displaystyle -\frac{1}{256\unicode[STIX]{x1D702}^{10}}\left[9-58\unicode[STIX]{x1D702}^{2}+4\unicode[STIX]{x1D702}^{4}+\frac{1}{\tilde{\unicode[STIX]{x1D706}}}(18-116\unicode[STIX]{x1D702}^{2}+8\unicode[STIX]{x1D702}^{4})\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\frac{1}{\tilde{\unicode[STIX]{x1D706}}^{2}}(9+122\unicode[STIX]{x1D702}^{2}+68\unicode[STIX]{x1D702}^{4})\right]\text{e}^{-1/(2\unicode[STIX]{x1D702}^{2})}.\end{eqnarray}$$

Of course, letting $\tilde{\unicode[STIX]{x1D706}}$ tend to infinity in (3.2)–(3.4) retrieves the asymptotics given in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013). Observe that for (3.2) to be an asymptotic expansion, $\tilde{\unicode[STIX]{x1D706}}$ cannot vanish (inspect $P_{1}^{\ast }$ in (3.3) and $P_{2}^{\ast }$ in (3.4)). In practice, the expansion (3.2) is valid for $\sqrt{\unicode[STIX]{x1D6FD}}\gg \text{Max}(1,\tilde{\unicode[STIX]{x1D706}}^{-1})$ . On the basis of (3.2)–(3.4), an asymptotic expansion is derived for the dimensionless drag force from (2.22):

(3.5) $$\begin{eqnarray}F\sim \frac{4}{\unicode[STIX]{x1D6FD}^{1/2}}-\left[\frac{6(\tilde{\unicode[STIX]{x1D706}}+1)}{\tilde{\unicode[STIX]{x1D706}}}\right]\frac{1}{\unicode[STIX]{x1D6FD}}+\frac{3}{2}\left[\frac{17+\tilde{\unicode[STIX]{x1D706}}(2+\tilde{\unicode[STIX]{x1D706}})}{\tilde{\unicode[STIX]{x1D706}}^{2}}\right]\frac{1}{\unicode[STIX]{x1D6FD}^{3/2}}\quad \text{for }\sqrt{\unicode[STIX]{x1D6FD}}\gg \text{Max}(1,\tilde{\unicode[STIX]{x1D706}}^{-1}).\end{eqnarray}$$

Let us introduce the dimensionless permeability $\overline{k}=24k/R=24k^{\ast }/(eR)>0$ which compares the porous slab typical pore size $r_{p}=\sqrt{k^{\ast }}$ with the length $\sqrt{eR}$ . From our second assumption (2.23), it turns out that $k^{\ast }\ll eR$ , i.e. that $\overline{k}\ll 24$ . The requirement $\sqrt{\unicode[STIX]{x1D6FD}}\gg \text{Max}(1,\tilde{\unicode[STIX]{x1D706}}^{-1})$ then becomes $\sqrt{\overline{k}}\gg \text{Max}(\unicode[STIX]{x1D6FF}/R,l/R)$ . These constraints tell us to what extent the ratio $k^{\ast }/(eR)$ must be large enough for given parameters $(\unicode[STIX]{x1D6FF}/R,l/R)$ . Using (2.22) and (3.5) provides the following approximation of the dimensional force $F_{s}$ experienced by a sphere translating normal and very close to a permeable (porous) slab with upper surface subject to the Saffman condition:

(3.6) $$\begin{eqnarray}\displaystyle F_{s} & {\sim} & \displaystyle 6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}(V_{p}-V_{0})R\left\{\frac{4}{(\overline{k})^{1/2}}-\frac{6(\tilde{\unicode[STIX]{x1D706}}+1)}{\tilde{\unicode[STIX]{x1D706}}\overline{k}}\left(\frac{\unicode[STIX]{x1D6FF}}{R}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\frac{3}{2(\overline{k})^{3/2}}\left[\frac{17+\tilde{\unicode[STIX]{x1D706}}(2+\tilde{\unicode[STIX]{x1D706}})}{\tilde{\unicode[STIX]{x1D706}}^{2}}\right]\left(\frac{\unicode[STIX]{x1D6FF}}{R}\right)^{2}\right\}\quad \text{for }\sqrt{\overline{k}}\gg \text{Max}(\unicode[STIX]{x1D6FF}/R,l/R).\end{eqnarray}$$

According to (3.2)–(3.6), the first approximations of the lubrication pressure and of the resulting force exerted on a sphere close to the permeable $(k>0)$ porous slab are independent of the slip parameter $\tilde{\unicode[STIX]{x1D706}}$ . For a large permeability this was not a priori expected on physical grounds. The Saffman condition is found to affect the lubrication pressure and dimensional force only at the next order.

Table 1. Normalized force $F$ versus the length $L$ for the no-slip permeable case (i.e. as in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013)) and $\unicode[STIX]{x1D6FD}=1,100$ . For $L<\infty$ the lubrication pressure is calculated numerically by solving (2.20) for the Saffman boundary condition $(\unicode[STIX]{x1D6FF}_{BJ}=0)$ and $\tilde{\unicode[STIX]{x1D706}}\rightarrow \infty$ . The last column for $L=\infty$ is the analytical lubrication pressure given in appendix A.

The range of validity of the asymptotic expansion (3.5) is evaluated by comparing in table 2 its predictions against the numerical results of three-digit accuracy . Clearly, the asymptotic expansion does not actually hold for $\tilde{\unicode[STIX]{x1D706}}=1$ when $\unicode[STIX]{x1D6FD}=10,50$ . (This is actually still the case when $\tilde{\unicode[STIX]{x1D706}}=0.1$ for $\unicode[STIX]{x1D6FD}=1000$ with $F=0.0930$ and $F_{a}=0.1421$ and also for $\unicode[STIX]{x1D6FD}=5000$ with $F=0.0474$ and $F_{a}=0.0507$ . In contrast, for $\tilde{\unicode[STIX]{x1D706}}=0.1$ and $\unicode[STIX]{x1D6FD}=10\,000$ the closer values $F=0.0348$ and $F_{a}=0.0360$ are obtained.) However, as seen in table 2, the asymptotic expansion provides already for $\unicode[STIX]{x1D6FD}=50$ and $\tilde{\unicode[STIX]{x1D706}}\geqslant 10$ the same three-digit accuracy and may therefore be seen as a handy formula in the range $\unicode[STIX]{x1D6FD}\geqslant 50$ for this domain of the slip parameter $\tilde{\unicode[STIX]{x1D706}}$ . It also nicely matches the numerics for $\tilde{\unicode[STIX]{x1D706}}\geqslant 1$ when $\unicode[STIX]{x1D6FD}\geqslant 100$ .

Table 2. Computed normalized force in the Saffman case, $F$ , for a large range of values of the dimensionless permeability, $\unicode[STIX]{x1D6FD}$ , and asymptotic estimate, $F_{a}$ , for either $\unicode[STIX]{x1D6FD}\gg 1$ (§ 3.1) or $\unicode[STIX]{x1D6FD}\ll 1$ (§ 3.2), both for various values of the parameter $\tilde{\unicode[STIX]{x1D706}}$ (see (2.19)). The asymptotic values $F_{a}$ for $\unicode[STIX]{x1D6FD}\leqslant 1$ have been obtained using the small $\unicode[STIX]{x1D6FD}$ approximation given in § 3.2. Each line with a $\ast$ symbol reports the closed-form (order $O(\unicode[STIX]{x1D6FD})$ ) approximation in small $\unicode[STIX]{x1D6FD}$ obtained by retaining only the first two terms in (3.10).

3.2 Asymptotic expansion for $\unicode[STIX]{x1D6FD}\ll 1$

For a vanishing porous slab permeance $k=k^{\ast }/e$ , the dimensionless permeability $\unicode[STIX]{x1D6FD}$ might be small compared with unity even for a very small normalized sphere–slab gap $\unicode[STIX]{x1D6FF}/R$ (see (2.23)). More precisely, since $\unicode[STIX]{x1D6FD}=\overline{k}(R/\unicode[STIX]{x1D6FF})^{2}$ the condition $\unicode[STIX]{x1D6FD}\ll 1$ is ensured for $\sqrt{\overline{k}}\ll \unicode[STIX]{x1D6FF}/R$ . In such circumstances a regular expansion of the problem (2.20) for small  $\unicode[STIX]{x1D6FD}$ is worked out, as also done in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013) for the no-slip case. Omitting the details, the pressure $P$ is approximated as

(3.7) $$\begin{eqnarray}P(s)=P_{0}(s)+\unicode[STIX]{x1D6FD}P_{1}(s)+\unicode[STIX]{x1D6FD}^{2}P_{2}(s)+O(\unicode[STIX]{x1D6FD}^{3})\quad \text{for }\unicode[STIX]{x1D6FD}\ll 1,\end{eqnarray}$$

with the zeroth-order pressure $P_{0}$ and first-order pressure $P_{1}$ derived in closed form (again using the Maple software), that is

(3.8) $$\begin{eqnarray}P_{0}(s)=\frac{3\tilde{\unicode[STIX]{x1D706}}^{2}}{256}\log \left[\frac{\tilde{\unicode[STIX]{x1D706}}(1+s^{2})}{4+\tilde{\unicode[STIX]{x1D706}}(1+s^{2})}\right]+\frac{2+3\tilde{\unicode[STIX]{x1D706}}(1+s^{2})}{64(1+s^{2})^{2}}\end{eqnarray}$$

and a much more complicated form for $P_{1}$ given in appendix C. Finally, the second-order pressure $P_{2}$ was calculated numerically using a finite-difference scheme (as in § 2.3) to solve the following problem (from the extended Reynolds equation (2.20)):

(3.9a,b ) $$\begin{eqnarray}\frac{1}{s}\frac{\text{d}}{\text{d}s}\left\{\left[\frac{(4+\tilde{\unicode[STIX]{x1D706}}H)H^{3}}{1+\tilde{\unicode[STIX]{x1D706}}H}\right]s\frac{\text{d}P_{2}}{\text{d}s}\right\}=P_{1},\quad \left(\frac{\text{d}P_{2}}{\text{d}s}\right)_{s=0}=P_{2}(\infty )=0.\end{eqnarray}$$

Following (3.7) and applying (2.22), the expansion for the normalized force $F$ reads

(3.10) $$\begin{eqnarray}F=F_{0}(\tilde{\unicode[STIX]{x1D706}})+\unicode[STIX]{x1D6FD}F_{1}(\tilde{\unicode[STIX]{x1D706}})+\unicode[STIX]{x1D6FD}^{2}F_{2}(\tilde{\unicode[STIX]{x1D706}})+O(\unicode[STIX]{x1D6FD}^{3})\quad \text{for }\unicode[STIX]{x1D6FD}\ll 1.\end{eqnarray}$$

Although the analytical expression for the pressure $P_{1}$ is very involved (see appendix C), it is nevertheless possible to find a simple expression for $F_{1}(\tilde{\unicode[STIX]{x1D706}})$ . The results are

(3.11a,b ) $$\begin{eqnarray}\displaystyle F_{0}(\tilde{\unicode[STIX]{x1D706}})=\frac{1}{8}(2-3\tilde{\unicode[STIX]{x1D706}})+\frac{3\tilde{\unicode[STIX]{x1D706}}}{32}(4+\tilde{\unicode[STIX]{x1D706}})X,\quad X=\log \left(\frac{4+\tilde{\unicode[STIX]{x1D706}}}{\tilde{\unicode[STIX]{x1D706}}}\right), & & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle F_{1}(\tilde{\unicode[STIX]{x1D706}})=-\frac{1}{384}-\frac{3\tilde{\unicode[STIX]{x1D706}}}{256}-\frac{3\tilde{\unicode[STIX]{x1D706}}^{2}(4-X)}{512}+\frac{3\tilde{\unicode[STIX]{x1D706}}^{3}X}{2048}[1+3X]+\frac{9\tilde{\unicode[STIX]{x1D706}}^{4}X^{2}}{8192}. & & \displaystyle\end{eqnarray}$$

By inspection of (3.11)–(3.12) and as shown in figure 2 it is observed that both coefficients $F_{0}$ and $F_{1}$ remain of the order of unity in the entire range $\tilde{\unicode[STIX]{x1D706}}\geqslant 0$ .

It should be noted that $F_{0}(\tilde{\unicode[STIX]{x1D706}})$ is actually the normalized force experienced by a sphere located close to (i.e. in the lubrication regime $\unicode[STIX]{x1D6FF}\ll R$ ) an impermeable surface on which the usual Navier slip condition holds. In this case the result (3.11) is recovered from Vinogradova (Reference Vinogradova1995) lubrication analysis (see equations (3.18)–(3.19) there), by letting in her notation $R_{2}\rightarrow \infty ,b\rightarrow 0$ and $k\rightarrow \infty$ with $kb=l$ being the slip length. We used Maple computer algebra software to check this limit. In principle, we should also retrieve the function $F_{0}(\tilde{\unicode[STIX]{x1D706}})$ from Goren (Reference Goren1973) lubrication results (his equation (17)), converting his notation with $K\!n_{s}=l/R,h=\unicode[STIX]{x1D6FF}$ and $K\!n_{w}\rightarrow 0$ . (We could not however succeed in retrieving his result and suspect a misprint in his equation (17) for two reasons: first, there are two analogous terms $4ST$ and $2ST$ in the same parentheses and, second, his results from table 1 are not recovered from his formula (17).)

As $\tilde{\unicode[STIX]{x1D706}}$ becomes large the no-slip boundary condition applies and the zeroth-order pressure $P_{0}$ tends to the function obtained in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013) while $F_{0}(\tilde{\unicode[STIX]{x1D706}})\sim 1$ and $F_{1}(\tilde{\unicode[STIX]{x1D706}})\sim -1/24$ . We also recover the results by Goren (Reference Goren1979) for a very thin membrane.

Like $P_{0}$ and $P_{1}$ , the pressure $P_{2}$ solution of (3.9) is of order unity whatever the value of $\tilde{\unicode[STIX]{x1D706}}\geqslant 0$ . Therefore, and in contrast to the previous case of large permeability $\unicode[STIX]{x1D6FD}$ , no restriction bearing on the slip parameter $\tilde{\unicode[STIX]{x1D706}}$ appears here for the expansions (3.7) and (3.10) to be valid. Again invoking (2.22) eventually provides the following expansion for the dimensionless force:

(3.13) $$\begin{eqnarray}\displaystyle F_{s} & = & \displaystyle 6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}(V_{p}-V_{0})R\left\{F_{0}(\tilde{\unicode[STIX]{x1D706}})\left(\frac{R}{\unicode[STIX]{x1D6FF}}\right)+F_{1}(\tilde{\unicode[STIX]{x1D706}})\overline{k}\left(\frac{R}{\unicode[STIX]{x1D6FF}}\right)^{3}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,F_{2}(\tilde{\unicode[STIX]{x1D706}})\overline{k}^{2}\left(\frac{R}{\unicode[STIX]{x1D6FF}}\right)^{5}+O\left(\overline{k}^{3}\left(\frac{R}{\unicode[STIX]{x1D6FF}}\right)^{7}\right)\right\}\quad \text{for }\unicode[STIX]{x1D6FD}\ll 1.\end{eqnarray}$$

As for the case of large permeability, the asymptotic estimate (3.10) of the normalized force $F$ for $\unicode[STIX]{x1D6FD}\ll 1$ was tested against the numerical predictions. The results displayed in table 2 show the approximation (3.10) both at $O(\unicode[STIX]{x1D6FD}^{2})$ ‘first order’ and $O(\unicode[STIX]{x1D6FD}^{3})$ ‘second order’. From table 2 the agreement between the three-digit numerics and both (first-order and second-order) asymptotic approximations for the normalized force turns out to be amazingly good (that is, within $O(10^{-4})$ accuracy) up to $\unicode[STIX]{x1D6FD}=0.1$ in the entire range $\tilde{\unicode[STIX]{x1D706}}\geqslant 1$ . Note also that the first-order approximation, obtained by retaining the first two terms on the right-hand side of (3.10), may be considered as a handy formula for practical use.

Figure 2. Coefficients $F_{0}$  (○) and $F_{1}$  (●) of the normalized force, see (3.11) and (3.12), versus $\tilde{\unicode[STIX]{x1D706}}$ for low dimensionless permeability, $\unicode[STIX]{x1D6FD}\ll 1$ .

3.3 Numerical results

In practice the permeability parameter $\unicode[STIX]{x1D6FD}$ can take values in the entire range $\unicode[STIX]{x1D6FD}\geqslant 0$ . For $\unicode[STIX]{x1D6FD}$ neither small nor large the lubrication pressure $P$ and resulting normalized force $F$ have to be calculated numerically. This is achieved in the present work, with a three-digit accuracy, using the numerical method described in § 2.3.

Figure 3 displays the force $F$ versus $\unicode[STIX]{x1D6FD}$ in the range $[0.01,1000]$ for the Saffman slip case $\tilde{\unicode[STIX]{x1D706}}=0.1,1,10,100$ and also for the Ramon no-slip case (associated with the value $\tilde{\unicode[STIX]{x1D706}}=\infty$ ). The slip parameter is also written $\tilde{\unicode[STIX]{x1D706}}=\unicode[STIX]{x1D6FF}/l$ with $l$ the porous slab slip length. The curves in figure 3 reveal that $F$ is deeply sensitive to the ratio $\unicode[STIX]{x1D6FF}/l$ in a wide range (say here in the range $[0.01,100]$ ) of the permeability $\unicode[STIX]{x1D6FD}$ . Not surprisingly, at given values of $\unicode[STIX]{x1D6FD}$ and of the small sphere–slab gap $\unicode[STIX]{x1D6FF}$ , the normalized force decreases as the slip increases (i.e. as $l$ increases). Moreover, at given slip parameter $\tilde{\unicode[STIX]{x1D706}}$ the force quickly reaches for $\unicode[STIX]{x1D6FD}\leqslant O(1)$ its limit value for the impermeable $(\unicode[STIX]{x1D6FD}=0)$ slab case. This limit value agrees with the one obtained by Vinogradova (Reference Vinogradova1995) for the impermeable slipping wall case. Of course, as shown by the leading order in expansion (3.5), all curves collapse as $\unicode[STIX]{x1D6FD}$ becomes large enough. Moreover, it is observed in figure 3 that curves for $\tilde{\unicode[STIX]{x1D706}}\geqslant 1$ almost coincide as soon as $\unicode[STIX]{x1D6FD}\geqslant 1000$ .

Figure 3. Normalized force $F$ versus the permeability parameter $\unicode[STIX]{x1D6FD}$ for $\tilde{\unicode[STIX]{x1D706}}=0.1$  (○),  $\tilde{\unicode[STIX]{x1D706}}=1$ (▫), $\tilde{\unicode[STIX]{x1D706}}=10$ (▵) and $\tilde{\unicode[STIX]{x1D706}}=100$ (solid curve). For comparison purposes, a few values obtained for the no-slip Ramon case (associated with $\tilde{\unicode[STIX]{x1D706}}=\infty$ ) are also indicated by $\ast$ symbols.

The lubrication pressure $P(s)$ is also displayed in figure 4 versus the normalized radial coordinate in the gap, $s$ , for $\unicode[STIX]{x1D6FD}=1,100$ and the same values of $\tilde{\unicode[STIX]{x1D706}}$ as in figure 3. As shown for $\unicode[STIX]{x1D6FD}=1$ in figure 4(a), the pressure is strongly sensitive to the value of $\tilde{\unicode[STIX]{x1D706}}$ but is mainly exerted in the narrow domain $s\leqslant 1.5$ close to the symmetry axis. In contrast, figure 4(b) reveals that for $\unicode[STIX]{x1D6FD}=100$ the pressure becomes less sensitive to the parameter $\tilde{\unicode[STIX]{x1D706}}$ but is mainly distributed in the extended range $s\leqslant 3$ .

Figure 4. Lubrication pressure $P$ versus $s$ for $\tilde{\unicode[STIX]{x1D706}}=0.1$  (○),  $\tilde{\unicode[STIX]{x1D706}}=1$ (▫), $\tilde{\unicode[STIX]{x1D706}}=10$ (▵), $\tilde{\unicode[STIX]{x1D706}}=100$ (solid curve) and also for the $\tilde{\unicode[STIX]{x1D706}}=\infty$ case (no curve but * symbols). (a) Case $\unicode[STIX]{x1D6FD}=1$ . (b) Case $\unicode[STIX]{x1D6FD}=100$ .

4.1 Time-dependent problem and numerical treatment

As displayed in figure 1, the sphere centre is $O^{\prime }$ with $\boldsymbol{O}\boldsymbol{O}^{\prime }=(R+\unicode[STIX]{x1D6FF})\boldsymbol{e}_{z}$ . The sphere–slab gap $\unicode[STIX]{x1D6FF}$ now depends on time $t$ . The sphere density $\unicode[STIX]{x1D70C}_{s}$ is uniform, so that the sphere settles without rotating at the translational velocity $-V_{p}\boldsymbol{e}_{z}$ with $V_{p}=-\text{d}\unicode[STIX]{x1D6FF}/\text{d}t$ . Since the fluid has a uniform density $\unicode[STIX]{x1D70C}$ the sphere buoyancy-corrected mass is $m=4\unicode[STIX]{x03C0}(\unicode[STIX]{x1D70C}_{s}-\unicode[STIX]{x1D70C})R^{3}/3$ . Accounting for gravity, the flow field perturbed by the sphere is $(\boldsymbol{u},p+\unicode[STIX]{x1D70C}\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{x})$ in the fluid domain ${\mathcal{D}}$ and $(\boldsymbol{u}^{\prime },p^{\prime }+\unicode[STIX]{x1D70C}\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{x})$ in the slab ${\mathcal{S}}$ with $(\boldsymbol{u},p)$ and $(\boldsymbol{u}^{\prime },p^{\prime })$ still governed by (2.2)–(2.5) and the lubrication approximation (2.10). The total (hydrodynamic plus gravity) force exerted on the sphere with volume ${\mathcal{V}}$ is $(F_{s}+\unicode[STIX]{x1D70C}{\mathcal{V}}g)\boldsymbol{e}_{z}$ where $F_{s}$ is related to the normalized force $F$ by (2.22) (in which one sets $V_{0}=0$ ). Neglecting the sphere inertia, which corresponds to either a liquid or a gas at low Stokes number, the time-dependent sphere ‘location’ $\unicode[STIX]{x1D6FF}(t)$ is the solution to the following problem:

(4.1a,b ) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D6FF}}{\text{d}t}=-\frac{mg\unicode[STIX]{x1D6FF}}{6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}R^{2}F},\quad \unicode[STIX]{x1D6FF}(t=0)=\unicode[STIX]{x1D6FF}_{0},\end{eqnarray}$$

with $\unicode[STIX]{x1D6FF}_{0}$ the prescribed initial sphere–slab gap.

As outlined after (2.22), the normalized force $F$ solely depends upon $\unicode[STIX]{x1D6FF}/R$ and the dimensionless parameters $(\unicode[STIX]{x1D6FD},\tilde{\unicode[STIX]{x1D706}})$ . From the definitions (3.1), the time-dependent normalized sphere ‘location’ $\overline{\unicode[STIX]{x1D6FF}}=\unicode[STIX]{x1D6FF}/R$ is thus determined by the time scale $\unicode[STIX]{x1D70F}_{g}=\sqrt{24}\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}R^{2}/(mg)$ as in Ramon et al. (Reference Ramon, Huppert, Lister and Stone2013), the given initial sphere location $\overline{\unicode[STIX]{x1D6FF}}_{0}=\unicode[STIX]{x1D6FF}_{0}/R$ and the auxiliary parameters $(k/R,l/R)$ . For convenience, the normalized time $\overline{t}=t/\unicode[STIX]{x1D70F}_{g}$ is introduced. From (4.1), we obtain the ratio $\overline{\unicode[STIX]{x1D6FF}}$ versus $\overline{t}$ , for given quantities $(\overline{\unicode[STIX]{x1D6FF}}_{0},k/R,l/R)$ , by solving the problem

(4.2a,b ) $$\begin{eqnarray}\frac{\text{d}\overline{\unicode[STIX]{x1D6FF}}}{\text{d}\overline{t}}=-\left[\frac{(2/3)^{1/2}}{F}\right]\overline{\unicode[STIX]{x1D6FF}},\quad \overline{\unicode[STIX]{x1D6FF}}(\overline{t}=0)=\overline{\unicode[STIX]{x1D6FF}}_{0}.\end{eqnarray}$$

Except for a no-slip and impermeable slab (see below), this task requires a numerical procedure. In the present work we resort to a predictor–corrector marching-in-time algorithm (except for the first two time steps for which a second-order Runge–Kutta scheme is used) to solve (4.2). This procedure is benchmarked against the analytical solution available for a no-slip and impermeable slab (i.e. for the case of a usual no-slip solid boundary wall $\unicode[STIX]{x1D6F4}_{u}$ ). Indeed, in such circumstances $k=\unicode[STIX]{x1D6FD}=0$ while $\tilde{\unicode[STIX]{x1D706}}\rightarrow \infty$ so that (recall (3.10) and our discussion after (3.12)) $F=1$ whatever $\overline{\unicode[STIX]{x1D6FF}}$ . Analytically solving (4.2) thus becomes possible and yields the analytical solution $\overline{\unicode[STIX]{x1D6FF}}_{a}=\overline{\unicode[STIX]{x1D6FF}}_{0}\text{e}^{-\sqrt{2/3}\overline{t}}$ .

As shown in table 3, the analytical solution $\overline{\unicode[STIX]{x1D6FF}}_{a}$ is accurately retrieved by the numerics.

Table 3. Analytical and numerical values of normalized sphere–slab gap $\overline{\unicode[STIX]{x1D6FF}}=\unicode[STIX]{x1D6FF}/R$ at different dimensionless times $\bar{t}=1,2,3,5,10$ in the impermeable $(\unicode[STIX]{x1D6FD}=0)$ and no-slip $(\overline{\unicode[STIX]{x1D706}}=\infty )$ case. Here $\overline{\unicode[STIX]{x1D6FF}}_{0}=0.2$ and (4.2) is solved numerically by the marching-in-time algorithm taking for the dimensionless force $F$ either (i) its analytical value $F=1;$ or (ii) its numerical Ramon value (as in table 1); or (iii) its numerical value calculated for the Saffman condition with $\tilde{\unicode[STIX]{x1D706}}=\infty$ .

We compute the sphere trajectory $\overline{\unicode[STIX]{x1D6FF}}(\overline{t})$ for prescribed sphere initial location $\overline{\unicode[STIX]{x1D6FF}}_{0}$ and parameters $(k/R,l/R)$ . Both parameters $\unicode[STIX]{x1D6FD}$ and $\tilde{\unicode[STIX]{x1D706}}$ evolve on the trajectory with, as the sphere approaches the slab, an increase of $\unicode[STIX]{x1D6FD}$ and a decrease of $\tilde{\unicode[STIX]{x1D706}}$ . The product $C=\unicode[STIX]{x1D6FD}\tilde{\unicode[STIX]{x1D706}}^{2}/24=(k/R)/(l/R)^{2}\geqslant 0$ is however constant on each trajectory with a value entirely prescribed by the parameters $(k/R,l/R)$ (i.e. indepently of the initial location $\overline{\unicode[STIX]{x1D6FF}}_{0}$ ).

For a permeable slab (that is, for a positive slab permeance ratio $k/R$ ), note that $\unicode[STIX]{x1D6FD}\rightarrow \infty$ as the sphere approaches the slab. Moreover, for a sufficiently small sphere–slab normalized gap $\overline{\unicode[STIX]{x1D6FF}}$ then $\sqrt{\unicode[STIX]{x1D6FD}}$ becomes large enough for (3.5) to hold with $F\sim 4\unicode[STIX]{x1D6FD}^{-1/2}=(2/3)^{1/2}(k/R)^{-1/2}\overline{\unicode[STIX]{x1D6FF}}$ and the sphere normalized velocity $[\text{d}\overline{\unicode[STIX]{x1D6FF}}/\text{d}\overline{t}]\boldsymbol{e}_{z}$ therefore adopts the time-independent value $-V_{pc}\boldsymbol{e}_{z}$ with $V_{pc}=\sqrt{k/R}$ . As found in § 3.1, the sphere velocity $V_{p}$ tends to $V_{pc}$ for $\sqrt{\overline{k}}\gg \text{Max}(\overline{\unicode[STIX]{x1D6FF}},l/R)$ with $\overline{k}=24k/R$ .

4.2 Settling trajectories

The sphere trajectory sensitivity to the parameters $k/R$ and $l/R$ is illustrated in figure 5 which plots $\overline{\unicode[STIX]{x1D6FF}}$ versus $\overline{t}$ for the initial location $\overline{\unicode[STIX]{x1D6FF}}_{0}=0.2$ and various values of the permeability parameter $k/R=0.001,0.01$ and of the slip parameter $l/R=0,0.01,0.1$ .

Figure 5. Normalized sphere–slab gap $\overline{\unicode[STIX]{x1D6FF}}$ versus dimensionless time $\overline{t}$ for a no-slip and impermeable $(k/R=0)$ slab (dashed curve) and for a permeable slab, with either $k/R=0.001$ (a) or $k/R=0.01$ (b). The permeable slab is either no-slip with $l/R=0$ (*) or slipping with the Saffman boundary condition with normalized slip length $l/R=0.01$ (▫) and $l/R=0.1$ (▵). Here $\overline{\unicode[STIX]{x1D6FF}}_{0}=0.2$ and the fictitious straight trajectory $\unicode[STIX]{x1D6EC}$ (see § 4.1) of a sphere settling at the constant velocity $-V_{pc}\boldsymbol{e}_{z}$ is also displayed (straight dash-dotted line).

As expected, the sphere settles faster towards a porous slab than towards a no-slip impermeable (dashed line) plane wall. In addition, for a given value of $k/R$ the sphere approaches faster the slap as the slip parameter $l/R$ increases. Trajectories for a slipping slab either for $l/R=0.01$ or $l/R=0.1$ differ either weakly or moderately, respectively, from the one prevailing in the no-slip case (see Ramon et al. Reference Ramon, Huppert, Lister and Stone2013). Moreover, these differences decrease in magnitude as $k/R$ increases.

The slopes of trajectories for slip and no-slip porous slabs approach each other as $\overline{\unicode[STIX]{x1D6FF}}$ decreases. This can be observed by tracking on a trajectory for given $l>0$ the slope ratio (see (4.2)) $s_{r}(\overline{\unicode[STIX]{x1D6FF}})=F(\infty ,\unicode[STIX]{x1D6FD})/F(\tilde{\unicode[STIX]{x1D706}},\unicode[STIX]{x1D6FD})>1$ versus the pair $(\tilde{\unicode[STIX]{x1D706}},\unicode[STIX]{x1D6FD})$ which evolves as $\overline{\unicode[STIX]{x1D6FF}}$ drops. This can be done using table 4 and then figure 3.

Table 4. Values of $\overline{k}^{1/2},\unicode[STIX]{x1D6FD}$ and $\tilde{\unicode[STIX]{x1D706}}$ for different parameters $k/R,l/R$ and normalized gap  $\overline{\unicode[STIX]{x1D6FF}}$ . Here $\unicode[STIX]{x1D6FD}_{n}$ and $\tilde{\unicode[STIX]{x1D706}}_{n}$ are obtained for $\overline{\unicode[STIX]{x1D6FF}}=10^{-n}$ with $n=1,2,3$ .

Hence, for the $l/R=0.1$ trajectory in figure 5(a) note that, from table 4, $(\tilde{\unicode[STIX]{x1D706}},\unicode[STIX]{x1D6FD})=(\tilde{\unicode[STIX]{x1D706}}_{1},\unicode[STIX]{x1D6FD}_{1})=(1,2.4)$ for $\overline{\unicode[STIX]{x1D6FF}}=0.1$ and $(\tilde{\unicode[STIX]{x1D706}},\unicode[STIX]{x1D6FD})=(\tilde{\unicode[STIX]{x1D706}}_{2},\unicode[STIX]{x1D6FD}_{2})=(0.1,240)$ for $\overline{\unicode[STIX]{x1D6FF}}=0.01$ . From inspection of figure 3, it thus follows that $s_{t}(\overline{\unicode[STIX]{x1D6FF}}_{1})>s_{t}(\overline{\unicode[STIX]{x1D6FF}}_{2})$ .

As predicted at the end of § 4.1, each trajectory slope tends to the time-independent limit $-V_{pc}$ for $\overline{\unicode[STIX]{x1D6FF}}$ sufficiently small in the following sense: $\sqrt{\overline{k}}\gg \text{Max}(\overline{\unicode[STIX]{x1D6FF}},l/R)$ . This trend is checked by also plotting in figure 5(a,b) the fictitious straight trajectory $\unicode[STIX]{x1D6EC}$ of a sphere that would settle with the previous time-independent velocity $V_{pc}=\sqrt{k/R}$ . For  $\overline{\unicode[STIX]{x1D6FF}}\sim 0.01$ the slopes of $\unicode[STIX]{x1D6EC}$ and the trajectories for a permeable slab become close for $k/R=0.001$ and very close for $k/R=0.01$ , respectively. This is because, as shown in table 4, the condition $\sqrt{\overline{k}}\geqslant \text{Max}(0.01,l/R)$ or $\sqrt{\overline{k}}\gg \text{Max}(0.01,l/R)$ is met for $k/R=0.001$ or $k/R=0.01$ , respectively.