1. Introduction
The presence of surfactants can significantly impact the flow in the vicinity of fluid interfaces. For bounded fluid interfaces, the convection of interface-bound surfactant molecules can lead to their stacking up at stagnation points, which in turn results in Marangoni stress, effectively rendering the interface incompressible through stress directed along the negative concentration gradient (Manikantan & Squires Reference Manikantan and Squires2020). Well-known examples are the formation of a stagnant surfactant film in front of an obstacle piercing the surface of a flowing liquid (Merson & Quinn Reference Merson and Quinn1965; Scott Reference Scott1982; Harper Reference Harper1992) and the reduced mobility of drops or bubbles moving through a liquid due to an increased surfactant concentration on the downstream portion of their interface (Savic Reference Savic1953; Levich Reference Levich1962; Davis & Acrivos Reference Davis and Acrivos1966; Sadhal & Johnson Reference Sadhal and Johnson1983).
Special attention has been given to liquid flow over ‘superhydrophobic’ surfaces with gas-filled (or liquid-filled) cavities due to their potential for drag reduction (Rothstein Reference Rothstein2010; Schönecker, Baier & Hardt Reference Schönecker, Baier and Hardt2014; Lee, Choi & Kim Reference Lee, Choi and Kim2016). The influence of surfactants on such flows has been observed experimentally (Kim & Hidrovo Reference Kim and Hidrovo2012; Bolognesi, Cottin-Bizonne & Pirat 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; Li et al. Reference Li, Li, Tan, Wang, Huang, Xiang, Lv and Duan2020) and analysed theoretically (Gaddam et al. Reference Gaddam, Agrawal, Joshi and Thompson2018; Landel et al. Reference Landel, Peaudecerf, Temprano-Coleto, Gibou, Goldstein and Luzzatto-Fegiz2020; Baier & Hardt Reference Baier and Hardt2021; Temprano-Coleto et al. Reference Temprano-Coleto, Smith, Peaudecerf, Landel, Gibou and Luzzatto-Fegiz2021). Also in this case, surfactant molecules stack up at stagnation points and can severely impede the drag reduction compared with expectations based on a surfactant-free situation. In the extreme case of large surfactant concentrations or low enough shear rates, Marangoni stress at the interface can completely balance the viscous shear stress, rendering the interface immobile (Peaudecerf et al. Reference Peaudecerf, Landel, Goldstein and Luzzatto-Fegiz2017; Baier & Hardt Reference Baier and Hardt2021). One such situation where this can occur is homogeneous shear flow over a flat gas–liquid interface embedded in an entirely flat surface, resulting in a constant shear rate on the interface. However, when the interface experiences an inhomogeneous shear stress, for example on a curved gas–liquid interface protruding into the channel (Song et al. Reference Song, Song, Hu, Du, Du, Choi and Rothstein2018) or sidewalls imposing an inhomogeneous shear field (Li et al. Reference Li, Li, Tan, Wang, Huang, Xiang, Lv and Duan2020), recirculation zones can occur on the interface at sufficiently high surfactant concentrations. A similar phenomenon is observed when vesicles attached to a flat interface are subjected to shear flow, where recirculation becomes observable in the lipid bilayer enclosing the vesicle (Woodhouse & Goldstein Reference Woodhouse and Goldstein2012; Honerkamp-Smith et al. Reference Honerkamp-Smith, Woodhouse, Kantsler and Goldstein2013). In reality, flat gas–liquid interfaces are more an exception than a rule. Gas dissolution in the liquid or the pressure drop in a channel with superhydrophobic walls deforms gas–liquid interfaces. It was shown experimentally that the interface curvature influences the flow over bubble mattresses very significantly (Steinberger et al. Reference Steinberger, Cottin-Bizonne, Kleimann and Charlaix2007; Karatay et al. Reference Karatay, Haase, Visser, Sun, Lohse, Tsai and Lammertink2013). In contrast, how the curvature of a surfactant-laden interface influences the flow along the interface has remained largely unexplored.
We aim here at a qualitative understanding of shear-driven flow along a long narrow gas-filled groove in a planar surface, similar to the experimental set-up employed by Song et al. (Reference Song, Song, Hu, Du, Du, Choi and Rothstein2018), by modelling the surfactants as an incompressible non-viscous surface fluid. We thus assume that even small variations in surface concentrations impose such large Marangoni stresses on the bulk fluid that the surface flow is rendered virtually incompressible (large Marangoni number), while viscous stresses within the surface fluid are insignificant compared with viscous stresses in the bulk (small Boussinesq number) (Barentin et al. Reference Barentin, Muller, Ybert, Joanny and Di Meglio2000; Elfring, Leal & Squires Reference Elfring, Leal and Squires2016; Manikantan & Squires Reference Manikantan and Squires2020). This simplification allows for an analytical solution of the flow field for small deflections of the gas–liquid interface away from the planar surface. For a shear-free gas–liquid interface without surfactants, the analogous situation of pressure-driven or shear-induced flow over surfaces with grooves containing gas pockets with curved menisci has been analysed numerically (Ng & Wang Reference Ng and Wang2011; Li, Alame & Mahesh Reference Li, Alame and Mahesh2017; Ageev, Golubkina & Osiptsov Reference Ageev, Golubkina and Osiptsov2018; Alinovi & Bottaro Reference Alinovi and Bottaro2018), analytically (Sbragaglia & Prosperetti Reference Sbragaglia and Prosperetti2007; Crowdy Reference Crowdy2010, Reference Crowdy2015, Reference Crowdy2016; Asmolov, Nizkaya & Vinogradova Reference Asmolov, Nizkaya and Vinogradova2018; Kirk Reference Kirk2018) and experimentally (Karatay et al. Reference Karatay, Haase, Visser, Sun, Lohse, Tsai and Lammertink2013; Kim & Park Reference Kim and Park2019). Similarly, thermocapillary flow along such surfaces due to a thermal gradient along the grooves (Baier, Steffes & Hardt Reference Baier, Steffes and Hardt2010) has received considerable attention in the case of curved menisci (Kirk et al. Reference Kirk, Karamanis, Crowdy and Hodes2020; Yariv & Crowdy Reference Yariv and Crowdy2020; Yariv & Kirk Reference Yariv and Kirk2021). Both of these situations bear some similarities to the one involving an incompressible surface fluid, and similar techniques can be employed in their solution.
2. Mathematical model
The configuration under investigation is sketched in figure 1, showing the view along the flow direction in (a) and a side view in (b). A Newtonian liquid of viscosity $\mu$ is driven by application of a shear stress
$\tau _\infty$ along the
$z$-axis far above a planar surface containing a single long, narrow gas-filled groove of width
$2a$ and length
$2b$, with
$a \ll b$, oriented in the same direction as the applied shear stress. The region of interest is an
$x$–
$y$-plane perpendicular to the flow direction not too close to the ends of the groove at
$|z| = b$, such that the flow is described sufficiently well by a unidirectional velocity field
$w(x,y)$ along the groove, unaffected by its finite length. In this region, the deflection
$y=h(x)$ of the gas–liquid interface
${\mathcal {S}}$ can be assumed to have the form of a circular arc, independent of the
$z$-coordinate along the groove, and we assume the flow to be slow enough that the shape of the gas–liquid interface is unaffected by viscous stresses acting on it. Furthermore, the viscosity of the gas in the cavity is assumed to be small enough that viscous stresses acting on the interface from the gas are negligible compared with stresses exerted by the liquid.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_fig1.png?pub-status=live)
Figure 1. (a) Sketch of the geometry with a liquid above a surface with a long, narrow gas-filled cavity of width $2a$ and length
$2b$. The gas–liquid interface,
${\mathcal {S}}$, is assumed to have the shape of a circular arc and is laden with an insoluble surfactant. Application of a shear stress
$\tau _\infty$ far away from the surface and in a direction normal to the
$x$–
$y$-plane drives a Couette flow along the groove. (b) Side view of the geometry. The region of interest is sufficiently far away from the ends of the groove, such that the deflection
$h(x)$ of the gas–liquid interface can be considered independent of
$z$.
The gas–liquid interface is assumed to be covered by an incompressible, inviscid surface fluid, insoluble in the liquid. Since the flow is unidirectional, the mass conservation equation for the surface fluid becomes
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn1.png?pub-status=live)
where the integral is along the part ${\mathcal {S}}$ of the gas–liquid interface cut by the
$x$–
$y$-plane containing the liquid domain
$\varOmega$ of interest. In the limit of vanishing Reynolds number, the local momentum conservation for the liquid together with the local form of the mass conservation equation for the surface fluid can be obtained by minimizing the dissipation rate (Batchelor Reference Batchelor2000; Kim & Karrila Reference Kim and Karrila2005) subject to the constraint (2.1),
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn2.png?pub-status=live)
where $c$ is a Lagrange multiplier. Expanding
$\delta \varPhi = \varPhi [w+\delta w] - \varPhi [w]$ to first order in the variation
$\delta w$, which is assumed to vanish on all boundaries of
$\varOmega$ except for
${\mathcal {S}}$, and setting it equal to zero, leads to
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn3.png?pub-status=live)
and, with the normal vector $\boldsymbol {n}$ pointing into the liquid domain,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn4.png?pub-status=live)
The former of these is the Stokes equation for unidirectional shear-driven flow, while the latter describes a constant shear stress at the interface, reminiscent of a situation with thermal Marangoni flow along a groove with a constant temperature gradient along its surface (Baier et al. Reference Baier, Steffes and Hardt2010). However, while Marangoni flow is induced towards regions of higher surface tension only, here, the boundary condition (2.1) requires a backflow on parts of the interface. This is particularly striking when the interface is flat; in this case the shear rate dictated by the far-field condition,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn5.png?pub-status=live)
extends all the way to the planar surface, since the velocity field and Lagrange multiplier
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn6.png?pub-status=live)
solve the Laplace equation (2.3) with boundary conditions (2.4), (2.5) and (2.1). Thus, the velocity vanishes at the gas–liquid interface and the Lagrange multiplier $c_0$ is the gradient of the surface pressure within the incompressible surface fluid. In the following, we will use this velocity field as the starting point of a perturbation expansion in the dimensionless deflection
$\varepsilon = h(0)/a$ at the centre of the gas–liquid interface.
We remark that the boundary condition (2.4) is compatible with a constant Marangoni stress opposing the main flow in the limit of large Marangoni number when surface compressibility is small. In this case, the surfactant concentration remains virtually constant on the entire interface and even small gradients in surfactant concentration lead to large stresses opposing a compression of the surfactant layer. This is briefly explored in Appendix A.
2.1. Dimensionless formulation
Using the length scale $a$ and the velocity scale
$u_0=a\tau _\infty /\mu$ allows us to introduce dimensionless coordinates
$(X,Y)=(x/a,y/a)$ and a dimensionless velocity
$W(X,Y)=w(aX,aY)/u_0$ in
$\varOmega$. The integral boundary condition (2.1) then reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn7.png?pub-status=live)
and the no-slip condition on the planar surface $\bar {{\mathcal {L}}}$, condition (2.4) on the gas–liquid interface
${\mathcal {S}}$ and the far-field condition (2.5) become
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn8.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn9.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn10.png?pub-status=live)
where $C=\mu c/\tau _\infty$ is the dimensionless Lagrange multiplier and
$\widetilde {\boldsymbol {\nabla }}$ is the gradient in the dimensionless coordinates
$(X,Y)$. Since
$W$ solves the Laplace equation (2.3),
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn11.png?pub-status=live)
shear-driven Stokes flow can be written as the imaginary part of an analytic function
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn12.png?pub-status=live)
by introducing the complex coordinate $Z=X+{\rm i}Y$. Denoting the complex derivative of
$f$ as
$f'$, and writing
$\mathrm {Re}[\,{\cdot }\,]$ and
$\mathrm {Im}[\,{\cdot }\,]$ for real and imaginary parts, the Cauchy–Riemann conditions,
$\mathrm {Re}[f']=\partial _X V= \partial _Y W$ and
$\mathrm {Im}[f']=-\partial _Y V=\partial _X W$, allow rewriting of the boundary conditions (2.8)–(2.10) as conditions for
$f$ and its derivative.
2.2. Domain perturbation
We parameterize the deflection $y=h(x)$ of the gas–liquid interface
${\mathcal {S}}$ by its maximal dimensionless deflection
$\varepsilon = h(0)/a$ at its centre. The radius of curvature of its circular arc then becomes
$r=(a^{2}+h(0)^{2})/(2h(0))=a(\varepsilon ^{-1}+\varepsilon )/2$, counted as positive for
$\varepsilon >0$, that is, when the deflection is into the upper half-plane. Up to second order in
$\varepsilon$ the deflection then becomes
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn13.png?pub-status=live)
where $\sigma (\varepsilon )=\{1 \text { for } \varepsilon >0; -1 \text { for } \varepsilon <0; 0 \text { for } \varepsilon =0 \}$ is the sign function. Thus, the dimensionless deflection
$H(X)=h(aX)/a$ has the expansion
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn14.png?pub-status=live)
For small values of the dimensionless deflection $\varepsilon$, the velocity field
$W$ can be found using a domain perturbation method (Leal Reference Leal2007), by projecting the boundary conditions (2.9) and (2.7) onto
$\mathcal {L}$, the projection of
${\mathcal {S}}$ into the real axis (cf. figure 1a). For this, we use a regular perturbation expansion in
$\varepsilon$ for
$f$,
$W$ and
$C$,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn15.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn16.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn17.png?pub-status=live)
where we have used the fact that, for a flat interface, $f_0(Z)=Z$,
$W_0(X,Y) = Y$ and
$C_0=1$. The boundary condition (2.9) projected onto
${\mathcal {L}}$ then becomes, up to second order in
$\varepsilon$,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn18.png?pub-status=live)
where the Laplace equation (2.11) was used in the second equality. Similarly, for the projected integral boundary condition (2.7) we obtain
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn19.png?pub-status=live)
Note that the measure for the arclength containing the square root does not contribute to the expansion to order $\varepsilon ^{2}$ in (2.19).
2.3. Solution
At each order in $\varepsilon$ we successively seek solutions for
$f_i(Z)$ and
$W_i(X,Y)$ in the whole upper half-plane
$\varOmega _0=\{Z\,|\,\mathrm {Im}[Z]\geq 0\}$, obeying the boundary conditions (2.8) on the no-slip surface
$\bar {\mathcal {L}}$, (2.10) in the far field and the projections (2.18) and (2.19) on
${\mathcal {L}}$.
2.3.1. Order
$\varepsilon ^{0}$
For the flat surface we already established
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn20.png?pub-status=live)
2.3.2. Order
$\varepsilon ^{1}$
Using (2.8), (2.18) and (2.10), $W_1(X,Y)$ obeys the boundary conditions
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn21.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn22.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn23.png?pub-status=live)
We recognize this as the boundary conditions for the flow above an infinite plane, driven by a constant shear rate $C_1$ along a strip of constant width on an otherwise no-slip surface. The well-known solution is (Philip Reference Philip1972)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn24.png?pub-status=live)
Note that $f_1(Z)$ and hence
$W_1(Z)$ vanish far from the surface. While the far-field boundary condition explicitly only forces the shear rate to vanish, this behaviour is expected, as the momentum flux occurs between the surface of the groove (where momentum is introduced) and the surrounding no-slip surface. Inserting
$W_1$ into (2.19), we obtain for the Lagrange multiplier
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn25.png?pub-status=live)
2.3.3. Order
$\varepsilon ^{2}$
The boundary conditions for $W_2$ inferred from (2.8), (2.18) and (2.10) become
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn26.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn27.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn28.png?pub-status=live)
It is instructive to rewrite these as conditions on the real and imaginary parts of the analytic function $f_2'$, using the fact that
$\partial _X W_2 = 0$ on
$\bar {\mathcal {L}}$, together with the Cauchy–Riemann conditions,
$H_1(X)=1-X^{2}$, and (2.24)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn29.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn30.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn31.png?pub-status=live)
We recognize this as a mixed boundary value problem for $f'_2(Z)$ on the real line, which can be converted into a Riemann–Hilbert problem for which the methods of solution are well developed (Gakhov Reference Gakhov1966; Lawrentjew & Schabat Reference Lawrentjew and Schabat1967; Muskhelishvili Reference Muskhelishvili2008). In particular, using the Keldysh–Sedov formula (Gakhov Reference Gakhov1966, § 46.3; see also Appendix B) we obtain
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn32.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn33.png?pub-status=live)
The line integrals above can be performed by standard techniques (England Reference England2003; Muskhelishvili Reference Muskhelishvili2013; Gogolin Reference Gogolin2014), and the term proportional to $A$ is a solution obeying homogeneous boundary conditions (i.e.
$\mathrm {Im}[f'(X)]=0$ on
$\bar {\mathcal {L}}$,
$\mathrm {Re}[f'(X)]=0$ on
${\mathcal {L}}$,
$f'(X+{\rm i}\infty )=0$). The function
$f_2$ can be obtained by taking the antiderivative of (2.33),
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn34.png?pub-status=live)
where we have taken into account that, just as for $W_1$, the velocity field
$W_2$ vanishes in the far field, requiring
$A=0$. The Lagrange multiplier
$C_2$ is obtained from (2.19) as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn35.png?pub-status=live)
This concludes our determination of the velocity field as a perturbation series up to second order in the dimensionless deflection $\varepsilon$. As an aside, we note that the solution (2.24) at order
$\varepsilon ^{1}$ also appears as the last part of the solution (2.34) at order
$\varepsilon ^{2}$, corresponding to a constant shear rate
$C_2$ instead of
$C_1$ on
${\mathcal {L}}$.
2.4. Branch cuts and choice of analytic functions
The principal branches of the analytic functions $\sqrt {Z}$ and
$\log (Z)$ both have branch cuts on the negative real axis. By the nature of the solution method, some of the functions in our solution are discontinuous across the section
$-1< X<1$ of the real axis. However, for
$\varepsilon <0$ the functions sought after must be analytic across this line segment. We thus replace
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn36.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn37.png?pub-status=live)
which have branch cuts along the real axis except on $-1< X<1$, agree with their original expressions in the first quadrant and thus solve the Laplace equation in the region of interest while obeying the same boundary conditions when approaching the real line from above. We will assume these replacements to have been carried out in the expressions (2.24) and (2.34) for
$f_1$ and
$f_2$ without further mentioning them. In order to illustrate the resulting expressions, we have plotted their imaginary parts, i.e. their contribution to the velocity, in figure 2. In order to simplify the presentation we restrict the plot of
$\mathrm {Im} [f_2(Z)]$ to negative values, since the corresponding region (shown in grey) is of no interest for the problem at hand. Above the real axis both functions have a very similar appearance, while marked differences are apparent for
$Y<0$. This is an indication that our approximation may be better for positive deflection,
$\varepsilon > 0$, than for negative deflection.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_fig2.png?pub-status=live)
Figure 2. Contour plot of (a) $\mathrm {Im} [f_1(Z)]$ and (b)
$\mathrm {Im} [f_2(Z)]$. In the grey region
$\mathrm {Im} [f_2(Z)]>0$; since here
$X>1$ and
$Y<0$, this region lies outside the considered liquid domain. The interval between contour lines is 0.05.
2.5. Finite element calculation
For comparison with the analytical results, numerical calculations were performed using the commercial finite element solver COMSOL Multiphysics (version 5.6, COMSOL AB, Stockholm, Sweden), employing its ’Coefficient form PDE’ interface. The symmetry of the problem with respect to reflection at the $Y$-axis allows us to restrict the calculation to the region with
$X\ge 0$ of the liquid domain of figure 1(a). Specifically, the Laplace equation, (2.11), was solved in a quadratic domain
$0< X,Y< D$ of width and height
$D = 25$ with a circular-arc section, corresponding to the deflected gas–liquid interface, added below (or removed above) the
$X$-axis at
$X<1$. On the circular arc the integral conservation equation (2.7) is prescribed as a constraint, while a Dirichlet condition,
${W=0}$, enforces the no-slip condition on the rest of the bottom surface. A constant shear rate
$\partial _Y W = 1$ is applied on the top surface at
$Y=D$, and a vanishing shear rate,
$\partial _X W = 0$, is assumed on the left and right edges at
$X=0$ and
$D$, corresponding to a symmetry condition. The domain is discretized using quadratic Lagrange elements on a triangular mesh with cells of size
$h_B=0.05$ away from the surface and
$h_S=h_B/5$ on the circular arc, with a maximal element growth rate of 1.01. It was verified that in the range
$-0.3 < \varepsilon < 0.4$ the velocity at the centre of the interface at
$(X,Y)=(0,\varepsilon )$ changes by less than 0.01 % when quadrupling the domain size
$D$ and by less than 0.4 % when halving the element size
$h_B$ and
$h_S$, indicating that the results are independent of the grid, and the influence of the finite extent of the domain plays no significant role.
3. Results and discussion
Velocity fields $W(X,Y)=\mathrm {Im}[f(X+{\rm i}Y)]$ are shown in figure 3 for
$\varepsilon$ varying between
$-$0.2 and 0.3 in steps of 0.1. Due to the symmetry of the problem, only half of the cavity and its immediate vicinity is shown. As can be seen, for
$\varepsilon >0$ the velocity at the centre of the interface near
$X=0$ is positive and becomes negative towards the edge of the cavity close to
$X=1$. By contrast, for
$\varepsilon <0$ a backflow is induced at the centre of the cavity, while the flow velocity is positive towards the edge of the cavity. It is also apparent that the velocity on the interface is relatively small, and when leaving the interface into the fluid domain is quickly dwarfed by the increasing velocity due to the constant shear rate applied in the far field. As expected, compared with a pure Couette flow, the velocity profile in the vicinity of the groove attains slightly larger values for negative deflection and slightly smaller values for positive deflection. Due to the smallness of the interface velocity the situation is not much different from the case of a solid protrusion into the channel. Note that without an incompressible surface fluid the dimensionless velocity at the centre of a flat (
$\varepsilon =0$) shear-free interface is 1 (Philip Reference Philip1972), and thus the corresponding velocity field in the vicinity of the groove becomes markedly different.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_fig3.png?pub-status=live)
Figure 3. Contour plot of velocity $W(X,Y)$ for (a)
$\varepsilon =0.1$, (b)
$\varepsilon =0.2$, (c)
$\varepsilon =0.3$, (d)
$\varepsilon =-0.1$, (e)
$\varepsilon =-0.2$. The interval between contour lines is 0.025 and the contour
$W(X,Y)=0$ is shown as a dashed white line.
The velocity on the interface is more clearly shown in figure 4(a), where it is plotted for the same values of $\varepsilon$ as in figure 3. Here, it becomes particularly transparent that, with an incompressible surface fluid, the velocity only reaches a few per cent of the values that would be reached on a shear-free interface. At the same time, the corresponding velocities from the numerical calculations using the finite element discretization are shown as dashed lines. It is evident that the analytical solution agrees well with the numerical calculation in the chosen range of
$\varepsilon$, with the largest deviations occurring close to the centre of the cavity. The fact that backflow occurs on the interface, with flow in opposite directions close to the centre of the groove and close to its edges, is evidently a prerequisite for the mass conservation of the incompressible surface fluid. Its direction of transport at the centre reflects the fact that a deflection of the interface above or below the reference plane leads to a respectively increased or decreased shear rate on the interface compared with the planar surface, see figure 4(c).
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_fig4.png?pub-status=live)
Figure 4. (a) Velocity on the interface, $W(X,\varepsilon H_1(X))$ (solid lines) and the corresponding numerical results (dashed lines) for
$\varepsilon =-0.2$,
$-$0.1, 0, 0.1 0.2, 0.3. The arrow indicates increasing values of
$\varepsilon$. (b) Velocity
$W(0,\varepsilon )$ on the interface at
$X=0$ (blue solid line), projected velocity
$\tilde W_S(0)$ ((3.1), yellow dashed line) and corresponding numerical values (green circles). (c) Shear rate,
$\partial _Y W(0,\varepsilon )$, on the interface at
$X=0$ (blue solid line),
$C(\varepsilon )$ (yellow dashed line) and corresponding numerical values (green circles).
To further assess the quality of the analytical solution, the velocity close to the centre of the cavity is plotted as a function of $\varepsilon$ in figure 4(b) (blue line) together with corresponding values from the numerical calculations (green circles). Note that the quality of the approximation is not symmetric in
$\varepsilon$ and deviations become lager more quickly for negative deflections than for positive ones. This expected behaviour was already alluded to above in the discussion of the functions
$f_1(Z)$ and
$f_2(Z)$ in § 2.4.
Another estimate for the range of validity of the expansion can be obtained by investigating the expansion of $W_S(X) = W(X,\varepsilon H_1(X))$ to second order in
$\varepsilon$ as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn38.png?pub-status=live)
Since similar expansions are performed during the projection of boundary conditions in (2.18) and (2.19), the difference between $\tilde W_S(0)$ and
$W(0,\varepsilon )$ should indicate when the expansion becomes inadequate. Remarkably,
$\tilde W_S(0)$, shown as a yellow dashed line in figure 4(b), closely traces the numerical values obtained and is thus a faithful indicator of the range of validity of our expansion.
Equation (3.1) thus is an excellent approximation for determining the extremal velocity encountered at the midpoint of the channel in the interval $-0.2\lesssim \varepsilon \lesssim 0.3$. It can also be used to assess the extremum occurring towards the edges of the groove. When considering (3.1) to first order in
$\varepsilon$ only, the location of this extremum lies at
$X=\sqrt {1-C_1^{2}/4}\approx 0.905$, while for the full equation (3.1) the location
$X \approx 0.9 \pm 0.02$ of this extremum is only weakly dependent on
$\varepsilon$ in the considered interval of
$\varepsilon$ and thus
$\tilde {W}_S(0.9)$ remains within 1 % of the value obtained from (3.1) at the location of the extremum. Similarly, the position where the interfacial velocity changes sign varies little with
$\varepsilon$ and according to (3.1) lies between
$X=0.54$ and 0.51 in the considered interval of
$\varepsilon$, tightly straddling the position
$X=\sqrt {1-C_1^{2}}\approx 0.529$ of the zero of (3.1) to first order in
$\varepsilon$.
Another value of interest is the shear rate on the interface. In our expansion this can be obtained directly from $\partial _Y W(0,\varepsilon )=\mathrm {Re}[f'({\rm i}\varepsilon )]$. At the same time, from (2.9), the shear rate at the interface is encoded in the Lagrange multiplier
$C(\varepsilon )$. Again, any deviation between the two is an indication for the quality of our approximation. The curves for both these quantities are shown as the respective blue solid and yellow dashed lines in figure 4(c), together with the numerically obtained values as green circles. Note that, in the numerical calculations, the shear is constant on the entire interface. As can be seen, the shear rate at the centre of the interface is less sensitive to the approximation than the velocity. However, for the shear rate the discrepancies are expected to become largest towards the edges of the groove, and this is indeed reflected in our solution (not shown).
4. Conclusion and outlook
We have presented an analytical solution for shear-driven liquid flow along a single bounded gas-filled groove, embedded in an otherwise planar surface, when the gas–liquid interface protrudes slightly above or below the planar surface and is covered by an incompressible surface fluid. The flow velocities displayed in figures 3 and 4(a) show that backflow occurs close to the edges of the groove when the interface protrudes above the planar surface and close to its centre when it deflects below the plane. Evidently, such regions with backflow are mandatory since on a bounded groove as much surface fluid is transported along the flow direction as against it.
Qualitatively, the same behaviour of the interface flow predicted for an incompressible surface fluid was observed experimentally at a single gas-filled groove in the presence of large concentrations of surfactant by Song et al. (Reference Song, Song, Hu, Du, Du, Choi and Rothstein2018). We take this as evidence for the adequacy of our model for describing such situations. Nevertheless, there are some discrepancies, as experimentally the backflow on the interface is not always as prominent as in our prediction. This may be due to interfacial concentrations not being large enough for the surfactant film to become fully incompressible or due to surfactant dissolving in the liquid. Moreover, experiments were performed using pressure-driven Poiseuille flow in a relatively shallow channel instead of pure shear flow, affecting the details of the boundary conditions. While a more complex model for surfactants taking into account adsorption/desorption kinetics and an equation of state or an effective surface viscosity at large surfactant concentrations could be incorporated in a more complete model, it would be difficult to capture such details in an analytical description. Despite the mentioned shortcomings, we hope that our analytic model will be a valuable reference for describing flow with large interfacial concentrations of surfactants.
Flow over a superhydrophobic surface in the Cassie state is often characterized by reporting an apparent slip length (Rothstein Reference Rothstein2010; Lee et al. Reference Lee, Choi and Kim2016). It should therefore be of interest to extend the present study from a single groove to a parallel array of gas-filled grooves in order to investigate the impact of large surfactant concentrations on the observable slip length. Moreover, as mentioned in the introduction, thermal Marangoni flow along a grooved surface can also be approximately described by a constant shear stress along the grooves (Baier et al. Reference Baier, Steffes and Hardt2010). The similarity to the flow over an incompressible surface fluid promises some synergy between the investigations of both situations.
More generally, flow over arbitrary gas- or liquid-filled patches (e.g. circular holes) covered by an incompressible surface fluid, with interfaces protruding above or below a flat surface is of interest for various designs of superhydrophobic surfaces. In this respect, an extension to a slightly compressible surface fluid, possibly exhibiting surface viscosity, may be of interest for a more complete picture. Such configurations may also be relevant in other situations such as flow over vesicles attached to a solid wall (Woodhouse & Goldstein Reference Woodhouse and Goldstein2012; Honerkamp-Smith et al. Reference Honerkamp-Smith, Woodhouse, Kantsler and Goldstein2013).
On a higher level of abstraction, we hypothesize that the results reported in this paper may hint at a quite general class of fluid dynamic phenomena: that the flow along a surfactant-covered liquid surface is very sensitive to the surface deformation. In our case, the surface flow is suppressed on a flat surface, while a flow emerges on a deformed surface. It is conceivable that the flow field itself deforms a liquid surface, or that the surface deformation is controlled by an external parameter (such as pressure), which in turn would influence the flow pattern on the surface.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Surface incompressibility
In § 2 the interface was modelled as containing an incompressible surface fluid. In this appendix we briefly discuss the compatibility of this simple model with momentum conservation and transport of a surfactant species at the interface in the limit of large Marangoni number.
The steady-state interfacial stress balance in the Boussinesq–Scriven model reads (Edwards, Brenner & Wasan Reference Edwards, Brenner and Wasan1991; Slattery, Sagis & Oh Reference Slattery, Sagis and Oh2007; Manikantan & Squires Reference Manikantan and Squires2020)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn39.png?pub-status=live)
where $\boldsymbol {I}_s = \boldsymbol {I}-\boldsymbol {n}\boldsymbol {n}$ is the interface projection operator and
$\boldsymbol {\nabla }_s = \boldsymbol {I}_s\boldsymbol {\cdot }\boldsymbol {\nabla }$ the interface gradient. Also,
$p^{(+)}$,
$p^{(-)}$,
$\boldsymbol {\tau }^{(+)}$ and
$\boldsymbol {\tau }^{(-)}$ are the pressures and viscous stress tensors in the fluid on the side of the interface the normal vector
$\boldsymbol {n}$ points to and away from, respectively, with
$\boldsymbol {\tau }^{(\,{\cdot }\,)}=\mu ^{(\,{\cdot }\,)}(\boldsymbol {\nabla }\boldsymbol {u}^{(\,{\cdot }\,)}+(\boldsymbol {\nabla }\boldsymbol {u}^{(\,{\cdot }\,)})^{\rm T})$. On the interface it is assumed that the velocities of both fluid phases are identical,
$\boldsymbol {u}=\boldsymbol {u}^{(+)}=\boldsymbol {u}^{(-)}$. Here,
$\varPi (\varGamma ) = \gamma _0 - \gamma (\varGamma )$ is the surface pressure with
$\gamma (\varGamma )$ being the interfacial tension of an interface with surfactant concentration
$\varGamma$ and
$\gamma _0$ the interfacial tension of the clean interface. The interfacial rheology is captured in the Boussinesq–Scriven stress
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn40.png?pub-status=live)
where $\mu _s$ and
$\kappa _s$ are the surface shear and the surface dilatational viscosity, respectively.
In our model with a gas–liquid interface, we assume that the viscous stress $\boldsymbol {\tau }^{(-)}$ on the gas side is negligible compared with the viscous stress
$\boldsymbol {\tau }^{(+)}$ on the liquid side. For small Boussinesq numbers
${{Bq}}_\mu =\mu _s/(\mu a)$ and
${{Bq}}_\kappa =\kappa _s/(\mu a)$, the intrinsic surface stresses can be neglected compared with the stresses exerted by the adjacent fluid. In our case the
$z$-component of (A1) along the grove then reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn41.png?pub-status=live)
where in the last step we have assumed that the surfactant concentration has the form $\varGamma =\varGamma _0(1+\delta \tilde {\varGamma })$ with
$\delta \tilde {\varGamma } \ll 1$. It is useful to introduce the Marangoni modulus (Manikantan & Squires Reference Manikantan and Squires2020)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn42.png?pub-status=live)
as a measure of the interfacial elasticity or the amount of work needed for compressing an interface with surfactants. As in § 2.1 we use the length scale $a$ to introduce dimensionless coordinates, here writing
$\tilde {Z}=z/a$ for the
$z$-coordinate, and the velocity scale
$u_0=a\tau _\infty /\mu$. The non-dimensional form of (A3) is approximately
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn43.png?pub-status=live)
where we have introduced the Marangoni number
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn44.png?pub-status=live)
Consequently, for a constant shear stress along the groove, the surfactant gradient becomes arbitrarily small for large ${{Ma}}$, with
$\delta \tilde {\varGamma }$ remaining small for not too large groove lengths. Since the left-hand side of (A5) is constant in our model, this equation corresponds to (2.4) (or (2.9)).
The surfactant flux at the interface in the $z$-direction along the grove is governed by convection and diffusion,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn45.png?pub-status=live)
where we have substituted $\delta \tilde {\varGamma }$ according to (A5) and in the last step taken the limit of large Marangoni number, assuming that
$\tilde {Z}=z/a$ is not too large and the Péclet number
$Pe=au_0/D$ is not too small. Additionally, in our approximation the transverse velocity on the groove vanishes, which is consistent with a negligible gradient in the surfactant concentration in this direction. Thus, for large
${{Ma}}$ the surfactant concentration can be considered as constant,
$\varGamma _0$, on the interface, with
$\varPi$ taking the role of the pressure in the momentum equation for the incompressible surface fluid. Integrating (A7) across the width of the interface then leads to (2.1) of our model. Finally, the pressure difference between the phases is governed by the mean curvature
$H$ of the interface,
$2H=-\boldsymbol {\nabla }_s\boldsymbol {\cdot }\boldsymbol {n}$, such that
$p^{(+)} - p^{(-)} = 2 \gamma H$. As for large
${{Ma}}$ the surface tension
$\gamma$ stays nearly constant on the interface, this is consistent with the circular-arc cross-section of the interface assumed in our model.
Typical values of the shear rates employed in experiments are $\dot {\gamma }=\tau _\infty /\mu =0.1 \ldots 1\,\text {s}^{-1}$, for water with
$\mu \simeq 1$ mPa s, and a typical length scale is
$a = 0.1 \ldots 1$ mm. With the Marangoni modulus scaling as
$E_0 \simeq k_B T \varGamma _0$ at a surfactant concentration
$\varGamma _0 \simeq 0.01 \ldots 1 \text { nm}^{-2}$, typical values for the Marangoni number lie in the range
${{Ma}} \simeq k_B T \varGamma _0 / (a \tau _\infty ) \simeq 4\times (10^{1} \ldots 10^{5})$. With surface diffusion coefficients in the range
$D \simeq (0.1\ldots 1)\times 10^{-9}\, \mbox{m}^{2} \,\mbox{s}^{-1}$ the corresponding Péclet numbers are in the range
${{Pe}} = au_0/D\simeq a^{2} \dot {\gamma }/D\simeq 1 \ldots 10^{4}$. It is thus expected that in many experimental scenarios the assumption of an incompressible surfactant phase is well justified.
Appendix B. Outline of a derivation of the Keldysh–Sedov formula
In § 2.3 a holomorphic function was obtained in the upper half-plane, obeying certain boundary conditions on the real line, by using the Keldysh–Sedov formula. Here, we give a brief sketch of how to obtain this expression based on the behaviour of holomorphic functions at cuts on the real line. Standard references for these techniques are Muskhelishvili (Reference Muskhelishvili2008) or Gakhov (Reference Gakhov1966) while England (Reference England2003) gives an excellent brief introduction.
Consider a function $\varPhi (Z)$ that is holomorphic in the complex plane with the exception of possible cuts on intervals located on the real line. For such a function it is useful to define the limit when approaching a point
$X$ on the real line from above or below
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn46.png?pub-status=live)
When $\varPhi (Z)$ has purely real values in an interval on the real line, the Schwartz reflection principle,
$\overline{\varPhi (Z)}=\varPhi (\bar {Z})$, where the overbar denotes complex conjugation, leads to
$\overline {\varPhi ^{+}(X)}=\varPhi ^{-}(X)$. As an example, the function
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn47.png?pub-status=live)
has a cut on the interval $[-1,1]$ of the real line and obeys
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn48.png?pub-status=live)
According to the Sokhotski–Plemelj theorem (Gogolin Reference Gogolin2014, § 1.4.1)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn49.png?pub-status=live)
where $\mathcal {P}$ denotes the principal value of the integral. Hence, the function
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn50.png?pub-status=live)
is another example that obeys the jump condition
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn51.png?pub-status=live)
The conditions (B3) and (B6) can be used to obtain a holomorphic function $g(Z)$ in the upper half-plane that obeys boundary conditions as in (2.29) and (2.30)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn52.png?pub-status=live)
Introducing $\psi (Z) = R(Z)g(Z)$ and using (B3), these can be converted to
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn53.png?pub-status=live)
which has the form of (B6). Finally, using (B5), we obtain
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20220929194441093-0311:S0022112022007753:S0022112022007753_eqn54.png?pub-status=live)
Apart from the solution of the homogeneous problem, this is the Keldysh–Sedov formula of (2.32). Further generalizations and details on the conditions for applicability can be found in the standard references listed above.