3.1 Introducing the Marangoni concentration $k^{\ast }\equiv k\,Ma$ for small concentrations

The key assumption we propose is that the normalized interfacial surfactant concentration $\unicode[STIX]{x1D6E4}$ is sufficiently small, such that (2.6) and (2.17) can be linearized. The same assumption was made by Harper (Reference Harper2004) for the study of air bubbles rising in contaminated water. Hence, we obtain kinetics congruent with the Henry isotherm (Chang & Franses Reference Chang and Franses1995), namely

(3.1)$$\begin{eqnarray}\displaystyle S(c_{I},\unicode[STIX]{x1D6E4})\approx Bi(kc_{I}-\unicode[STIX]{x1D6E4}), & & \displaystyle\end{eqnarray}$$

(3.2)

In many realistic situations where surfactants are not artificially added, we indeed expect to have low effective bulk concentrations, i.e. $k\ll 1$, which generally lead to small interfacial concentrations $\unicode[STIX]{x1D6E4}$. The interfacial concentration is usually away from saturation, i.e. $\unicode[STIX]{x1D6E4}\ll 1$, because the maximum packing concentration $\hat{\unicode[STIX]{x1D6E4}}_{m}$ used in surfactant models is in fact based on geometrical arguments (Rosen & Kunjappu Reference Rosen and Kunjappu2012) or on achieving a good fit to experimental data based on a specific kinetic model (Chang & Franses Reference Chang and Franses1995). Hence, $\hat{\unicode[STIX]{x1D6E4}}_{m}$ is usually not attained even when the bulk concentration is at the critical micellar concentration. We also have $A\lesssim 1$ for common surfactants (see Prosser & Franses Reference Prosser and Franses2001). We discuss further the relevance of our assumption $\unicode[STIX]{x1D6E4}\ll 1$ in the context of applications in § 7.

We take advantage of the linearization of (2.6) and (2.17) to propose a parameter reduction in our problem, by introducing the following rescaled effective Marangoni concentration and surface concentration:

(3.3a,b)$$\begin{eqnarray}k^{\ast }\equiv Ma\,k,\quad \unicode[STIX]{x1D6E4}^{\ast }\equiv Ma\,\unicode[STIX]{x1D6E4}.\end{eqnarray}$$

Substituting $k=k^{\ast }/Ma$ and $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}^{\ast }/Ma$ into (2.5), (2.14) and (2.13), with the Henry kinetics (3.1) and (3.2), we obtain a set of equations where $k$ and $Ma$ have been combined to form $k^{\ast }$, thereby reducing by one the number of dimensionless groups. This can be verified by examining the updated version of the complete set of equations, (2.2)–(2.13), which becomes

(3.4)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(3.5)$$\begin{eqnarray}\displaystyle & \displaystyle Re\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\boldsymbol{u})=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(3.6)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}c)=\frac{1}{Pe}\unicode[STIX]{x1D6FB}^{2}c, & \displaystyle\end{eqnarray}$$

(3.7)

with boundary conditions

(3.8)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}=\mathbf{0}\quad \text{on all solid surfaces}, & \displaystyle\end{eqnarray}$$

(3.9)$$\begin{eqnarray}\displaystyle & \displaystyle v=0\quad \text{on the interface}, & \displaystyle\end{eqnarray}$$

(3.10)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}=0\quad \text{on all solid surfaces}, & \displaystyle\end{eqnarray}$$

(3.11)

(3.12)$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D712}k^{\ast }}{Pe}\left.\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}\right|_{I}=Bi(k^{\ast }c_{I}-\unicode[STIX]{x1D6E4}^{\ast })\quad \text{on the interface}, & & \displaystyle\end{eqnarray}$$

(3.13)

such that the three quantities $k$, $\unicode[STIX]{x1D6E4}$ and $Ma$ have been replaced by the dimensionless number $k^{\ast }$ and the variable $\unicode[STIX]{x1D6E4}^{\ast }$.

3.2 Scaling theory for surfactant dynamics

To make further progress in modelling the shear rate $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y|_{I}$ and velocity $u_{I}$, we perform a scale analysis on the equations in our problem, starting with rearranging (3.12), which expresses continuity of surfactant fluxes between the bulk and the interface

(3.14)$$\begin{eqnarray}\left.\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}\right|_{I}=\frac{Bi\,Pe}{\unicode[STIX]{x1D712}k^{\ast }}(k^{\ast }c_{I}-\unicode[STIX]{x1D6E4}^{\ast }).\end{eqnarray}$$

For steady flows, adsorption and desorption fluxes between the bulk and the interface are in balance overall, implying

(3.15)$$\begin{eqnarray}\int _{-g/2}^{g/2}\left.\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}\right|_{I}\text{d}x=\frac{Bi\,Pe}{\unicode[STIX]{x1D712}k^{\ast }}\int _{-g/2}^{g/2}(k^{\ast }c_{I}-\unicode[STIX]{x1D6E4}^{\ast })\,\text{d}x=0,\end{eqnarray}$$

such that, by the mean value theorem, there is a point on the interface where $\unicode[STIX]{x2202}c/\unicode[STIX]{x2202}y|_{I}=0$ and $k^{\ast }c_{I}=\unicode[STIX]{x1D6E4}^{\ast }$. With a flow in the positive $x$-direction, interfacial surfactant $\unicode[STIX]{x1D6E4}^{\ast }$ is advected downstream, such that the beginning of the interface has a lower surfactant concentration, implying that $\unicode[STIX]{x1D6E4}^{\ast }<k^{\ast }c_{I}$, and that an adsorption flux exists from the bulk onto the beginning of the interface, such that $\unicode[STIX]{x2202}c/\unicode[STIX]{x2202}y|_{I}>0$ there, as illustrated in figure 2(b). By the same argument, near the end of the interface, a higher surfactant concentration $\unicode[STIX]{x1D6E4}^{\ast }>k^{\ast }c_{I}$ leads to desorption from the interface into the bulk, implying $\unicode[STIX]{x2202}c/\unicode[STIX]{x2202}y|_{I}<0$. Therefore, somewhere along the interface, we must have $\unicode[STIX]{x2202}c/\unicode[STIX]{x2202}y|_{I}=0$. We designate by $x_{0}$ this location where the kinetics flux $S=0$, as depicted in figure 2(b).

In addition, at the beginning of the interface, $c_{I}$ is less than the bulk concentration, i.e. $c_{I}<1$ with our non-dimensionalization, whereas towards the end of the interface, where surfactants accumulate, $c_{I}>1$. This means that, at a specific location along the interface, the concentration near the interface is equal to the background bulk concentration, that is $c_{I}=1$. Taking $c_{I}\sim 1$ along the interface, we then find that (3.14) implies that the interfacial concentration scales as $\unicode[STIX]{x1D6E4}^{\ast }\sim k^{\ast }$.

Next, assuming that the variations of $c_{I}$ and $\unicode[STIX]{x1D6E4}$ scale in the same way for the adsorption region, $-g/2<x<x_{0}$, and the desorption region, $x_{0}<x<g/2$, we have

(3.16a,b)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}^{\ast }\sim k^{\ast }\mp \unicode[STIX]{x0394}\unicode[STIX]{x1D6E4}^{\ast },\quad c_{I}\sim 1\mp \unicode[STIX]{x0394}c_{I}, & & \displaystyle\end{eqnarray}$$

for the adsorption ($-$) and desorption ($+$) regions, respectively (see figure 2b). The quantities $\unicode[STIX]{x0394}\unicode[STIX]{x1D6E4}^{\ast }$ and $\unicode[STIX]{x0394}c_{I}$ are the characteristic variations of $\unicode[STIX]{x1D6E4}^{\ast }$ and $c_{I}$, respectively. We must have $\unicode[STIX]{x0394}\unicode[STIX]{x1D6E4}^{\ast }>k^{\ast }\unicode[STIX]{x0394}c_{I}>0$ to satisfy the direction of the kinetic flux, as described above.

From the relation between Marangoni stress and surfactant gradient (3.2), we also have

(3.17)$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D6E4}^{\ast }\sim g\unicode[STIX]{x1D6FE}_{Ma},\end{eqnarray}$$

where

(3.18)

is the average shear rate induced by Marangoni stresses along the interface, such that $\unicode[STIX]{x1D6FE}_{Ma}=0$ corresponds to free slip at the interface and $\unicode[STIX]{x1D6FE}_{Ma}=1$ corresponds to a no-slip interface. Then, a scale analysis of (3.14) gives

