2.1. Dimensional modelling

We consider the steady propagation of a two-dimensional bubble inside a channel of height $2\hat {H}$ (see figure 1). We orient the $\hat {x}$ -axis and the $\hat {z}$ -axis along and perpendicular to the lower wall, respectively. We define $\boldsymbol{i}$ and $\boldsymbol{k}$ as the unit vectors in the $\hat {x}$ - and $\hat {z}$ -directions, respectively. We assume that buoyancy effects are negligible, so the flow is symmetric across the centreline of the channel and we can restrict our attention to the region $0\leqslant \hat {z} \leqslant \hat {H}$ . We define the bubble’s surface in this region to be at $\hat {z}=\hat {h}(\hat {x})$ . The normal to the liquid surface pointing into the bubble and the tangent to the bubble’s surface are denoted by $\boldsymbol{n}$ , and $\boldsymbol{t}$ , respectively, and are given by

(2.1) \begin{align} \boldsymbol{n} & = \dfrac {-\hat {h}_{\hat {x}} \boldsymbol{i}+\boldsymbol{k}}{\sqrt {1+\hat {h}_{\hat {x}}^2}}, \quad \quad \boldsymbol{t} = \dfrac {\boldsymbol{i}+\hat {h}_{\hat {x}}\boldsymbol{k}}{\sqrt {1+\hat {h}_{\hat {x}}^2}}, \end{align}

where the subscript variable means partial differentiation with respect to that variable.

We assume the motion of the bubble is sufficiently slow that the flow is in the Stokes regime and we move into a frame of reference in which the bubble is stationary and the walls travel at a velocity $-\hat {U}_b\boldsymbol{i}$ , where $\hat {U}_b$ is the bubble propagation speed. We denote the liquid velocity as $\hat {\boldsymbol{u}}=(\hat {u},\hat {v})$ and the pressure as $\hat {p}$ . The motion of the liquid is then governed by the Stokes equations

(2.2a) \begin{align} \hat {\boldsymbol{\nabla }}\boldsymbol{\cdot }\hat {\boldsymbol{u}} &=0, \quad \quad\quad\;\;\; \text{in}\quad \hat {\varOmega }, \end{align}

(2.2b) \begin{align} \hat {\boldsymbol{\nabla \!}} \hat {p} &=\hat {\mu }\hat {{\nabla} }^2\hat {\boldsymbol{u}} \quad\quad\, \text{in}\quad \hat {\varOmega }, \end{align}

where $\hat {\mu }$ is the constant liquid viscosity, $\hat {\varOmega }$ denotes the liquid region and $\hat {\boldsymbol{\nabla }}=(\partial /\partial \hat {x},\partial /\partial \hat {z} )$ is the two-dimensional gradient operator.

The surfactant concentration on the surface of the bubble, $\hat {\varGamma }$ , is governed by the advection–diffusion–reaction equation (Stone Reference Stone1990)

(2.3) \begin{align} \frac {{\rm d}}{\textrm{d} \hat {s}}\big (\hat {\varGamma \!} u_{\!{S}} \big ) &= \hat {D}\frac {{d}^2\hat {\varGamma }}{\textrm{d}\hat {s}^2} +\hat {k}_1 \hat {C}-\hat {k}_2 \hat {\varGamma }, \quad \quad \text{on} \quad \hat {z}=\hat {h}(\hat {x}), \end{align}

where $\hat {s}$ is arclength, $\hat {D}$ is the surface diffusion coefficient, $\hat {C}$ is the bulk concentration and

(2.4) \begin{equation} {u \!}_{{S}}=\bigl .\hat {\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{t}\bigr |_{\hat {z}=\hat {h}} \end{equation}

is the tangential surface velocity. We have supplemented (2.3) with a linear reaction term, with rate constants $\hat {k}_1$ and $\hat {k}_2$ , because the surfactant concentrations in the bulk and on the surface are assumed to be close to equilibrium. There are numerous numerical studies that include nonlinear reaction kinetics, such as Fujioka & Grotberg (Reference Fujioka and Grotberg2005) and Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019). In general one would solve a coupled advection–diffusion equation for $\hat {C}$ but, as mentioned above, we focus on the ‘bulk equilibrium’ limit in which the bulk concentration does not vary significantly, so we take $\hat {C}$ to be a known constant. We further assume that the surfactant quickly adsorbs or desorbs onto the surface of the bubble and thus the surfactant concentration is close to equilibrium. This implies that surfactant cannot accumulate and thus the bubble cannot rigidify, a phenomenon seen in surfactant systems without fast reactions (see e.g. Park Reference Park1992).

On the Hele-Shaw cell wall, we supply the no-slip boundary condition

(2.5a) \begin{align} \hat {\boldsymbol{u}} &= -\hat {U}_b\boldsymbol{i} \quad \quad \text{on} \quad \hat {z}=0. \end{align}

On the bubble’s surface, we supply a kinematic condition and normal and tangential stress balances as follows:

