1. Introduction
Superhydrophobic surfaces (SHSs) offer a promising avenue for decreasing drag in laminar and turbulent flows (Lee et al. Reference Lee, Choi and Kim2016; Park et al. Reference Park, Choi and Kim2021). Hydrophobic chemistry and microscopic topography on the SHSs entrap gas bubbles that lubricate the flow and decrease the wall-average shear stress compared with solid walls. This design presents numerous opportunities for industrial and environmental applications, including enhanced cooling (Lam et al. Reference Lam, Hodes and Enright2015) and reduced emissions (Xu et al. Reference Xu, Grabowski, Yu, Kerezyte, Lee, Pfeifer and Kim2020). Consequently, SHSs have attracted significant attention in the academic literature (Park et al. Reference Park, Sun and Kim2014; Schönecker et al. Reference Schönecker, Baier and Hardt2014; Türk et al. Reference Türk, Daschiel, Stroh, Hasegawa and Frohnapfel2014; Cheng et al. Reference Cheng, Xu and Sui2015; Seo & Mani Reference Seo and Mani2018; Rastegari & Akhavan Reference Rastegari and Akhavan2019; Kirk et al. Reference Kirk, Karamanis, Crowdy and Hodes2020; Landel et al. Reference Landel, Peaudecerf, Temprano-Coleto, Gibou, Goldstein and Luzzatto-Fegiz2020; Tomlinson et al. Reference Tomlinson, Peaudecerf, Temprano-Coleto, Gibou, Luzzatto-Fegiz, Jensen and Landel2023b ). However, recent investigations have highlighted potential challenges. Interfacial displacement and deflection (Ng & Wang Reference Ng and Wang2009; Teo & Khoo Reference Teo and Khoo2010), gas-phase flow dynamics (Game et al. Reference Game, Hodes, Keaveny and Papageorgiou2017) and heat- and surfactant-induced Marangoni stresses (Peaudecerf et al. Reference Peaudecerf, Landel, Goldstein and Luzzatto-Fegiz2017; Kirk et al. Reference Kirk, Karamanis, Crowdy and Hodes2020) have been shown to impede the anticipated drag reduction. Consequently, in certain applications, SHSs may not outperform solid walls. Motivated by these questions, this study investigates drag reduction in a surfactant-contaminated flow confined between a SHS and a solid wall. Specifically, we focus on the computationally challenging regime of weak bulk molecular diffusion, a regime which has remained largely unexplored in prior literature.
Surfactants have been measured in both artificial (Hourlier-Fargette et al. Reference Hourlier-Fargette, Dervaux, Antkowiak and Neukirch2018; Temprano-Coleto et al. Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023) and natural (Pereira et al. Reference Pereira, Ashton, Sabbaghzadeh, Shutler and Upstill-Goddard2018; Frossard et al. Reference Frossard2019) environments. When a liquid flow becomes contaminated with soluble surfactant, these chemical compounds are transported throughout the flow, adsorbing onto interfaces and desorbing into the bulk. The surfactant distribution within the bulk is regulated by the bulk Péclet number (
$\textit {Pe}$
), forming concentration boundary layers when
$\textit {Pe}$
is large. As the flow carries surfactant towards the downstream stagnation point of an interface, it accumulates, forming a concentration gradient. Depending on factors such as the flow properties, surfactant characteristics and geometry, the resulting adverse Marangoni force can cause the theoretically shear-free interface to behave like a no-slip (or partially no-slip) wall (Peaudecerf et al. Reference Peaudecerf, Landel, Goldstein and Luzzatto-Fegiz2017). Experimental support for this mechanism can be found in studies by Kim & Hidrovo (Reference Kim and Hidrovo2012), Bolognesi et al. (Reference Bolognesi, Cottin-Bizonne and Pirat2014), Schäffel et al. (Reference Schäffel, Koynov, Vollmer, Butt and Schönecker2016), Peaudecerf et al. (Reference Peaudecerf, Landel, Goldstein and Luzzatto-Fegiz2017), Song et al. (Reference Song, Song, Hu, Du, Du, Choi and Rothstein2018) and Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023).
Peaudecerf et al. (Reference Peaudecerf, Landel, Goldstein and Luzzatto-Fegiz2017) and Landel et al. (Reference Landel, Peaudecerf, Temprano-Coleto, Gibou, Goldstein and Luzzatto-Fegiz2020) conducted numerical simulations of a two-dimensional (2-D) channel flow featuring streamwise-periodic SHSs using COMSOL, employing surfactant properties that are characteristic of sodium dodecylsulphate (SDS). Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023) extended this methodology to encompass a three-dimensional channel flow featuring streamwise- and spanwise-periodic SHSs, while Sundin & Bagheri (Reference Sundin and Bagheri2022) investigated 2-D channel flow laden with surfactants over streamwise-periodic liquid-infused surfaces. In all of these configurations, numerical simulations of the Navier–Stokes equations and advection–diffusion equations for bulk and interfacial surfactant incurred substantial computational costs, particularly in the high-Péclet-number regime, where bulk-concentration boundary layers can be extremely thin. To address these challenges of high Péclet number values, we develop a numerical method based on Chebyshev collocation. This numerical method clusters nodes around the bulk-concentration boundary layers, reducing computational demand compared with numerical simulations. We complement numerical simulations with asymptotic analysis that reveals the underlying flow structures and offers direct predictions of flow properties.
The distinction between strictly insoluble and very weakly soluble surfactant can be subtle. In the former case, there is no exchange between bulk and interface and the amount of surfactant present on a SHS plastron must be prescribed as an input parameter. In the latter case, adsorption and desorption processes control bulk and interface exchange and the amount of surfactant present on a SHS plastron is determined by coupled transport processes on the interface and in the bulk. In many practical applications, diffusive effects are weak, and might not be expected to influence the manner in which surfactant compromises drag reduction. However, the surfactant flux between bulk and interface is mediated by bulk diffusion, enabling it to play a central role in determining drag reduction, even at very high Péclet numbers. In such circumstances, weak bulk diffusion can be the cause of weak apparent solubility, even when adsorption and desorption processes are very fast. Furthermore, while bulk diffusion cannot be neglected, interfacial diffusion can be, but then adsorption and desorption must be accommodated by weak stretching or compression of the interface. Taking these factors together, we will see how structures such as the classical stagnant cap, of some prescribed size when surfactant is strictly insoluble, must be replaced by a more mobile structure that we term a quasi-stagnant cap.
Scaling theories have been developed to analyse the slip length over SHSs and liquid-infused surfaces with streamwise-periodic gratings (Landel et al. Reference Landel, Peaudecerf, Temprano-Coleto, Gibou, Goldstein and Luzzatto-Fegiz2020; Sundin & Bagheri Reference Sundin and Bagheri2022) and SHSs with streamwise- and spanwise-periodic gratings (Temprano-Coleto et al. Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023). These theories typically assume small surfactant concentrations and uniform shear stresses at the interface. For common surfactants such as SDS, Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023) demonstrated a significant slip length when the grating length exceeds a modified depletion length and a mobilisation length. The modified depletion length depends on the depletion length (Manikantan & Squires Reference Manikantan and Squires2020), the height of the channel and the ratio of interfacial and bulk diffusivities. It characterises the interfacial length above which interfacial diffusion is weak compared to bulk–interface exchanges. The mobilisation length depends on the surfactant concentration and surfactant properties through the Marangoni, Damköhler and Biot numbers, and the SHS geometry (Temprano-Coleto et al. Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023). The mobilisation length being generally larger than the modified depletion length, it is the key length scale above which an interface can display significant slip. Landel et al. (Reference Landel, Peaudecerf, Temprano-Coleto, Gibou, Goldstein and Luzzatto-Fegiz2020), Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023) and Sundin & Bagheri (Reference Sundin and Bagheri2022) used numerical simulations to calibrate the empirical coefficients linked to the bulk-concentration boundary-layer thickness. Unlike these prior studies, our asymptotic approach to the problem bypasses the need for empirical coefficient fitting, extending the slip-length scalings identified in previous studies (Landel et al. Reference Landel, Peaudecerf, Temprano-Coleto, Gibou, Goldstein and Luzzatto-Fegiz2020; Temprano-Coleto et al. Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023) and uncovering new drag reduction regimes for surfactant-contaminated SHSs.
Tomlinson et al. (Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a
) developed an asymptotic theory for laminar channel flow featuring streamwise- and spanwise-periodic grooves and then investigated the impact of spatio–temporal fluctuations in surfactant concentration over the SHS (Tomlinson et al. Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2024). Their theory assumed that the channel is long and bulk diffusion is strong enough to eliminate cross-channel concentration gradients. They derived a spatiallyone-dimensional (1-D) long-wave model that accommodated non-uniform shear stresses at the interface and enabled an analysis of the soluble stagnant-cap regime at low bulk Péclet numbers. The numerical simulations performed at moderate bulk Péclet number by Sundin & Bagheri (Reference Sundin and Bagheri2022) show adsorption at the upstream end of the stagnant cap, but they did not offer a scaling theory for this regime. Similar stagnant-cap scenarios, where there is no adsorption at the upstream end of the interface, were studied by Baier & Hardt (Reference Baier and Hardt2021) and Mayer & Crowdy (Reference Mayer and Crowdy2022) in the context of insoluble surfactant, with a linear and nonlinear equation of state, respectively. These studies were extended to include protrusion at the liquid–gas interface: Baier & Hardt (Reference Baier and Hardt2022) found recirculating interfacial flows for weak protrusions, which Rodriguez-Broadbent & Crowdy (Reference Rodriguez-Broadbent and Crowdy2023) extended to larger protrusions, deriving an effective slip length. Crowdy et al. (Reference Crowdy, Curran and Papageorgiou2023) investigated a linear extensional flow between two fluids. Assuming that
$\textit {Pe}$
is zero, they showed that solubility and strong exchange with the bulk can reduce the length of the stagnant cap and make the interface less immobile. We show how ‘remobilisation’ due to strong bulk–surface exchange can also arise when
$Pe$
is large, revealing the coupling between the bulk-concentration boundary layer and the underlying interface, which is particularly intricate when Marangoni effects are weak.
Tomlinson et al. (Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a
) delineated three key asymptotic regions of parameter space: one dominated by Marangoni effects, with minimal drag reduction; and others dominated by advection and diffusion, with significant drag reduction. They also identified when their 1-D long-wave model first breaks down due to 2-D effects. In this paper, we employ a combination of asymptotic and numerical methods to investigate the impact of soluble surfactant on a laminar pressure-driven channel flow confined between a streamwise-periodic SHS and a solid wall. We assume that the period of the SHS is significantly longer than the height of the channel and that the bulk and surface Péclet numbers are large. To address the numerical challenges of these high-Péclet-number flows, we develop asymptotic solutions for weak diffusion, which may offer insights into a broader range of surfactant-contaminated flows than flows over SHSs. For example, numerical simulations of flows involving soluble-surfactant-contaminated drops and bubbles exhibit velocity and surfactant distributions reminiscent of stagnant caps at large but finite
$\textit {Pe}$
(Oguz & Sadhal Reference Oguz and Sadhal1988; Tasoglu et al. Reference Tasoglu, Demirci and Muradoglu2008). However, the structure of these stagnant caps has only been described in the insoluble limit (Harper Reference Harper2004; Palaparthi et al. Reference Palaparthi, Papageorgiou and Maldarelli2006).
Additionally, we construct asymptotic expressions for practical quantities such as the slip length and drag reduction. These asymptotic expressions not only aid in understanding the physics of surfactant-contaminated flows over SHSs but also facilitate the prediction of the slip length and drag reduction in experiments.
In § 2, we formulate the 2-D problem. We non-dimensionalise the fluid and surfactant equations in the long-wave limit and derive 1-D and 2-D long-wave models for bulk and interfacial surfactant behaviour. The calculations underpinning drag reduction predictions are laid out in appendices. In § 3, we present our main asymptotic and numerical findings. We analyse the drag reduction and compare the predicted flow and surfactant fields with results obtained from COMSOL simulations and experiments. In § 4, we illustrate the uses of our study and discuss the key asymptotic results that may be generalisable to broader surfactant-contaminated flows.
2. Formulation
2.1. Governing equations of the 2-D model
Figure 1. A 2-D laminar pressure-driven channel flow transporting soluble surfactant confined between a streamwise-periodic SHS and a solid wall. We position the origin of the Cartesian coordinate system,
$(\hat {x}, \, \hat {y})$
, at the centre of the liquid–gas interface. Each periodic cell has channel height
$2\hat {H}$
and period length
$2\hat {P}$
. The SHS, characterised by gas fraction
$\phi$
, has an interface region of length
$2\phi \hat {P}$
and a solid region of length
$2(1 - \phi )\hat {P}$
; the area above these regions defines subdomains
$\hat {\mathcal{D}}_1$
and
$\hat {\mathcal{D}}_2$
, respectively, as specified in (2.1).
We explore a steady 2-D laminar pressure-driven channel flow contaminated with soluble surfactant, confined between a streamwise-periodic SHS and a solid wall (see figure 1). The streamwise and wall-normal directions are represented by
$\hat {x}$
and
$\hat {y}$
coordinates respectively, where hats indicate dimensional quantities. Assuming the fluid to be incompressible and Newtonian, we define the velocity field
$\hat {\boldsymbol {u}}=(\hat {u}(\hat {x},\,\hat {y}),\,\hat {v}(\hat {x},\,\hat {y}))$
, pressure field
$\hat {p}=\hat {p}(\hat {x},\,\hat {y})$
, bulk surfactant concentration
$\hat {c}=\hat {c}(\hat {x},\,\hat {y})$
and interfacial surfactant concentration
$\hat {\Gamma }=\hat {\Gamma }(\hat {x})$
. Due to the periodicity of the SHS in the streamwise direction, our analysis focuses on a single periodic cell with dimensions
$2\hat {P}$
in length and
$2\hat {H}$
in height. At
$\hat {y}=0$
, there is a solid ridge of length
$2(1-\phi )\hat {P}$
and a liquid–gas interface, or plastron, of length
$2\phi \hat {P}$
, where
$\phi$
is the gas fraction. The interface is assumed to be flat. There is a solid boundary at
$\hat {y}= 2\hat {H}$
. The periodic domain is partitioned into two subdomains,
\begin{align} \hat {\mathcal{D}}_1 &= \{\hat {x}\in [-\phi \hat {P},\, \phi \hat {P}]\} \times \{\hat {y}\in [0, \,2 \hat {H}]\}, \end{align}
\begin{align} \hat {\mathcal{D}}_2 &= \{\hat {x}\in [\phi \hat {P}, \, (2 - \phi )\hat {P}]\} \times \{\hat {y}\in [0,\, 2 \hat {H}]\}, \end{align}
as illustrated in figure 1.
A comprehensive discussion of the equations governing fluid and surfactant behaviour in a three-dimensional geometry featuring streamwise- and spanwise-periodic SHSs was presented in Tomlinson et al. (Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a
). Here, we provide a succinct overview of the model for a streamwise-periodic SHS in two dimensions. Within domains
$\hat {\mathcal{D}}_1$
and
$\hat {\mathcal{D}}_2$
, we have Stokes equations with an advection–diffusion equation for bulk surfactant:
\begin{align} \hat {\boldsymbol {\nabla }} \cdot \hat {\boldsymbol {u}} = 0, \quad \hat {\mu } \hat {\nabla }^2 \hat {\boldsymbol {u}}- \hat {\boldsymbol {\nabla }} \hat {p} = \boldsymbol {0}, \quad \hat {D} \hat {\nabla }^2 \hat {c} - \hat {\boldsymbol {u}}\cdot \hat {\boldsymbol {\nabla }} \hat {c} = 0, \end{align}
where
$\hat {\mu }$
is the dynamic viscosity and
$\hat {D}$
is the bulk diffusivity. Linearising the equation of state, we assume
$\hat {\sigma }_{\hat {x}} = - \hat {A}\hat {\Gamma }_{\hat {x}}$
, where
$\hat {\sigma }$
is the surface tension and
$\hat {A}$
is the surface activity. For
$\hat {x} \in [-\phi \hat {P}, \, \phi \hat {P}]$
and
$\hat {y}=0$
, we impose the tangential stress balance, no penetration, linear adsorption–desorption kinetics and an advection–diffusion equation for interfacial surfactant:
\begin{align} \hat {\mu } \hat {u}_{\hat {y}} - \hat {A} \hat {\Gamma }_{\hat {x}} = 0, \quad \hat {v}=0, \quad \hat {D} \hat {c}_{\hat {y}} - \hat {K}_a \hat {c} + \hat {K}_d \hat {\Gamma } =0,\nonumber \\ \quad \hat {D}_I \hat {\Gamma }_{\hat {x}\hat {x}} + \hat {K}_a \hat {c} - \hat {K}_d \hat {\Gamma } -(\hat {u} \hat {\Gamma })_{\hat {x}} = 0, \end{align}
where
$\hat {D}_I$
is the interfacial diffusivity,
$\hat {K}_a$
is the adsorption rate and
$\hat {K}_d$
is the desorption rate. At
$\hat {x} = \pm \phi \hat {P}$
and
$\hat {y}=0$
, there is no interfacial flux of surfactant:
\begin{align} \hat {u}\hat {\Gamma } - \hat {D}_I \hat {\Gamma }_{\hat {x}} = 0. \end{align}
For
$\hat {x} \in [\phi \hat {P}, \, (2-\phi )\hat {P}]$
and
$\hat {y}=0$
, and likewise for
$\hat {x} \in [-\phi \hat {P}, \, (2 - \phi )\hat {P}]$
and
$\hat {y}=2\hat {H}$
, we impose no slip, no penetration and no bulk flux of surfactant:
\begin{align} \hat {u} = 0, \quad \hat {v} = 0, \quad \hat {c}_{\hat {y}} =0. \end{align}
For
$\hat {\boldsymbol {q}} \equiv (\hat {\boldsymbol {u}}, \,\hat {p}_{\hat {x}},\,\hat {c})$
, continuity and periodicity conditions on bulk quantities are given by
\begin{align} \hat {\boldsymbol {q}}((\phi \hat {P})^-,\,\hat {y}) = \hat {\boldsymbol {q}}((\phi \hat {P})^+,\,\hat {y}), \quad \hat {\boldsymbol {q}}(-\phi \hat {P},\,\hat {y}) = \hat {\boldsymbol {q}}((2- \phi )\hat {P},\,\hat {y}), \end{align}
where the superscript
$-$
(
$+$
) indicates evaluation in
$\hat {\mathcal{D}}_1$
(
$\hat {\mathcal{D}}_2$
). The pressure field is continuous across
$x=\pm \phi \hat {P}$
but is not periodic. Within
$\hat {\mathcal{D}}_1$
, the bulk surfactant equations and boundary conditions, (2.2c
), (2.3c
) and (2.5c
), can be combined to give
\begin{align} \frac {\text{d}}{\text{d}\hat {x}} \int _{0}^{2 \hat {H}} (\hat {u}\hat {c} - \hat {D} \hat {c}_{\hat {x}}) \, \text{d}\hat {y} - (\hat {K}_d \hat {\Gamma } - \hat {K}_a \hat {c}(\hat {x},\, 0)) = 0, \end{align}
which, combined with the interfacial surfactant equation (2.3d
), leads to an expression for the total flux of surfactant,
$\hat {K}$
, which must be uniform:
\begin{align} \int _{0}^{2 \hat {H}} (\hat {u}\hat {c} - \hat {D} \hat {c}_{\hat {x}}) \, \text{d}\hat {y} + (\hat {u}(\hat {x},\, 0)\hat {\Gamma } - \hat {D}_I \hat {\Gamma }_{\hat {x}}) = \hat {K}. \end{align}
Integrating the interfacial surfactant equation (2.3d ) over the plastron and utilising the no-flux condition (2.4) yields a constraint on net adsorption and desorption:
\begin{align} \int _{ - \phi \hat {P}}^{\phi \hat {P}} (\hat {K}_d \hat {\Gamma } - \hat {K}_a \hat {c}(\hat {x},\, 0)) \, \text{d}\hat {x} = 0. \end{align}
Similarly, within
$\hat {\mathcal{D}}_2$
, (2.3c
) and (2.5c
) can be combined to give
\begin{align} \int _{0}^{2 \hat {H}} (\hat {u}\hat {c} - \hat {D} \hat {c}_{\hat {x}}) \, \text{d}\hat {y} = \hat {K}. \end{align}
From (2.2a
), (2.3b
) and (2.5b
), it follows that the volume flux of fluid in
$\hat {\mathcal{D}}_1$
and
$\hat {\mathcal{D}}_2$
,
$\hat {Q}$
, is also uniform, where
\begin{align} \hat {Q} = \int _{0}^{2 \hat {H}} \hat {u} \, \text{d}\hat {y}. \end{align}
The fluxes
$\hat {Q}$
and
$\hat {K}$
are prescribed in this model.
The flow is driven in the streamwise direction by a cross-channel-averaged pressure drop per period,
$\Delta \hat {p} \equiv \langle \hat {p}\rangle (-\phi \hat {P}) - \langle \hat {p}\rangle ((2-\phi )\hat {P}) \gt 0$
, where
$\langle \cdot \rangle \equiv \int _{0}^{2 \hat {H}} \cdot \, \text{d}\hat {y} / (2\hat {H})$
. We define the normalised drag reduction
\begin{align} {DR} = \frac {\Delta \hat {p}_I - \Delta \hat {p}}{\Delta \hat {p}_I - \Delta \hat {p}_U}, \end{align}
where
$\Delta \hat {p} = \Delta \hat {p}_I$
when the interface is immobilised due to surfactant (
${DR}=0$
) and
$\Delta \hat {p} = \Delta \hat {p}_U$
when the interface is unaffected by surfactant (
${DR}=1$
). We aim to determine the dependence of
$DR$
on the dimensional parameters of the problem (
$\phi$
,
$\hat {P}$
,
$\hat {H}$
,
$\hat {\mu }$
,
$\hat {D}$
,
$\hat {A}$
,
$\hat {D}_I$
,
$\hat {K}_a$
,
$\hat {K}_d$
,
$\hat {K}$
and
$\hat {Q}$
) in the singular high-Péclet-number regime.
2.2. Non-dimensionalisation
We non-dimensionalise the governing equations (2.1)–(2.12) using
$\hat {Q}/\hat {H}$
for the velocity scale,
$\hat {K}/\hat {Q}$
for the bulk concentration scale and
$\hat {K}_a \hat {K}/(\hat {K}_d\hat {Q})$
for the interfacial concentration scale. We write
\begin{align} x = \frac {\hat {x}}{\hat {P}}, \quad y = \frac {\hat {y}}{\epsilon \hat {P}}, \quad u = \frac {\hat {u}}{\hat {Q}/\hat {H}}, \quad v = \frac {\hat {v}}{\epsilon \hat {Q}/\hat {H}}, \nonumber \\ p = \frac {\hat {p}}{\hat {\mu } \hat {Q}/(\epsilon \hat {H}^2)}, \quad c = \frac {\hat {c}}{\hat {K}/\hat {Q}}, \quad \Gamma = \frac {\hat {\Gamma }}{\hat {K}_a\hat {K}/(\hat {K}_d\hat {Q})}, \end{align}
where
$\epsilon \equiv \hat {H}/ \hat {P}$
is the slenderness parameter. The subdomains (2.1) become
\begin{align} \mathcal{D}_1 &= \{x\in [-\phi, \,\phi ]\} \times \{y\in [0, \,2]\}, \end{align}
\begin{align} \mathcal{D}_2 &= \{x\in [\phi, \,2 - \phi ]\} \times \{y\in [0,\, 2]\}. \end{align}
In
$\mathcal{D}_1$
and
$\mathcal{D}_2$
, the bulk equations (2.2) become
\begin{align} u_x + v_{y} = 0, \quad \epsilon ^2 u_{xx} + u_{yy} - p_x = 0,\nonumber \\ \epsilon ^4 v_{xx} + \epsilon ^2 v_{yy} - p_y = 0, \quad (\epsilon ^2 c_{xx} + c_{yy})/{\textit {Pe}} - \epsilon u c_x - \epsilon v c_y = 0, \end{align}
with
${\textit {Pe}} = \hat {Q}/\hat {D}$
the bulk Péclet number. For
$x\in [-\phi, \, \phi ]$
and
$y=0$
, the interface conditions (2.3) give
\begin{align} u_y - \epsilon {\textit {Ma}} \Gamma _{x} = 0, \quad v =0, \quad c_y - {\textit {Da}} (c -\Gamma ) =0, \nonumber \\ \quad \epsilon ^2 \Gamma _{xx}/{\textit {Pe}}_I + {\textit {Bi}} ( c - \Gamma ) - \epsilon (u \Gamma )_{x} = 0, \end{align}
with
${\textit {Ma}} = \hat {A}\hat {K}_a \hat {K}\hat {H}/(\hat {\mu }\hat {K}_d\hat {Q}^2)$
the Marangoni number,
${\textit {Da}}= \hat {K}_a \hat {H}/\hat {D}$
the Damköhler number,
${\textit {Pe}}_I= \hat {Q}/\hat {D}_I$
the interfacial Péclet number and
${\textit {Bi}} = \hat {K}_d \hat {H}^2/\hat {Q}$
the Biot number. At
$x = \pm \phi$
and
$y=0$
, the no-flux condition (2.4) becomes
\begin{align} u \Gamma - \epsilon \Gamma _{x}/{\textit {Pe}}_I = 0. \end{align}
For
$x\in [\phi, \, 2 - \phi ]$
and
$y=0$
(and for
$x\in [-\phi, \, 2-\phi ]$
and
$y=2$
), the solid boundary conditions (2.5) give
\begin{align} u = 0, \quad v = 0, \quad \ c_{y} =0. \end{align}
For
$\boldsymbol {q} = (\boldsymbol {u}, \, p_x, \,c)$
, the matching conditions (2.6) become
\begin{align} \boldsymbol {q}(\phi ^-,\,y) = \boldsymbol {q}(\phi ^+,\,y), \quad \boldsymbol {q}(-\phi, \,y) = \boldsymbol {q}(2- \phi, \,y). \end{align}
For
$x\in [-\phi, \, 2 - \phi ]$
, the surfactant-flux conditions (2.7)–(2.10) yield
\begin{align} \frac {\text{d}}{\text{d}x} \int _{0}^{2} \left ( u c - \frac {\epsilon c_x}{{\textit {Pe}}} \right ) \, \text{d} y - \frac {{\textit {Da}}}{\epsilon {\textit {Pe}}}(\Gamma - c(x,\, 0)) &= 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}
\begin{align} \int _{0}^{2} \left ( u c - \frac {\epsilon c_x}{{\textit {Pe}}}\right ) \, \text{d}y + \frac {{\textit {Da}}}{{\textit {Bi}} {\textit {Pe}}} \left ( u(x,\, 0) \Gamma - \frac {\epsilon \Gamma _{x}}{{\textit {Pe}}_I}\right ) &= 1 \quad \text{in} \quad \mathcal{D}_1, \end{align}
\begin{align} \int _{ - \phi }^{\phi } ( \Gamma - c(x,\, 0) ) \, \text{d}x &= 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}
\begin{align} \int _{0}^{2} \left (u c - \frac {\epsilon c_x}{{\textit {Pe}}} \right ) \, \text{d}y &= 1 \quad \text{in} \quad \mathcal{D}_2, \\[2pt] \nonumber \end{align}
and in
$\mathcal{D}_1$
and
$\mathcal{D}_2$
, the volume-flux condition (2.11) gives
\begin{align} \int _{0}^{2} u \, \text{d}y = 1. \end{align}
Finally, the drag reduction (2.12) becomes
\begin{align} {DR} = \frac {\Delta p_I - \Delta p}{\Delta p_I - \Delta p_U}, \end{align}
where
$\Delta {p} \equiv \langle p\rangle (-\phi ) - \langle p\rangle (2-\phi )$
and
$\langle \cdot \rangle \equiv \int _{0}^{2} \cdot \, \text{d} y / 2$
.
In the limit
$\epsilon \ll 1$
, we proceed to solve the leading-order boundary-value problem and flux constraints defined in (2.15)–(2.21) for
$u$
,
$v$
,
$p$
,
$c$
and
$\Gamma$
using asymptotic and numerical methods. Subsequently, we evaluate the drag reduction (2.22) and investigate its dependence on the seven dimensionless groups (
$\epsilon$
,
$\textit {Pe}$
,
$\textit {Ma}$
,
${\textit {Pe}}_I$
,
$\textit {Bi}$
,
$\textit {Da}$
and
$\phi$
) that characterise the geometry, flow, liquid and surfactant. Moreover, full solutions of (2.15)–(2.21) obtained as COMSOL simulations are presented in § 3.
2.3. The long-wave limit of the 2-D model at high Péclet numbers
We investigate distinguished limits of (2.15)–(2.21) to uncover different physical balances. We assume that
$1/{\textit {Pe}} = O(\epsilon )$
,
$1/{\textit {Pe}}_I = O(\epsilon )$
,
${\textit {Bi}} = O(\epsilon )$
,
${\textit {Da}} = O(1)$
and
${\textit {Ma}} = O(1/\epsilon )$
as
$\epsilon \rightarrow 0$
, denoting this as the ‘2-D long-wave model’ in which streamwise and cross-channel concentration gradients are comparable at leading order. We rescale
$1/{\textit {Pe}} = \epsilon / \mathscr{P}$
,
$1/{\textit {Pe}}_I = \epsilon /\mathscr{P}_I$
,
${\textit {Bi}} = \epsilon \mathscr{B}$
and
${\textit {Ma}} = \mathscr{M}/\epsilon$
, assuming that
$\mathscr{P}$
,
$\mathscr{P}_I$
,
$\mathscr{B}$
and
$\mathscr{M}$
remain
$O(1)$
as
$\epsilon \rightarrow 0$
. We substitute the expansions
\begin{align} (u,\,v, \,p,\, c,\, \Gamma ) = (u_0,\,v_0, \,p_0,\, c_0, \,\Gamma _0) + \epsilon ^2 (u_1,\,v_1, \,p_1,\, c_1, \,\Gamma _1) + \cdots \end{align}
into the governing equations (2.15)–(2.22) and take the leading-order approximation.
In
$\mathcal{D}_1$
and
$\mathcal{D}_2$
, flow is driven by a streamwise pressure gradient, and wall-normal diffusion is comparable to advection. The bulk equations (2.15) become
\begin{align} u_{0x} + v_{0y} = 0, \quad u_{0yy} - p_{0x} = 0, \quad p_{0y} = 0, \quad c_{0yy}/\mathscr{P} - u_0 c_{0x} - v_0 c_{0y} = 0. \end{align}
For
$x\in [-\phi, \, \phi ]$
and
$y=0$
, the interfacial flow is inhibited by the surfactant gradient and exchange is comparable to advection. The interface conditions (2.16) yield
\begin{align} u_{0y} - \mathscr{M} \Gamma _{0x} = 0, \quad v_0 =0, \quad c_{0y} - {\textit {Da}} (c_0 -\Gamma _0) =0, \nonumber \\ \mathscr{B} ( c_0 - \Gamma _0) - (u_0 \Gamma _0)_{x} = 0. \end{align}
For
$x\in [\phi, \, 2 - \phi ]$
and
$y=0$
(and for
$x\in [-\phi, \, 2-\phi ]$
and
$y=2$
), the solid wall boundary conditions (2.18) become
\begin{align} u_0 = 0, \quad v_0 = 0, \quad \ c_{0y} =0. \end{align}
Short 2-D Stokes-flow regions arise at the junctions between
$\mathcal{D}_1$
and
$\mathcal{D}_2$
. In the present long-wave theory, in which surface diffusion appears at the next order, we impose the no-flux condition (2.17) at
$x = \pm \phi$
and
$y=0$
:
\begin{align} u_0 \Gamma _0 = 0. \end{align}
Continuity of bulk concentration (2.19) between
$\mathcal{D}_1$
and
$\mathcal{D}_2$
requires
\begin{align} c_0(\phi ^-,\,y) = c_0(\phi ^+,\,y), \quad c_0(-\phi, \,y) = c_0(2- \phi, \,y). \end{align}
In § 3, comparison with COMSOL simulations will allow us to investigate the effect of these Stokes-flow regions. Streamwise surfactant transport in the bulk and interface is dominated by advection and exchange, so (2.20a ) becomes
\begin{align} \frac {\text{d}}{\text{d}x} \int _{0}^{2} u_0 c_0 \, \text{d} y - \frac {Da}{\mathscr{P}}(\Gamma _0 - c_0(x,\, 0)) &= 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}
\begin{align} \int _{0}^{2} u_0 c_0 \, \text{d}y + \frac {{\textit {Da}}}{\mathscr{B}\mathscr{P}} \left ( u_0(x,\, 0) \Gamma _0 \right ) &= 1 \quad \text{in} \quad \mathcal{D}_1, \end{align}
\begin{align} \int _{ - \phi }^{\phi } ( \Gamma _0 - c_0(x,\, 0) ) \, \text{d}x &= 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}
\begin{align} \int _{0}^{2} u_0 c_0 \, \text{d}y &= 1 \quad \text{in} \quad \mathcal{D}_2. \end{align}
In
$\mathcal{D}_1$
and
$\mathcal{D}_2$
, the volume-flux condition (2.21) gives
\begin{align} \int _{0}^{2} u_0 \, \text{d}y = 1. \end{align}
The leading-order drag reduction (2.22) is given by
\begin{align} {DR}_0 = \frac {\Delta p_I - \Delta p_0}{\Delta p_I - \Delta p_U}. \end{align}
The leading-order pressure field simplifies to
$p_0 = p_0(x)$
using the wall-normal equation (2.24c
). The streamwise velocity field is driven by the pressure gradient
$p_{0x}$
and is inhibited by the surfactant gradient
$\Gamma _{0x}$
. Using linear superposition, we can write
\begin{align} u_0 = \tilde {U} p_{0x} + \mathscr{M} \bar {U}\Gamma _{0x} \quad \text{in} \quad \mathcal{D}_1, \quad u_0 = \breve {U} p_{0x} \quad \text{in} \quad \mathcal{D}_2, \end{align}
where
$\tilde {U} \equiv y^2/2 - 2$
,
$\bar {U} \equiv y - 2$
and
$\breve {U} \equiv y^2/2 - y$
are velocity contributions satisfying the following boundary-value problems:
\begin{align} \tilde {U}_{yy} = 1, \quad \text{subject to} \quad \tilde {U}_y(0)=0, \quad \tilde {U}(2)=0; \\[-9pt] \nonumber \end{align}
\begin{align} \bar {U}_{yy} = 0, \quad \text{subject to} \quad \bar {U}_y(0)= 1, \quad \bar {U}(2) = 0; \\[-9pt] \nonumber \end{align}
\begin{align} \breve {U}_{yy} = 1, \quad \text{subject to} \quad \breve {U}(0)= 0, \ \quad \breve {U}(2)=0. \\[-9pt] \nonumber \end{align}
Substituting (2.32) into (2.30), the volume flux constraint gives
\begin{align} \tilde {Q} p_{0x} + \mathscr{M} \bar {Q} \Gamma _{0x} =1 \quad \text{in} \quad \mathcal{D}_1, \quad \breve {Q} p_{0x} = 1 \quad \text{in} \quad \mathcal{D}_2, \end{align}
where
$\tilde {Q} \equiv -8/3$
,
$\bar {Q} \equiv -2$
and
$\breve {Q} \equiv -2/3$
. Using the continuity equation (2.24a
), the wall-normal velocity field is forced by
$p_{0xx}$
and
$\Gamma _{0xx}$
, with velocity contributions that complement those given in (2.32)–(2.34):
\begin{align} v_0 = \tilde {V} p_{0xx} + \mathscr{M} \bar {V} \Gamma _{0xx} \quad \text{in} \quad \mathcal{D}_1, \quad v_0 = 0 \quad \text{in} \quad \mathcal{D}_2, \end{align}
where
$\tilde {V} \equiv -y^3/6 + 2 y$
and
$\bar {V} \equiv -y^2/2 + 2 y$
(so that
$\tilde {V}_{y} = -\tilde {U}$
,
$\tilde {V}(0)=0$
,
$\bar {V}_{y} = -\bar {U}$
,
$\bar {V}(0)= 0$
). The no-penetration condition (2.26b
) is satisfied, as
$v_0(2) = \tilde {V}(2) p_{0xx} + \mathscr{M}\bar {V}(2)\Gamma _{0xx} = -\tilde {Q} p_{0xx} - \mathscr{M}\bar {Q}\Gamma _{0xx} = 0$
, using the volume-flux condition (2.36a
).
Table 1. A summary of the transport coefficients in the 2-D long-wave model, (2.38)–(2.42), with their definition and physical interpretation. The transport coefficients
$\alpha$
,
$\beta$
,
$\gamma$
,
$\delta$
and
$\nu$
are written in terms of
$\textit {Pe}$
,
$\textit {Da}$
,
$\textit {Bi}$
,
$\textit {Ma}$
,
${\textit {Pe}}_I$
and
$\epsilon$
defined in § 2.2. Compared with the three-dimensional transport coefficients in Tomlinson et al. (Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a
),
$\alpha _{{3D}}\propto \epsilon ^2\alpha$
,
$\beta _{{3D}}\propto \beta$
,
$\gamma _{{3D}}\propto \gamma$
,
$\delta _{{3D}}\propto \epsilon ^2\delta$
and
$\nu _{{3D}}\propto \epsilon ^2\nu$
.
We substitute (2.32) and (2.37) into the bulk and interfacial surfactant equations (2.24)–(2.29) to derive the 2-D long-wave model. To facilitate its numerical computation, we retain
$O(\epsilon ^2)$
bulk and interfacial streamwise diffusion operators to make the problem more elliptic, simplifying the continuity conditions and smoothing the flow between
$\mathcal{D}_1$
and
$\mathcal{D}_2$
, later seeking the limit in which these effects are weak. The bulk surfactant equation (2.24c
) gives
\begin{align} \alpha (\epsilon ^2 c_{0xx} + c_{0yy}) - \bigg (\frac {3}{4} - \frac {3y^2}{16}\bigg ) & c_{0x} \nonumber \\ - \frac {\gamma }{\beta } \bigg (-\frac {1}{2} + y - \frac {3y^2}{8}\bigg )& \Gamma _{0x} c_{0x} - \frac {\gamma }{\beta } \bigg (\frac {y}{2} - \frac {y^2}{2} + \frac {y^3}{8}\bigg ) \Gamma _{0xx} c_{0y} = 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}
\begin{align} &\alpha (\epsilon ^2 c_{0xx} + c_{0yy}) - \bigg (\frac {3y}{2} - \frac {3y^2}{4}\bigg ) c_{0x} = 0 \quad \text{in} \quad \mathcal{D}_2, \\[6pt] \nonumber \end{align}
where
$\alpha$
,
$\beta$
and
$\gamma$
are given in table 1. This table defines the five primary parameter combinations that we use below and gives their physical interpretation. For
$x\in [-\phi, \, \phi ]$
and
$y=0$
, the boundary condition specifying continuity of surfactant flux and the interfacial surfactant equation (2.25c
,
d
) become
\begin{align} c_{0y} - \frac {\nu }{\alpha }(c_{0} - \Gamma _0) = 0, \quad \nu ( c_{0} - \Gamma _0) - \frac {3}{4}\beta \Gamma _{0x} + \frac {1}{2}\gamma (\Gamma _{0x}\Gamma _0)_x + \epsilon ^2 \delta \Gamma _{0xx} = 0, \end{align}
where
$\nu$
and
$\delta$
are given in table 1. At
$x= \pm \phi$
and
$y=0$
, no flux of interfacial surfactant (2.27) gives
\begin{align} \frac {3}{4}\beta \Gamma _{0} - \frac {1}{2}\gamma \Gamma _{0x}\Gamma _0 - \epsilon ^2 \delta \Gamma _{0x} = 0. \end{align}
For
$x\in [\phi, \, 2 - \phi ]$
and
$y=0$
(and for
$x\in [-\phi, \, 2-\phi ]$
and
$y=2$
), no flux of bulk surfactant (2.26c
) becomes
\begin{align} c_{0y} =0. \end{align}
The surfactant-flux conditions (2.29a , b , d ) become
\begin{align} \frac {\text{d}}{\text{d}x} \int _{0}^{2} \left (\bigg (\frac {3}{4} - \frac {3y^2}{16}\bigg ) c_{0} + \frac {\gamma }{\beta } \bigg (-\frac {1}{2} + y - \frac {3y^2}{8}\bigg ) \Gamma _{0x} c_{0} - \epsilon ^2 \alpha c_{0x} \right )\text{d} y \nonumber \\ - \nu (\Gamma _0 - c_0(x,\, 0)) &= 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}
\begin{align} \int _{0}^{2} \left (\bigg (\frac {3}{4} - \frac {3y^2}{16}\bigg ) c_{0} + \frac {\gamma }{\beta } \bigg (-\frac {1}{2} + y - \frac {3y^2}{8}\bigg ) \Gamma _{0x} c_{0} - \epsilon ^2 \alpha c_{0x} \right )\text{d}y \nonumber \\ + \frac {3}{4}\beta \Gamma _{0} - \frac {1}{2}\gamma \Gamma _{0x}\Gamma _0 - \epsilon ^2 \delta \Gamma _{0x} &= 1 \quad \text{in} \quad \mathcal{D}_1, \end{align}
\begin{align} \int _{0}^{2} \left ( \bigg (\frac {3y}{2} - \frac {3y^2}{4}\bigg ) c_0 - \epsilon ^2 \alpha c_{0x} \right ) \, \text{d}y &= 1 \quad \text{in} \quad \mathcal{D}_2, \\[6pt] \nonumber \end{align}
and in
$\mathcal{D}_1 \cup \mathcal{D}_2$
, the volume-flux condition (2.30) becomes
\begin{align} \tilde {Q} p_{0x} + (\gamma / \beta ) \bar {Q} \Gamma _{0x} = 1 \quad \text{in} \quad \mathcal{D}_1, \quad \breve {Q} p_{0x} = 1 \quad \text{in} \quad \mathcal{D}_2. \end{align}
To summarise, the seven-parameter Stokes-flow system (2.15)–(2.21) involving five variables
$(\hat {u}, \hat {v}, \, \hat {p}, \, \hat {c}, \, \hat {\Gamma })$
has been reduced to the seven-parameter long-wave system (2.38)–(2.42) involving two variables
$(\hat {c}, \, \hat {\Gamma })$
. Further, by neglecting streamwise diffusion in the bulk,
$-\epsilon ^2 \alpha c_{0x}$
, and at the interface,
$-\epsilon ^2 \delta \Gamma _{0x}$
, the system reduces to a five-parameter system (
$\alpha$
,
$\beta$
,
$\gamma$
,
$\nu$
,
$\phi$
). This long-wave model disregards short (but passive) inner regions associated with the small parameters
$\epsilon$
and
$\delta$
.
The numerical solution of (2.38)–(2.42) (Appendix A) is generally less expensive (in terms of memory and runtime) than COMSOL simulations. The numerical method iteratively solves linearised bulk- and interfacial-concentration problems until convergence is achieved. Convergence in the nonlinear problem can be slow, necessitating small tolerances. Accurate resolution of singularities at
$x=\pm \phi$
is essential to achieving convergence; accordingly, we employ domain decomposition and Chebyshev collocation techniques that cluster nodes around
$y=0,\,2$
and
$x = -\phi, \, \phi, \,2-\phi$
. Once
$c_0$
and
$\Gamma _0$
are determined,
$u_0$
and
$v_0$
can be calculated using (2.32) and (2.37), respectively.
Following Tomlinson et al. (Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a
), we express the leading-order drag reduction (
${DR}_0$
) in terms of the transport parameters and geometry (
$\beta$
,
$\gamma$
and
$\phi$
) as follows. When the interface is unaffected by surfactant,
$\Delta p_U = -2\phi /\tilde {Q} - 2(1-\phi )/\breve {Q}$
. When the interface is immobilised by surfactant,
$\Delta p_I = -2/\breve {Q}$
. Between these two limits,
$\Delta p_0 = -2\phi /\tilde {Q} - 2(1-\phi )/\breve {Q} + \gamma \bar {Q} \Delta \Gamma _0 / (\beta \tilde {Q})$
, where
$\Delta \Gamma _0 = \Gamma _0(\phi ) - \Gamma _0(-\phi )$
. Substituting these expressions into (2.31) gives the leading-order drag reduction as
\begin{align} {DR}_0 = 1 - \frac {\gamma \Delta \Gamma _0}{3\phi \beta }. \end{align}
When evaluated using COMSOL, we write
${DR}_0 = {DR}_{NS}$
. An alternative measure used commonly in the SHS literature is the effective slip length,
$\lambda _0$
, which can be evaluated from
${DR}_0$
using (Tomlinson et al. Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a
)
\begin{align} \lambda _0 = \frac {{DR}_0(\Delta p_I - \Delta p_U)}{\Delta p_I[\Delta p_U {DR}_0 + \Delta p_I(1 - {DR}_0)]}. \end{align}
Before presenting the results of the 2-D long-wave model, it is helpful to recall key features of the problem at low Péclet numbers, when a 1-D description applies.
Figure 2. Maps of (
$\alpha$
,
$\gamma$
)-parameter space for
$\epsilon \ll 1$
, using parameters given in table 1, showing the magnitude of drag reduction (
$DR_0$
) across different asymptotic regions. Existing predictions from the 1-D long-wave model of Tomlinson et al. (Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a
) are combined with predictions derived in Appendix C from the 2-D long-wave model for (a) strong exchange and (b) weak exchange. The Marangoni-dominated region (
$M$
) incorporates subregions
$M^{{1D}}$
,
$M^{{1D}}_E$
,
$M^{{2D}}$
and
$M^{{2D}}_E$
, for which
${DR}_0 \ll 1$
; the advection-dominated region (
$A$
) incorporates subregions
$A^{{1D}}$
,
$A^{{1D}}_E$
,
$A^{{2D}}$
and
$A^{{2D}}_E$
, for which
$1-{DR}_0 \ll 1$
; and the diffusion-dominated region (
$D$
) incorporates subregions
$D^{{1D}}$
and
$D^{{1D}}_E$
, for which
$1-{DR}_0 \ll 1$
. Here it is assumed that
$\alpha \sim \delta$
and
$\beta \sim 1$
.
2.4. The 1-D long-wave model
Figure 2 shows maps of
$(\alpha, \,\gamma )$
-parameter space, illustrating how the present 2-D model relates to existing 1-D results. In the strong-cross-channel-diffusion limit, for which
$\alpha \gg 1$
(with
$\epsilon ^2\alpha =O(1)$
and
$\epsilon ^2\delta =O(1)$
), cross-channel concentration gradients are small. Expanding the bulk and interfacial concentration fields
\begin{align} c_0 = \bar {c}_0(x) + \bar {c}_1(x,\,y)/\alpha + \cdots, \quad \Gamma _0 = \bar {\Gamma }_0(x) + \bar {\Gamma }_1(x)/\alpha + \cdots, \end{align}
equations (2.38)–(2.42) reduce to the ‘1-D long-wave model’, given by
\begin{align} \bar {c}_{0x} - 2 \epsilon ^2\alpha \bar {c}_{0xx} - \nu (\bar {\Gamma }_0 - \bar {c}_0) &= 0 \quad \text{in} \ \mathcal{D}_1, \end{align}
\begin{align} \frac {3}{4} \beta \bar {\Gamma }_{0x} - \frac {1}{2} \gamma (\bar {\Gamma }_{0x} \bar {\Gamma }_0)_x - \epsilon ^2\delta \bar {\Gamma }_{0xx} + \nu (\bar {\Gamma }_0 - \bar {c}_0) &= 0 \quad \text{in} \ \mathcal{D}_1, \end{align}
\begin{align} \bar {c}_0 - 2 \epsilon ^2\alpha \bar {c}_{0x} +\frac {3}{4} \beta \bar {\Gamma }_0 - \frac {1}{2} \gamma \bar {\Gamma }_{0x} \bar {\Gamma }_0 - \epsilon ^2\delta \bar {\Gamma }_{0x} &= 1 \quad \text{in} \ \mathcal{D}_1, \end{align}
\begin{align} \bar {c}_0 - 2 \epsilon ^2\alpha \bar {c}_{0x} &= 1 \quad \text{in} \ \mathcal{D}_2. \\[6pt] \nonumber \end{align}
The system in (2.47) is solved subject to continuity of bulk surfactant, continuity of bulk-surfactant flux and no flux of interfacial surfactant:
\begin{align} \bar {c}_0(\phi ^-) = \bar {c}_0(\phi ^+), \\[-33pt] \nonumber \end{align}
\begin{align} \bar {c}_0(-\phi ) = \bar {c}_0(2-\phi ), \\[-33pt] \nonumber \end{align}
\begin{align} \bar {c}_0(\phi ^-) - 2\epsilon ^2\alpha \bar {c}_{0x}(\phi ^-) = \bar {c}_0(\phi ^+) - 2\epsilon ^2\alpha \bar {c}_{0x}(\phi ^+), \\[-33pt] \nonumber \end{align}
\begin{align} \bar {c}_0(-\phi ) - 2\epsilon ^2\alpha \bar {c}_{0x}(-\phi ) = \bar {c}_0(2-\phi ) - 2\epsilon ^2\alpha \bar {c}_{0x}(2-\phi ), \\[-33pt] \nonumber \end{align}
\begin{align} \frac {3}{4} \beta \bar {\Gamma }_0(\pm \phi ) - \frac {1}{2} \gamma \bar {\Gamma }_0(\pm \phi ) \bar {\Gamma }_{0x}(\pm \phi ) - \epsilon ^2\delta \bar {\Gamma }_{0x}(\pm \phi ) = 0. \end{align}
The system of ordinary differential equations in (2.47)–(2.48) can be solved numerically using the method outlined in Appendix A of Tomlinson et al. (Reference Tomlinson, Peaudecerf, Temprano-Coleto, Gibou, Luzzatto-Fegiz, Jensen and Landel2023). The drag reduction can be evaluated using (2.44). As illustrated in figure 2, the 1-D long-wave model (2.47)–(2.48) exhibits multiple asymptotic regimes in the
$(\alpha, \,\gamma )$
-parameter space, which provides useful context for the results to follow. We highlight regions dominated by Marangoni effects (
$M$
), interfacial advection (
$A$
) and interfacial diffusion (
$D$
); a subscript
$E$
denotes cases in which
$\Gamma _0(x)$
and
$c_0(x,0)$
can be out of equilibrium. We summarise key asymptotic results from Tomlinson et al. (Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a
) that apply in the strong-exchange (
$\bar {c}_0 \approx \bar {\Gamma }_0$
, large
$\nu$
) and weak-exchange (
$\bar {c}_0 \neq \bar {\Gamma }_0$
, small
$\nu$
) limits in Appendix B, in regions labelled with superscript
$1D$
. We now proceed to investigate the emergence of 2-D flow structures as bulk diffusion weakens, i.e. for
$\alpha$
small (
$\ll 1$
) or at intermediate values
$\sim 1$
.
3. Results
In § 3.1, we investigate the leading-order drag reduction (
${DR}_0$
), interfacial and bulk surfactant concentrations (
$\Gamma _0$
and
$c_0$
), as well as the streamwise and wall-normal velocity components (
$u_0$
and
$v_0$
), for strong bulk–surface exchange (
$\nu = 100$
). We vary the bulk diffusion
$(\alpha )$
and surfactant strength (
$\gamma$
), exploring both strong- and weak-cross-channel-diffusion regimes by solving the 2-D long-wave model (2.38)–(2.42). We focus primarily on the weak-cross-channel-diffusion regime (
$\alpha = \delta \ll 1$
), as the strong-cross-channel-diffusion regime (§ 2.4) was addressed in Tomlinson et al. (Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a
). For simplicity, we constrain the partition coefficient (
$\beta$
) and interfacial diffusion (
$\delta$
), assuming that
$\alpha = \delta$
,
$\beta = 1$
and
$\phi =0.5$
throughout. However, asymptotic solutions are derived in Appendix C for generic variables
$\alpha$
,
$\beta$
,
$\delta$
and
$\phi$
, which we use to validate numerical simulations and identify dominant physical mechanisms. In § 3.2, we consider the weak-exchange limit (
$\nu = 0.01$
). As illustrated in figure 2, in both strong- and weak-exchange limits, we identify regions of parameter space dominated by Marangoni effects (
$M$
), advection (
$A$
) and diffusion (
$D$
). Additionally, we compare numerical solutions of the 2-D long-wave model (2.38)–(2.42) with numerical solutions of the Stokes and advection–diffusion equations (2.15)–(2.21) in the weak-cross-channel-diffusion limit.
3.1. Strong exchange
Figure 3 summarises the surfactant concentration profiles and drag reduction in the strong-exchange limit, computed using the 2-D long-wave model, (2.38)–(2.42). For large values of
$\gamma$
, with strong Marangoni effect, the interface is almost immobilised (figure 3
a,b), the
$\Gamma _0$
profile is almost linear and the bulk surfactant concentration transitions from a 2-D (figure 3
a) to a 1-D (figure 3
b) distribution as
$\alpha$
increases. Surfactant adsorbs onto the interface across the upper half of the plastron and desorbs from it in the lower half. In contrast, for small
$\gamma$
with weaker Marangoni effects, interfacial surfactant accumulates near the downstream end of the plastron, lengthening the adsorption region and compressing the desorption region (figure 3
f).
Figure 3. The leading-order drag reduction (
${DR}_0$
), bulk (
$c_0$
) and interfacial (
$\Gamma _0$
) surfactant concentration fields for
$\beta = 1$
,
$\epsilon = 0.1$
,
$\nu = 100$
and
$\phi = 0.5$
, computed using (2.38)–(2.42) when bulk–surface exchange is strong. In the Marangoni-dominated (
$M$
) region, the SHS is no-slip (
${DR}_0\ll 1$
), and in the advection-dominated (
$A$
) and diffusion-dominated (
$D$
) regions, the interface is shear-free (
${DR}_0\approx 1$
). Fields
$c_0$
,
$\Gamma _0$
,
$\langle c_0 \rangle$
and
$c_0(x, \, 0)$
are plotted in (a)
$M^{{2D}}$
and (b)
$M^{{1D}}$
. (c) Contours of
${DR}_0$
, (d) plots of
${DR}_0$
for different
$\gamma$
and (e) plots of
${DR}_0$
for different
$\alpha$
, where (c,d,e) are compared with asymptotic predictions (B1b
) in
$M^{{1D}}$
, (B2b
) in
$D^{{1D}}$
, (C22b
) in
$M^{{2D}}$
and (3.2b
) in
$A^{{2D}}$
. The dashed line in (e) represents the largest
$\gamma$
for which
$\Gamma _0(-\phi )=0$
when
$\alpha \ll 1$
. Fields
$c_0$
,
$\Gamma _0$
,
$\langle c_0 \rangle$
and
$c_0(x, \, 0)$
are plotted in (f)
$A^{{2D}}$
(notice the quasi-stagnant-cap profile) and (g)
$D^{{1D}}$
, where
$x=x_0$
is plotted using (C55).
Figure 3(c–e) illustrates the variation of
${DR}_0$
with
$\alpha$
and
$\gamma$
. In the present strong-exchange limit, drag reduction decreases as
$\alpha$
decreases and the bulk concentration becomes more 2-D. Reduced bulk diffusion promotes the formation of bulk-concentration boundary layers, while slowing fluxes between the bulk and interface. We characterise this transition by comparing the established large-
$\alpha$
asymptotes (in regions
$M^{{1D}}$
and
$D^{{1D}}$
of the
$(\alpha, \gamma )$
plane; see (B1), (B2)) with new small-
$\alpha$
limits (derived below) in regions
$M^{{2D}}$
and
$A^{{2D}}$
. We use these limits to characterise the immobilisation of the interface as
$\gamma$
increases, highlighting in particular the transition from region
$A^{{2D}}$
to region
$M^{{2D}}$
for small
$\alpha$
. The central asymptotic subregions
$M^{{1D}}$
and
$A^{{1D}}$
for
$1\ll \alpha \ll 1/\epsilon ^2$
in figure 2(a), which formally provide a bridge between the 1-D and 2-D long-wave models, are smoothed out by transitional effects associated with the nominally small geometric parameter
$\epsilon$
taking the value 0.1 in figure 3. We now discuss these new regions, starting with
$M^{{2D}}$
.
Figure 4. Schematics of the asymptotic structure of the bulk-concentration boundary layer. Weak diffusion ensures that
$c_0$
is approximately uniform in the core of the channel, varying primarily in a thin concentration boundary layer near the SHS. Blue (pink) regions illustrate regions where surfactant is drawn from (released into) the bulk onto (from) the interface.
$(\textit {a})$
The bulk-concentration boundary layer when Marangoni effects are strong (region
$M^{{2D}}$
) and the surfactant distribution almost immobilises the entire interface.
$(\textit {b})$
The bulk-concentration boundary layer when Marangoni effects are weak (region
$A^{{2D}}$
), creating a slip region with low surfactant concentration upstream of a quasi-stagnant region in which interfacial surfactant accumulates. Details of each asymptotic region are provided in Appendix C.
Region
$M^{{2D}}$
(figures 2
a and 3
a) emerges when
$\alpha \ll 1$
and
$\gamma$
is large. In this region, Marangoni effects dominate, with weak bulk diffusion resulting in a bulk-concentration boundary layer near
$y=0$
. We provide a scaling argument for
${DR}_0$
as follows, with a detailed derivation given in Appendix C.1 and a schematic of the asymptotic structure in figure 4(a). At the interface,
$\Gamma _0$
is in equilibrium with
$c_0$
, and the interface is nearly immobile, exhibiting a linear surfactant distribution with a slope of size
$1/\gamma$
. The bulk concentration is close to unity, and
$\Gamma _0\lt 1$
(
$\Gamma _0\gt 1$
) in the upstream (downstream) half of the plastron, leading to adsorption (desorption). Perturbations to the bulk-concentration boundary layer are driven by the linear surface concentration distribution, giving the boundary layer a self-similar structure. The thickness of the boundary layer is
$\alpha ^{1/3} \ll 1$
(characteristic of a shear flow), resulting in fluxes of size
$\alpha c_{0y} \sim \alpha ^{2/3}$
onto and off the interface. Adsorption and desorption are accommodated by weak stretching and compression of the interface, generating smaller contributions to
$\Gamma _0$
of size
$\alpha ^{2/3}/\gamma ^2$
and to
${DR}_0$
of size
$\alpha ^{2/3}/\gamma$
, using (2.44). The constraint
$\alpha \ll 1$
ensures that the boundary layer is thin and we require
$\gamma \gg 1$
and
$\alpha ^{2/3} \ll \gamma$
for the correction to
$\Gamma _0$
to be small (although a more precise condition will emerge below). The leading-order surfactant distribution and drag reduction are given by
\begin{align} c_0(x, \, 0) \approx \Gamma _0 \approx 1 + \frac {3\beta }{2\gamma }\left (x - \frac {\phi }{5}\right ), \quad {DR}_0 \approx \frac {m_1\alpha ^{2/3}\phi ^{5/3}}{\gamma }, \end{align}
where the coefficient
$m_1 \approx 0.79$
is given explicitly in Appendix C.1. The range of validity of the leading-order solution in (3.1) is extended from
$\gamma \gg 1$
to
$\gamma \sim 1$
in (C22), where the extension must be evaluated numerically. Therefore, we give the simpler formula above, but compare (C22) with numerical results from the 2-D long-wave model at small
$\alpha$
in figure 3(c–e), which capture numerical results successfully. In particular, decreasing
$\gamma$
leads to a steepening of the interfacial surfactant gradient until
$\Gamma _0$
approaches zero at the upstream contact line. This transition provides a lower bound on region
$M^{{2D}}$
, namely
$\Gamma _0(-\phi )=0$
when
$\gamma = 9\phi \beta /5$
(see vertical dashed lines in figure 3
e). Hence, we plot asymptotic predictions up to this limit in figure 3(e).
When
$\alpha \ll 1$
and
$\gamma$
is sufficiently small, we transition into region
$A^{{2D}}$
(figures 2
a and 3
f), where advection dominates Marangoni effects at the interface, while interacting with a strongly coupled bulk-concentration boundary layer. We distinguish two primary regions along the interface: (i) a slip region,
$x\in [-\phi, \,x_0)$
(for some
$x_0$
), at the upstream end of the interface, where Marangoni effects are weak and surfactant adsorbs from the bulk; and (ii) a quasi-stagnant region,
$x\in (x_0,\,\phi ]$
, at the downstream end of the interface, where surfactant gradients decelerate the interface through Marangoni effects without completely immobilising the interface (in contrast with a classic stagnant cap), while desorbing the accumulated surfactant. We now provide a scaling argument for
${DR}_0$
that is supported by a detailed derivation in Appendix C.2; a schematic of the asymptotic structure is given in figure 4(b). In the slip region, of length
$x_0+\phi =O(1)$
,
$c_0$
varies by
$O(1)$
across a boundary-layer thickness
$y\sim \alpha ^{1/2}$
, yielding a flux per unit length of size
$\alpha c_{0y}\sim \alpha ^{1/2}$
onto the interface. This allows
$\Gamma _0$
to grow along the slip region from zero to
$O(\alpha ^{1/2})$
. Although
$c_0$
is effectively zero at the interface at leading order, it is coupled to
$\Gamma _0$
through an
$O(\alpha ^{1/2})$
correction. The slip region therefore delivers an interfacial flux of size
$O(\alpha ^{1/2})$
at
$x\approx x_0$
, which is carried through short regions across which the interface decelerates rapidly at the leading edge of the quasi-stagnant region (i.e. at
$x\approx x_0$
). The sizes of the short regions are indicated in figure 4(b). Across the quasi-stagnant region, of length
$L=\phi -x_0$
(to be determined),
$\Gamma _0$
is approximately linear with a slope of size
$1/\gamma$
, setting
$c_0$
and
$\Gamma _0$
to be of size
$L/\gamma$
. In the bulk, diffusion balances shear-flow transport over a vertical length
$y \sim (\alpha L)^{1/3}$
and therefore
$\alpha c_{0y} \sim (\alpha L)^{2/3}/\gamma$
. Integrated over
$L$
, this accommodates the
$O(\alpha ^{1/2})$
interfacial flux, resulting in
$L \sim \gamma ^{3/5}/\alpha ^{1/10}$
. Consequently, the drag reduction arising from the surfactant adsorbed and desorbed in the slip and quasi-stagnant regions is proportional to
${DR}_0 \sim \gamma \Delta \Gamma \sim L$
.
In Appendix C.2, we present a more detailed analysis showing that
\begin{align} c_0(x, \, 0) \approx \Gamma _0 \approx \begin{cases} \dfrac {4(\alpha (x+\phi ))^{1/2}}{\sqrt {3\pi }\beta } \quad \text{for} \quad x\in [-\phi, \, x_0), \\ \dfrac {3\beta }{2\gamma }(x-x_0) \quad \text{for} \quad x\in (x_0,\,\phi ], \end{cases} \quad {DR}_0 = 1 - \dfrac {a_1 \gamma ^{3/5}}{\alpha ^{1/10}},\nonumber \\ { u_0(x, \, 0) \approx \begin{cases}3/4 \quad \text{for} \quad x\in [-\phi, x_0), \\ \dfrac {3}{4} \left [G\left (\dfrac {x-x_0}{\gamma \alpha ^{1/2}}\right )\right ]^{-1} \text{for} \quad x\approx x_0, \\[6pt] \dfrac {3\alpha ^{2/3} C^{\prime}(0)}{5\beta } \Bigg ((x-x_0)^{2/3} - \dfrac {(\phi -x_0)^{5/3}}{x-x_0}\Bigg ) \quad \text{for} \ x\in (x_0,\,\phi ]. \end{cases} } \\[-6pt] \nonumber \end{align}
The coefficient
$a_1 = (-5\lambda /(6C^{\prime}(0)))^{3/5}/(2\phi )$
in (3.2b
) where
$\lambda = (16(x_0+\phi )/ (3\pi \beta ^2))^{1/2}$
;
$C^{\prime}(0) \approx -0.92$
in (3.2c
) is given explicitly in Appendix C.2; the function
$G$
in (3.2
$c$
) is given implicitly in (C41). The length of the quasi-stagnant cap is
\begin{align} L=\phi -x_0= \frac {2 \phi a_1 \gamma ^{3/5}}{\alpha ^{1/10}}. \end{align}
In (3.2),
$x\in [-\phi, x_0)$
,
$x\approx x_0$
and
$x\in (x_0,\phi ]$
are shorthand for the slip, deceleration and quasi-stagnant regions shown in figure 4(b); the leading-order approximations in (3.2) match together as explained in Appendix C. Equation (3.2c
) highlights the role of interfacial compression in the quasi-stagnant region, accommodating desorption.
Asymptotes for
${DR}_0$
,
$c_0(x,\,0)$
and
$\Gamma _0$
in figure 3(c–f) are evaluated numerically for a given
$\alpha$
and
$\gamma$
and capture the behaviour exhibited by the numerical simulations of the 2-D long-wave model. Increasing
$\gamma$
from low values causes lengthening of the quasi-stagnant region, pushing
$x_0$
towards the upstream contact line and reducing
$DR_0$
. Figure 3(e) shows how (3.2b
) applies almost until
$x_0$
approaches
$-\phi$
. For small
$\alpha$
,
$c_0(-\phi, \, 0) \approx \Gamma _0(-\phi ) \approx 0$
for
$\gamma \ll 1$
, causing the
$c_0$
field to develop a singular first derivative at the upstream contact point. This concentration-gradient singularity significantly impacts the numerical solution (in terms of both the required resolution and run time) of the 2-D long-wave model and COMSOL simulations (below).
Figure 5. The leading-order streamwise velocity field (
$u_0$
) and surfactant concentration fields (
$c_0$
and
$\Gamma _0$
) for
$\beta = 1$
,
$\nu = 100$
,
$\epsilon = 0.1$
and
$\phi = 0.5$
evaluated using (2.38)–(2.42) (lines) and COMSOL simulations (2.15)–(2.21) (symbols), when bulk–surface exchange is strong. (a,c) Plots of
$c_0$
and
$u_0$
, respectively, for
$\alpha = 0.1$
and
$\gamma = 1$
, when the flow is in the Marangoni-dominated region with weak cross-channel diffusion (
$M^{{2D}}$
). (b,d) Plots of
$c_0$
and
$u_0$
, respectively, for
$\alpha = 0.1$
and
$\gamma = 0.1$
, when the flow is in the advection-dominated region with weak cross-channel diffusion (
$A^{{2D}}$
), where
$x=x_0$
is plotted using (C55).
Further insight into the
$M^{{2D}}$
and
$A^{{2D}}$
solutions is provided in figure 5, which compares
$u_0$
,
$c_0$
and
$\Gamma _0$
computed using the 2-D long-wave model, (2.38)–(2.42), with the Stokes and advection–diffusion equations in COMSOL, (2.15)–(2.21), with the same parameters as the examples shown in figure 3(a, f). For the solution in the
$M^{{2D}}$
region,
$u_0$
exhibits a slight deviation from no slip at the liquid–gas interface (figure 5
c), generating a weak vertical flow towards the interface. The concentration profiles in the two simulations match well (figure 5
a), with the leading-order drag reduction predicted from the 2-D long-wave model,
${DR}_0 = 0.106$
, in close agreement with the COMSOL simulations,
${DR}_{NS} = 0.109$
. For the solution in the region
$A^{{2D}}$
, depicted in figure 5(b,d), the liquid–gas interface is nearly shear-free and
$u_0\approx 3/4$
in the slip region, where the interfacial surfactant concentration is low, before falling towards zero across the quasi-stagnant region, where the gradient of
$\Gamma _0$
is large. The asymptotic predictions (3.2
c) capture respectively the slight decrease of the slip velocity in the deceleration region and then the rapid fall towards zero across the quasi-stagnant region; the shift parameter in (C41) was chosen as
$X_0 = 0.5$
to minimise the error against the COMSOL simulations.
Figure 5 illustrates the distinction between a classical stagnant cap of a strictly insoluble surfactant, for which
$u_0$
vanishes abruptly in the cap region where
$\Gamma _0$
has a large gradient, and the present quasi-stagnant-cap structure, where desorption of the soluble surfactant in the bulk boundary layer is mediated by molecular diffusion, and which leads to a small but non-zero slip velocity (3.2
$c$
) across the quasi-stagnant-cap region. The predictions for the drag reduction and concentration field from the 2-D long-wave model and numerical simulations in COMSOL agree, with
${DR}_0 = 0.745$
and
${DR}_{NS} = 0.747$
. Both methods capture weak interfacial stretching (allowing adsorption) in the slip region (
$u_{0x}\gt 0$
for
$x\in [-\phi, \, x_0]$
in figure 5
d) and stronger compression (allowing desorption) in the quasi-stagnant region (
$u_{0x}\lt 0$
for
$x\in [x_0,\,\phi ]$
in figure 5
d). There are discrepancies in
$u_0$
at both the upstream and downstream contact lines. At the upstream contact line, the 2-D long-wave theory does not capture the Stokes-flow region, where
$u_0$
transitions sharply from zero to the slip velocity. At the downstream contact line, the retention of weak streamwise diffusion in the 2-D long-wave model leads to a small but non-zero value of
$u_0$
, ensuring zero surfactant flux. In contrast, the small-
$\alpha$
asymptotic solution, which neglects streamwise diffusion, correctly predicts that the slip velocity vanishes at the contact line. Since the small-
$\alpha$
solution represents a limiting case of the 2-D long-wave model as streamwise diffusion vanishes, the two solutions align in the limit
$\epsilon \to 0$
, where the 2-D long-wave model’s prediction for
$u_0(\phi, 0)$
also tends to zero at the downstream edge. Further features of the surface velocity profile at the tip of the quasi-stagnant region, which explain the dramatic thickening of the boundary layer in figures 3(f) and 4(b), are discussed in Appendix C.2.
3.2. Weak exchange
Figure 6. The leading-order drag reduction (
${DR}_0$
), bulk (
$c_0$
) and interfacial (
$\Gamma _0$
) surfactant concentration fields for
$\beta = 1$
,
$\epsilon = 0.1$
,
$\nu = 0.01$
and
$\phi = 0.5$
, computed using (2.38)–(2.42) when bulk–surface exchange is weak. In the Marangoni-dominated (
$M_E$
) region, the SHS is mostly no slip (
${DR}_0\ll 1$
), and in the advection-dominated (
$A_E$
) and diffusion-dominated (
$D_E$
) regions, the interface is mostly shear-free (
${DR}_0\approx 1$
). Fields
$c_0$
,
$\Gamma _0$
,
$\langle c_0 \rangle$
and
$c_0(x, \, 0)$
are plotted in (a)
$M^{{2D}}_E$
and (b)
$M^{{1D}}_E$
. (c) Contours of
${DR}_0$
, (d) plots of
${DR}_0$
for different
$\gamma$
and (e) plots of
${DR}_0$
for different
$\alpha$
, where (c,d,e) are compared to (B4b
) in
$M^{{1D}}_E$
, (B5b
) in
$D^{{1D}}_E$
, (B4b
) in
$M^{{2D}}_E$
and (3.4b
) in
$A^{{2D}}_E$
. The dashed line in (e) represents the largest
$\gamma$
for which
$\Gamma _0(-\phi )=0$
when
$\alpha \ll 1$
. Fields
$c_0$
,
$\Gamma _0$
,
$\langle c_0 \rangle$
and
$c_0(x, \, 0)$
are plotted in (f)
$A^{{2D}}_E$
(notice the classical stagnant-cap profile) and (g)
$D^{{1D}}_E$
(where streamwise diffusion has smoothed the stagnant cap), where
$x=x_0$
is plotted using (C62).
When bulk–surface exchange is weak (
$\nu = 0.01$
), 2-D effects are again prevalent in bulk-concentration profiles at small
$\alpha$
(figure 6). However, the impact of bulk boundary layers on the drag reduction is more modest than in the strong-exchange limit, in that they have less of an effect on the leading-order drag reduction. Solutions at large
$\gamma$
again have approximately linear interfacial surfactant profiles that immobilise the interface (figure 6
a,b). These profiles are almost decoupled from the bulk. They act as a weak sink/source combination driving adsorption/desorption. Solutions with weak Marangoni effects (small
$\gamma$
) show the formation of a more classical stagnant cap at small
$\alpha$
(figure 6
f), which is smoothed by streamwise diffusion as
$\alpha =\delta$
increases (figure 6
g).
Figure 6(c–e) illustrates the dependency of
${DR}_0$
on
$\alpha$
and
$\gamma$
, evaluated against existing limits
$M_E^{{1D}}$
and
$D_E^{{1D}}$
(see (B4), (B5)) for large
$\alpha$
and introducing new limits
$M_E^{{2D}}$
and
$A_E^{{2D}}$
at small
$\alpha$
. Importantly, weak coupling between the bulk and the interface suppresses the dependence of
$DR_0$
on
$\alpha$
in the limit of very small bulk diffusivity, making the drag-reduction contours horizontal in figure 6(c) for
$\alpha \to 0$
.
The following scaling argument applies to region
$A_E^{2D}$
(figures 2
b and 6
f), where Marangoni effects and bulk diffusion are weak. Here, the interfacial concentration is approximated by the classical distribution for nearly insoluble surfactant, which is close to zero along the interface except at the downstream end where it forms a stagnant cap of length
$L$
and slope
$1/\gamma$
. Thus,
$\Gamma _0$
is of magnitude
$L/\gamma$
, and its integral over the whole interface is
$O(L^2/\gamma )$
, while the bulk concentration remains close to unity. The surfactant flux-balance constraint (2.29c
), which matches net absorption to net desorption over the interface, therefore requires
$L^2/\gamma =O(1)$
, implying that
$\Gamma _0$
is
$O(\gamma ^{-1/2})$
and
$DR_0$
is
$O(\gamma ^{1/2})$
. This scaling argument, which mirrors the 1-D case, is supported by asymptotic calculations in Appendix C.3, which yield
\begin{align} c_0 \approx 1, \quad \Gamma _0 \approx \begin{cases} \ 0 \quad \text{for} \quad x\in [-\phi, \, x_0], \\ \ \dfrac {3\beta }{2\gamma }(x-x_0) \quad \text{for} \quad x\in [x_0,\,\phi ], \end{cases} \quad{DR}_0 \approx 1 - \left (\dfrac {2\gamma }{3\phi \beta }\right )^{1/2}, \end{align}
where
$\phi - x_0 = (8 \phi \gamma /(3\beta ))^{1/2}$
is the length of the stagnant cap. Bulk diffusion has a higher-order effect, influencing exchange via adsorption and desorption, and hence weak stretching and compression of the interface. In figure 6(c–f), dashed black asymptotes for
$c_0$
,
$\Gamma _0$
and
${DR}_0$
match those results from the 2-D long-wave model, almost up to
$\gamma = 3\beta \phi /2$
in figure 6(e) (see vertical dashed lines), which denotes the conditions where the stagnant cap reaches
$x=-\phi$
and the slip region of this nonlinear distribution vanishes. The disappearance of the stagnant-cap distribution with increasing
$\gamma$
is linked to a rapid decrease in
$DR_0$
, owing to the loss of the slip region. The same rapid transition is recovered if we approach region
$A^{2D}$
from region
$M_E^{2D}$
, using (B4
b), i.e. for decreasing
$\gamma$
.
Figure 7. The leading-order streamwise velocity field (
$u_0$
) and surfactant concentration fields (
$c_0$
and
$\Gamma _0$
) for
$\beta = 1$
,
$\nu = 0.01$
,
$\epsilon = 0.1$
and
$\phi = 0.5$
evaluated using (2.38)–(2.42) (lines) and COMSOL simulations (2.15)–(2.21) (symbols) when bulk–surface exchange is weak. (a,c) Plots of
$c_0$
and
$u_0$
, respectively, for
$\alpha = 0.1$
and
$\gamma = 1$
, when the flow is in the Marangoni-dominated region with weak cross-channel diffusion (
$M^{{2D}}_E$
). (b,d) Plots of
$c_0$
and
$u_0$
, respectively, for
$\alpha = 0.1$
and
$\gamma = 0.1$
, when the flow is in the advection-dominated region with weak cross-channel diffusion (
$A^{{2D}}_E$
) and
$x = x_0$
is plotted using (C62).
A comparison between the solutions evaluated using the 2-D long-wave model, (2.38)–(2.42), and COMSOL simulations is given in figure 7, focusing on the examples discussed in figure 6(a,f). Figure 7(a,c) shows
$c_0$
,
$\Gamma _0$
and
$u_0$
in region
$M_E^{{2D}}$
. The slip velocity is approximately zero for both symbols and lines, as the interfacial surfactant gradient effectively immobilises the interface. The 2-D long-wave model predicts
${DR}_0 = 0.027$
, while the numerical simulations yield
${DR}_{NS} = 0.056$
(due to the 4 % difference in
$\Delta \Gamma _0$
shown in figure 7
a). Figure 7(b,d) depicts stagnant-cap solutions in region
$A_E^{{2D}}$
. The upstream end of the liquid–gas interface (
$x\in [-\phi, \, x_0]$
) is shear-free, resulting in a slip velocity
$u_0 = 3/4$
(substituting
$\Gamma _{0x} = 0$
into (2.32a
), (2.36a
),
$u_0 = \tilde {U} p_{0x}$
, where
$\tilde {U}(0) = - 2$
and
$p_{0x} = - 3 / 8$
). Conversely, the downstream end of the liquid–gas interface (
$x\in [x_0,\, \phi ]$
) does not exhibit slip, leading to a near-zero streamwise velocity at the downstream stagnation point. The transitional and downstream regions induce a strong wall-normal flow that drives bulk surfactant into the core of the channel. We find
${DR}_0 = 0.650$
and
${DR}_{NS} = 0.662$
, demonstrating agreement between the predictions of the 2-D long-wave model and COMSOL simulations, validating the model’s accuracy under a wide range of conditions.
3.3. Comparison with experiments
To facilitate comparison with experiments, in table 2, we convert
${DR}_0$
to slip lengths (
$\lambda _e$
) using the non-dimensionalisation in Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023). Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023) proposed that for microchannel applications,
$\lambda _e \sim (\phi /\epsilon )^2 / L_m^2$
(using our notation), where
$\phi /\epsilon = \phi \hat {P}/\hat {H}$
is half the length of the liquid–gas interface divided by the channel height,
$L_m^2 = (\hat {n}\hat {R}\hat {T} \hat {K}_a^2 \hat {c}_0)/(\hat {D}\hat {\mu }\hat {K}_d^2)$
is the squared mobilisation length and
$\lambda _e \sim (\phi / \epsilon ) \lambda _0$
in (2.45). The slip-length formulae in Landel et al. (Reference Landel, Peaudecerf, Temprano-Coleto, Gibou, Goldstein and Luzzatto-Fegiz2020) and Sundin & Bagheri (Reference Sundin and Bagheri2022) have the same dependencies in this limit. We find that, in all regimes, the slip length shares the same dependence on
$1/L_m^{2}$
, with quadratic dependence on
$\phi /\epsilon$
only when
$1 \ll \alpha \ll 1/\epsilon ^{2}$
and
$\gamma \gg 1$
(region
$M^{{1D}}$
). However, when
$\alpha \ll 1$
and
$\gamma \gg \alpha ^{1/6}$
,
$\lambda _e$
has a stronger (8/3 power) dependence on
$\phi /\epsilon$
(region
$M^{{2D}}$
), and when
$\alpha \gg 1/\epsilon ^2$
and
$\gamma \gg \epsilon ^2 \alpha$
, we find a weaker (linear) dependence (region
$M^{{1D}}$
). The remaining asymptotic solutions for
${DR}_0$
are summarised in figure 2. When
$1-{DR}_0 \ll 1$
, the slip length simplifies to
$\lambda _e = (\Delta p_I - \Delta p_U) / (\epsilon \Delta p_I \Delta p_U)$
(
$A^{{1D}}$
and
$A^{{2D}}$
). Hence, outside of the specific microchannel applications considered in Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023), whose predictions agree with those in region
$M^{{1D}}$
(as discussed below), we have found new parameters in regions
$M^{{2D}}$
that influence surfactant impairment of SHS drag reduction. As diffusion weakens, the deviation between the drag-reduction predictions in
$M^{{1D}}$
and
$M^{{2D}}$
will increase.
Table 2. Summary of the leading-order drag reduction
$({DR}_0)$
scaling in regions
$M^{{2D}}$
,
$M^{{1D}}$
and
$M^{{1D}}_E$
and the corresponding leading-order slip length (
$\lambda _e$
) scaling using the non-dimensionalisation in Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023), in terms of the mobilisation length
$L_m$
. Hats denote dimensional quantities. The drag reduction is converted to the slip length using (2.45) for
${DR}_0 \ll 1$
and
$\lambda _e \sim (\phi /\epsilon ) \lambda _0$
, where
$\hat {L}_m^2 = (\hat {n}\hat {R}\hat {T} \hat {H}^2 \hat {K}_a^2 \hat {c}_0)/(\hat {D}\hat {\mu }\hat {K}_d^2)$
,
$\phi /\epsilon = \phi \hat {P}/\hat {H}$
,
$L_m = \hat {L}_m/\hat {H}$
and
$b = \hat {Q}/\hat {D}$
.
For laboratory experimental studies documented in the literature, we can use our 2-D long-wave model to estimate surfactant effects via the drag reduction, employing parameters typical of microchannel applications in laminar flows. Estimates for these parameters, based on the surfactant SDS, using models and experimental data from Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023), are given in table 3. The remaining parameters
$\hat {P}$
,
$\hat {Q}$
,
$\hat {H}$
and
$\hat {K}$
vary across experiments. Peaudecerf et al. (Reference Peaudecerf, Landel, Goldstein and Luzzatto-Fegiz2017) use
$\alpha \approx 1$
and
$4.7\times 10^{1} \lessapprox \gamma \lessapprox 2.9\times 10^4$
(at the
$M^{{1D}}/M^{{2D}}$
boundary; figure 2) and Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023) use
$\alpha \approx 25$
and
$10^{-1} \lessapprox \gamma \lessapprox 3.7\times 10^{2}$
(at the
$A^{{1D}}/M^{{1D}}$
boundary; figure 2). Qualitatively, all of these studies report surfactant effects, which is consistent with their location in the non-dimensional parameter space. All of these studies operate within the strong-exchange regime for the values of
$\hat {K}_a$
and
$\hat {K}_d$
given in table 3. However, as
$\hat {K}_a$
and
$\hat {K}_d$
are merely estimates based on fitting to microchannel experiments, weak exchange could be achieved, for instance, when
$\hat{K}_a = 3.4 \times 10^{-8}\,{\textrm{m}\, \textrm{s}}^{-1}$
and
$\hat{K}_d = 7.5\times 10^{-5}\,{\textrm{s}}^{-1}$
. Quantitatively, we compare predicted values of slip velocity using our theoretical model with experimental results in table 4. Using the experimental parameters from Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023), where
$\alpha \approx 25$
and
$\epsilon \approx 0.002$
, we confirm
$1\ll \alpha \ll 1/\epsilon ^2$
and they are within the regime
$M^{{1D}}$
. From (B1
b), we calculate
${DR}_0$
, convert to
$\lambda _e$
using table 2 and then to
$u_{Ic}/u_{Ic}^{\textrm {clean}}$
following Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023). Agreement is observed across
$\phi /\epsilon = 253$
, 418 and 587, which is remarkable as our theory does not include any fitting parameters. The deviation from Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023) at
$\phi /\epsilon = 762$
could be attributed to experimental error, which is consistent with the expected behaviour of increasing slip with grating length.
Table 3. A summary of the parameters in the dimensional problem, (2.1)–(2.12), with their values based on the surfactant SDS, models and experimental data from Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023).
Table 4. A comparison between the slip velocity normalised by the clean (surfactant-free) value,
$u_{Ic}/u_{Ic}^{\textrm {clean}}$
, predicted using the current model (2.38)–(2.42), and experimental data from Peaudecerf et al. (Reference Peaudecerf, Landel, Goldstein and Luzzatto-Fegiz2017) and Temprano-Coleto et al. (Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2023).The parameter
$\phi /\epsilon$
is half the length of the liquid–gas interface divided by the channel height.
4. Discussion
In this study, we develop a long-wave theory to analyse the behaviour of a 2-D laminar pressure-driven channel flow contaminated with soluble surfactant. The channel is bounded by a SHS with streamwise-periodic grooves and a solid wall. We linearise the equation of state and adsorption–desorption kinetics, solving Stokes and advection–diffusion equations in the long-wave limit. Our investigation focuses primarily on regimes where concentration gradients in the streamwise and cross-channel directions are comparable. By numerically solving the 2-D long-wave model, we delineate asymptotic regions within the parameter space. Under conditions of strong cross-channel diffusion, our findings align with the 1-D results in Tomlinson et al. (Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a ). Conversely, when cross-channel diffusion is weak, we reveal new regions of the parameter space, where the drag reduction exhibits a non-trivial dependence on the thickness of the bulk-concentration boundary layer and surfactant strength (discussed below). We derive asymptotic solutions for the boundary-layer problem, validating them against both the numerical solution of the 2-D long-wave model (2.38)–(2.42) and the full Stokes and advection–diffusion equations (2.15)–(2.21) in COMSOL. These complementary methodologies, ranging from asymptotic solutions to long-wave numerical simulations, give physical insight for applications where numerical simulations can be computationally prohibitive.
We have explored how the interfacial concentration (
$\Gamma _0$
), bulk concentration (
$c_0$
) and leading-order drag reduction (
${DR}_0$
) are influenced by bulk diffusion (
$\alpha$
), surface advection (
$\beta$
), surfactant strength (
$\gamma$
), interfacial diffusion (
$\delta$
) and bulk–surface exchange (
$\nu$
) (figures 3 and 6); the dimensionless parameters
$\alpha$
,
$\beta$
,
$\gamma$
,
$\delta$
and
$\nu$
are given in table 1. In cases where cross-channel diffusion dominates, we recover the Marangoni-dominated (
$M^{{1D}}$
and
$M^{{1D}}_E$
), advection-dominated (
$A^{{1D}}$
and
$A^{{1D}}_E$
) and diffusion-dominated (
$D^{{1D}}$
and
$D^{{1D}}_E$
) regions identified in Tomlinson et al. (Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a
) (see figure 2). However, under weak cross-channel diffusion, we reveal new subregions dominated by Marangoni effects (
$M^{{2D}}$
and
$M^{{2D}}_E$
) and advection (
$A^{{2D}}$
and
$A^{{2D}}_E$
). In these regions, a concentration boundary layer forms in the bulk, influencing the drag reduction when bulk–surface exchange is strong. When Marangoni effects dominate, the interface is immobilised and the surfactant distribution along it is approximately linear. The asymptotic solutions for
${DR}_0$
are summarised in figure 2.
When both bulk diffusion and Marangoni effects are weak, part of the interface is shear-free and the surfactant distribution forms either what we call a quasi-stagnant cap (when bulk–interface exchange is strong) or a classical stagnant cap (when exchange is weak). The classical stagnant cap for an insoluble surfactant-contaminated SHS has been described in Baier & Hardt (Reference Baier and Hardt2021) and Mayer & Crowdy (Reference Mayer and Crowdy2022). The quasi-stagnant cap has received much less attention and has a particularly intricate structure, illustrated in figure 4(b). Surface stretching and compression accommodate weak adsorption from, and desorption to, the bulk across thin concentration boundary layers, ‘remobilising’ the interface (analogous to a mechanism described by Crowdy et al. (Reference Crowdy, Curran and Papageorgiou2023) at low
$\textit {Pe}$
– and so without the involvement of bulk boundary layers – for a linear extensional flow). The upstream ‘slip’ region has a bulk-concentration boundary layer of thickness
$O(\alpha ^{1/2})$
, which is to be expected above an almost fully mobile interface. Weak stretching draws surfactant from the bulk to the interface via diffusive adsorption. The slip region transitions abruptly to a quasi-stagnant region, which grows to a thickness
$O(\gamma ^{1/5}\alpha ^{3/10})$
. This differs from a classical stagnant cap (of a fully insoluble surfactant) in having weak surface compression at a rate that balances diffusive desorption. The transition from the mobile to the quasi-immobile interface takes place across two nested regions at the tip of the quasi-stagnant cap: a short deceleration region, where the surface velocity falls abruptly, displacing the bulk boundary layer upwards towards the core; and a slightly longer transition region, across which shear in the boundary layer balances weakening surface advection (figure 4
$b$
). Shear then dominates in the bulk-concentration boundary layer along the quasi-stagnant region.
It is likely that the physical balances arising in these regions may emerge in other flow configurations involving soluble surfactant transport near confined interfaces at high Péclet numbers, such as the cap forming at the rear of a rising drop or bubble. For example, computations by Oguz & Sadhal (Reference Oguz and Sadhal1988) and Tasoglu et al. (Reference Tasoglu, Demirci and Muradoglu2008) revealed quasi-stagnant caps at large but finite
$\textit {Pe}$
. Interfacial flux balances between diffusion-limited adsorption (slip) and desorption (stagnant region) were given by Harper (Reference Harper2004) and Palaparthi et al. (Reference Palaparthi, Papageorgiou and Maldarelli2006), treating the size of the cap as a parameter. Here, we determine the size of the cap by solving the bulk-concentration boundary layer, which leads (for example) to the interfacial surfactant concentration being of size
$\beta /(\gamma ^{2/5}\alpha ^{1/10})$
in the
$A^{{2D}}$
regime.
In summary, our study compares asymptotic and numerical solutions for a 2-D laminar pressure-driven channel flow, confined by a streamwise-periodic SHS and a solid wall, and contaminated with soluble surfactant. While numerical solutions demand significant computational resources, our asymptotic solutions provide a cost-effective alternative, particularly in regimes at high Péclet numbers where accurate numerical solutions of very thin boundary layers are computationally very demanding. These asymptotic solutions offer accurate predictions devoid of empirical fitting coefficients and provide physical insights into the mechanisms governing drag reduction across a large part of the parameter space, including in particular high-Péclet-number flow regimes.
Funding
We acknowledge support from CBET–EPSRC (EPSRC Ref. EP/T030739/1, NSF no. 2054894), as well as partial support from ARO MURI W911NF-17-1-0306. F. T-C. acknowledges support from a Distinguished Postdoctoral Fellowship from the Andlinger Center for Energy and the Environment. For the purpose of open access, the authors have applied a Creative Commons Attribution (CCBY) licence to any Author Accepted Manuscript version arising.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical solutions for weak cross-channel diffusion
The numerical solution to (2.38)–(2.42) can be divided into two subroutines, employing Chebyshev collocation techniques as outlined by Trefethen (Reference Trefethen2000). First, given
$\Gamma _0^{\textit{old}}$
, we evaluate
$c_0^{\textit{new}}$
in
$\mathcal{D}_1$
and
$\mathcal{D}_2$
. The subdomains of the periodic cell are mapped to the 2-D canonical Chebyshev collocation domain,
$\mathcal{D}_n = \{\xi \in [-1, \, 1]\} \times \{\eta \in [-1, \, 1]\}$
, using the transformations
\begin{align} (\xi _1, \, \xi _2, \, \eta ) = \left ( (x + \phi )/\phi - 1, \, (x-\phi )/(1-\phi ) - 1, \, y - 1\right ). \end{align}
We then discretise using
$N = (N_\xi + 1)\times (N_\eta + 1)$
nodes in each subdomain, at points
\begin{align} (\xi _{1,\,i}, \, \xi _{2,\,i}, \, \eta _j) = (\cos (i \pi / N_\xi ), \, \cos (i \pi / N_\xi ), \, \cos (j \pi / N_\eta )), \end{align}
where
$i = 0, \, 1, \, \ldots, \, N_\xi$
and
$j = 0, \, 1, \, \ldots, \, N_\eta$
. Continuity boundary conditions (2.28) are enforced on
$c_0^{\textit{new}}$
at
$i=0$
and
$N_\xi$
, supplemented with continuity of bulk diffusive surfactant flux between subdomains:
\begin{align} c_{0, \, 1}^{\textit{new}} = c_{0, \, 2}^{\textit{new}}, \quad \frac {\partial c_{0, \, 1}^{\textit{new}}}{\partial \xi _1} = \frac {\partial c_{0, \, 2}^{\textit{new}}}{\partial \xi _2}. \end{align}
Boundary conditions (2.39)–(2.41) are enforced at
$j = 0$
and
$N_\eta$
, ensuring continuity of surfactant flux at the liquid–gas interface, no flux of surfactant at the solid wall in
$\mathcal{D}_1$
and no flux of surfactant at the solid walls in
$\mathcal{D}_2$
. Discretising the bulk surfactant equation (2.38) in conservative form in domains
$\mathcal{D}_1$
and
$\mathcal{D}_2$
and incorporating the discretised matching conditions (represented by
$\boldsymbol {{M}}_{12}$
and
$\boldsymbol {{M}}_{21}$
) gives
\begin{align} \begin{pmatrix} \boldsymbol {{A}}_1 & \boldsymbol {{M}}_{12} \\ \boldsymbol {{M}}_{21} & \boldsymbol {{A}}_{2} \end{pmatrix} \begin{pmatrix} \boldsymbol {c}_{0, \, 1}^{\textit{new}} \\[5pt] \boldsymbol {c}_{0, \, 2}^{\textit{new}} \end{pmatrix} = \begin{pmatrix} \boldsymbol {f}_{1} \\ \boldsymbol {0} \end{pmatrix}, \end{align}
where
$\boldsymbol {f}_{1}$
is the forcing due to the interfacial surfactant concentration
$\Gamma _0^{\textit{old}}$
in (2.39a
).
Second, given a
$c_0^{\textit{new}}$
, we evaluate
$\Gamma _0^{\textit{new}}$
for
$x\in [-\phi, \, \phi ]$
and
$y=0$
. The liquid–gas interface is mapped to the 1-D canonical Chebyshev collocation domain,
$\mathcal{D}_i = \{\xi \in [-1, \, 1]\}$
, using
$\xi _1$
in (A1), discretised using
$N_\xi + 1$
nodes at points
$\xi _{1,\,i}$
in (A2). Discretising the interfacial surfactant equation (2.39
b) for
$x\in [-\phi, \, \phi ]$
and
$y=0$
and incorporating the discretised no-flux conditions (2.40) results in
\begin{align} \boldsymbol {{B}}(\boldsymbol {\Gamma }_{0}^{\textit{old}}) \boldsymbol {\Gamma }_{0}^{\textit{new}} = \boldsymbol {g}, \end{align}
where
$\boldsymbol {g}$
is the forcing due to the bulk surfactant concentration
$c_0^{\textit{new}}$
in (2.39b
).
Therefore, to evaluate
${DR}_0$
for a given
$\alpha$
,
$\beta$
,
$\gamma$
,
$\delta$
,
$\epsilon$
,
$\nu$
and
$\phi$
, we choose an initial guess based on the 1-D long-wave model discussed in § 2.4. The initial guess can be substituted into the forcing
$\boldsymbol {f}_1$
in the bulk-concentration subroutine and the linear system inverted to calculate
$\boldsymbol {c}_0^{\textit{new}}$
. The bulk concentration can then be substituted into the forcing
$\boldsymbol {g}$
in the interfacial concentration subroutine and the linear system inverted to calculate
$\boldsymbol {\Gamma }_{0}^{\textit{new}}$
. We then set
$\boldsymbol {\Gamma }_{0}^{\textit{old}}=\boldsymbol {\Gamma }_{0}^{\textit{new}}$
and repeat the above steps. This procedure of updating
$\boldsymbol {c}_0$
and
$\boldsymbol {\Gamma }_{0}$
is then continued until convergence is achieved for both the bulk and interfacial concentration fields.
Appendix B. Key results of the 1-D long-wave model
Here we outline key limits of the 1-D long-wave model, (2.47) and (2.48), exploiting analogous findings as in Tomlinson et al. (Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a
). For strong exchange, with
$\alpha \gg 1$
and
$\gamma \gg \max (1, \, \epsilon ^2 \alpha )$
, a region dominated by Marangoni effects (
$M^{{1D}}$
; figure 2
a) is characterised by the immobilisation of the liquid–gas interface. The surfactant distribution and drag reduction are given by
\begin{align} \bar {c}_0 \approx \bar {\Gamma }_0 \approx 1 + \frac {3\beta }{2\gamma }\left (x - \frac {\phi (E+1)}{(E-1)}\right ), \quad {DR}_0 \approx \frac {2}{\gamma }\left (\epsilon ^2(2\alpha +\delta ) + \frac {\phi (E+1)}{(E-1)}\right ), \end{align}
where
$E\equiv \exp ((1-\phi )/(\epsilon ^2\alpha ))$
. Thus,
${DR}_0 \approx 2\phi /\gamma$
for
$1\ll \alpha \ll 1/\epsilon ^2$
and
${DR}_0\approx 2\epsilon ^2(2\alpha + \delta )/\gamma$
for
$\alpha \gg 1/\epsilon ^2$
. Region
$M^{{1D}}$
transitions into a diffusion-dominated region (
$D^{{1D}}$
; figure 2
a) when
$\alpha \gg 1/\epsilon ^2$
and
$\gamma \sim \epsilon ^2\alpha$
, and into an advection-dominated region (
$A^{{1D}}$
; figure 2
a) when
$1 \ll \alpha \ll 1/\epsilon ^2$
and
$\gamma \sim 1$
. For
$\gamma \ll \epsilon ^2 \alpha$
and
$\alpha \gg 1/\epsilon ^2$
(region
$D^{{1D}}$
),
\begin{align} \bar {c}_0 \approx \bar {\Gamma }_0 \approx \frac {4\alpha +2\delta (1 -\phi )}{\alpha (4+3\beta \phi ) +2\delta (1-\phi )}, \ \ {DR}_0 \approx 1 - \frac {\gamma (1-\phi )}{\epsilon ^2(\alpha (4+3\beta \phi ) +2\delta (1-\phi ))}, \end{align}
and for
$\gamma \ll 1$
and
$1 \ll \alpha \ll 1/\epsilon ^2$
(region
$A^{{1D}}$
),
\begin{align} \bar {c}_0 \approx \bar {\Gamma }_0 \approx \frac {4}{3\beta + 4} + \frac {3\beta }{3\beta +4}\exp \left (\frac {(4+3\beta )(x-\phi )}{4\epsilon ^2(2\alpha + \delta )}\right ), \quad {DR}_0 \approx 1 - \frac {\gamma }{\phi (3\beta +4)}. \end{align}
The onset of 2-D effects in the 1-D long-wave model was estimated in Tomlinson et al. (Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a
) to arise via the emergence of shear dispersion, where
$\alpha \sim 1$
and
$\gamma \ll 1$
or
$\alpha \sim 1/\gamma$
and
$\gamma \gg 1$
. We refine such estimates by solving the 2-D long-wave model, which reveals regions
$M^{{2D}}$
and
$A^{{2D}}$
shown in figure 2(a).
A similar picture emerges when bulk–surface exchange is weak. A region dominated by Marangoni effects (
$M^{{1D}}_{E}$
; figure 2
b) arises for
$\delta \gg 1$
and
$\gamma \gg \max (1,\,\epsilon ^2 \alpha )$
. Here,
$\bar {c}_0$
decouples from
$\bar {\Gamma }_0$
and is approximately equal to its background value (i.e.
$\bar {c}_0 \approx 1$
). The surface concentration
$\bar {\Gamma }_0$
varies around this value to fulfil the net adsorption–desorption condition, (2.29), decreasing the drag reduction compared with the strong-exchange problem:
\begin{align} \bar {c}_0 \approx 1, \quad \bar {\Gamma }_0 \approx 1 + \frac {3\beta x}{2\gamma }+O(\epsilon ^2), \quad {DR}_0 \approx \frac {2\epsilon ^2 \delta }{\gamma }. \end{align}
The region
$M^{{1D}}_{E}$
transitions to a diffusion-dominated region (
$D^{{1D}}_{E}$
; figure 2
b) when
$\gamma \sim \epsilon ^2 \delta$
and
$\delta \gg 1/\epsilon ^2$
and an advection-dominated region (
$A^{{1D}}_{E}$
; figure 2
b) when
$\gamma \sim 1$
and
$1\ll \delta \ll 1/\epsilon ^2$
. For
$\gamma \ll \epsilon ^2 \delta$
and
$\delta \gg 1/\epsilon ^2$
(region
$D_E^{{1D}}$
),
\begin{align} \bar {c}_0 \approx 1, \quad \bar {\Gamma }_0 \approx 1 + \frac {3\beta x}{4\epsilon ^2 \delta }, \quad {DR}_0 \approx 1 - \frac {\gamma }{2\epsilon ^2\delta }, \end{align}
and for
$\gamma \ll 1$
and
$1\ll \delta \ll 1/\epsilon ^2$
(region
$A_E^{{1D}}$
),
\begin{align} \bar {c}_0 \approx 1, \quad \bar {\Gamma }_0 \approx \begin{cases} \ 0 \quad \text{for} \quad x\in [-\phi, \, x_0], \\ \ \dfrac {3\beta }{2\gamma }(x-x_0) \quad \text{for} \quad x\in [x_0,\,\phi ], \end{cases}\quad {DR}_0 \approx 1 - \left (\frac {2\gamma }{3\phi \beta }\right )^{1/2}, \end{align}
where
$x_0 = \phi - (8 \phi \gamma /(3\beta ))^{1/2}$
gives the length of the stagnant cap.
Appendix C. Asymptotic solutions for weak cross-channel diffusion
We derive asymptotic solutions to the 2-D long-wave model, (2.38)–(2.42), when bulk diffusive effects are weak and confined to a boundary layer along
$y=0$
. In this limit, streamwise diffusion terms do not appear at leading order, which means that we do not resolve inner regions around
$x = -\phi$
,
$\phi$
and
$2-\phi$
. First, we consider the strong-exchange regime, where Marangoni effects (Appendix C.1; figure 4
$a$
) and advection (Appendix C.2; figure 4
$b$
) dominate. Second, we consider the weak-exchange regime (Appendix C.3). Throughout, we assume
$\beta =O(1)$
and
$\phi = O(1)$
.
Assuming that
$c_0 \approx 1$
in the core of the channel, we anticipate that the transport equations (2.38) for the bulk surfactant involve a diffusive boundary layer near
$y=0$
. The diffusive boundary layer is made up of components over the interface and ridge, which are governed by
\begin{align} \alpha c_{0yy} - \left (\frac {3}{4}+\frac {\gamma }{\beta }\left (-\frac {1}{2} + y\right )\Gamma _{0x}\right )c_{0x} - \frac {\gamma }{\beta }\left (\frac {y}{2} - \frac {y^2}{2}\right )\Gamma _{0xx} c_{0y} &=0 \quad \text{in} \quad \mathcal{D}_1, \end{align}
\begin{align} \alpha c_{0yy} - {3y} c_{0x}/2 &=0 \quad \text{in} \quad \mathcal{D}_2, \\[6pt] \nonumber \end{align}
where we have assumed
$y\ll 1$
and have neglected second-order (third-order) terms in
$O(y^2)$
(
$O(y^3)$
) in the streamwise (wall-normal) velocity. For
$x\in [-\phi, \, \phi ]$
and
$y=0$
, the bulk–interface flux condition and the interfacial transport equation for the surfactant in (2.39) give
\begin{align} \alpha c_{0y} = \nu (c_0 - \Gamma _0) = \beta \bigg [\left (\frac {3}{4}-\frac {\gamma }{2\beta }\Gamma _{0x}\right ) \Gamma _0 \bigg ]_x. \end{align}
At
$x=\pm \phi$
and
$y=0$
, the no-surfactant-flux conditions at the stagnation points (2.40) become
\begin{align} \left (\frac {3}{4}-\frac {\gamma }{2\beta }\Gamma _{0x}\right ) \Gamma _0 =0. \end{align}
For
$x\in [\phi, \, 2 - \phi ]$
and
$y=0$
, the no-flux condition at the solid wall (2.41) gives
\begin{align} c_{0y} =0. \end{align}
For
$x\in [-\phi, \, 2 - \phi ]$
and
$y\rightarrow 1$
(shorthand for
$y$
values in the core of the channel, far outside the boundary layer along
$y=0$
), the core condition is
\begin{align} c_0 \rightarrow 1. \end{align}
The boundary-layer equations (C1)–(C5) can be integrated across the boundary layer to derive excess-surfactant-flux conditions. The total flux constraint (2.42b , c ) is satisfied at leading order in the core of the channel, implying that
\begin{align} \!\!\!\beta \left (\frac {3}{4}-\frac {\gamma }{2\beta }\Gamma _{0x}\right ) \Gamma _0 + \int _{0}^1 \!\!\left (\frac {3}{4}+\frac {\gamma }{\beta }\left (-\frac {1}{2} + y\right )\Gamma _{0x}\right )\left (c_0-1\right ) \, \text{d}y &= 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}
\begin{align} \int _{-\phi }^\phi c_{0y} \, \text{d} x = \int _{-\phi }^\phi (c_0-\Gamma _0) \, \text{d} x &= 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}
\begin{align} \int _{0}^1 \frac {3}{2}y\left (c_0-1\right ) \, \text{d}y &= 0 \quad \text{in} \quad \mathcal{D}_2, \\[6pt] \nonumber \end{align}
provided
$c_0\rightarrow 1$
sufficiently quickly where the boundary layer meets the core flow. Equation (C6b
) is derived by integrating (C2) and using (C3). The condition in (C6a
) indicates that the extra surfactant flux carried by the flow at the interface balances the flux deficit in the boundary layer (where
$c_0\leqslant 1$
), which happens over the upstream part of the interface. Conversely, over the downstream part of the interface, the reduced interfacial surfactant flux is compensated by a surplus of bulk surfactant in the boundary layer (where
$c_0\geqslant 1$
). This can be seen in the boundary-layer profiles at small
$\alpha$
values, as depicted in figures 3(a,f) and 6(a,f). The condition in (C6b
) specifies that the total amounts of surfactant adsorbed and desorbed along the interface must match, with continuity between the bulk diffusive flux (first integral) and the adsorption–desorption flux (second integral). Motivated by the simulations and scaling arguments provided in § 3.1, we propose the asymptotic structures for regions
$M^{{2D}}$
and
$A^{{2D}}$
outlined in figure 4 and derive asymptotic solutions of the surfactant distributions in the following subsections for the different regions in the parameter space where boundary layers exist.
C.1.
Strong Marangoni effect and strong bulk–surface exchange: region
$M^{{2D}}$
We first consider the boundary-layer equations, (C1)–(C6), in the limit where bulk diffusion is weak and bulk–surface exchange is strong, i.e.
$\alpha \ll 1$
and
$\nu$
large. Suppose the interfacial surfactant gradient leads to a significant adverse Marangoni effect, immobilising the interface. Surfactant adsorbs onto the interface along its upstream part and desorbs along the downstream part. Motivated by the scaling arguments given in § 3.1, we set
\begin{align} \Gamma _0 = M_1 + \frac {3 \beta x}{2 \gamma } - \alpha ^{2/3} G(x) + \cdots, \quad c_0 = C(x, \, Y), \end{align}
where
$y = \alpha ^{1/3} Y$
balances advection and diffusion in the bulk-surfactant equation (C1) and
$M_1$
is a constant to be determined. Substituting (C7) into (C1)–(C6), (C1) becomes, at leading order,
\begin{align} C_{YY} - 3Y C_x/2 =0. \end{align}
For
$x\in [-\phi, \, \phi ]$
and
$Y=0$
, we enforce
$\Gamma _0 = c_0$
as bulk–surface exchange is strong, so that (C2) becomes
\begin{align} C_{Y} - \frac {\gamma }{2}\left [G_x\left (\frac {3\beta x}{2\gamma } + M_1 \right )\right ]_x= 0. \end{align}
At
$x=\pm \phi$
and
$Y=0$
, the no-flux conditions (C3) at the stagnation points are
\begin{align} G_x =0, \end{align}
and for
$x\in [\phi, \, 2 - \phi ]$
and
$Y=0$
, the no-flux condition (C4) at the solid wall is
\begin{align} C_{Y} =0. \end{align}
In the bulk of the periodic domain, for
$x\in [-\phi, \, 2 - \phi ]$
and
$Y\rightarrow \infty$
, the core condition (C5) becomes
\begin{align} C \rightarrow 1. \end{align}
The excess-surfactant-flux conditions (C6) reduce to
\begin{align} \int _{0}^\infty \frac {3 Y}{2 }\left (C-1\right ) \, \text{d}Y + \frac {G_x}{2}\left (\frac {3\beta x}{2} + \gamma M_1 \right ) &= 0 \quad \text{in} \quad \mathcal{D}_1, \end{align}
\begin{align} \int _{0}^\infty \frac {3Y}{2}\left (C-1\right ) \, \text{d}Y &= 0 \quad \text{in} \quad \mathcal{D}_2. \end{align}
We construct a similarity solution to (C8)–(C13), assuming that the length of the interface is long compared with the boundary-layer thickness, and that there is no interaction between adjacent plastrons. The bulk concentration is split into two contributions that grow from the leading edge of the plastron:
\begin{align} C = 1 - \left (1 - M_1 + \frac {3\beta \phi }{2\gamma } \right )F(\eta ) + \frac {3\beta }{2\gamma }(x+\phi )H(\eta ) \end{align}
for some functions
$F$
and
$H$
depending on the similarity variable
$\eta = Y/(x+\phi )^{1/3}$
. The function
$F$
captures the response of the bulk field to the sudden drop in concentration along
$Y=0$
at
$x=-\phi$
;
$H$
captures the response to the linear profile of
$\Gamma _0$
in
$x\gt -\phi$
. Substituting (C14) into (C8)–(C12) leads to the boundary-value problems
\begin{align} F_{\eta \eta } + \eta ^2 F_\eta / 2 = 0, \quad F(0)=1, \quad F(\infty ) =0, \\[-6pt] \nonumber \end{align}
\begin{align} H_{\eta \eta } + \eta ^2 H_\eta / 2 - 3 \eta H/2 = 0, \quad H(0)=1, \quad H(\infty ) =0, \\[12pt] \nonumber \end{align}
which have solutions
\begin{align} F = \frac {\Gamma \left(\dfrac {1}{3}, \, \eta ^3/6 \right)}{\Gamma \left(\dfrac {1}{3} \right)}, \quad H = \frac { (4 + \eta ^3) \Gamma \left(\dfrac {1}{3}, \, \eta ^3/6 \right) -6^{2/3} \eta \exp (-\eta ^3/6) }{4 \Gamma \left(\dfrac {1}{3} \right)}. \end{align}
Here
$\Gamma (z)$
is the gamma function and
$\Gamma (s,\,z)$
is the upper incomplete gamma function. We then use the surfactant-flux conditions (C13) to calculate the response of the interfacial concentration field to the bulk field (C14):
\begin{align} \frac {\gamma G_x}{2}\left (\frac {3\beta x}{2 \gamma } + M_1\right ) = \,& \frac {3}{2} (x+\phi )^{2/3}\left (1 - M_1 + \frac {3\beta \phi }{2\gamma }\right )\int _{0}^\infty \eta F \, \text{d}\eta \nonumber \\ & - \frac {9\beta }{4\gamma }(x+\phi )^{5/3} \int _{0}^\infty \eta H \, \text{d}\eta . \\[1pt] \nonumber \end{align}
The total flux of surfactant must be continuous between domains. The no-flux condition (C10) at
$x=\phi$
requires that
\begin{align} M_1 = 1 + \frac {3\beta \phi }{2\gamma }\left (\int _{0}^\infty \eta F \, \text{d}\eta - 2 \int _{0}^\infty \eta H \, \text{d}\eta \right )\bigg /\int _{0}^\infty \eta F \, \text{d}\eta = 1 - \frac {3\beta \phi }{10\gamma }. \end{align}
Integrating the total-flux condition (C18) across the plastron reveals that
\begin{align} \Delta G &= \int _{-\phi }^{\phi }\bigg \{\frac {27\beta \phi }{10\gamma } (x+\phi )^{2/3}\int _{0}^\infty \eta F \, \text{d}\eta \notag\\ & \quad - \frac {9\beta }{4\gamma }(x+\phi )^{5/3} \int _{0}^\infty \eta H \, \text{d}\eta \bigg \}\bigg /\left [\frac {\gamma }{2}\left (1 + \frac {3\beta }{2 \gamma }\left (x - \frac {\phi }{5}\right )\right )\right ]\,\text{d}x. \end{align}
Hence, in region
$M^{{2D}}$
, substituting (C7) into (2.44), the leading-order interfacial concentration and drag reduction are given by
\begin{align} \Gamma _0 = 1 + \frac {3\beta }{2\gamma }\left (x - \frac {\phi }{5}\right ) + \cdots, \quad {DR}_0 = \frac {\gamma \alpha ^{2/3}\Delta G}{3\phi \beta }. \end{align}
The leading-order interfacial surfactant distribution
$\Gamma _0$
in region
$M^{{2D}}$
is a linear profile with mean value close to 1, as found numerically (see figure 5
a). This approximation requires
$\Gamma _0\gt 0$
at
$x=-\phi$
, so that
$\gamma \gt 9\beta \phi /5$
. For
$\gamma \gg 1$
, (C21) reduces to
\begin{align} \Delta G= \frac {81\beta (2\phi )^{8/3}}{50\gamma ^2}\int _{0}^\infty \eta F \, \text{d}\eta - \frac {27\beta (2\phi )^{8/3}}{16\gamma ^2} \int _{0}^\infty \eta H \, \text{d}\eta \equiv \frac {\Delta J}{\gamma ^2}, \\[-12pt] \nonumber \end{align}
\begin{align} \Gamma _0 = 1 + \frac {3\beta }{2\gamma }\left (x - \frac {\phi }{5}\right ) + \cdots, \quad {DR}_0 = \frac {\alpha ^{2/3}\Delta J(\beta, \, \phi )}{3\phi \beta \gamma } \equiv \frac {m_1 \alpha ^{2/3}\phi ^{5/3}}{\gamma }, \\[6pt] \nonumber \end{align}
where
$m_1 = 81 (3/2)^{2/3}/(50 {\Gamma }({1}/{3})) \approx 0.79$
, as in (3.1). This expression for the leading-order drag reduction (C23b
) assumes strong Marangoni effects (
$\gamma \gg 1$
), weak diffusion (
$\alpha \ll 1$
) and strong bulk–surface exchange, and shows how drag reduction in this limit arises from differences between the two self-similar components
$F$
and
$H$
of the bulk-concentration field. As mentioned previously, this asymptotic prediction agrees with the numerical simulations in the
$M^{{2D}}$
region shown in figure 3(a,c–e).
C.2.
Weak Marangoni effect and strong bulk–surface exchange: region
$A^{{2D}}$
We now consider the boundary-layer equations (C1)–(C6) when
$\gamma \ll 1$
and
$\alpha \ll 1$
. Weak Marangoni effects imply that the interface cannot support any stress over most of its length. Consequently, there exists a slip flow along most of the plastron, characterised by a slip velocity in (C1a
) close to
$3/4$
. A short transition region, mediated by surface diffusion, exists between a no-slip and a slip flow at the upstream contact line, but this is disregarded here. Surfactant adsorbs onto the interface in the slip region (in
$-\phi \lt x\lt x_0$
, say) and accumulates (and desorbs) in a region at the downstream end of the plastron (in
$x_0\lt x\lt \phi$
). Across this ‘quasi-stagnant’ region, the surface velocity smoothly decreases to zero at the downstream contact line. Motivated by simulations and the scaling arguments given in § 3.1, we propose the asymptotic structure illustrated in figure 4(b), in which ‘deceleration’ and ‘transition’ regions are nested at the tip of the quasi-stagnant region over length scales indicated. The formal requirement that the deceleration, transition and quasi-stagnant layers are nested (
$\gamma \alpha ^{1/2}\ll \gamma ^{3/4}\alpha ^{1/8}\ll \gamma ^{3/5}/\alpha ^{1/10} \lesssim 1$
) is
$\gamma \ll \alpha ^{-3/2}$
. We now discuss these regions in turn.
C.2.1.
The slip region
$(-\phi \lt x\lt x_0)$
In the slip region, surfactant adsorbs onto the interface, causing the interface to carry a spatially growing component of the overall flux. Consequently, the bulk advective flux falls with
$x$
. Here, we set
\begin{align} \Gamma _0=\alpha ^{1/2}G(x), \quad c_0=C(x,Y), \end{align}
where
$y=\alpha ^{1/2}Y$
balances advection and diffusion in the bulk-transport equation (C1). When we substitute (C24) into (C1)–(C6), the bulk-transport equation (C1) in
$\mathcal{D}_1$
becomes
\begin{align} C_{YY} - 3 C_x/4 =0. \end{align}
At the interface
$\alpha ^{1/2} G(x) = C(x,0)$
in the strong-exchange limit and the flux condition (C2) reduces to
\begin{align} C_{Y} (x,0)= 3\beta G_x/4. \end{align}
The no-flux condition at the upstream stagnation point (C3) is
\begin{align} G(-\phi ) =0. \end{align}
For
$x\in [-\phi, \, x_0]$
and
$Y\rightarrow \infty$
, the core condition (C5) is
\begin{align} C \rightarrow 1. \end{align}
We construct a similarity solution to (C25)–(C28) by expanding the surfactant fields for
$\alpha \ll 1$
as
\begin{align} C(x,Y)=C_0(\eta )+(\alpha (x+\phi ))^{1/2}C_1(\eta )+\cdots, \quad G(x)=(x+\phi )^{1/2}C_1(0), \end{align}
with the similarity variable given by
$\eta =Y/(\phi +x)^{1/2}$
. Substituting (C29) into (C25)–(C28) yields the boundary-value problems
\begin{align} C_0^{\prime\prime}+3\eta C_0^{\prime}/8=0, \quad C_0(0)=0, \quad C_0(\infty ) = 1, \\[-24pt] \nonumber \end{align}
\begin{align} C_1^{\prime\prime}+3\eta C_1^{\prime}/8-3C_1/8=0, \quad 3\beta C_1(0) = 8C_0^{\prime}(0), \quad C_1(\infty )=0. \\[6pt] \nonumber \end{align}
These are solved by
\begin{align} C_0=\text{erf}(\sqrt {3}\eta /4), \quad C_1=4\text{ierfc}(\sqrt {3}\eta /4)/(\sqrt {3}\beta ), \end{align}
where
$\text{erf}(z)$
is the error function,
$\text{erfc}(z)$
is the complementary error function and
$\mathrm{ierfc}(z) \equiv {\textrm{e}}^{-z^2}/\sqrt {\pi }-z\,\text{erfc}(z)$
is the integrated complementary error function. We find that the interfacial concentration field in the slip region is given by
\begin{align} \Gamma _0=\frac {4}{\beta } \left (\frac {\alpha (x+\phi )}{3 \pi }\right )^{1/2}+\cdots, \end{align}
contributing to (3.2
$a$
). The interfacial surfactant distribution
$\Gamma _0$
in region
$A^{{2D}}$
is a nonlinear boundary-layer profile scaling as
$\alpha ^{1/2}x^{1/2}$
in the slip region of the interface (up to
$x_0$
), as found numerically (see figure 5
b). There is a large vertical concentration gradient across the boundary layer. Adsorption of surfactant onto the interface is accommodated by the horizontal gradient of the advective flux
$\beta u_0 \Gamma _0 \sim 3\beta \Gamma _0/4$
. Marangoni effects remain subdominant, so the slip flow is uniform to this order, but the concentration rises with
$x$
to accommodate the adsorbed material. The total flux adsorbed across the slip region, up to the location
$x=x_0$
(to be determined), is given by
\begin{align} \int _{-\phi }^{x_0} c_{0y}(x,0)\,\text{d}x= \left (\frac {9\alpha (x_0+\phi )}{3 \pi }\right )^{1/2}+\cdots, \end{align}
which is subdominant to the flux in the core. The slip region terminates abruptly, meeting a very short deceleration region near
$x=x_0$
nested at the start of the quasi-stagnant region (figure 4
b).
C.2.2.
Deceleration region
$( x-x_0 =O(\gamma \alpha ^{1/2}))$
In the deceleration region (figure 4
b),
$u_0(x,\,0)$
varies by
$O(1)$
over a short horizontal scale
$\Delta x$
and
$\Gamma _0\sim \alpha ^{1/2}$
(where
$\sim$
denotes ‘scales like’), based on continuity with the slip region. Thus,
$\Delta x\sim \gamma \alpha ^{1/2}$
for
$\gamma \Gamma _{0x}$
to vary by
$O(1)$
in (C2). Bulk diffusion balances advection by
$u_0$
over a vertical length scale
$y\sim (\alpha \Delta x)^{1/2}\sim \gamma ^{1/2}\alpha ^{3/4}$
using (C1). Diffusion therefore influences the concentration field near the interface, while advection dominates transport in the upper part of the boundary layer. At the downstream limit of this region,
$u_0$
falls towards zero, with
$\Gamma _{0x}\approx 3\beta /(2\gamma )$
. The deceleration region is too short for there to be appreciable desorption so flux conservation requires that
$u_0\sim \alpha ^{1/2}/\Gamma _0\sim \alpha ^{1/2} \gamma /(x-x_0)$
at the outlet of the region. We write the interfacial flux as
\begin{align} u_0(x,\,0) \Gamma _0 =\frac {3}{\beta }\left (\frac {\alpha (x_0+ \phi )}{3\pi }\right )^{1/2} \equiv \frac {3 \alpha ^{1/2}}{4} \lambda, \quad \text{where} \quad \lambda = \frac {4}{\beta }\left (\frac {x_0+\phi }{3\pi }\right )^{1/2}, \end{align}
as
$u_0 = 3/4 + \cdots$
and
$\Gamma _0$
is (C33) in the slip region. We then rescale, using
\begin{align} x=x_0+\gamma \alpha ^{1/2} X, \ \ \ y=\gamma ^{1/2}\alpha ^{3/4} Y, \ \ \ \Gamma _0=\alpha ^{1/2}\lambda G(X), \ \ \ c_0=\alpha ^{1/2}\lambda C(X,Y). \end{align}
Substituting (C36) into (C1)–(C6), the bulk-surfactant equation (C1) in
$\mathcal{D}_1$
is given by
\begin{align} C_{YY} - \left (\frac {3}{4}-\frac {\lambda }{2\beta } G_X\right ) C_X - \frac {\lambda Y}{2 \beta }G_{XX} C_Y =0, \end{align}
neglecting terms of
$O(\alpha ^{3/4})$
. We have
$G(X) = C(X,0)$
in the strong-exchange limit and the interfacial-surfactant equation (C2) simplifies to
\begin{align} \left [G\left (\frac {3}{4} - \frac {\lambda }{2\beta } G_X \right )\right ]_X = 0. \end{align}
Integrating (C38) and matching it to the slip region, where
$G \rightarrow 1$
for
$X\rightarrow -\infty$
, gives
\begin{align} \frac {3}{4}G - \frac {\lambda }{2\beta } G G_X = \frac {3}{4}, \end{align}
which can be solved directly to give
\begin{align} G+\log (G-1)=(3\beta /(2\lambda ))(X-X_0) \end{align}
for some
$X_0$
(as in Jensen & Halpern (Reference Jensen and Halpern1998), for example). Substituting
$G$
into (C35a
), we compare
$u_0(x,\,0)$
with the COMSOL simulations in figure 5(d), showing how
$u_0$
falls steeply from
$3/4$
. The core condition (C5) becomes
\begin{align} C \rightarrow \infty \quad \mathrm{as}\quad Y\rightarrow \infty . \end{align}
Equations (C37)–(C41) constitute a nonlinear boundary-layer problem that describes the abrupt ejection of bulk concentration towards the core of the channel in figure 3(f). We only consider the downstream limit of (C37)–(C41) here. The upstream condition on
$C$
is provided from the near-wall component of the concentration in the slip region, which grows in amplitude with
$Y$
in (C32). The downstream limit has self-similar form:
\begin{align} G\approx \frac {3\beta }{2\lambda }(X-X_0)+\cdots, \quad C\approx \frac {3\beta }{2\lambda }(X-X_0)\hat {C}(\eta )+\cdots, \end{align}
where the similarity variable is given by
$\eta ={Y}/(X-X_0)$
for some
$X_0$
and
$\hat {C}$
satisfies
$\hat {C}_{\eta \eta }=(\lambda /2\beta )\hat {C}$
, which has solutions
$\exp [\pm \sqrt {\lambda /2\beta }\eta ]$
. Concentration contours are constant along lines
$Y=\sqrt {2\beta /\lambda }(X-X_0)$
, representing ejection of bulk surfactant away from the interface. Expressed in terms of the original variables using (C36), the downstream limit (C42) of this region yields
\begin{align} \Gamma _0 \approx \frac {3\beta }{2\gamma }(x-x_0),\quad u_0(x,\,0) \approx \frac {\lambda \alpha ^{1/2} \gamma }{2\beta (x-x_0)} \quad \mathrm{for}\quad \gamma \alpha ^{1/2}\ll x-x_0\ll \gamma ^{3/4}\alpha ^{1/8}, \end{align}
which delivers an interfacial advective flux
$3\lambda \alpha ^{1/2}/4$
into a very short transition region near
$x=x_0$
. Surface diffusion is a subdominant effect in this region provided
$\epsilon ^2 \ll \gamma /\alpha ^{1/2}$
. The leading-order interfacial concentration profile in the deceleration region, (C43a
), contributes to the second part of (3.2a
) (for
$x_0\leqslant x$
) and (C43b
) explains the steep fall in surface velocity in figure 5(d) at the end of the slip region; the incident boundary layer of thickness
$\alpha ^{1/2}$
is thickened to
$\alpha ^{1/2}/u_0$
as
$u_0$
diminishes (see the arrow in figure 3
f). The surface velocity in this regime, computed using (C35) with
$\Gamma _0$
evaluated using (C39), is shown in figure 5(d).
C.2.3.
Transition region
$(x- x_0=O(\gamma ^{3/4}\alpha ^{1/8}))$
In the transition region, there is a balance in the bulk between vertical diffusion, advection by
$u_0$
(which is
$O(1)$
) and advection by shear (which is affected by viscous and Marangoni effects). Hence, over a horizontal length scale
$\Delta x$
,
$\alpha /y^2\sim u_0/\Delta x\sim y/\Delta x$
; thus
$y\sim u_0\sim (\alpha \Delta x)^{1/3}$
. Matching
$u_0$
to the deceleration region requires
$u_0\sim \alpha ^{1/2}\gamma /\Delta x$
, implying
$\Delta x\sim \gamma ^{3/4}\alpha ^{1/8}$
and
$y\sim u_0\sim \gamma ^{1/4}\alpha ^{3/8}$
. The interfacial concentration is linear at leading order (
$\Gamma _0\approx 3\beta (x-x_0)/(2\gamma ))$
with a correction of size
$u_0 \Delta x/\gamma \sim \alpha ^{1/2}$
. The bulk concentration is set by the dominant interfacial concentration and is
$O(\Delta x/\gamma )\sim \alpha ^{1/8}/\gamma ^{1/4}$
. The region is again too short for there to be appreciable desorption, so the interfacial flux is still
$O(\alpha ^{1/2})$
. To provide a bridge between the deceleration region and the quasi-stagnant region further downstream, we reintroduce shear
$3Y/2$
to the bulk transport equation, using the scaling
\begin{align} x=x_0+\alpha ^{1/8}\gamma ^{3/4}X, \quad y=\alpha ^{3/8}\gamma ^{1/4}Y, \nonumber \\ \Gamma _0=\frac {3\beta }{2\gamma }(x-x_0)+\alpha ^{1/2}G(X), \quad c_0=\frac {\alpha ^{1/8}}{\gamma ^{1/4}}C(X,Y). \end{align}
Substituting (C44) into (C1)–(C6), the bulk transport equation (C1) in
$\mathcal{D}_1$
simplifies to
\begin{align} C_{YY} + \left (\frac {G_X}{2\beta } -\frac {3Y}{2} \right )C_X- \frac {G_{XX}}{2\beta } Y C_Y = 0. \end{align}
Again
$C(X,0)=G(X)$
and the interfacial transport equation (C2) becomes
\begin{align} - \left [\frac {3\beta X}{4}G_X\right ]_X = 0. \end{align}
Using (C46), there is no exchange of surfactant with the bulk to leading order. Integrating (C46) in
$X$
and matching to the deceleration region,
$XG_X \rightarrow \lambda$
as
$X\rightarrow -\infty$
. This can be integrated in
$X$
to give
$G=-\lambda \log X+A_1$
for some constant
$A_1$
, which we substitute into (C45), to obtain
\begin{align} C_{YY}-\frac {\lambda }{2\beta } \left (\frac {1}{X}C_X -\frac {YC_Y}{X^2}\right ) -\frac {3}{2}YC_X = 0. \end{align}
This linear boundary-layer problem bridges the scalings
$Y\sim X$
for
$X \ll 1$
and
$Y\sim X^{1/3}$
for
$X\gg 1$
, but we do not solve the full problem here. Expressed in terms of original variables using (C44), the downstream limit of the transition region demonstrates that
\begin{align} \Gamma _0\approx \frac {3\beta (x-x_0)}{2\gamma }, \quad u_0(x,\,0)\approx \frac {\lambda \alpha ^{1/2}\gamma }{2\beta (x-x_0)}\quad \mathrm{for} \quad \gamma ^{3/4}\alpha ^{1/8}\ll x-x_0\ll 1. \end{align}
Note that the interfacial concentration profile in this subregion, (C48
a), remains linear and exactly the same as in the deceleration region, see (C48
a), by continuity. Over this part of the interface, the interfacial flux is dominated by advection:
$u_0(x,0)\Gamma _0\approx 3\lambda \alpha ^{1/2}/4$
. This flux is delivered to the rest of the interface, in the last region that we call the quasi-stagnant region for
$x_0 \lt x \lt \phi$
. The drag reduction can only be determined once the interfacial concentration profile across the whole interface is determined.
C.2.4.
The quasi-stagnant region (
$x_0\lt x\lt \phi$
)
In the quasi-stagnant region, compression of the almost immobile interface promotes desorption, returning surfactant to the bulk. Therefore, the interfacial flux of surfactant decreases and the bulk flux increases with
$x$
. Following the scaling arguments and numerical simulations outlined in § 3.1, we write
\begin{align} \Gamma _0=\frac {3\beta }{2\gamma }(x-x_0)+\cdots, \quad c_0=\frac {3\beta }{2\gamma }(x-x_0)C(\eta ) + \cdots, \end{align}
where the similarity variable
$\eta =y/(\alpha (x-x_0))^{1/3}$
is obtained by balancing advection and diffusion in the bulk-transport equation (C1). Substituting (C49) into (C1)–(C6), the bulk-surfactant equation (C1) in
$\mathcal{D}_1$
simplifies to
\begin{align} C_{\eta \eta } + \eta ^2 C_{\eta }/2 - 3\eta C/2 = 0, \quad C(0)=1, \quad C(\infty ) = 0. \end{align}
Equation (C50) is the same as (C16), whose solution is given in (C17
b). At the interface
$x\in [x_0, \, \phi ]$
and
$Y=0$
, the interfacial-surfactant equation (C2) becomes
\begin{align} \beta \left [ u_0(x,\,0)(x-x_0) \right ]_x = (\alpha (x-x_0))^{2/3}C^{\prime}(0). \end{align}
Integrating (C51) and applying the no-flux condition
$u_0(\phi, \,0)=0$
gives
\begin{align} \beta u_0(x,\,0) = - \frac {3\alpha ^{2/3} C^{\prime}(0)}{5} \left (\frac {(\phi -x_0)^{5/3}}{x-x_0}-(x-x_0)^{2/3}\right ), \end{align}
where
$C^{\prime}(0) = - 3^{5/3}/(2^{4/3}{\Gamma }({1}/{3}))\approx - 0.92$
. This profile successfully captures the COMSOL prediction in figure 5(d). In figure 5(d), the small discrepancy between the asymptotic results (C35), (C52) and the 2-D long-wave model arises from the retention of surface diffusion in its numerical approximation. Matching (C52) to the surface velocity in the transition region (C48) yields
\begin{align} 5\lambda \gamma = -6 \alpha ^{1/6} C^{\prime}(0) (\phi -x_0)^{5/3}. \end{align}
The constants
$\lambda$
and
$x_0$
have to be evaluated numerically for given
$\alpha$
,
$\beta$
,
$\gamma$
and
$\phi$
using (C35) and (C53). Thus, in the quasi-stagnant region of
$A^{{2D}}$
, the interfacial concentration is given by
\begin{align} \Gamma _0 \approx \frac {3\beta }{2\gamma }\left (x - x_0\right ) = \frac {3\beta }{2\gamma }\left (x - \phi + \left (\frac {- 5\lambda \gamma }{6\alpha ^{1/6}C^{\prime}(0)}\right )^{3/5} \right ), \end{align}
using (C53). As anticipated, the leading-order interfacial surfactant concentration (C54) obtained in the quasi-stagnant-cap region is linear and equal to the profiles obtained in the two previous short subregions: the deceleration (Appendix C.2.2) and the transition regions (Appendix C.2.3). This linear profile agrees with our numerical results in region
$A^{{2D}}$
for the quasi-stagnant part of the interface, i.e. for
$x_0\lt x\lt \phi$
, as depicted in figure 5(b). Only now can we substitute (C54) into (2.44), to obtain the leading-order drag reduction in region
$A^{{2D}}$
, which is given by
\begin{align} {DR}_0 = 1 - \frac {1}{2\phi }\left (\frac {-5\lambda }{6 C^{\prime}(0)} \right )^{3/5}\frac {\gamma ^{3/5}}{\alpha ^{1/10}} \equiv 1 - \frac {a_1 \gamma ^{3/5}}{\alpha ^{1/10}}, \end{align}
as in (3.2
b). Equations (C54) and (C55) show how bulk diffusion, captured in the non-dimensional coefficient
$\alpha ^{1/10}$
, has a weak but non-zero influence on the interfacial surfactant profile that determines
$DR_0$
. The solution is valid as long as the layers are nested (
$\gamma \lesssim \alpha ^{-3/2}$
) and
$x_0 \gt -\phi$
, such that
$1-{DR}_0 \lt 1$
. In order to determine
$x_0$
,
$\lambda$
and
$a_1$
, we solve (C35) and (C53) numerically to get
$\lambda$
and
$x_0$
for given
$\alpha$
,
$\beta$
,
$\gamma$
and
$\phi$
, from which we can evaluate
$a_1 = (-5\lambda /(6C^{\prime}(0)))^{3/5}/(2\phi )$
and therefore
${DR}_0$
.
In summary, to acquire the solution in the quasi-stagnant region, where the surfactant gradient decelerates the flow towards the downstream stagnation point, (C54)–(C55), we have resolved a nested asymptotic structure which includes deceleration (Appendix C.2.2) and transition (Appendix C.2.3) regions. The nested asymptotic structure has been matched to a slip region at the upstream end of the plastron (Appendix C.2.1), where the surfactant gradient has a weak effect on the flow and slips from the upstream stagnation point. These four asymptotic regions are detailed in figure 4(b), which we used to describe and explain the numerical solution that we see in figure 3(f). Using the expressions for
$u_0(x,\,0)$
and
$\lambda$
in the quasi-stagnant and slip regions, (C52) and (C35
b), respectively, the surfactant desorbed across the quasi-stagnant region is given by
\begin{align} \int _{x_0}^\phi c_{0y}(x,\,0)\,\text{d}x = - \left (\frac {9 \alpha (x_0+\phi )}{3\pi }\right )^{1/2}+\cdots, \end{align}
which matches the total amount of surfactant adsorbed in the slip region in (C34), as required by the conservation of surfactant.
C.3.
Weak Marangoni effect and weak bulk–surface exchange: region
$A^{{2D}}_E$
In region
$A^{{2D}}_E$
, we simplify the boundary-layer equations (C1)–(C6), where bulk diffusion and bulk–surface exchange are weak compared with Marangoni effects and advection, i.e.
$\alpha \ll 1$
and small
$\nu$
. In this regime, there is a one-way coupling of the bulk concentration onto the interfacial concentration, with
$c_0 \approx 1$
everywhere, and the interfacial surfactant distribution is in the classical stagnant-cap regime (He et al. Reference He, Maldarelli and Dagan1991), with a streamwise-averaged concentration imposed by the bulk concentration that is approximately equal to 1. Here, the surfactant gradient generated by the flow renders the downstream end of the interface no slip and the upstream end of the interface shear-free. We expand the surfactant field as follows:
\begin{align} c_0 = 1 + \alpha C_1(x,\,y) + \cdots, \quad \Gamma _0 = G_0(x) + \alpha G_{1}(x) + \cdots, \end{align}
and substitute (C57) into (C1)–(C6). In the interfacial surfactant equation (C2), there is no flux of bulk surfactant onto the interface, such that
$c_{0y}=0$
. For
$x\in [-\phi, \, \phi ]$
and
$y=0$
, we have
\begin{align} \frac {3\beta }{4} G_{0} - \frac {\gamma }{2} G_{0x}G_{0} = 0, \end{align}
at leading order, and where we used the no flux of interfacial surfactant condition (C3). The surfactant flux condition (C6b
) in
$\mathcal{D}_1$
with
$c_0 \sim 1$
reduces to
\begin{align} \int _{-\phi }^\phi G_{0} \,\text{d}x = 2\phi . \end{align}
A piecewise linear solution to the nonlinear ordinary differential equation (C58) exists for
$2 \phi \gt 4\gamma / 3\beta$
, such that in
$A^{{2D}}_E$
the interfacial-surfactant distribution is given by
\begin{align} G_{0} = \begin{cases} 0\quad \text{for} \quad -\phi \lt x \lt x_0, \\ \frac {3\beta }{2\gamma } (x - x_0) \quad \text{for} \quad x_0 \lt x \lt \phi . \end{cases} \end{align}
This piecewise linear stagnant-cap profile for the interfacial surfactant concentration agrees with our numerical results shown in figure 7(b). We evaluate
$x_0$
using the no-net-flux condition (C59), such that
\begin{align} x_0 = \phi - \left (\frac {8 \phi \gamma }{3\beta }\right )^{1/2}, \end{align}
where
$x_0\gt -\phi$
to ensure existence of the no-surfactant part of the profile in (C60). Hence, for
$2\phi \gt 4\gamma / 3\beta$
, the leading-order drag reduction calculated by substituting (C60) into (2.44) is given by
\begin{align} {DR}_0 = 1 - \left (\frac {2 \gamma }{3 \phi \beta }\right )^{1/2}. \end{align}
This solution for the
$A^{{2D}}_E$
region, already given in (3.4c
), is equivalent to that in region
$A^{{1D}}_E$
in the strong-cross-channel-diffusion problem (Tomlinson et al. Reference Tomlinson, Gibou, Luzzatto-Fegiz, Temprano-Coleto, Jensen and Landel2023a
). It agrees with numerical results shown in figure 6.
Open access