(3.19)$$\begin{eqnarray}\left.\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}\right|_{I}\sim \frac{\unicode[STIX]{x0394}c_{I}}{\unicode[STIX]{x1D6FF}}\sim \frac{Bi\,Pe}{\unicode[STIX]{x1D712}k^{\ast }}(g\unicode[STIX]{x1D6FE}_{Ma}-k^{\ast }\unicode[STIX]{x0394}c_{I}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}$ is the typical thickness of the diffusive layer of bulk surfactant. To estimate $\unicode[STIX]{x1D6FF}$, we can use the bulk advection–diffusion equation (2.4). At high Péclet numbers, $Pe\gg 1$, the diffusive layer of surfactant forms a thin boundary layer. As explained in detail in appendix B, there are two main asymptotic regimes depending on whether there is slip or not at the interface. For large slip and small interfacial shear rate, $\unicode[STIX]{x1D6FE}_{Ma}\ll 1$, we can show that the boundary layer thickness scales as (see appendix B)

(3.20)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D6FF}}{g}=\unicode[STIX]{x1D6FF}_{0,1}(1+\unicode[STIX]{x1D6FF}_{1,1}g^{2}Pe)^{-1/2}\quad \text{for }g\lesssim 1, & \displaystyle\end{eqnarray}$$

(3.21)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D6FF}}{g}=\unicode[STIX]{x1D6FF}_{0,2}(1+\unicode[STIX]{x1D6FF}_{1,2}gPe)^{-1/2}\quad \text{for }g\gtrsim 1, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}_{0,1}$, $\unicode[STIX]{x1D6FF}_{1,1}$, $\unicode[STIX]{x1D6FF}_{0,2}$, $\unicode[STIX]{x1D6FF}_{1,2}$ are empirical parameters which need to be determined. We note that the scaling $\unicode[STIX]{x1D6FF}\sim Pe^{-1/2}$ at large Péclet numbers corresponds to having a uniform velocity in the diffusive boundary layer, consistently with the case $\unicode[STIX]{x1D6FE}_{Ma}\ll 1$.

For negligible slip at the interface and $\unicode[STIX]{x1D6FE}_{Ma}\sim 1$, we obtain

(3.22)$$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FF}}{g}=\unicode[STIX]{x1D6FF}_{0,3}(1+\unicode[STIX]{x1D6FF}_{1,3}g^{2}Pe)^{-1/3},\end{eqnarray}$$

for any $g>0$, and with $\unicode[STIX]{x1D6FF}_{0,3}$ and $\unicode[STIX]{x1D6FF}_{1,3}$ two empirical parameters which need to be determined. This corresponds to the Lévêque regime (Lévêque Reference Lévêque1928; Landel et al. Reference Landel, Thomas, McEvoy and Dalziel2016), giving a power law $\unicode[STIX]{x1D6FF}\sim Pe^{-1/3}$ at large Péclet numbers owing to a linear shear rate profile in the diffusive boundary layer. The scalings (3.20)–(3.22) assume that: (i) the variation of the bulk concentration along the interface is sufficiently smooth; (ii) the boundary layer is not confined vertically, i.e. $\unicode[STIX]{x1D6FF}\lesssim 1$; and (iii) the diffusive boundary layers between consecutive interfaces are independent. As we will discuss in § 6, our scaling prediction remains accurate even for confined diffusive boundary layers $\unicode[STIX]{x1D6FF}\sim 1$.

With $\unicode[STIX]{x1D6FF}$ assumed known in terms of $g$ and $Pe$, we rearrange (3.19) to solve for $\unicode[STIX]{x0394}c_{I}$,

(3.23)$$\begin{eqnarray}\unicode[STIX]{x0394}c_{I}\sim \unicode[STIX]{x1D6FE}_{Ma}\frac{\displaystyle \frac{Bi\,Pe}{\unicode[STIX]{x1D712}k^{\ast }}\,g\,\unicode[STIX]{x1D6FF}}{1+\displaystyle \frac{Bi\,Pe}{\unicode[STIX]{x1D712}}\unicode[STIX]{x1D6FF}},\end{eqnarray}$$

such that, dividing by $\unicode[STIX]{x1D6FF}$, we obtain a scaling that relates the kinetics flux to the shear

(3.24)$$\begin{eqnarray}\left.\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}\right|_{I}\sim \frac{\unicode[STIX]{x0394}c_{I}}{\unicode[STIX]{x1D6FF}}\sim \unicode[STIX]{x1D6FE}_{Ma}\frac{\displaystyle \frac{Bi\,Pe}{\unicode[STIX]{x1D712}k^{\ast }}g}{1+\displaystyle \frac{Bi\,Pe}{\unicode[STIX]{x1D712}}\unicode[STIX]{x1D6FF}}.\end{eqnarray}$$

3.3 Scaling for the interfacial velocity and for the slip length

We now seek a scaling expression for $u_{I}$. We integrate the interfacial advection–diffusion equation (3.7) from the upstream stagnation point $x=-g/2$ to $x_{0}$. We find

(3.25)

where we used the no-slip boundary condition (3.8) at $x=-g/2$ for the left-hand side, the no-flux boundary condition (3.11) at $x=-g/2$ for the first term on the right-hand side, as well as the continuity of flux condition (3.12) for the last term. To write the right-hand side in terms of $\unicode[STIX]{x1D6FE}_{Ma}$, note that $\unicode[STIX]{x1D6E4}^{\ast }|_{x_{0}}\sim k^{\ast }$ and $\text{d}\unicode[STIX]{x1D6E4}^{\ast }/\text{d}x|_{x_{0}}\sim \unicode[STIX]{x1D6FE}_{Ma}$. For the last term, we use (3.24) to scale the integral

(3.26)

Substituting into (

) we obtain a scaling relation between interfacial velocity and shear. Introducing empirical prefactors (to be determined) ahead of each term, we write

(3.27)$$\begin{eqnarray}u_{I}|_{x_{0}}=\frac{2}{a_{1}}\frac{1}{k^{\ast }}\left(\frac{1}{Pe_{I}}+a_{2}\frac{g^{2}Bi}{1+\displaystyle \frac{BiPe}{\unicode[STIX]{x1D712}}\unicode[STIX]{x1D6FF}}\right)\unicode[STIX]{x1D6FE}_{Ma},\end{eqnarray}$$

where $a_{1},a_{2}$ are empirical parameters; the choice of writing $2/a_{1}$ for the overall prefactor leads to a more convenient expression for the results later in § 6.1.

As also noted in the previous section, scaling expressions for the boundary layer thickness $\unicode[STIX]{x1D6FF}$ are given by (3.20), (3.21) or (3.22), which depend on $g$, $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D6FE}_{Ma}$. Therefore, our scaling is a nonlinear function of the Marangoni shear rate. However, in the comparison of our model with numerical simulations (see § 6), we find that the nonlinear dependence of $\unicode[STIX]{x1D6FE}_{Ma}$ with $\unicode[STIX]{x1D6FF}$ is actually weak. Consequently, we adopt only (3.22) in our model. This proves to be a good approximation and allows us to regard $\unicode[STIX]{x1D6FF}$ as independent from $\unicode[STIX]{x1D6FE}_{Ma}$.

Furthermore, we note that a first-order linear expansion of the concentrations $c$ and $\unicode[STIX]{x1D6E4}$ near $x_{0}$ predicts $a_{2}=1/8$ since $x_{0}=0$ (that is, $x_{0}$ is at the mid-gap location) due to the balance of desorption and adsorption fluxes along the interface.

A characteristic scale for the slip length near $x_{0}$, which corresponds to the mid-gap of the interface under our assumptions, is therefore simply

(3.28)$$\begin{eqnarray}\unicode[STIX]{x1D706}_{x_{0}}=\frac{u_{I}|_{x_{0}}}{\unicode[STIX]{x1D6FE}_{Ma}}=\frac{2}{a_{1}}\frac{1}{k^{\ast }}\left(\frac{1}{Pe_{I}}+a_{2}\frac{g^{2}Bi}{1+\displaystyle \frac{BiPe}{\unicode[STIX]{x1D712}}\unicode[STIX]{x1D6FF}}\right).\end{eqnarray}$$

This scaling prediction shows that the local slip length $\unicode[STIX]{x1D706}_{x_{0}}$ depends strongly on the Marangoni concentration $k^{\ast }=k\,Ma$ and the normalized gap length $g$. It is intuitive that increasing the gap length tends to increase the slip length, since it would reduce the concentration gradient at the interface and thus the opposing Marangoni stress. In contrast, increasing the effective bulk surfactant concentration $k$ or the Marangoni number tends to reduce the slip length, as expected. We also find increasing the bulk or interfacial Péclet numbers, $Pe$ or $Pe_{I}$, reduces $\unicode[STIX]{x1D706}_{x_{0}}$. Increasing the Biot or $\unicode[STIX]{x1D712}$ numbers has a positive effect on the slip length.

However, we note that (3.28) is only a local measure of the characteristic slip length near the middle of the interface, where $S(x_{0})=0$. In order to have an effective or global slip length over the entire SHS which takes into account all interfaces and solid ridges, we also need to model the channel flow over the SHS. In the next section, we analyse the remaining governing equations for the flow, i.e. the continuity and Navier–Stokes equations (3.4) and (3.5), to study how the flow is affected by a SHS with a surfactant-induced Marangoni stress over the interfaces.

4.1 Stokes flow model for SHS channels with surfactant contamination