(2.5b) \begin{align} \hat {\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{n} &= 0 \quad \quad\quad \quad\quad \quad \text{on} \quad \hat {z}=\hat {h}(\hat {x}), \end{align}

(2.5c) \begin{align} \boldsymbol{n}\boldsymbol{\cdot }\hat {\boldsymbol{\sigma }}\boldsymbol{\cdot }\boldsymbol{n} &= -\hat {p}_b+\hat {\gamma } \hat {\kappa } \quad \quad\; \text{on} \quad \hat {z}=\hat {h}(\hat {x}), \end{align}

(2.5d) \begin{align} \boldsymbol{t}\boldsymbol{\cdot }\hat {\boldsymbol{\sigma }}\boldsymbol{\cdot }\boldsymbol{n} &= \frac {\textrm{d}\hat {\gamma }}{\textrm{d}\hat {s}} \quad \quad\quad \quad\quad \text{on} \quad \hat {z}=\hat {h}(\hat {x}), \end{align}

where $\hat {\gamma }$ is the (no longer constant) surface tension, $\hat {\kappa }$ is the curvature, $\hat {p}_b$ is the constant pressure inside the bubble and $\hat {\boldsymbol{\sigma }}$ is the viscous stress tensor, given by

(2.5e) \begin{align} \hat {\boldsymbol{\sigma }}=-\hat {p}\boldsymbol{I}+\hat {\mu }\big (\hat {\boldsymbol{\nabla }}\hat {\boldsymbol{u}} +{\hat {\boldsymbol{\nabla }}\hat {\boldsymbol{u}}}^{\text{T}}\big ), \end{align}

where $\boldsymbol{I}$ denotes the two-dimensional identity tensor and the superscript $\text{T}$ denotes the transpose.

We close our model with an equation of state, which relates the surface tension, $\hat {\gamma }$ , to the surfactant surface concentration, $\hat {\varGamma }$ . Since we are assuming that the surfactant concentrations in the bulk and on the surface are close to equilibrium, we also supply a linear equation of state, i.e.

(2.6) \begin{align} \hat {\gamma } = \hat {\gamma }_0 + \frac {\textrm{d}\hat {\gamma }}{\textrm{d}\hat {\varGamma }}\bigg |_{\hat {\varGamma }=\hat {\varGamma }_0 }(\hat {\varGamma }-\hat {\varGamma }_0), \end{align}

where $\hat {\varGamma }_0 = \hat {k}_1\hat {C}/\hat {k}_2$ is the equilibrium concentration of surfactant and $\hat {\gamma }_0$ is the surface tension at equilibrium.

2.2. Non-dimensionalisation

We non-dimensionalise the system (2.2)–(2.6) as follows (in which dimensionless variables are denoted without hats):

(2.7) \begin{align} (\hat {x},\hat {z},\hat {s},\hat {h}(\hat {x}))&=\hat {H}(x,z,s,h(x)), \quad \quad \hat {\boldsymbol{u}}= \hat {U}_b\boldsymbol{u}, \quad \quad \hat {\gamma }=\hat {\gamma }_0 \gamma , \nonumber \\ (\hat {p},\hat {p}_b)&=\frac {\hat {\gamma }_0}{\hat {H}}(p,p_b), \quad \quad\quad \quad\;\, \hat {\kappa }=\frac {\kappa }{\hat {H}}, \quad \quad \hat {\varGamma } = \hat {\varGamma }_0 \varGamma . \end{align}

The dimensionless governing equations are given by

(2.8) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{u} &= 0, \quad\quad \boldsymbol{\nabla \!}p \quad \quad= {\textit{Ca}} {\nabla} ^2\boldsymbol{u} \quad\quad \text{in }0\lt z\lt h(x). \end{align}

The dimensionless versions of the boundary conditions (2.5) are given by

(2.9a) \begin{align} \boldsymbol{u} &= -\boldsymbol{i} \quad\quad \text{on} \quad z=0, \end{align}

and, on the bubble’s surface $z=h(x)$ ,

(2.9b) \begin{align} \frac {{\rm d}}{\textrm{d} s}\big (\varGamma\! u_{\!{S}}\big ) &= \frac {1}{{\textit{Pe}}}\frac {\rm{d}^2\varGamma }{\textrm{d} s^2} +k (1-\varGamma ),\\[-12pt]\nonumber \end{align}

(2.9c) \begin{align} v&=uh_x,\\[-12pt]\nonumber \end{align}

(2.9d) \begin{align} p_b-p-\frac {2{\textit{Ca}}}{1+h_x^2}\big (\big (1-h_x^2\big )u_x +h_x(u_z +v_x)\big ) &= \frac {(1-{{\textit{Ma}}}(\varGamma -1) )h_{xx}} {\big (1+h_x^2\big )^{3/2}}, \end{align}

(2.9e) \begin{align} \frac {1}{\big (1+h_x^2\big )^{1/2}}\big (-4h_x u_x +\big (1-h_x^2\big )(u_z+v_x)\big ) &= -{\frac {{\textit{Ma}}}{{\textit{Ca}}}}\, \varGamma _x, \end{align}

where

(2.10) \begin{align} &\frac {\rm{d}}{\textrm{d} s} =\frac {1}{\sqrt {1+h_x^2}}\frac {\rm{d}}{\textrm{d} x} \quad\quad \text{and} \quad\quad u_{\!{S}} = \left .\frac {u+h_x\! v}{\sqrt {1+h_x^2}}\right |_{z=h(x)}. \end{align}

Figure 2.

Schematic of the front and rear of a bubble, showing the three regions of interest in each: regions 1, front- and rear-cap regions, regions 2, thin film regions and regions 3 transition regions.

The model (2.8)–(2.9) contains the dimensionless parameters

(2.11) \begin{align} {\textit{Ca}} &= \frac {\hat {\mu}\hat {U}_b}{\hat {\gamma }_0}, &{\textit{Pe}} &= \frac {\hat {H}\hat {U}_b}{\hat {D}}, & {{\textit{Ma}}}& {= -\frac {\hat {\varGamma }_0}{\hat {\gamma }_0}\,\frac {\textrm{d}\hat {\gamma }}{\textrm{d}\hat {\varGamma }},} & k& = \frac {\hat {k}_2 \hat {H}}{\hat {U}_b}, \end{align}

namely the capillary number, the Péclet number, the Marangoni number and the dimensionless reaction constant, respectively. We consider slow flow in which ${\textit{Ca}}\ll 1$ . Furthermore, we assume that ${\textit{Pe}}\gg 1$ , and hence we take the limit ${\textit{Pe}}\rightarrow \infty$ , so (2.9b ) reduces to an advection–reaction equation for $\varGamma$ . Finally, we also consider the regime where $k\gg 1$ so that the surface and bulk surfactants are approximately in equilibrium. We observe from (2.9b ) that $\varGamma \rightarrow 1$ as $k\rightarrow \infty$ and therefore define $\tau =k{\textit{Ca}}^{1/3}(\varGamma -1)$ . Here, the scaling of $\tau$ with ${\textit{Ca}}^{1/3}$ is designed to achieve dominant balances in the transition region between the capillary-static meniscus and the thin film on the channel wall (marked Region 3 in figure 2), as we will show below. We then identify a distinguished limit in which, although the surface concentration of surfactant is almost constant, the Marangoni stress at the free surface is retained in the model at leading order. To this end, we define the dimensionless elasticity parameter

(2.12) \begin{equation} E=\frac {{\textit{Ma}}}{{\textit{Ca}}k}=-\frac {\hat {\varGamma }_0}{\hat {\mu}\hat {H}\hat {k}_2}\,\frac {\textrm{d}\hat {\gamma }}{\textrm{d}\hat {\varGamma }}, \end{equation}

and take $E=O(1)$ while letting $k\rightarrow \infty$ . The surfactant conservation (2.9b ) and boundary conditions (2.9d )–(2.9e ) thus reduce to

(2.13a) \begin{align} {\textit{Ca}}^{1/3}\frac {\textrm{d}u_{\!{S}}}{\textrm{d}s} &= \tau, \end{align}

(2.13b) \begin{align} p_b-p-\frac {2{\textit{Ca}}}{1+h_x^2}\big (\big (1-h_x^2\big )u_x +h_x(u_z +v_x)\big ) &= \frac {(1-{\textit{ECa}}^{2/3}\tau )h_{xx}}{\big (1+h_x^2\big )^{3/2}}, \end{align}

(2.13c) \begin{align} \frac {{\textit{Ca}}^{1/3}}{\big (1+h_x^2\big )^{1/2}}\big (-4h_x u_x +\big (1-h_x^2\big )(u_z+v_x)\big ) &= -E\tau _x, \end{align}

all on $z=h(x)$ .

Following the above simplifying assumptions, the model now contains just two dimensionless parameters: ${\textit{Ca}}$ , which is assumed to be small, and $E$ , which we take to be $O(1)$ . Waters & Grotberg (Reference Waters and Grotberg2002) study a similar system in which $E\ll 1$ and the Marangoni stress is therefore negligible at leading order. The boundary conditions (2.13) represent the simplest model that contains both the Marangoni stress due to surface gradients in surfactant concentration (2.13c ) and the depletion of adsorbed surfactant by surface dilatation (2.13a ). The resulting system can also be interpreted as a model for surface viscosity (Gounley et al. Reference Gounley, Boedec, Jaeger and Leonetti2016).

2.3. Summary of assumptions

Here, we summarise the assumptions that have been made in this section.

1. (i) The bulk surfactant concentration $\hat {C}$ is sufficiently large that it does not vary significantly due to interaction with the surface, and we can treat it as effectively constant.

2. (ii) The surfactant is highly soluble, so the surface and bulk surfactant concentrations are approximately in equilibrium, i.e. $k\gg 1$ .

3. (iii) Surface diffusion of surfactant is negligible, so ${\textit{Pe}}\ll 1$ .

Following these assumptions, our model for a highly soluble surfactant is given by (2.8) with the boundary conditions (2.9a ), (2.9c ) and (2.13). In addition, we will assume shortly that the bubble propagates sufficiently slowly for the capillary number ${\textit{Ca}}$ to be small.

3.1. Regions and scalings

We consider the small- ${\textit{Ca}}$ limit and perform a perturbation expansion in powers of ${\textit{Ca}}^{1/3}$ , following the analysis presented by Park & Homsy (Reference Park and Homsy1984), for the similar problem of a surfactant-free bubble. As shown by Park & Homsy, in the small- ${\textit{Ca}}$ limit the problem splits into three regions of interest (see figure 2).

1. (i) The front-cap region. Here, the free surface is capillary static and hence, to leading order, is a circular cap.

2. (ii) The thin film region. Through a lubrication analysis, we find there is a thin liquid film of constant thickness between the channel wall and the bubble.

3. (iii) The transition region. Here, both viscous and capillary forces are important, which allows us to smoothly transition from the circular cap to the thin film region.

We include for completeness the analysis of the front-cap and thin film regions, which follows directly from Park & Homsy. It is in the transition region where the effect of the surfactant is included, and thus the equations deviate from those found by Park & Homsy.

3.2. Region 1: front-cap region

In this region we expand our variables in powers of ${\textit{Ca}}^{1/3}$ : $p(x,z)\sim p_0(x,z) +{\textit{Ca}}^{1/3} p_1(x,z) +{\cdots}$ and so forth. At leading order, the equations of motion (2.8) become

(3.1a,b) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_0 &= 0, \quad\quad \boldsymbol{\nabla \!}p_0 =\boldsymbol{0}, \end{align}

while the boundary conditions (2.13a ) and (2.13c ) are both satisfied identically by $\tau _0=0$ . Equation (3.1b ) implies that $p_0=\text{constant}$ , and the normal stress balances at the bubble’s surface (2.13b ) at leading order reads

(3.2) \begin{align} \Delta p_0 &= \frac {h_0^{\prime\prime}}{(1+h_0^{\prime 2})^{3/2}} \quad\quad \text{on} \quad z=h_0(x), \end{align}

where $\Delta p_0 = p_b - p_0$ is the leading-order difference between the constant pressure inside the bubble, $p_b$ , and the fluid pressure.

To find the leading-order shape of the bubble, we impose the conditions $h_0(x)\rightarrow 1$ and $h_0^{\prime}(x)\rightarrow \infty$ at the front tip of the bubble (which we define to be at $x=0$ ), and $h_0(x_c)=h_0^{\prime}{}(x_c)=0$ , where $x_c\lt 0$ is the a priori unknown location of the point where the leading-order meniscus encounters the cell wall. We thus find that $\Delta p_0=1$ and $x_c=-1$ , and the leading-order shape of the meniscus a circular cap of radius 1 (Park & Homsy Reference Park and Homsy1984), given by

(3.3) \begin{align} h_0(x) &= 1-\sqrt {1-(x+1)^2}, \quad\quad \text{for}\quad x\in (-1,0). \end{align}

In (1.2) we require knowledge of the pressure drop across the meniscus up to $O({\textit{Ca}}^{2/3})$ . Following Park & Homsy (Reference Park and Homsy1984), we find that $p_1,\,h_1\equiv 0$ ; then at $O({\textit{Ca}}^{2/3})$ we find that $p_2$ is a constant, denoted by $-\beta _1$ , and

(3.4) \begin{align} \beta _1 = \frac {h_2^{\prime\prime}}{(1+h_0^{\prime 2})^{3/2}}, \end{align}

which can be solved to give

(3.5) \begin{align} h_2(x)=\beta _1\sqrt {1-(x+1)^2}, \end{align}

where $\beta _1$ is the a priori unknown $O({\textit{Ca}}^{2/3})$ pressure correction.

3.3. Region 2: thin film region

In the thin film region, we rescale the variables as

(3.6) \begin{align} x& = \tilde {x}, \quad\quad z= {\textit{Ca}}^{2/3}\tilde {z}, \quad\quad h= {\textit{Ca}}^{2/3}\tilde {h}, \quad\quad p-p_b={\textit{Ca}}^{2/3}\tilde {p},\nonumber \\ u&=\tilde {u}, \quad\quad v= {\textit{Ca}}^{2/3}\tilde {v} \quad\quad\, \tau ={\textit{Ca}}^{1/3} \tilde {\tau }. \end{align}

Again we expand the variables in powers of ${\textit{Ca}}^{1/3}$ . At leading order, the equations of motion (2.8) become

(3.7a–c) \begin{align} \tilde {u}_{0\tilde {x}}+\tilde {v}_{0\tilde {z}} &= 0, \quad\quad \tilde {u}_{0\tilde {z}\tilde {z}} =0, \quad\quad \tilde {p}_{0\tilde {z}} =0. \end{align}

The kinematic boundary conditions (2.9a ) and (2.9c ) become

(3.8a) \begin{align} \tilde {u}_0=-1,\quad \tilde {v}_0&=0 \quad \quad\quad \quad \text{on} \quad \tilde {z}=0, \end{align}

(3.8b) \begin{align} \tilde {v}_0&=\tilde {u}_0\tilde {h}_{0\tilde {x}} \quad \quad\, \text{on} \quad \tilde {z}=\tilde {h}_0(\tilde {x}), \end{align}

and the remaining boundary conditions (2.13a )–(2.13c ) reduce to

(3.9a–c) \begin{align} \tilde {\tau }_0&=-\tilde {u}_{0x}, & \tilde {p}_0&=\tilde {h}_{0xx}, & \tilde {u}_{0\tilde {z}}&=0 \quad \quad \text{on} \quad \tilde {z}=\tilde {h}_0(\tilde {x}). \end{align}

Hence, we find that $\tilde {\tau }_0\equiv 0$ , $\tilde {u}_0\equiv -1$ , $\tilde {p}_0\equiv 0$ and $\tilde {h}_0(\tilde {x}) = \text{constant}$ . Thus, region 2 has a constant (a priori unknown) film thickness

(3.10) \begin{equation} \tilde {h}_0=\frac {\hat {h}_\infty }{\hat {H}\!{\textit{Ca}}^{2/3}}, \end{equation}

where $\hat {h}_\infty$ is the dimensional film thickness indicated in figure 1.

Next, we examine the transition region, which allows us to smoothly transition from the constant film region to the circular cap at the front of the bubble.

3.4. Region 3: transition region

In the transition region, we shift the origin to $x=-1$ and rescale the variables as

(3.11) \begin{align} \begin{aligned} x+1 &= {\textit{Ca}}^{1/3}X, & z&= {\textit{Ca}}^{2/3}Z,& h &= {\textit{Ca}}^{2/3}H, & p=P,\\ u&=U,& v&= {\textit{Ca}}^{1/3}V, &\tau &=G. \end{aligned} \end{align}

We again expand the variables in powers of ${\textit{Ca}}^{1/3}$ .

The leading-order equations of motion (2.8) become

(3.12a–c) \begin{align} U_{0X}+V_{0Z}&=0, & P_{0X}&=U_{0ZZ}, & P_{0Z}&=0 \quad \quad \text{in }\quad 0\lt Z\lt H_0(X). \end{align}

The boundary conditions (2.9a ) and (2.9c ) become

(3.13a) \begin{align} U_0=-1,\quad V_0&=0 \quad \quad\quad \quad\;\; \text{on} \quad Z=0, \end{align}

(3.13b) \begin{align} V_0 &=U_0 H_{0X} \quad \quad \text{on} \quad Z=H_0(X), \end{align}

and the remaining boundary conditions (2.13) become

(3.14a–c) \begin{align} G_0&=-U_{S0}^{\prime}, & -P_0&=H_0^{\prime\prime}, & U_{0Z}&= -{\textit{EG}}_0^{\prime} \quad \quad \text{on} \quad Z=H_0(X), \end{align}

where $U_{S0}(X) = U_0(X,H_0(X))$ .

Using (3.12)–(3.14), we find that $U_0$ is given by

(3.15) \begin{align} U_0(X,Z) =-1 +\frac {1}{2}H_0^{\prime\prime\prime}(X) (2H_0(X) Z-Z^2) - {\textit{EG}}_0^{\prime}(X) Z. \end{align}

By integrating (3.15) across the liquid layer (from $Z=0$ to $Z=H_0(X)$ ) we find that the flux of liquid in the $X$ -direction is given by

(3.16) \begin{align} Q= \frac {1}{3}H_0^{\prime\prime\prime} H_0^3 - \frac {1}{2}{\textit{EG}}_0^{\prime} H_0^2 - H_0. \end{align}

By conservation of mass this quantity must equal the flux of liquid in the thin film region, where $Q=-\tilde {h}_0$ . We thus obtain an equation for the height of the bubble’s surface in the transition region , namely

(3.17) \begin{align} H_0^{\prime\prime\prime} = \frac {3(H_0-\tilde {h}_0)}{H_0^3} +\frac {3{\textit{EG}}_0^{\prime}}{2H_0}. \end{align}

Next, using (3.15) and (3.17), we find that (3.14a ) becomes

(3.18) \begin{align} \frac {{d}}{\textrm{d}X}\left (\frac {1}{4}{\textit{EG}}_0^{\prime}H_0 +\frac {3\tilde {h}_0}{2H_0}\right )&=G_0. \end{align}

Equations (3.17) and (3.18) form a closed system for the film profile, $H_0(X)$ , and the perturbation to the surface concentration of surfactant, $G_0(X)$ , in the transition region. In addition, we enforce the matching conditions

(3.19a) \begin{align} H_0(X) &\rightarrow \tilde {h}_0, \quad \quad\quad\quad\;\;\;\, G_0(X) \rightarrow 0 \quad \quad \text{as}\quad X\rightarrow -\infty , \end{align}

(3.19b) \begin{align} H_0(X)&\sim \frac {1}{2}X^2 +\beta _1, \quad \quad G_0(X) \rightarrow 0 \quad \quad \text{as}\quad X\rightarrow \infty . \end{align}

In principle, the solution of the system (3.17)–(3.18) subject to the far-field behaviour (3.19) determines the a priori unknown constants $\tilde {h}_0$ and $\beta _1$ along with $H_0$ and $G_0$ .

In the next section, we analyse the problem (3.17)–(3.19).

4.1. Normalisation

We begin by normalising the equations (3.17) and (3.18) by scaling the variables as

(4.1) \begin{align} &\xi = \frac {X+\mathcal{S}}{\tilde {h}_0}, \quad \quad \eta (\xi ) = \frac {H_0(X)}{\tilde {h}_0}, \quad \quad g(\xi ) = \tilde {h}_0G_0(X). \end{align}

The (3.17) and (3.18) are translation invariant, so we introduce an arbitrary shift $\mathcal{S}$ to simplify the forthcoming analysis. Under these scalings, (3.17) and (3.18) become

(4.2a) \begin{align} \eta ''' &= \frac {3(\eta -1)}{\eta ^3} +\frac {3\mathcal{E}\!g'}{2\eta }, \end{align}

(4.2b) \begin{align} \frac {\mathcal{E}}{4} (\eta g')'&=g+\frac {3\eta '}{2\eta ^2}, \end{align}

where

(4.3) \begin{equation} \mathcal{E}=\frac {E}{\tilde {h}_0}. \end{equation}

The boundary conditions (3.19) imply that

(4.4a) \begin{align} \eta (\xi ) \rightarrow 1,\quad\quad &g(\xi ) \rightarrow 0 \quad\quad \text{as}\quad \xi \rightarrow -\infty , \end{align}

(4.4b) \begin{align} &g(\xi ) \rightarrow 0 \quad\quad \text{as}\quad \xi \rightarrow \infty . \end{align}

The solution of (4.2a ) can be shown to behave quadratically as $\xi \rightarrow \infty$ so

(4.5) \begin{align} \eta (\xi ) \sim \frac {1}{2}a \xi ^2 +b\xi +c \quad\quad \text{as}\quad \xi \rightarrow \infty , \end{align}

where $a$ , $b$ and $c$ are constants. Notice that the coefficients are not uniquely determined due to the arbitrary choice of origin for $\xi$ . However, the translation-invariant groups $a$ and $ac -( {1}/{2})b^2$ are uniquely determined. By comparison with (3.19), we see that they are related to the a priori unknown constants $\tilde {h}_0$ and $\beta _1$ by

(4.6) \begin{align} a&=\tilde {h}_0, \quad\quad ac-\frac {1}{2}b^2 = \beta _1. \end{align}

The solution strategy for the problem (4.2)–(4.4) is explained in the following subsection. Once $\eta$ and $g$ have been computed for a given value of $E$ , the surface velocity of the thin film is calculated from (3.15), giving

(4.7) \begin{align} U_{\!S} = \frac {1}{2}-\frac {\mathcal{E}}{4}\eta g' -\frac {3}{2\eta }. \end{align}

4.2. Solution

We solve (4.2) numerically by shooting from $\xi \rightarrow -\infty$ . Linearising (4.2) about the far-field behaviour (4.4a ), we find that

(4.8) \begin{align} \eta (\xi ) &\sim 1 + \sum _{n=1}^5 A_n \textrm{e}^{\lambda _n \xi },& g(\xi ) &\sim \sum _{n=1}^5 B_n \textrm{e}^{\lambda _n \xi } \quad \quad \text{as} \quad \xi \rightarrow -\infty , \end{align}

where the $\lambda _n$ are roots of the quintic polynomial

(4.9) \begin{align} \lambda ^5 - \frac {4}{\mathcal{E}}\lambda ^3 - 12\lambda ^2 +\frac {12}{\mathcal{E}}=0. \end{align}

This equation has two real and positive roots (which we label $\lambda _1$ and $\lambda _2$ ), one real and negative (labelled $\lambda _3$ ), and a complex conjugate pair with negative real part (labelled $\lambda _c$ and $\overline {\lambda }_c$ ). We require the solution to decay as $\xi \rightarrow -\infty$ , so only the positive eigenvalues are permitted. Hence, the decaying linearised far-field behaviour is given by

(4.10a) \begin{align} \eta (\xi ) &\sim 1 + A_1 \textrm{e}^{\lambda _1 \xi } + A_2 \textrm{e}^{\lambda _2 \xi } \quad \quad \text{as} \quad \xi \rightarrow -\infty , \end{align}

(4.10b) \begin{align} g(\xi ) &\sim \frac {\lambda _1\big (\lambda _1^3-12\big )A_1}{6}\,\textrm{e}^{\lambda _1\xi } +\frac {\lambda _2\big (\lambda _2^3-12\big )A_2}{6}\,\textrm{e}^{\lambda _2\xi } \quad \quad \text{as} \quad \xi \rightarrow -\infty , \end{align}

where $A_1$ and $A_2$ are a priori unknown constants. Due to the translation invariance we may (e.g.) set $A_1 = \pm 1$ by choice of $\mathcal{S}$ in (4.1). (Although translation allows us to set the coefficient of one exponential to have magnitude 1, we do not know its sign in advance; however, we always find that in the front transition region $A_1\gt 0$ .) We then determine $A_2$ via the shooting method to ensure our solution satisfies $g(\xi )\rightarrow 0$ as $\xi \rightarrow \infty$ .

For each value of $\mathcal{E}$ , we use the above shooting method to solve for $\eta$ and $g$ , then read off the coefficients $\{a,b,c\}$ in the quadratic behaviour (4.5) as $\xi \rightarrow \infty$ . We then use (4.3) and (4.6) to determine ${a=}\tilde {h}_0$ and $\beta _1$ parametrically as functions of $E$ . However, the shooting problem can become delicate for small or large values of $\mathcal{E}$ . In the next two subsections, we present asymptotic results for these two limits.

4.3. Small- $\mathcal{E}$ limit

In the limit where $\mathcal{E}$ is small we expand

(4.11) \begin{align} \eta &\sim \eta _{0}+\mathcal{E}\eta _{1}+{\cdots} , \quad \quad g\sim g_{0}+\mathcal{E}g_{1} +{\cdots} . \end{align}

This regime is similar to that studied by Waters & Grotberg (Reference Waters and Grotberg2002) for a surfactant-laden liquid plug and by Ginley & Radke (Reference Ginley and Radke1989) who considered a bubble in a capillary tube.

At $O(1)$ in (4.2), we find that

(4.12a) \begin{align} \eta _0^{\prime\prime\prime} &= \frac {3(\eta _{0}-1)}{\eta _{0}^3}, \end{align}

(4.12b) \begin{align} g_{0} &=-\frac {3\eta _{0}^{\prime}}{2\eta _{0}^2}. \end{align}

The decoupled (4.12a ) for $\eta _{0}$ is the same Landau–Levich equation used by Bretherton (Reference Bretherton1961) to determine the shape of a surfactant-free bubble in the transition region. The solution for $\eta _0$ is uniquely determined, up to an arbitrary translation, and the corresponding leading-order surfactant concentration profile $g_0(\xi )$ , given by (4.12b ), is plotted in figure 3. Although the limit $\mathcal{E}\rightarrow 0$ is singular, removing the highest derivative in (4.2b ), we see that $g_0$ tends to zero in the far field, as required, and no boundary-layer behaviour is produced, as also found by Ginley & Radke (Reference Ginley and Radke1989) and Waters & Grotberg (Reference Waters and Grotberg2002). The coefficients in the quadratic behaviour

(4.13a) \begin{align} \eta _{0}(\xi )&\sim \frac {1}{2}a_0 \xi ^2 +b_0 \xi +c_0 \quad \quad \text{as}\quad \xi \rightarrow \infty , \end{align}

are also determined uniquely and, in particular, we have $a_0\approx 1.337$ and $a_0c_0-({1}/{2})b_0^2\approx 3.88$ , as found by Bretherton (Reference Bretherton1961).

Figure 3.

The surfactant concentration in the front transition region, $g(\xi )$ , with $\mathcal{E}\rightarrow 0$ (black), $\mathcal{E}=1$ (blue), $\mathcal{E}=4$ (red) and $\mathcal{E}\rightarrow \infty$ (purple).

To find the correction to $\tilde {h}_0$ and $\beta _1$ due to the effect of surfactants, we proceed to first order in (4.2a ) to obtain the equation

(4.14) \begin{align} \eta _{1}^{\prime\prime\prime} = \frac {12(3-2\eta _{0})\eta _{1} + 9(2\eta _{0}^{\prime 2}-\eta _{0}\eta _0^{\prime\prime})}{4\eta _{0}^4}, \end{align}

for the correction to the thin film height. We solve (4.14) in the same fashion as (4.12a ) by shooting from $\xi \rightarrow -\infty$ . Again the solution to (4.14) behaves quadratically for large positive $\xi$

(4.15a) \begin{align} \eta _{01} &\sim \frac {1}{2}a_1 \xi ^2 +b_1 \xi +c_1 \quad \quad \text{as}\quad \xi \rightarrow \infty , \end{align}

where the constants $a_1,\,b_1,\,c_1$ are in principle determined (up to an arbitrary translation) by the solution of (4.14). In particular we find that $a_1\approx -0.0146$ and $a_0c_1+c_0a_1-b_0b_1\approx 0.58$ .

Finally, we obtain the small- $E$ expansions for ${a=}\tilde {h}_0$ and $\beta _1$ from (4.6). Note that the definition (4.3) of $\mathcal{E}$ involves $\tilde {h}_0$ , so we have to manipulate the expansions to remove the dependence on $\tilde {h}_0$ to get

(4.16a) \begin{align} {a}&\sim a_0+\frac {a_1}{a_0}\,E\approx 1.337-0.011E&&\text{as}\quad E\rightarrow 0, \end{align}

(4.16b) \begin{align} \beta _1&\sim a_0c_0-\frac {1}{2}b_0^2+\frac {a_0c_1+c_0a_1-b_0b_1}{a_0}\,E\approx 3.88+0.43E&&\text{as}\quad E\rightarrow 0. \end{align}

We note that (4.16a ) and (4.16b ) differ from (26) in Waters & Grotberg (Reference Waters and Grotberg2002) because they include a factor of 3 in their ${\textit{Ca}}$ , which induces a factor of $3^{1/3}$ in the definition of $E$ . We also note that Waters & Grotberg’s expression for the pressure drop is twice ours, due to the cylindrical instead of two-dimensional geometry that they study.

4.4. Large- $\mathcal{E}$ limit

In the other extreme where $\mathcal{E}$ is large, we expand

(4.17) \begin{align} \eta &\sim \eta _{0}+\frac {1}{\mathcal{E}}\,\eta _{1}+{\cdots} , \quad \quad g\sim \frac {1}{\mathcal{E}}\,g_{1} +{\cdots} . \end{align}

At $O(1)$ in (4.2) we find

(4.18a) \begin{align} \eta _0^{\prime\prime\prime} &= \frac {3(\eta _{0}-1)}{\eta _{0}^3}+\frac {3g_{1}^{\prime}}{2\eta _{0}}, \end{align}

(4.18b) \begin{align} \frac {1}{4}(g_{1}^{\prime}\eta _{0})' &= \frac {3\eta _{0}^{\prime}}{2\eta _{0}^2}. \end{align}

We can integrate (4.18b ) and substitute into (4.18a ) to obtain

(4.19a) \begin{align} \eta _0^{\prime\prime\prime}&=\frac {12(\eta _{0}-1)}{\eta _{0}^3}, \end{align}

(4.19b) \begin{align} g_{1}&=\frac {1}{2}\,\eta _{0}\eta _0^{\prime\prime}-\frac {1}{4} (\eta _{0}^{\prime})^2-\frac {1}{2}\,\beta _{10}, \end{align}

where $\beta _{10}$ is the leading-order approximation for the coefficient $\beta _1$ .

Once again, (4.19a ) is the Landau–Levich equation and it is similar to the surfactant-free equation (4.12a ) found by Bretherton (Reference Bretherton1961) except with an additional factor of 4 in the numerator. This additional factor of 4 induces an increase in the thin film height and correction to the pressure drop by a factor of $4^{2/3}$ , i.e.

(4.20a) \begin{align} {a} &\rightarrow 4^{2/3}\boldsymbol{\cdot }1.337\approx 3.369 \quad \quad\quad \;\;\, \text{as} \quad E\rightarrow \infty , \end{align}

(4.20b) \begin{align} \beta _1 &\rightarrow \beta _{10}\approx 4^{2/3}\boldsymbol{\cdot }3.88\approx 9.78 \quad \quad \text{as} \quad E\rightarrow \infty . \end{align}

In this limit, we find that the surface velocity (4.7) is given by $U_{\!S}\equiv -1$ , which corresponds to the bubble interface travelling at the same velocity as the walls of the cell.

These results reproduce the large-Marangoni-number limit reported in previous studies (see, for example, Ratulowski & Chang Reference Ratulowski and Chang1990; Park Reference Park1992; Stebe & Barthes-Biesel Reference Stebe and Barthes-Biesel1995; Shen et al. Reference Shen, Gleason, McKinley and Stone2002). However, we also evaluate the correction to the surfactant concentration, given by (4.19b ) and plotted in figure 5, where it is evident that $g_1$ does not satisfy the far-field condition $g_1(\xi )\rightarrow 0$ as $\xi \rightarrow -\infty$ . This apparent inconsistency can be resolved by examining an outer region in which

(4.21) \begin{align} \xi &=\mathcal{E}^{1/2}\varXi , & g(\xi )&=\mathcal{E}^{-1}\varPsi (\varXi ), & \eta (\xi )&\sim 1+\text{exponentially small terms}, \end{align}

so that (4.2b ) is transformed to

(4.22) \begin{equation} \varPsi ''=4\varPsi , \end{equation}

up to exponentially small corrections. By matching with (4.19b ) we thus obtain the leading-order outer solution

(4.23) \begin{equation} \varPsi _0(\varXi )=-\frac {1}{2}\,\beta _{10}\textrm{e}^{2\varXi }. \end{equation}

4.5. Results

In figure 3, we plot the correction from equilibrium to the surfactant concentration, $g$ , in the front transition region when $\mathcal{E}\rightarrow 0$ , $\mathcal{E}=1$ , $\mathcal{E}=4$ and $\mathcal{E}\rightarrow \infty$ . We use the arbitrary shift $\mathcal{S}$ introduced in § 4.1 to align the peaks of the concentration profiles. In the limit $\mathcal{E}\rightarrow \infty$ , $g$ vanishes across the entire domain, but in all other cases, we observe that $g\lt 0$ and so the surfactant concentration is everywhere below equilibrium in the front transition region. Similar concentration profiles were observed by Stebe & Barthes-Biesel (Reference Stebe and Barthes-Biesel1995) in a system with an elevated bulk concentration. In figure 4 we plot the film height in the transition region and, for all values of $\mathcal{E}$ , we observe similar profiles to those found by Bretherton (Reference Bretherton1961) and Park & Homsy (Reference Park and Homsy1984) for a surfactant-free bubble.

Figure 4.

The free-surface profile $\eta (\xi )$ in the front transition region, with $\mathcal{E}\rightarrow 0$ (black), $\mathcal{E}=1$ (blue), $\mathcal{E}=4$ (red).

Figure 5.

The leading-order (a) perturbation to the surfactant concentration, $g_1(\xi )$ , and (b) free-surface profile, $\eta _0(\xi )$ , in the front transition region in the limit $\mathcal{E}\rightarrow \infty$ .

We plot the lowest-order perturbation to the surfactant concentration, $g_1(\xi )$ , and the film height, $\eta _0(\xi )$ , in the limit $\mathcal{E}\rightarrow \infty$ in figure 5. The leading-order solution evidently does not satisfy the downstream boundary condition $g_1(\xi )\rightarrow 0$ as $\xi \rightarrow -\infty$ , implying that there must be a boundary layer at infinity, as explained in § 4.4. In figure 6 we plot the leading-order surface velocity, $U_{\!S}$ , for $\mathcal{E}\rightarrow 0$ , $\mathcal{E}=1$ , $\mathcal{E}=4$ and $\mathcal{E}\rightarrow \infty$ . We observe that, for finite $\mathcal{E}$ , there is a stagnation point (in the frame of the bubble) within the transition region. Its location is close to the minimum point of $g$ in figure 3, because the flow directed outwards advects surfactant away from the stagnation point. The presence of a stagnation point along the front of the bubble is a prevalent feature of gas bubbles in Hele-Shaw cells or capillary tubes, even in systems with more complicated surfactant dynamics and non-zero Reynolds numbers (Fujioka & Grotberg Reference Fujioka and Grotberg2005; Zheng, Fujioka & Grotberg Reference Zheng, Fujioka and Grotberg2007).

Figure 6.

The surface velocity, $U_S(\xi )$ , in the front transition region with $\mathcal{E}\rightarrow 0$ (black), $\mathcal{E}=1$ (blue), $\mathcal{E}=4$ (red) and $\mathcal{E}\rightarrow \infty$ (purple).

Figure 7.

The normalised thin film height $a$ versus elasticity parameter $E$ . The solid curve is from the numerical solution of (4.2), and dashed curves are the asymptotic predictions: (4.16a ) for small $E$ and (4.20a ) for large $E$ . (a) A log–linear plot to show the full range of $E$ . (b) The solution for $0\leqslant E\leqslant 5.77$ .

The normalised height of the thin film, $a$ , is plotted as a function of the elasticity parameter $E$ in figure 7. As $E\rightarrow 0$ , surfactant effects become negligible and the thin film height approaches Bretherton’s result of $1.337$ for a surfactant-free bubble (Bretherton Reference Bretherton1961). At the other extreme, when $E$ is large, $a$ approaches 3.369, which is larger by a factor of $4^{2/3}$ , as expected. Interestingly, (4.16a ) predicts a decrease in the thin film height for small $E$ (see figure 7 b), however, figure 7(b) shows that the asymptotic result (4.16a ) (also obtained by Waters & Grotberg Reference Waters and Grotberg2002) quickly becomes redundant, and the normalised thin film height, $a$ , increases with $E$ thereafter.

The correction to the pressure drop across the front meniscus, $\beta _1$ , is plotted as a function of $E$ in figure 8. Again when $E$ is small we recover the Bretherton (Reference Bretherton1961) result that $\beta _1\approx 3.88$ . We observe that $\beta _1$ is a monotonic increasing function and when $E$ is large $\beta _1$ approaches 9.78, in agreement with (4.20b ). In the numerical simulation of (4.2), accurate convergence for the value of $\beta _1$ could not be achieved for values of $E\lesssim 0.{5}$ due to the sensitivity of the numerical shooting method, caused by the singular nature of the system (4.2) as $E\rightarrow 0$ . In general, it is harder to compute the value of $\beta _1$ than $a$ at small $E$ because a significantly larger value of $\xi$ is needed to robustly extract the value of $ac-b^2/2$ from the quadratic function (4.5) than to determine $a$ . The numerical approach is thus useful provided $E\gtrsim O(1)$ , while the asymptotic approximation (4.16b ) is useful when $E$ is small, and we are reassured by figure 8 that there is at least a small overlap region where they approximately agree.

Figure 8.

The correction to the pressure drop across the front meniscus, $\beta _1$ , versus elasticity parameter $E$ . The solid curve is from the numerical solution of (4.2), and dashed curves are the asymptotic predictions: (4.16b ) for small $E$ and (4.20b ) for large $E$ .

5.1. Regions

As for the front meniscus, in the small- ${\textit{Ca}}$ limit the problem at the rear of the bubble splits into three regions of interest (see figure 2). In particular, for the rear-cap region, we can follow the same analysis as in § 3.2 and find that the leading-order shape and surfactant concentration are given by

(5.1a) \begin{align} h_0(x) &= 1-\sqrt {1-(l-1 +x)^2}, \end{align}

(5.1b) \begin{align} h_2(x) &= \beta _2\sqrt {(x+l)(2-l -x)}, \end{align}

(5.1c) \begin{align} \varGamma _0(x)&=1, \end{align}

where $l$ is the dimensionless length of the bubble, and $\beta _2$ is the a priori unknown $O({\textit{Ca}}^{2/3})$ correction to the pressure drop across the rear meniscus.

As in § 3.3, in the thin film region (region 2), the film height $\tilde {h}_0$ is constant, and now in principle known from the solution for the front meniscus (see figure 7). In the next section we analyse the equation in the rear transition region in a similar manner to § 4.

6.1. Normalisation

We again normalise by scaling the variables as

(6.1) \begin{align} &\xi = \frac {X+\mathcal{S}}{\tilde {h}_0}, \quad\quad \eta (\xi ) = \frac {H_0(X)}{\tilde {h}_0}, \quad \quad g(\xi ) = \tilde {h}_0 G_0(X), \end{align}

where $\mathcal{S}$ is an arbitrary shift of our coordinates. We obtain exactly the same (4.2) as the front transition region, i.e.

(6.2a) \begin{align} \eta '''&= \frac {3(\eta -1)}{\eta ^3} + \frac {3\mathcal{E}g'}{2\eta }, \end{align}

(6.2b) \begin{align} \frac {\mathcal{E}}{4}(g'\eta )'&=g+\frac {3\eta '}{2\eta ^2}. \end{align}

We solve (6.2) numerically now by shooting from $\xi \rightarrow \infty$ . We find that the decaying linearised solution is given by

(6.3a) \begin{align} \eta &\sim 1\pm \textrm{e}^{\lambda _3 \xi } + S \textrm{e}^{\lambda _R\xi }\cos (\lambda _I \xi -q), \end{align}

(6.3b) \begin{align} g&\sim \pm \frac {6 \lambda _3}{\mathcal{E} \lambda _3^2 -4}\textrm{e}^{\lambda _3 \xi } -S \! \varLambda _c(\lambda _R,\lambda _I) \textrm{e}^{\lambda _R\xi }\cos (\lambda _I \xi -q)-S\! \varLambda _s(\lambda _R,\lambda _I) \textrm{e}^{\lambda _R\xi }\sin (\lambda _I \xi -q), \end{align}

as $\xi \rightarrow \infty$ . Here, $S$ and $q$ are a priori unknown shooting parameters, $\lambda _3$ is the real negative solution of (4.9) and $\lambda _c= \lambda _R +\textrm{i}\lambda _I$ is the complex root with negative real part. The coefficients are given by

(6.4a) \begin{align} \varLambda _c(\lambda _R,\lambda _I) = \frac {6 \lambda _R \left(-4+\mathcal{E}\left(\lambda _R^2+\lambda _I^2\right)\right)}{16+8\mathcal{E}\left(\lambda _I^2-\lambda _R^2\right)+\mathcal{E}^2 \left(\lambda _R^2+\lambda _I^2\right)^2}, \end{align}

(6.4b) \begin{align} \varLambda _s(\lambda _R,\lambda _I) = \frac {6 \lambda _I \left(4+\mathcal{E}\left(\lambda _R^2+\lambda _I^2\right)\right)}{16+8\mathcal{E}\left(\lambda _I^2-\lambda _R^2\right)+\mathcal{E}^2 \left(\lambda _R^2+\lambda _I^2\right)^2}. \end{align}

Note again that the $\pm$ occurs in (6.3) because, although translation allows us to set the coefficient of the exponential to be of magnitude 1, we do not know its sign in advance. Finally, we now have that $\eta _0$ behaves quadratically for large negative $\xi$ , i.e.

(6.5) \begin{align} \eta &\sim \frac {1}{2}A \xi ^2 +B\xi +C \quad \quad \text{as}\quad \xi \rightarrow -\infty . \end{align}

We solve (6.2) for each value of $\mathcal{E}=E/\tilde {h}_0(E)$ by applying a shooting method with two unknown parameters $S$ and $q$ , which are fixed by ensuring $g(\xi ) \rightarrow 0$ and $\eta ''(\xi )\rightarrow \tilde {h}_0(E)$ as $\xi \rightarrow -\infty$ , where $\tilde {h}_0(E)$ is as shown in figure 7. The first condition corresponds to matching the surfactant concentration in the thin film to the equilibrium concentration in the rear cap (see § 5.1), and the second ensures that the thin film height at the rear meniscus matches the height of the thin film deposited at the front meniscus. Following the matching procedure laid out in § 4.1, we then obtain the $O({\textit{Ca}}^{2/3})$ correction to the pressure drop across the rear meniscus as

(6.6) \begin{align} \beta _2 = AC-\frac {1}{2}B^2. \end{align}

This two-parameter shooting problem can be extremely sensitive, so we examine the limiting cases using asymptotic analysis.

6.2. Small- $\mathcal{E}$ limit

In the extreme where $\mathcal{E}$ is small we expand

(6.7a) \begin{align} \eta &\sim \eta _{0}+\mathcal{E}\eta _{1}+{\cdots} , \end{align}

(6.7b) \begin{align} g&\sim g_{0}+\mathcal{E}g_{1} +{\cdots} . \end{align}

Then at $O(1)$ in (6.2) we again find that the equations reduce to (4.12). We note that (6.2b ) is singular in the limit $\mathcal{E}\rightarrow 0$ ; however, for the same reasons as presented in § 4.3 there is no boundary-layer behaviour and the solution of (4.12b ) satisfies all the relevant boundary conditions.

Again, the (4.12a ) for $\eta _{0}$ decouples and is just the usual Landau–Levich equation obtained for a surfactant-free bubble. At first order in (6.2) we again obtain (4.14) for the correction to the bubble’s surface. We solve (4.14) by shooting from $\xi \rightarrow \infty$ . Following the matching methodology laid out in § 4.1, we thus find that the $O({\textit{Ca}}^{2/3})$ correction to the pressure drop is given by

(6.8) \begin{align} \beta _2&\sim -1.13-0.73E \quad \quad \text{as}\quad E\rightarrow 0. \end{align}

The leading term in (6.8) is Bretherton’s classical result for the rear meniscus of a surfactant-free bubble (Bretherton Reference Bretherton1961), and the second term is the first correction due to the presence of surfactant.

6.3. Large- $E$ limit

In this limit we follow the same methodology as in § 4.4 to obtain

(6.9) \begin{align} \eta _0^{\prime\prime\prime}=\frac {12(\eta _{0}-1)}{\eta _{0}^3}. \end{align}

Matching with the rear-cap solutions (5.1) we find that the correction to the pressure drop is then given by

(6.10) \begin{align} \beta _2 \sim -4^{2/3}\boldsymbol{\cdot }1.13 \approx -2.85 \quad \quad \text{as}\quad E\rightarrow \infty . \end{align}

Again this is a factor of $4^{2/3}$ larger than the original Bretherton (Reference Bretherton1961) result. This extends the large-Marangoni-number limit reported in many studies (see, for example, Ratulowski & Chang Reference Ratulowski and Chang1990; Park Reference Park1992) to the rear meniscus.

6.4. Results

In figure 9 we show example solutions for $E=1.36$ ( $\mathcal{E}=1)$ and $E=5.76$ ( $\mathcal{E}=4)$ . In figure 9(b) we observe that the surfactant concentration can be both above and below the equilibrium concentration in the rear transition region, in contrast to the front transition region, where the concentration is always below equilibrium. For these specific solutions, we find that $\beta _2 \approx -1.50$ for $E=1$ and $\beta _2\approx -2.33$ for $E=4$ , which are greater in magnitude than the pressure drop $\beta _2\approx -1.13$ for a surfactant-free bubble found by Bretherton (Reference Bretherton1961). In figure 9(c) we plot the corresponding surface velocities, $U_{\!S}$ for the same values of $E$ . Similarly to the front meniscus, we observe a stagnation point (in the frame of the bubble) in the transition region. However, here the flow is directed into the stagnation point, resulting in a local increase in the surfactant concentration.

Figure 9.

(a) The surface profile, $\eta$ , (b) the surfactant concentration, $g_1$ , and (c) the surface velocity, $U_{\!S}$ , in the rear transition region with $\mathcal{E}=1$ , (black) $\mathcal{E}=4$ (red).

7.1. Force balance

We are now in a position to include the effect of surfactants in the models presented by Booth et al. (Reference Booth, Griffiths and Howell2023, Reference Booth, Griffiths and Howell2025a ,Reference Booth, Wu, Griffiths, Howell, Nunes and Stone b ) and Wu et al. (Reference Wu, Booth, Griffiths, Howell, Nunes and Stone2024) for the motion of an approximately circular bubble in a Hele-Shaw cell moving due to a uniform background flow $\hat {U}_{\!{f}}\boldsymbol{i}$ . Booth et al. (Reference Booth, Griffiths and Howell2023) find that the dimensionless velocity $\boldsymbol{U \!}_b$ of such a bubble is determined by the force balance

(7.1) \begin{equation} \frac {\boldsymbol{U \!}_b}{|\boldsymbol{U \!}_b|^{1/3}} =\frac {\delta }{\pi }\oint _{\partial \varOmega } -p\boldsymbol{n}\, \textrm{d}s, \end{equation}

where $\partial \varOmega$ is the bubble’s surface as viewed from above (see figure 10) and the Bretherton parameter is defined by

(7.2) \begin{equation} \delta = \frac {3\sqrt {\pi }\varGamma (11/6)}{(\beta _1 - \beta _2)\varGamma (4/3)}\, \frac {{\textit{Ca}}_{\!{f}}^{1/3}}{\epsilon }. \end{equation}

Here, $\epsilon = \hat {H}/\hat {R}$ , where $\hat {R}$ is the radius of the bubble (measured from above), and ${\textit{Ca}}_{\!{f}} = \hat {\mu }\hat {U_f}/\hat {\gamma }$ is the capillary number based on the background flow speed, $\hat {U}_{\!{f}}$ , both of which are assumed to be small. In the distinguished limit where ${\textit{Ca}}_{\!{f}} = O(\epsilon ^3)$ as $\epsilon \rightarrow 0$ , so the viscous lubrication pressure balances the pressure drop across the menisci, the bubble is circular to leading order (Booth et al. Reference Booth, Griffiths and Howell2023). For a surfactant-free bubble, $\beta _1$ and $\beta _2$ are given by the values $\beta _1(0)\approx 3.88$ and $\beta _2(0)\approx -1.13$ originally calculated by Bretherton. This result is now easily generalised for a surfactant-laden bubble by using the expressions for $\beta _1(E)$ and $\beta _2(E)$ found in §§ 4 and 6, respectively. Crucially, we recall that the elasticity parameter $E$ , given by (2.11), is independent of the capillary number, so the values of $\beta _1$ and $\beta _2$ depend only on the given surfactant properties and concentration. Note that the values of $\beta _{1}$ and $\beta _2$ are both $O(1)$ for the entire range of values of $E\in [0,\infty ).$

Figure 10.

Plan view of a surfactant-laden bubble in a Hele-Shaw cell in a uniform background flow.

From (7.1), we find that the velocity of an isolated bubble in a uniform background flow is given by $\boldsymbol{U \!}_b = U_{\!{b}}\boldsymbol{i}$ , where

(7.3) \begin{align} \frac {U_b^{2/3}}{2-U_b} = \delta = \left (\frac {\beta _1(0)-\beta _2(0)}{\beta _1(E)-\beta _2(E)}\right )\delta _B, \end{align}

and we define the surfactant-free Bretherton parameter

(7.4) \begin{equation} \delta _B = \frac {3\sqrt {\pi }\varGamma (11/6)}{(\beta _1(0)-\beta _2(0))\varGamma (4/3)}\, \frac {{\textit{Ca}}_{\!{f}}^{1/3}}{\epsilon } \approx 1.12 \frac {{\textit{Ca}}_{\!{f}}^{1/3}}{\epsilon }. \end{equation}

7.2. Results

In figure 11 we plot $U_b$ versus $\delta _B$ , for a range of values of $E$ . Note that, if we plotted versus the Bretherton parameter, $\delta$ given by (7.2), then all the curves would collapse. Plotting $U_b$ versus $\delta _B$ allows us to analyse the effect of surfactant on the bubble velocity in comparison with a surfactant-free bubble experiencing the same flow conditions. We observe that the velocity of a surfactant-laden bubble ( $E\gt 0)$ at each $\delta _B$ is less than that of a surfactant-free bubble ( $E=0$ ) at the same value of $\delta _B$ . This trend continues as we increase $E$ , up to the limiting case $E\rightarrow \infty$ when $\delta = 4^{2/3}\delta _B$ , the maximum value that $\delta$ can take for a fixed $\delta _B$ . Hence, we always find that a surfactant-laden bubble travels more slowly than a surfactant-free bubble under the same flow conditions.

The form of (7.3) has the same structure as the expression found by Baué et al. (Reference Baué, Gans, Jannin, Reichert, Cantat and Jullien2025) (their (5.11)) for the velocity of a droplet in a highly soluble surfactant solution, in the limit as the droplet viscosity tends to zero. Note that their capillary number is calculated from the droplet velocity, whereas we use the capillary number ${\textit{Ca}}_{\!{f}}$ based on the background flow speed. One can make an analogy between their constant $1/K$ and the prefactor of $\delta _B$ in (7.4). Our model provides the dependence of this prefactor on the surfactant properties, which is missing in their work. Baué et al. (Reference Baué, Gans, Jannin, Reichert, Cantat and Jullien2025) found experimentally that, with highly soluble surfactants, the velocity of a bubble increases with its size, which is consistent with (7.3).

Figure 11.

The dimensionless bubble velocity, $U_b$ (7.3) as a function of the surfactant-free Bretherton parameter, $\delta _B$ (7.4) for a range of values of $E=0$ (black), $E=1.36$ ( $\mathcal{E}=1)$ (blue), $E=5.76$ ( $\mathcal{E}=4)$ (red), $E=\infty$ (purple), with $\beta _1$ and $\beta _2$ given by (4.6) and (6.6), respectively.