According to equation (3.2), interfacial surfactant concentration gradients can generate a Marangoni shear rate at the interface $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y|_{I}\approx \text{d}\unicode[STIX]{x1D6E4}^{\ast }/\text{d}x\geqslant 0$. In this section, we derive an expression for how interfacial stresses with arbitrary profile can affect the flow over a periodic SHS. The geometry follows the same schematic presented in figure 2. Such a periodic SHS arrangement was studied in detail by Lauga & Stone (Reference Lauga and Stone2003) and Teo & Khoo (Reference Teo and Khoo2009) for a shear-free interface, i.e. with $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y|_{I}=0$ along the interface, at low Reynolds number. Here we generalize their approach to also study the case where $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y|_{I}\geqslant 0$. We also assume $Re\ll 1$, such that (3.5) simplifies to the Stokes flow equation

(4.1)$$\begin{eqnarray}\unicode[STIX]{x1D735}p=\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}.\end{eqnarray}$$

Taking the curl of (4.1) and using the continuity equation (3.4), we find that the pressure field $p$ and the vorticity field, $\unicode[STIX]{x1D74E}=\unicode[STIX]{x1D735}\times \boldsymbol{u}$ are both solutions of Laplace’s equation. Using the superposition principle to solve Laplace’s equation for the vorticity, we decompose it as the sum of the two-dimensional Poiseuille flow component, which is a pressure-driven flow in a channel with full solid walls on both sides (denoted by a subscript $p$), and a deviating component (denoted by a subscript $d$), such that

(4.2)$$\begin{eqnarray}\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{p}+\unicode[STIX]{x1D714}_{d},\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{p}=y$. As the flow is incompressible, we can also use the streamfunction $\unicode[STIX]{x1D733}$, defined such that $\boldsymbol{u}=\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D733}$, and which is the solution of the biharmonic equation $\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D733}=\mathbf{0}$. Note that $\unicode[STIX]{x1D733}=(0,0,\unicode[STIX]{x1D6F9})$ for two-dimensional flows. The solution for the deviating component of the vorticity is obtained using separation of variables considering the periodicity of the flow with wavelength $L$. Integrating twice, we then obtain the deviation streamfunction (Lauga & Stone Reference Lauga and Stone2003; Teo & Khoo Reference Teo and Khoo2009). Noting that the mean pressure gradient imposed by the deviating field is zero, neglecting the constant of integration and using the no-flow boundary condition $v=0$ in (3.8) and (3.9), the deviating component of the streamfunction is

(4.3)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F9}_{d} & = & \displaystyle B\frac{y^{2}}{2}+Ey+\mathop{\sum }_{n=1}^{\infty }\{e_{n}[\cosh (k_{n}y)-\coth (k_{n})y\sinh (k_{n}y)]\nonumber\\ \displaystyle & & \displaystyle +\,d_{n}[\sinh (k_{n}y)-\tanh (k_{n})y\cosh (k_{n}y)]\}\cos (k_{n}x),\end{eqnarray}$$

where $k_{n}=2\unicode[STIX]{x03C0}n/L$, and $B$, $E$, $e_{n}$ and $d_{n}$ are unknowns to be determined using the other boundary conditions. The streamfunction for the Poiseuille component is

(4.4)$$\begin{eqnarray}\unicode[STIX]{x1D6F9}_{p}={\textstyle \frac{1}{6}}(3y-y^{3}).\end{eqnarray}$$

Up to this point, (4.3) and (4.4) are general solutions for any arrangement and geometry of a SHS in a two-dimensional Stokes flow channel: i.e. they are not limited to one-sided SHSs, symmetric patterns or a particular shear rate profile at the interface.

With our geometry, using the no-slip boundary condition on the solid wall side at $y=1$ for all $x$, we find $B=-E$ and

(4.5)$$\begin{eqnarray}g_{n}=\frac{e_{n}}{d_{n}}=-\frac{\sinh (k_{n})-k_{n}\cosh (k_{n})+k_{n}\tanh (k_{n})\sinh (k_{n})}{\cosh (k_{n})-k_{n}\sinh (k_{n})+k_{n}\coth (k_{n})\cosh (k_{n})}.\end{eqnarray}$$

Hence, the deviating streamfunction simplifies to

(4.6)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F9}_{d} & = & \displaystyle \left(-\frac{y^{2}}{2}+y\right)E+\mathop{\sum }_{n=1}^{\infty }d_{n}\{g_{n}[\cosh (k_{n}y)-\coth (k_{n})y\sinh (k_{n}y)]\nonumber\\ \displaystyle & & \displaystyle +\,\sinh (k_{n}y)-\tanh (k_{n})y\cosh (k_{n}y)\}\cos (k_{n}x).\end{eqnarray}$$

To determine the unknowns $E$ and $d_{n}$ for $n\geqslant 1$, we can use the no-slip boundary condition on the SHS side. At $y=-1$ for $g/2<|x|<L/2$, we have the condition

(4.7)$$\begin{eqnarray}0=2E+\mathop{\sum }_{n=1}^{\infty }d_{n}\unicode[STIX]{x1D6FC}_{n}\cos (k_{n}x),\end{eqnarray}$$

where

(4.8)$$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{n}=2k_{n}[\cosh (k_{n})-\tanh (k_{n})\sinh (k_{n})]-2\sinh (k_{n}).\end{eqnarray}$$

We then apply the last boundary condition on the interface, where we assume that there is an arbitrary shear rate profile $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y|_{I}(x)\geqslant 0$. Hence, we obtain the general condition

(4.9)$$\begin{eqnarray}0=\left.\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}\right|_{I}-1+E+\mathop{\sum }_{n=1}^{\infty }d_{n}\unicode[STIX]{x1D6FD}_{n}\cos (k_{n}x)\end{eqnarray}$$

for $|x|<g/2$, $y=-1$, and with

(4.10)$$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{n}=2k_{n}[g_{n}\coth (k_{n})\cosh (k_{n})-\tanh (k_{n})\sinh (k_{n})].\end{eqnarray}$$

To make further progress and obtain a relationship between the interfacial shear rate and the interfacial velocity, we now assume that the interfacial shear rate is uniform along the interface: $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y|_{I}=\unicode[STIX]{x1D6FE}_{Ma}$, where $0\leqslant \unicode[STIX]{x1D6FE}_{Ma}\leqslant 1$ corresponds to the interface-averaged surfactant-induced Marangoni shear rate, as introduced previously in (3.17). This assumption is consistent with having a uniform concentration gradient, following the linearized coupling condition (3.2). In the context of air bubbles rising in surfactant-contaminated water, this assumption is also consistent with the ‘uniformly retarded’ regime described for instance by Palaparthi et al. (Reference Palaparthi, Papageorgiou and Maldarelli2006). In the context of a SHS, a similar assumption was made by Schönecker & Hardt (Reference Schönecker and Hardt2013) to model viscous effects from a gas phase trapped inside the cavities of the SHS. This allowed them to decouple the flow above the interface from the flow in the cavity of the SHS. We discuss further the relevance of this assumption in applications in § 7. Hence, (4.9) becomes

(4.11)$$\begin{eqnarray}0=\unicode[STIX]{x1D6FE}_{Ma}-1+E+\mathop{\sum }_{n=1}^{\infty }d_{n}\unicode[STIX]{x1D6FD}_{n}\cos (k_{n}x).\end{eqnarray}$$

If $\unicode[STIX]{x1D6FE}_{Ma}=0$ in the equation above, the interface is stress free and the surfactant concentration gradient at the interface vanishes. The surface is completely immobilized if $\unicode[STIX]{x1D6FE}_{Ma}=1$, and the flow follows a channel Poiseuille flow.

Following Lauga & Stone (Reference Lauga and Stone2003), we can compute an approximation of the solution by truncating the series in equations (4.7) and (4.11) at $n=N-1$, multiplying (4.7) and (4.11) by $\cos (2\unicode[STIX]{x03C0}mr)$ for $m\in [0,N-1]$ (with $r=x/L$) and integrating them for $r\in (\unicode[STIX]{x1D719}/2,1/2)$ and $r\in (0,\unicode[STIX]{x1D719}/2)$, respectively, where $\unicode[STIX]{x1D719}$ is the gas fraction. Summing together the results for each $m$ in one single equation, we finally obtain a linear system of $N$ equations for the $N$ unknown coefficients $E$ and $d_{n}$ for $n\in [1,N-1]$, which we can solve numerically. The linear system in matrix form is, for $m\in [0,N-1]$ and $n\in [0,N-1]$,

(4.12)$$\begin{eqnarray}\unicode[STIX]{x1D63C}_{m,n}\boldsymbol{U}_{n}=\boldsymbol{B}_{m},\end{eqnarray}$$

with $U_{0}=E$ and $U_{n}=d_{n}$. The square matrix $\unicode[STIX]{x1D63C}_{m,n}$ has coefficients

(4.13)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{0,0}=1-\frac{\unicode[STIX]{x1D719}}{2}, & \displaystyle\end{eqnarray}$$

(4.14)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{0,n}=(\unicode[STIX]{x1D6FD}_{n}-\unicode[STIX]{x1D6FC}_{n})\frac{\sin (\unicode[STIX]{x03C0}n\unicode[STIX]{x1D719})}{2\unicode[STIX]{x03C0}n},\quad n>0 & \displaystyle\end{eqnarray}$$

(4.15)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{m,0}=-\frac{\sin (\unicode[STIX]{x03C0}m\unicode[STIX]{x1D719})}{2\unicode[STIX]{x03C0}m},\quad m>0 & \displaystyle\end{eqnarray}$$

(4.16)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{n,n}=\frac{\unicode[STIX]{x1D6FC}_{n}}{4}+(\unicode[STIX]{x1D6FD}_{n}-\unicode[STIX]{x1D6FC}_{n})\left(\frac{\unicode[STIX]{x1D719}}{4}+\frac{\sin (2\unicode[STIX]{x03C0}n\unicode[STIX]{x1D719})}{8\unicode[STIX]{x03C0}n}\right),\quad n>0 & \displaystyle\end{eqnarray}$$

(4.17)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{m,n}=(\unicode[STIX]{x1D6FD}_{n}-\unicode[STIX]{x1D6FC}_{n})\frac{1}{4\unicode[STIX]{x03C0}}\left(\frac{\sin (\unicode[STIX]{x03C0}(m+n)\unicode[STIX]{x1D719})}{m+n}+\frac{\sin (\unicode[STIX]{x03C0}(m-n)\unicode[STIX]{x1D719})}{m-n}\right),\quad m\neq n>0\quad & \displaystyle\end{eqnarray}$$

and the vector $\boldsymbol{B}_{m}$ has coefficients

(4.18)$$\begin{eqnarray}\displaystyle & \displaystyle B_{0}=(1-\unicode[STIX]{x1D6FE}_{Ma})\frac{\unicode[STIX]{x1D719}}{2}, & \displaystyle\end{eqnarray}$$

(4.19)$$\begin{eqnarray}\displaystyle & \displaystyle B_{m}=(1-\unicode[STIX]{x1D6FE}_{Ma})\frac{\sin (\unicode[STIX]{x03C0}m\unicode[STIX]{x1D719})}{2\unicode[STIX]{x03C0}m},\quad m>0. & \displaystyle\end{eqnarray}$$

Care must be taken at large $n$, where the system is not well conditioned, as pointed out by Teo & Khoo (Reference Teo and Khoo2009). We provide, as supplementary material, MATLAB routines solving the linear system (4.12).

4.2 Interfacial slip velocity

Once all the coefficients $E$ and $d_{n}$ are computed, the non-dimensional slip velocity at the interface $u_{I}$ can be determined to machine precision, depending on the size $N$ of the matrix $\unicode[STIX]{x1D63C}$, such that

(4.20)$$\begin{eqnarray}u_{I}=2E+\mathop{\sum }_{n=1}^{\infty }d_{n}\unicode[STIX]{x1D6FC}_{n}\cos (k_{n}x).\end{eqnarray}$$

Through its coefficients $E$ and $d_{n}$, (4.20) is a function of the uniform Marangoni interfacial shear rate $\unicode[STIX]{x1D6FE}_{Ma}$ and of the two non-dimensional geometrical parameters $g$ and $\unicode[STIX]{x1D719}$. Hence, we have

(4.21)$$\begin{eqnarray}u_{I}={\mathcal{G}}(\unicode[STIX]{x1D6FE}_{Ma},g,\unicode[STIX]{x1D719},x).\end{eqnarray}$$

The function ${\mathcal{G}}$ is only known implicitly through the solution of the linear system (4.12). In practice, it would be useful to obtain an explicit analytical solution, or at least a scaling expression for ${\mathcal{G}}$ which can give an approximate solution to the coupled surfactant–flow transport problem in combination with (3.27). In the linear system (4.12), we can factorize all the coefficients of $\boldsymbol{B}_{m}$ by $(1-\unicode[STIX]{x1D6FE}_{Ma})$. This means that $E$ and $d_{n}$ are proportional to $(1-\unicode[STIX]{x1D6FE}_{Ma})$ for all $n\geqslant 1$. Thus, the velocity at the interface is such that

(4.22)$$\begin{eqnarray}u_{I}=2(1-\unicode[STIX]{x1D6FE}_{Ma}){\mathcal{F}}(g,\unicode[STIX]{x1D719},x),\end{eqnarray}$$

where, again, ${\mathcal{F}}$ is an implicit function. Now, ${\mathcal{F}}$ is decoupled from the surfactant transport problem since it does not depend on $\unicode[STIX]{x1D6FE}_{Ma}$. It can thus be computed to arbitrary numerical precision for each couple of geometrical non-dimensional parameters $(g,\unicode[STIX]{x1D719})$ and for all $x$ by solving the linear system (4.12) in the surfactant-free case, i.e. setting $\unicode[STIX]{x1D6FE}_{Ma}=0$ in (4.18) and (4.19).

Figure 3.

Variation of the normalized mid-gap interfacial velocity, $u_{Ic}/(2(1-\unicode[STIX]{x1D6FE}_{Ma}))={\mathcal{F}}(g,\unicode[STIX]{x1D719},x=0)={\mathcal{F}}_{0}(g,\unicode[STIX]{x1D719})$ (see (4.22) and text), as a function of: $(a)$ the non-dimensional interfacial length $g$, for different non-dimensional gas fractions from $\unicode[STIX]{x1D719}=0.01$ to $0.99$ (shown with different colours, see legend); $(b)$ the gas fraction $\unicode[STIX]{x1D719}$, for various interfacial lengths $g$ (shown with different colours, see legend). The implicit function ${\mathcal{F}}_{0}(g,\unicode[STIX]{x1D719})$ has been computed by solving the linear system (4.12) with $N=500$, except in the more demanding cases of $0.01<\unicode[STIX]{x1D719}\leqslant 0.1$ ($N=2\,500$), $0.99\leqslant \unicode[STIX]{x1D719}\leqslant 1$ ($N=2\,500$) and $0\leqslant \unicode[STIX]{x1D719}\leqslant 0.01$ ($N=15\,000$). In (a) the black dotted line is plotted using (4.24) for $\unicode[STIX]{x1D719}\ll 1$ and $g\lesssim 1$, the black dot-dashed line is plotted using (C 24) for $g\lesssim 1$ and $\unicode[STIX]{x1D719}=0.99$ and the black dashed line is plotted using (4.25) for $\unicode[STIX]{x1D719}\ll 1$ and $g\gtrsim 1$. The black solid line corresponds to the maximum asymptotic value for $\unicode[STIX]{x1D719}\rightarrow 1$: $u_{I,c}\rightarrow u_{u}(y=-1)=2(1-\unicode[STIX]{x1D6FE}_{Ma})$ (see (4.27)). In (b), the black dashed line has been plotted using the asymptotic trend (4.28) for $g\gg 1$ (see also appendix C).

Based on this observation, ${\mathcal{F}}(g,\unicode[STIX]{x1D719},x)=u_{I}(x)/(2(1-\unicode[STIX]{x1D6FE}_{Ma}))$ is a normalized interfacial velocity. In figure 3(a), we plot on a log–log scale this normalized interfacial velocity at the middle of the gap, $x=0$:

(4.23)$$\begin{eqnarray}\frac{u_{Ic}}{2(1-\unicode[STIX]{x1D6FE}_{Ma})}={\mathcal{F}}(g,\unicode[STIX]{x1D719},x=0)={\mathcal{F}}_{0}(g,\unicode[STIX]{x1D719}),\end{eqnarray}$$

as a function of $g$ and for different $\unicode[STIX]{x1D719}$ (shown with different colours, see legend). For $\unicode[STIX]{x1D719}\ll 1$ and $g\lesssim 1$, the normalized interfacial velocity follows a linear asymptotic trend

(4.24)$$\begin{eqnarray}\frac{u_{Ic}}{2(1-\unicode[STIX]{x1D6FE}_{Ma})}\simeq \frac{g}{8},\end{eqnarray}$$

plotted with a black dotted line in figure 3(a). We can see that for $\unicode[STIX]{x1D719}=0.99$ the interfacial velocity still follows a linear scaling, although with a higher slope than in the asymptotic limit (4.24), as shown by the black dot-dashed line in figure 3(a), which was computed using (C 24). At large gap length, $g\gg 1$, and for low gas fraction, $\unicode[STIX]{x1D719}\ll 1$, the interfacial velocity collapses on the asymptotic plateau

(4.25)$$\begin{eqnarray}\frac{u_{Ic}}{2(1-\unicode[STIX]{x1D6FE}_{Ma})}\rightarrow \frac{1}{4},\end{eqnarray}$$

plotted with a black dashed line in figure 3(a). More details about the behaviour of the interfacial velocity $u_{Ic}$ with $\unicode[STIX]{x1D719}$ and $g$ and the two asymptotic limits (4.24) and (4.25) can be found in appendix C. The transition observed at $g\sim 1$ from a linear trend towards a plateau is due to the importance of the opposite wall at $y=1$ through viscous effects.

We note that the behaviour of $u_{I,c}/(2(1-\unicode[STIX]{x1D6FE}_{Ma}))$ is similar across all $g$ and for any $\unicode[STIX]{x1D719}$. This function goes from a linear behaviour for $g\lesssim 1$ to a plateau for $g\gtrsim 1$, and with simple asymptotics in the case $\unicode[STIX]{x1D719}\ll 1$. Most of the data in figure 3(a) follow these limiting regimes, suggesting that asymptotic results are sufficiently accurate in many applications.

This common behaviour of the interfacial velocity might also suggest that the velocity field follows a closed analytical form. However, we have not been able to demonstrate this theoretically from the biharmonic equation. As far as we are aware, the case of Stokes flow in a transverse channel with mixed boundary conditions changing twice (on one or both channel sides), which is reminiscent of the longitudinal-channel work of Philip (Reference Philip1972a), has not been shown to have a closed analytical form in the literature. It would be valuable to re-examine the present problem with conformal mapping tools similar to those used by Crowdy (Reference Crowdy2016, Reference Crowdy2017a).

Figure 3(b) plots curves of $u_{I,c}/(2(1-\unicode[STIX]{x1D6FE}_{Ma}))$ versus gas fraction $\unicode[STIX]{x1D719}$, with $g$ as a parameter. As the gas fraction $\unicode[STIX]{x1D719}$ increases towards $1$, the normalized interfacial velocity increases rapidly at any fixed $g$. In the limit $\unicode[STIX]{x1D719}\rightarrow 1$ we have

(4.26)$$\begin{eqnarray}\frac{u_{Ic}}{2(1-\unicode[STIX]{x1D6FE}_{Ma})}\rightarrow 1,\end{eqnarray}$$

which can be predicted from the velocity field with uniform boundary conditions at the top and bottom sides, that is $u(y=1)=0$ and $\text{d}u_{u}/\text{d}y(y=-1)=\unicode[STIX]{x1D6FE}_{Ma}$, respectively. The solution to the Stokes problem (4.1) with these uniform boundary conditions is independent of $x$,

(4.27)$$\begin{eqnarray}u_{u}={\textstyle \frac{1}{2}}(1-y^{2})+(1-\unicode[STIX]{x1D6FE}_{Ma})(1-y),\end{eqnarray}$$

and can be used to yield the limit of $u_{Ic}\rightarrow u_{u}(y=-1)$ for $\unicode[STIX]{x1D719}\rightarrow 1$. At $g\lesssim 1$, it is also clear from figure 3(b) that $u_{Ic}/(2(1-\unicode[STIX]{x1D6FE}_{Ma}))\rightarrow 1$ only for gas fractions very close to 1, i.e. in the limit $\unicode[STIX]{x1D719}\rightarrow 1$, as already observed in figure 3(a). This result confirms the range of validity of the first scaling (3.20) for the diffusive boundary layer thickness $\unicode[STIX]{x1D6FF}$. Then, we can show (see appendix C) that in the limit of large gap length, $g\gg 1$, the normalized interfacial velocity follows the asymptotic hyperbolic trend

(4.28)$$\begin{eqnarray}\frac{u_{Ic}}{2(1-\unicode[STIX]{x1D6FE}_{Ma})}\simeq \frac{1}{4-3\unicode[STIX]{x1D719}},\end{eqnarray}$$

plotted with a black dashed line in figure 3(b). The asymptotic result (4.28) is valid for any $\unicode[STIX]{x1D719}$. This result is consistent with (4.25) and (4.26).

4.3 Predictions of the interfacial shear rate, effective slip length and drag reduction

We now have two independent expressions relating the interfacial velocity $u_{I}$ and the Marangoni shear $\unicode[STIX]{x1D6FE}_{Ma}$. The scaling (3.27) was found based on the near-interface surfactant dynamics, whereas (4.22) was derived from a Stokes flow solution. Eliminating the interface velocity, we deduce a scaling expression for the average Marangoni shear rate,

(4.29)$$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{Ma}=a_{1}k^{\ast }{\mathcal{F}}_{0}(g,\unicode[STIX]{x1D719})\left(\frac{1}{Pe_{I}}+a_{2}\frac{g^{2}Bi}{1+{\displaystyle \frac{Bi\,Pe}{\unicode[STIX]{x1D712}}}\unicode[STIX]{x1D6FF}}+a_{1}k^{\ast }{\mathcal{F}}_{0}(g,\unicode[STIX]{x1D719})\right)^{-1},\end{eqnarray}$$

where $a_{1}$ and $a_{2}$ are the empirical parameters that were introduced in § 3.2. This predictive scaling depends only on the properties of the flow, fluid and surfactant through the non-dimensional numbers $k^{\ast }=k\,Ma$, $Pe_{I}$, $Bi$, $Pe$ and $\unicode[STIX]{x1D712}$, and on the two geometrical parameters $g$ and $\unicode[STIX]{x1D719}$. As noted earlier, this assumes a sufficiently small concentration of surfactant and a small Reynolds number in the flow, and the diffusive boundary layer thickness $\unicode[STIX]{x1D6FF}$ depends only weakly on $\unicode[STIX]{x1D6FE}_{Ma}$ following (3.20), (3.21) or (3.22). The parameters $a_{1}$, $a_{2}$, as well as $\unicode[STIX]{x1D6FF}_{0,i}$ and $\unicode[STIX]{x1D6FF}_{1,i}$ (with $i=1$, 2 or 3 for the scaling predictions (3.20), (3.21) or (3.22), respectively) for $\unicode[STIX]{x1D6FF}$, are determined empirically by fitting to our numerical simulations in § 6.

We can also compute a global effective slip length $\unicode[STIX]{x1D706}_{e}$ as defined by Lauga & Stone (Reference Lauga and Stone2003), which corresponds to the value $\unicode[STIX]{x1D706}_{e}$ such that an equivalent channel flow under the same pressure gradient, but with a uniform Navier slip boundary condition $u(y=-1)=\unicode[STIX]{x1D706}_{e}\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y|_{y=-1}$ replaces the mixed conditions of the SHS at the bottom boundary. We can show that the contribution of the effective slip length $\unicode[STIX]{x1D706}_{e}$ is such that the total volume flux in the channel is the sum of the background Poiseuille volume flux, $Q_{p}=2/3$, and the volume flux of the deviating flow,

(4.30)$$\begin{eqnarray}Q=Q_{p}+Q_{d}=\frac{2}{3}+\frac{2\unicode[STIX]{x1D706}_{e}}{\unicode[STIX]{x1D706}_{e}+2},\end{eqnarray}$$

where the maximum value for the deviating flux is $Q_{d}\rightarrow 2$, as $\unicode[STIX]{x1D706}_{e}\rightarrow \infty$. The effective slip length as a function of the deviating flux is

(4.31)$$\begin{eqnarray}\unicode[STIX]{x1D706}_{e}=\frac{2Q_{d}}{2-Q_{d}}.\end{eqnarray}$$

From equation (4.6), the deviating streamfunction is

(4.32)$$\begin{eqnarray}Q_{d}=\unicode[STIX]{x1D6F9}_{d}(y=1)-\unicode[STIX]{x1D6F9}_{d}(y=-1)=\frac{E}{2}-\left(-\frac{3E}{2}\right)=2E,\end{eqnarray}$$

and substituting into (4.31) yields

(4.33)$$\begin{eqnarray}\unicode[STIX]{x1D706}_{e}=\frac{2E}{1-E}.\end{eqnarray}$$

Following from the linearity of the governing equations, and of the boundary conditions, $E$ also scales linearly with $(1-\unicode[STIX]{x1D6FE}_{Ma})$. Accordingly, we can find the explicit dependence of $\unicode[STIX]{x1D706}_{e}$ with the Marangoni shear rate $\unicode[STIX]{x1D6FE}_{Ma}$,

(4.34)$$\begin{eqnarray}\unicode[STIX]{x1D706}_{e}=\frac{2(1-\unicode[STIX]{x1D6FE}_{Ma})E_{0}}{1-(1-\unicode[STIX]{x1D6FE}_{Ma})E_{0}},\end{eqnarray}$$

where $E_{0}$ is the first coefficient of the vector $\boldsymbol{U}_{n}$ (see (4.12)) in the surfactant-free case, i.e. $E_{0}=E$ for $\unicode[STIX]{x1D6FE}_{Ma}=0$, and $\unicode[STIX]{x1D6FE}_{Ma}$ is expressed by (4.29). As expected, $0\leqslant \unicode[STIX]{x1D706}_{e}<\infty$, since $0\leqslant \unicode[STIX]{x1D6FE}_{Ma}\leqslant 1$ and $0\leqslant E_{0}\leqslant 1$.

The corresponding drag reduction due to the presence of the SHS in our pressure-driven channel flows, inclusive of surfactant, can be computed as

(4.35)$$\begin{eqnarray}DR=1-\frac{C_{f}}{C_{f,p}}=1-\frac{\displaystyle \frac{\langle \hat{\unicode[STIX]{x1D70F}}_{s}\rangle }{\hat{\unicode[STIX]{x1D70C}}(\hat{Q}/2)^{2}}}{\displaystyle \frac{\langle \hat{\unicode[STIX]{x1D70F}}_{s}\rangle _{p}}{\hat{\unicode[STIX]{x1D70C}}(\hat{Q}_{p}/2)^{2}}},\end{eqnarray}$$

where $C_{f}=\langle \hat{\unicode[STIX]{x1D70F}}_{s}\rangle /(\hat{\unicode[STIX]{x1D70C}}(\hat{Q}/2)^{2})$ is the laminar friction coefficient for a pressure-driven flow through a SHS channel with surfactants and $C_{f,p}=\langle \hat{\unicode[STIX]{x1D70F}}_{s}\rangle _{p}/(\hat{\unicode[STIX]{x1D70C}}(\hat{Q}_{p}/2)^{2})$ is the laminar friction coefficient for the equivalent Poiseuille channel flow driven with the same pressure gradient and for the same channel height. The quantities $\langle \unicode[STIX]{x1D70F}_{s}\rangle$ and $\langle \unicode[STIX]{x1D70F}_{s}\rangle _{p}$ are the surface stresses averaged along both top and bottom surfaces for a SHS channel flow and a Poiseuille channel flow with the same geometry, respectively. Since, by construction, the pressure gradient is the same for the flow in the SHS channel and the Poiseuille channel flow, we have $\langle \hat{\unicode[STIX]{x1D70F}}_{s}\rangle =\langle \hat{\unicode[STIX]{x1D70F}}_{s}\rangle _{p}$. Then, using (4.30) and (4.31) into (4.35) we find

(4.36)$$\begin{eqnarray}DR=1-\left(1+\frac{3Q_{d}}{2}\right)^{-2}=1-\left(1+\frac{3\unicode[STIX]{x1D706}_{e}}{\unicode[STIX]{x1D706}_{e}+2}\right)^{-2}.\end{eqnarray}$$

The maximum possible drag reduction is $DR\rightarrow 15/16$ as $\unicode[STIX]{x1D706}_{e}\rightarrow \infty$. We can compute $\unicode[STIX]{x1D706}_{e}$ in (4.36) using (4.34) and (4.29). We also provide, as supplementary materials, MATLAB routines computing $\unicode[STIX]{x1D706}_{e}$, $DR$ and $\unicode[STIX]{x1D6FE}_{Ma}$ for any specified flow-related, surfactant or geometrical parameters.

6.1 Effective slip length

In figure 4, we compare our scaling predictions for the effective slip length $\unicode[STIX]{x1D706}_{e}^{theory}$ with the numerical results $\unicode[STIX]{x1D706}_{e}^{data}$. We compute $\unicode[STIX]{x1D706}_{e}^{theory}$ using (4.34), where $\unicode[STIX]{x1D6FE}_{Ma}$ follows (4.29) and the coefficients $E_{0}$ are computed by solving the linear problem (4.12) in the surfactant-free case for each couple of geometrical parameters $(g,\unicode[STIX]{x1D719})$. The empirical parameters $a_{1}$, $a_{2}$ in (4.29) and $\unicode[STIX]{x1D6FF}_{0,i}$, $\unicode[STIX]{x1D6FF}_{1,i}$, with $i=1$, 2 or 3 for $\unicode[STIX]{x1D6FF}$ (see equations (3.20)–(3.22)) can be determined using a least-squares fitting approach and the trust region reflective algorithm, as implemented in the package optimize.least_squares of Scipy (Jones et al. Reference Jones, Oliphant and Peterson2001).

First, we determine $\unicode[STIX]{x1D6FF}_{1,3}$ in $\unicode[STIX]{x1D6FF}$ by fitting a measure of the characteristic diffusive boundary layer thickness in our numerical simulations, calculated using (3.24), with the scaling model given in (3.22). We have only used the Lévêque scaling (3.22) for $\unicode[STIX]{x1D6FF}$ in (4.29). In our numerical simulations, the diffusive boundary layer mostly follows the Lévêque regime, which assumes a background linear shear flow, since the slip velocity $u_{I}$ is small. Moreover, as also noted earlier, the scaling model (4.29) for $\unicode[STIX]{x1D6FE}_{Ma}$ depends weakly on $\unicode[STIX]{x1D6FF}$. Hence, the choice of scaling for $\unicode[STIX]{x1D6FF}$, which can vary between (3.20), (3.21) or (3.22) depending on the geometry and the slip, does not appear to be critical. The fit gives

(6.1)$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{1,3}=0.0528,\end{eqnarray}$$

from the minimization of the sum of the squares of the relative distance of theory from data, i.e. $(\unicode[STIX]{x1D6FF}^{theory}-\unicode[STIX]{x1D6FF}^{data})^{2}/(\unicode[STIX]{x1D6FF}^{data})^{2}$. This prior independent determination of $\unicode[STIX]{x1D6FF}_{1,3}$ reduces the number of fitting parameters to three in (4.29): $a_{1}$, $a_{2}$ and $\unicode[STIX]{x1D6FF}_{0,3}$. This ensures a more accurate and robust fit for $\unicode[STIX]{x1D706}_{e}^{theory}$, less sensitive to the actual fitting technique used.

Then, using $\unicode[STIX]{x1D6FF}_{1,3}=0.0528$ in (4.29), we fit the effective slip length $\unicode[STIX]{x1D706}_{e}^{theory}$ given by (4.34) to $\unicode[STIX]{x1D706}_{e}^{data}$ computed via the deviating flux $Q_{d}$ using (4.31). Incidentally, computing $\unicode[STIX]{x1D706}_{e}^{data}$ using (4.31) gives an accurate and robust estimation of the effective slip length in our numerical simulations, as it relies solely on the integral quantity $Q_{d}=Q-Q_{p}$ (see (4.30)). Minimizing the sum of the squares of the absolute distance between $\unicode[STIX]{x1D706}_{e}^{theory}$ and $\unicode[STIX]{x1D706}_{e}^{data}$, we obtain

(6.2a-c)$$\begin{eqnarray}a_{1}=2.30,\quad a_{2}=0.319\quad \text{and}\quad \unicode[STIX]{x1D6FF}_{0,3}=1.68.\end{eqnarray}$$

As we can see in figure 4(a,b), the scaling predictions for $\unicode[STIX]{x1D706}_{e}^{theory}$, using the values for $a_{1}$, $a_{2}$, $\unicode[STIX]{x1D6FF}_{0,3}$ and $\unicode[STIX]{x1D6FF}_{1,3}$ stated in the previous paragraph, show an excellent agreement with $\unicode[STIX]{x1D706}_{e}^{data}$ over a very large range: $10^{-12}\lesssim \unicode[STIX]{x1D706}_{e}\lesssim 1$.

The nine data points at non-negligible Reynolds numbers, $1\leqslant Re\leqslant 10^{5}$ (identified with blue circles in figure 4a,b), also exhibit good agreement despite violating the low-Reynolds-number assumption made in our flow model (see § 4.1). As explained previously in § 5, although the full steady nonlinear Navier–Stokes equation (2.3) was used in the simulations, the flow remained in the laminar regime for all Reynolds numbers tested.

At large non-dimensional background concentrations, $1\leqslant k\leqslant 100$ (identified with green triangles in figure 4b), the scaling predictions underestimate slightly the slip length. This is due to the fact that the model assumes a low concentration of surfactant. However, the model still provides a practically useful prediction of the boundary condition at the interface, which can be effectively considered as no slip for all our simulations with $k\geqslant 1$. We also find that the maximum boundary layer thickness is $\unicode[STIX]{x1D6FF}=1.20$, which suggests that our scaling prediction is accurate even if the diffusive boundary layer is vertically confined.

We indicate in figure 4 (as well as in figures 5 and 6) data where the interface properties are strongly non-uniform, which are labelled by vermilion squares and orange diamonds. Qualitatively similar interface non-uniformities have been studied extensively in the context of air bubbles rising in surfactant-contaminated water (e.g. Bond & Newton Reference Bond and Newton1928; Frumkin & Levich Reference Frumkin and Levich1947; Levich Reference Levich1962; He, Maldarelli & Dagan Reference He, Maldarelli and Dagan1991), where they correspond to the ‘stagnant cap regime’. In this regime, an upstream part of the interface has a negligible surfactant gradient and can be considered as shear free ($\unicode[STIX]{x1D6FE}_{Ma}\rightarrow 0$), whilst the rest of the interface downstream has a large Marangoni shear ($\unicode[STIX]{x1D6FE}_{Ma}\rightarrow 1$), leading to a no-slip condition over a portion of the bubble known as the ‘stagnant cap’ (hereafter designated as SC). In the SC regime, advection of surfactant along the interface dominates relative to surfactant transport between the interface and the bulk. This makes possible highly non-uniform interfacial concentrations. Since transport between the interface and the bulk is mediated by both the diffusive boundary layer flux and the surfactant kinetics, the SC regime requires that advection along the interface must be large compared to either diffusive or kinetic fluxes (or both).

We briefly summarize here the bubble-flow analysis of Palaparthi et al. (Reference Palaparthi, Papageorgiou and Maldarelli2006), and translate it to SHS flow. For a bubble, the SC regime is found when the characteristic interfacial Péclet number is large, and either the adsorption–desorption kinetics flux $S$ is small, or the diffusive flux through the boundary layer is small compared with the interfacial advective flux. Denoting with a superscript ‘bubble’ the results of Palaparthi et al. (Reference Palaparthi, Papageorgiou and Maldarelli2006), they showed that this implies $Pe_{I}^{bubble}\gg 1$, and ${\mathcal{K}}_{I}=Bi^{bubble}(1+k)\ll 1$ or ${\mathcal{D}}_{I}=\unicode[STIX]{x1D712}^{bubble}(1+k)/(Pe^{bubble})^{1/2}\ll 1$. For a bubble, the characteristic length and velocity scales are the bubble radius and interfacial velocity in the surfactant-free case. In order to translate these canonical bubble results to SHSs, note that the bubble radius is analogous to the grating length ${\hat{g}}$. For the SHS, the characteristic velocity scale for these non-dimensional numbers is the mid-gap interfacial velocity in the surfactant-free case, namely $\hat{u} _{Ic}(\unicode[STIX]{x1D6FE}_{Ma}=0)$, which differs from the bulk characteristic velocity, such that $\hat{u} _{Ic}(\unicode[STIX]{x1D6FE}_{Ma}=0)=2{\mathcal{F}}_{0}\hat{U}$, according to (4.22). This contrasts slightly with contaminated air bubbles in water, where the characteristic interfacial velocity in the surfactant-free case scales as the far-field bulk velocity, owing to the absence of rigid no-slip walls. As shown in figure 3 and explained in detail in appendix C, we have ${\mathcal{F}}_{0}\sim 1$ for $g\gtrsim 1$ (as for bubbles) and ${\mathcal{F}}_{0}\sim g$ for $g\lesssim 1$.

Therefore, using our dimensionless group definitions of § 2, and using a ‘$g$’ subscript to characterize dimensionless groups where we use the length scale ${\hat{g}}$, rather than ${\hat{h}}$, we have $Pe_{I}^{bubble}\mapsto {\mathcal{F}}_{0}gPe_{I}={\mathcal{F}}_{0}Pe_{I,g}$ and ${\mathcal{K}}_{I}=Bi^{bubble}(1+k)\mapsto {\mathcal{K}}_{I,g}=Bi_{g}(1+k)/{\mathcal{F}}_{0}$, as well as ${\mathcal{D}}_{I}=\unicode[STIX]{x1D712}^{bubble}(1+k)/(Pe^{bubble})^{1/2}\mapsto {\mathcal{D}}_{I,g}=\unicode[STIX]{x1D712}_{g}(1+k)/({\mathcal{F}}_{0}Pe_{g})^{1/2}$.

Figure 5.

Comparison of the scaling predictions for the drag reduction $DR^{theory}$, computed using (4.36), with the numerical results from our simulations $DR^{data}$. We also show with red crosses the drag reduction $DR^{theory}$ estimated using a model neglecting the surfactant effect. Note that we plot in this graph only data for $DR\geqslant 10^{-4}$, to show more clearly our results in a range useful to applications. Equality between data and theory falls on the grey dashed line.

The ranges spanned by the quantities ${\mathcal{F}}_{0}Pe_{I,g}$, ${\mathcal{K}}_{I,g}$ and ${\mathcal{D}}_{I,g}$ are reported in table 1.

The distinction between the partial SC regime, where the SC fills only part of the interface, and the full SC regime, where the SC fills all the interface, is revealed by an inspection of the shear rate profiles along the interface (not shown here). In the partial SC regime, the shear rate increases abruptly from negligible values to $\unicode[STIX]{x1D6FE}_{Ma}\sim 1$ at a particular location along the interface. In the partial SC regime, the non-dimensional numbers in our simulations range approximately as: $2.5\times 10^{3}\leqslant {\mathcal{F}}_{0}Pe_{I,g}\leqslant 2.5\times 10^{5}$, $9.9\times 10^{-4}\leqslant {\mathcal{K}}_{I,g}\leqslant 0.4$ and $0.04\leqslant {\mathcal{D}}_{I,g}\leqslant 0.4$ (see also figure 8, appendix A, for the variations of these numbers across all our numerical simulations and for the different regimes, as well as supplementary table S1 for the value of each parameter for each simulation). In the full SC regime, the non-dimensional numbers range approximately as: $52\leqslant {\mathcal{F}}_{0}Pe_{I,g}\leqslant 2.5\times 10^{4}$, $2\times 10^{-2}\leqslant {\mathcal{K}}_{I,g}\leqslant 50$ and $4.0\times 10^{-5}\leqslant {\mathcal{D}}_{I,g}\leqslant 1.3$. The interfacial Péclet number is mostly higher in the partial SC regime than in the full SC regime, which is intuitively expected. We can see in figure 4(a) that the four data points in the partial SC regime (plotted with vermilion squares) are the only data points where the scaling predictions significantly underestimates the effective slip length with $\unicode[STIX]{x1D706}_{e}^{theory}\leqslant 5.2\times 10^{-5}$, whereas $\unicode[STIX]{x1D706}_{e}^{data}\geqslant 3.5\times 10^{-4}$. This discrepancy is due to the strong non-uniformity of the shear rate profile in the SC regime, not taken into account by our scaling model which is based on the assumption that the shear rate is approximately uniform along the interface (see (3.17)). The predictions $\unicode[STIX]{x1D706}_{e}^{theory}$ in the full SC regime, plotted with orange diamonds, are in reasonable agreement with the data $\unicode[STIX]{x1D706}_{e}^{data}$. We can see that $\unicode[STIX]{x1D706}_{e}^{theory}$ underestimates slightly the data, although by less than one order of magnitude for all our results in the full SC regime, with $0.25\leqslant \unicode[STIX]{x1D706}_{e}^{theory}/\unicode[STIX]{x1D706}_{e}^{data}\leqslant 2.4$.

The data plotted with black pluses in figure 4, i.e. not in the SC regime, are in a state analogous to the ‘uniformly retarded regime’ described by Palaparthi et al. (Reference Palaparthi, Papageorgiou and Maldarelli2006) in their study of air bubbles rising in contaminated water, where they make the case that this regime exists for ${\mathcal{K}}_{I,g}\sim 1$ and ${\mathcal{D}}_{I,g}\sim 1$. However, in our simulations, we find that the interfacial shear rate is in the ‘uniform’ regime, and thus satisfies our modelling assumption, over a range of ${\mathcal{K}}_{I,g}$ and ${\mathcal{D}}_{I,g}$ that spans several orders of magnitude, implying that the vast majority of the simulations satisfy our modelling assumptions. More specifically, we find that simulations in the ‘uniformly retarded regime’ have parameters that satisfy approximately $2.8\times 10^{-3}\leqslant {\mathcal{D}}_{I,g}\leqslant 4.4\times 10^{3}$ and $1.9\times 10^{-2}\leqslant {\mathcal{K}}_{I,g}\leqslant 3.2\times 10^{3}$. This is most likely due to the fact that some of our simulations are in an intermediate or transition regime between the SC regime and the uniformly retarded regime, and for which $\unicode[STIX]{x1D706}_{e}$ still follows our scaling prediction, although perhaps with slightly more scatter, as shown by some of the black pluses in figure 4.

In figure 4(c), we show the relative error between the scaling predictions $\unicode[STIX]{x1D706}_{e}^{theory}$ and the numerical results $\unicode[STIX]{x1D706}_{e}^{data}$ for the effective slip length. The error remains relatively small across all values of the average interfacial shear rate $\unicode[STIX]{x1D6FE}_{Ma}$. It is less than approximately 33 % for $\unicode[STIX]{x1D6FE}_{Ma}\leqslant 0.7$, except for the four simulations in the partial SC regime plotted with vermilion squares. The relative error is less than 1.7 for $0.7\leqslant \unicode[STIX]{x1D6FE}_{Ma}\leqslant 1$.

For comparison, we also show with red crosses in figure 4(c) the prediction from a surfactant-free model, which is obtained using (4.34) with $\unicode[STIX]{x1D6FE}_{Ma}=0$. Our model provides consistently better predictions than those that neglect the surfactant effect. In particular, the error made by neglecting surfactant effects becomes very large when the interfacial shear rate increases towards the Poiseuille value $\unicode[STIX]{x1D6FE}_{p}=1$. At low shear rate, $\unicode[STIX]{x1D6FE}_{Ma}\leqslant 0.3$ we can see that the two models have comparable (small) relative errors.

Overall, we find that our scaling model for $\unicode[STIX]{x1D706}_{e}$ provides excellent quantitative predictions across a large range of non-dimensional numbers, beyond the strict range of validity based on our modelling assumptions. Although our model predictions can underestimate the slip length in some cases (at large concentrations, and in the stagnant cap regime), our model remains practically useful as both theory and simulation yield negligible slip in those instances.

6.2 Drag reduction

We compare the drag reduction predicted by our theory ($DR^{theory}$) with the numerical results from our simulations ($DR^{data}$), as shown in figure 5. The value of $DR^{theory}$ is obtained from (4.36), where the corresponding values of $\unicode[STIX]{x1D706}_{e}^{theory}$ were shown earlier in figure 4. Similarly, $DR^{data}$ is calculated using $\unicode[STIX]{x1D706}_{e}^{data}$, whose values are also shown in figure 4. Using a log–log scale, we only plot data for $DR\geqslant 10^{-4}$, which correspond to the more meaningful range for practical applications. The predictions from our scaling model are in very good agreement with the numerical results. Data at even lower drag reductions (not shown here) still exhibit a very good agreement with our theoretical prediction.

In figure 5, we also plot, using red crosses, the drag reduction computed using a surfactant-free model. This is obtained by substituting the values for the surfactant-free $\unicode[STIX]{x1D706}_{e}^{theory}$ (plotted with red crosses in figure 4c) into (4.36). As may be expected, the surfactant-free theory almost always incorrectly predicts a larger drag reduction, with values often more than an order of magnitude larger than the actual ones. This clearly shows that the drag reduction potential of SHSs can be significantly overestimated in conditions where surfactants are important. This is consistent with the findings of Peaudecerf et al. (Reference Peaudecerf, Landel, Goldstein and Luzzatto-Fegiz2017), who showed that, for SHSs with rectangular longitudinal gratings, surfactant effects become important at very low concentrations, similar to background levels found in the environment. As may be expected, the few surfactant-free predictions in figure 5 that show better agreement with the numerical simulations correspond to lower values of $\unicode[STIX]{x1D6FE}_{Ma}$, when the surfactant-free predictions converge towards our model predictions (see figure 4c).

6.3 Interfacial shear rate

Figure 6.

Comparison of the scaling predictions for the average interfacial shear rate $\unicode[STIX]{x1D6FE}_{Ma}^{theory}$, computed using (4.29), with the numerical results from our simulations $\unicode[STIX]{x1D6FE}_{Ma}^{data}$, calculated by averaging the shear rate along the interface. The scaling predictions use empirical parameters: $a_{1}=2.30$ and $a_{2}=0.319$ for $\unicode[STIX]{x1D6FE}_{Ma}$, and $\unicode[STIX]{x1D6FF}_{0,3}=1.68$ and $\unicode[STIX]{x1D6FF}_{1,3}=0.0528$ for $\unicode[STIX]{x1D6FF}$ (see (3.22)), computed from the fit of $\unicode[STIX]{x1D706}_{e}$ (see § 6.1). In (a) we plot using a log–log scale $(1-\unicode[STIX]{x1D6FE}_{Ma})$ to reveal the behaviour at large shear rate, when $\unicode[STIX]{x1D6FE}_{Ma}^{data}\rightarrow \unicode[STIX]{x1D6FE}_{p}=1$. The predictions $\unicode[STIX]{x1D6FE}_{Ma}^{theory}$ for the four data points in the partial stagnant cap regime, plotted with vermilion squares, overestimate the data owing to the strong non-uniformity of the interfacial shear rate profile. In (b) we plot $(1-\unicode[STIX]{x1D6FE}_{Ma})$ over a larger range, revealing the error related to the singularity at the stagnation points. In (c) we plot $\unicode[STIX]{x1D6FE}_{Ma}$ to show more clearly the behaviour at intermediate shear rates, $0.7\leqslant \unicode[STIX]{x1D6FE}_{Ma}\leqslant 1$. In all plots, equality between data and theory would fall on the grey dashed line.

We compare in figure 6 the numerical results for the average interfacial shear rate $\unicode[STIX]{x1D6FE}_{Ma}^{data}$ with the theoretical predictions, $\unicode[STIX]{x1D6FE}_{Ma}^{theory}$ computed using (4.29) using the four empirical parameters optimized for $\unicode[STIX]{x1D706}_{e}$ in § 6.1: $a_{1}=2.30$ and $a_{2}=0.319$ for $\unicode[STIX]{x1D6FE}_{Ma}$, and $\unicode[STIX]{x1D6FF}_{0,3}=1.68$ and $\unicode[STIX]{x1D6FF}_{1,3}=0.0528$ for $\unicode[STIX]{x1D6FF}$ based on (3.22). The numerical results for $\unicode[STIX]{x1D6FE}_{Ma}^{data}$ have been computed by taking the spatial average of the interfacial shear rate in the interior of the interface $-g/2\leqslant x\leqslant g/2$.

In figure 6(a), we show $(1-\unicode[STIX]{x1D6FE}_{Ma})$ in a log–log plot to focus on the no-slip limit $\unicode[STIX]{x1D6FE}_{Ma}\rightarrow 1$. Over the limited range shown on this graph, we find good agreement between our scaling predictions and the data for all our numerical simulations where the interfacial shear rate is found to be approximately uniform along the interface (see the uniform regime, plotted with black plusses). Similar to $\unicode[STIX]{x1D706}_{e}$ shown in figure 4, we can see that the four data points in the partial SC regime (vermilion squares) with $(1-\unicode[STIX]{x1D6FE}_{Ma}^{data})\geqslant 0.4$ are the only ones where the predictions underestimate the data. As discussed earlier, this is due to the strong non-uniformity of the shear rate profile in the SC regime, in contradiction to the uniform assumption made in our model (see (3.17)). Nevertheless, the theory remains practically useful, as both model and simulation yield a shear that is essentially indistinguishable from that of a no-slip boundary.

In figure 6(b) (which is the inset of figure 6a), we plot $(1-\unicode[STIX]{x1D6FE}_{Ma}^{data})$ over the full range of values tested. As the average shear rate tends to the maximum Poiseuille value, $\unicode[STIX]{x1D6FE}_{Ma}^{data}\rightarrow \unicode[STIX]{x1D6FE}_{p}=1$ or equivalently $(1-\unicode[STIX]{x1D6FE}_{Ma}^{data})\rightarrow 0$, $(1-\unicode[STIX]{x1D6FE}_{Ma}^{theory})$ underestimates the data. The difference becomes significant for $(1-\unicode[STIX]{x1D6FE}_{Ma}^{theory})\lesssim 10^{-3}$. This is due to the singularity at the two stagnation points and the difficulty associated with resolving it numerically. The shear rate exhibits extreme variations very close to the stagnation points, whilst the shear rate remains flat in the interior of the interface with values very close to the Poiseuille shear rate. We note, however, that the effect of the singularity appears only in the limit $\unicode[STIX]{x1D6FE}_{Ma}^{data}\rightarrow \unicode[STIX]{x1D6FE}_{p}=1$, at values practically equivalent to a no-slip boundary condition at the interface.

This can also be seen in figure 6(c), where we plot $\unicode[STIX]{x1D6FE}_{Ma}$ directly, for $\unicode[STIX]{x1D6FE}_{Ma}\geqslant 0.7$. The scaling predictions consistently predict a no-slip boundary condition $\unicode[STIX]{x1D6FE}_{Ma}^{theory}\rightarrow 1$, as $\unicode[STIX]{x1D6FE}_{Ma}^{data}\rightarrow 1$. This shows that the actual error between $\unicode[STIX]{x1D6FE}_{Ma}^{theory}$ and $\unicode[STIX]{x1D6FE}_{Ma}^{data}$ is actually very small in this limit, where we find the simulations at large Reynolds numbers $Re\geqslant 1$ (blue circles), large non-dimensional concentrations $k\geqslant 1$ (green triangles) and in the full SC regime (orange diamonds). Predictions at intermediate values (shown by plusses in figure 6c), for $0.7\leqslant \unicode[STIX]{x1D6FE}_{Ma}^{theory}\leqslant 1$, show a good agreement with $\unicode[STIX]{x1D6FE}_{Ma}^{data}$ although with a slight overestimation.

Therefore, our scaling model also provides reasonable predictions across the whole range of interfacial shear rate values, even though the model has been fitted for $\unicode[STIX]{x1D706}_{e}$ and not for $\unicode[STIX]{x1D6FE}_{Ma}$. An agreement is found from intermediate to large values, provided the interface is not in a partial SC regime. Our scaling model remains accurate across a broad range of non-dimensional numbers (see table 1 and figure 8, appendix A) and in the full SC regime.