2.1. The Taylor bubble

The Taylor bubble in general comprises a capsular structure with a round cap connected to a region of approximately uniform thickness in its interior (figure 1). The control of the film thickness in the uniform region, $h_T$ , has received significant attention since it was considered experimentally by Fairbrother & Stubbs (Reference Fairbrother and Stubbs1935) and Taylor (Reference Taylor1961), and theoretically by Bretherton (Reference Bretherton1961). Bretherton (Reference Bretherton1961) shows that the relative size of $h_T$ depends crucially on the magnitude of the capillary number

(2.4) \begin{equation} C = \frac {\mu q_B}{\gamma d}, \end{equation}

representing the ratio of the size of viscous stresses to the size of capillary stresses on the scale of the capillary width $d$ . For small $C$ , Bretherton (Reference Bretherton1961) demonstrates an asymptotic structure defined by a frontal cap, with a leading-order circular cross-section dominated by surface tension, that is matched to the long interior of the inviscid bubble, where the exterior film is stagnant, through a region in which lubrication theory is applied. Analysis of this structure in the context of steady travelling-wave states yields the analytical expression for the film thickness in the approximately horizontal region

(2.5) \begin{equation} h_T \sim 1.3375\, \left ( \frac {\mu q_B}{\gamma d} \right )^{2/3} d \equiv 1.3375 \, C^{2/3} d \qquad (C \to 0), \end{equation}

implying a control by a combination of viscous, capillary and geometric parameters.

For moderate to large $C$ , viscous stresses play an important role at the bubble cap, and a full-Stokes resolution is necessary there. In general, the interior thickness of the Taylor bubble satisfies

(2.6) \begin{equation} h_T = B(C) d, \end{equation}

where $B(C)$ is a dimensionless function of $C$ only (figure 2), with the property that $B(C) \sim 1.3375 \, C^{2/3}$ as $C \to 0$ in conformity with (2.5). Full-Stokes numerical solutions (Reinelt & Saffman Reference Reinelt and Saffman1985) for two-dimensional Taylor bubbles have determined $B(C)$ over a broad range of $C$ for $C \gtrsim 10^{-2}$ (circular markers). For sufficiently large $C$ , capillary stresses ultimately become negligible compared with viscous stresses, and the interior thickness saturates towards a factor multiple of the capillary width, $h_T \sim 0.45\, d$ . We note that the analytical function

(2.7) \begin{equation} B(C) \approx \frac {1.3375}{2.95 + C^{-2/3}} \end{equation}

provides a good representation (black, solid curve in figure 2) that captures both the small- $C$ limiting result of Bretherton (2.5) (red, dashed) and the moderate- to large- $C$ values determined numerically by Reinelt & Saffman (Reference Reinelt and Saffman1985).

Mass conservation in the travelling-wave state implies that the constant translation speed of the bubble cap is given by (Bretherton Reference Bretherton1961)

(2.8) \begin{equation} U = \frac {q_B}{(1 - B)d}. \end{equation}

With $t$ representing the time since the injection is initiated, the leading-order position of the bubble nose can be characterised by

(2.9) \begin{equation} x_F(t) \sim Ut = \frac {q_Bt}{(1-B)d}, \end{equation}

or simply $x_F \sim (q_B/d) t$ in the limit of $B \to 0$ arising for $C \to 0$ .

Figure 2.

The empirical expression for the universal function $B(C) \equiv h_T/d$ defining the size of the interior thickness of the Taylor bubble $h_T$ to the capillary half-width, defined by (2.6). The numerical results of Reinelt & Saffman (Reference Reinelt and Saffman1985) are shown as blue circular markers. The small- $C$ result of Bretherton (Reference Bretherton1961) is shown as a red dashed line. The fitted analytical function (2.7) is shown as a solid black curve.

2.2. The necking zone

We define the necking zone as lying between the injection nozzle and the interior of the Taylor bubble, $0 \leqslant x \lesssim x_N(t)$ (figure 1), where $x_N(t)$ is a time-dependent characteristic scale of the developing necking disturbance to be predicted. Within the necking zone, the viscous film squeezes the Taylor bubble, instigating pinch-off.

In modelling the necking interface, we apply lubrication theory, based formally on the requirement that the interfacial gradient in this region is small

(2.10) \begin{equation} h_x \ll 1. \end{equation}

Satisfaction of lubrication theory in the necking zone is not obvious a priori because the magnitude of the interfacial gradient $h_x$ is dependent on the relative longitudinal and transverse length scales of the resulting solution, as characterised by the aspect ratio $\sim h_m(t)/x_N(t)$ , where $h_m(t)$ is the maximal film thickness

(2.11) \begin{equation} h_m(t) = \max _{x} \left ( h(x,t) \right )\!. \end{equation}

We follow the asymptotic heuristic of adopting lubrication theory and assess the self-consistency of its predictions in maintaining $h_x \ll 1$ within the necking zone a posteriori.

Lubrication theory is used extensively in modelling capillary flows. Examples include liquid film breakup in capillary tubes driven by the Rayleigh–Plateau instability (e.g. Hammond Reference Hammond1983; Gauglitz & Radke Reference Gauglitz and Radke1988; Rykner et al. Reference Rykner, Saikali, Bruneton, Mathieu and Nikolayev2024), two-phase fluid displacement in capillaries (e.g. Bretherton Reference Bretherton1961; Zhao et al. Reference Zhao, Pahlavan, Cueto-Felgueroso and Juanes2018; Pahlavan et al. Reference Pahlavan, Stone, McKinley and Juanes2019; Lu, Li & Gao Reference Lu, Li and Gao2023), coating (e.g. Landau & Levich Reference Landau and Levich1942; Eggers Reference Eggers2004, Reference Eggers2005; Snoeijer et al. Reference Snoeijer, Andreotti, Delon and Fermigier2007; Gao et al. Reference Gao, Li, Feng, Ding and Lu2016) and capillary-driven droplets (e.g. Hocking Reference Hocking1983; King & Bowen Reference King and Bowen2001; Savva & Kalliadasis Reference Savva and Kalliadasis2009; Kiradjiev et al. Reference Kiradjiev, Breward and Griffiths2019). The governing equations of lubrication theory are

(2.12a–b) \begin{gather} \mu u_{yy} = p_x, \qquad p_y = 0, \end{gather}

where $u(x,y,t)$ is the longitudinal velocity of the viscous fluid, $p(x,y,t)$ is the pressure field of the viscous fluid and we use subscripts to denote partial derivatives. With $\kappa = h_{\textit{xx}}$ representing the linearised interfacial curvature, and $\gamma$ the interfacial coefficient of surface tension, the interface is subject to the following conditions on the capillary wall, $y=d$ , and bubble interface, $y = d-h(x,t)$ :

(2.13a–c) \begin{gather} u(x,d,t) = 0, \qquad u_y(x,d-h,t) = 0, \qquad [p]_{y =(d-h)_{-}}^{y = (d-h)_+} = -\gamma \kappa , \end{gather}

representing conditions of no slip on the sides of the capillary, no stress at the interface with the inviscid bubble and the jump in capillary stress across the interface, respectively. Integrating (2.12b ) subject to the jump condition (2.13c ), we obtain the pressure in the film

(2.14) \begin{equation} p(x,t) = p_0 - \gamma \kappa (x,t), \end{equation}

where $p_0$ is an arbitrary constant reference pressure; owing to the assumed incompressibility of the fluids, $p_0$ will have no effect on the dynamics of the problem. Integrating (2.12a ) twice and applying (2.13a , b ), we obtain the velocity profile

(2.15) \begin{equation} u=\frac {p_x}{2\mu }(y-d)(y-d+2h), \end{equation}

and hence, with (2.13c ) and $\kappa = h_{\textit{xx}}$ , the volume flux (per unit width)

(2.16) \begin{equation} q(x,t) = \int _{d-h}^{d} u \; \textrm {d} {y} = \frac {\gamma }{3\mu } h^3 h_{\textit{xxx}}. \end{equation}

Substituting the above into the continuity equation of the fluid film, $h_t = -q_x$ , we obtain the governing nonlinear hyperdiffusion equation

(2.17) \begin{equation} h_t = -\frac {\gamma }{3\mu } \left (h^3 {h}_{\textit{xxx}}\right )_x. \end{equation}

Conditions (2.1) provide the two boundary conditions at the input nozzle

(2.18) \begin{equation} h(0,t) = h_0, \qquad q(0,t) = q_F. \end{equation}

We couple the necking film to the interior thickness of the Taylor bubble (2.6) by applying the matching condition

(2.19) \begin{equation} \lim _{x \to \infty } h = h_T, \end{equation}

which creates a connection between the dynamics of the necking zone and the interior thickness of the film within the Taylor bubble $h_T$ . For conformity with (2.19), we apply the initial condition

(2.20) \begin{equation} h(x,0) = h_T. \end{equation}

Equations (2.17)–(2.20) form a closed system describing the growth of the necking film $h(x,t)$ . At the critical time $t_*$ defined by

(2.21) \begin{equation} h_m(t_*) \equiv \max _{x} \left ( h(x,t_*) \right ) = d, \end{equation}

the viscous fluid spans the width of the channel, defining the pinch-off time of the bubble. In the confined planar configuration considered in this paper, we will find that the fluid film grows locally in the necking region independently on either side of the channel centreline until a point where the two fluid films meet and the bubble thickness becomes zero in finite time $t_*$ . Our focus will be on understanding the parametric control of the pinch-off time $t_*$ .

2.3. Dimensionless system

In order to maximally reduce the parametric dependence of the model (2.17)–(2.20), we define the following non-dimensional variables based on intrinsic scales:

(2.22) \begin{equation} x = \left ( \frac {\gamma h_T^4}{ \mu q_F } \right )^{1/3} \hat {x} , \qquad t = \left ( \frac {\gamma h_T^7}{\mu q_F^4} \right )^{1/3} \hat {t}, \qquad h = h_T \hat {h}. \end{equation}

Upon dropping hats, (2.17) becomes

(2.23) \begin{equation} h_t = -\frac {1}{3}\left (h^3h_{\textit{xxx}}\right )_x, \end{equation}

with dimensionless flux $q=({1}/{3})h^3h_{\textit{xxx}}$ . Conditions (2.18) and (2.19) become

(2.24a-c) \begin{gather} h(0,t) = \frac {H}{B}, \qquad q(0,t) = 1, \qquad \lim _{x\to \infty }h=1, \end{gather}

where $H = h_0/d$ and $B = h_T/d$ . The initial condition (2.20) becomes

(2.25) \begin{equation} h(x,0) = 1, \end{equation}

and the pinch-off criterion (2.21) is

(2.26) \begin{equation} \max _{x} \left ( h(x,t_*) \right ) = \frac {1}{B}. \end{equation}

The non-dimensionalisation has reduced the dependence of the solutions to two dimensionless numbers

(2.27) \begin{equation} B = \frac {h_T}{d}, \qquad H = \frac {h_0}{d}, \end{equation}

representing the ratio of the interior film thickness $h_T$ to the half-width of the capillary $d$ , and the ratio of the film thickness at the inlet to the half-width of the capillary, respectively. The number $B$ is correspondent with the quantity $B(C)$ given by the one-to-one function of $C$ represented by (2.7), and is thus a surrogate for the capillary number. Capillary numbers in microfluidic systems, for example, are characteristically small (e.g. Anna Reference Anna2016) with $C \lesssim 10^{-2}$ , for which (2.7) gives $B \lesssim 0.1$ . The second dimensionless number $H$ sets the level of confinement of the nozzle relative to the width of the capillary. The number is restricted to $0 \lt H\lt 1$ , with $H \ll 1$ representing a strongly confining nozzle.

Figure 3.

Example snapshots of the solution to the necking-zone system (2.23)–(2.25) with parameters $H=0.05$ and $B=0.05$ are shown in (a–d). The corresponding maximum height of the fluid film $h_m(t)$ (e) grows until the bubble thickness reaches zero and pinch-off occurs at a finite time $t_*$ , as shown in panel (d). The evolutions of the maximal film thickness $h_m(t)$ are plotted for $H=0.05$ (black) and $H=0.4$ (blue), each with $B=0.05$ . The pinch-off time in the former is indicated by the asterisk marker. For the case $H=0.4$ , the film thickness at the nozzle $H/B$ is indicated by a horizontal dashed line, illustrating its initial correspondence with the maximal thickness in that case.

Figure 4.

Numerical solutions to the time-dependent necking-zone model (2.23)–(2.25) (black) for $(a)$ $B=0.2$ and $H=0.4$ at $t=4,12,20,28,t_*\approx 37.09$ , $(c)$ $B=0.05$ and $H=0.4$ at $t=50,250,450,650,850$ , $t_*\approx 978.35$ and $(e)$ $B=0.05$ and $H=0.05$ at $t=50,250,450,650,$ $850,t_*\approx 1002.3$ . The red dashed line represents the parabola (3.4) at final time $ t_*$ . The centreline of the channel, representing the thickness of the film at which pinch-off occurs, $h_m(t_*) = 1/B$ , is indicated by a horizontal dashed grey line in the (a,c,e). The panels on the right plot the corresponding maximal film thickness $h_m$ as a function of time $ t$ , showing approach to a quasi-static solution with front position given by (3.14) (red dashed). The asymptotic prediction for $h_m(t)$ for the small- $H$ limit (3.21) is overlayed as blue crosses in the case $H=0.05$ .

Figure 5.

The time taken to pinch-off $ t _*(H,B)$ as a function of $B$ for $H = 0.05$ (black) and $H =0.4$ (blue), as predicted by numerical solutions to (2.23)–(2.25). The dashed red line represents the small- $B$ analytical prediction of (3.24).

3.1. Quasi-static theory

We formulate a simplified analytical theory based on utilising an asymptotic framework of quasi-static interfacial evolution. A theoretical approach of this kind has been applied previously, in particular, in the context of droplet spreading (e.g. Hocking & Rivers Reference Hocking and Rivers1982; Hocking Reference Hocking1983; King & Bowen Reference King and Bowen2001; Savva & Kalliadasis Reference Savva and Kalliadasis2009, Reference Savva and Kalliadasis2011; Vellingiri, Savva & Kalliadasis Reference Vellingiri, Savva and Kalliadasis2011; Savva & Kalliadasis Reference Savva and Kalliadasis2013; Kiradjiev et al. Reference Kiradjiev, Breward and Griffiths2019). The theory, applicable for small capillary number (equivalently $B \ll 1$ ), is based on a separation of the flow into two asymptotic zones: an outer region, wherein a fast diffusive time scale maintains the interface close to the shape of a static meniscus (a parabola for two-dimensional droplets; Kiradjiev et al. Reference Kiradjiev, Breward and Griffiths2019); and an inner zone localised near a frontal apparent contact-line position, wherein the flow is matched to the precursor film. The outer region evolves quasi-statically in response to a slow time scale associated with the evolution of the contact line. The leading-order equation of the outer region is the Young–Laplace equation describing a static meniscus subject to the condition that the interface connects to the contact-line position. The inner zone instead forms a travelling-wave state (exhibiting the small spatial oscillation in the interface of the kind we see in figure 3 a), the solution of which yields a matching condition (the Cox–Voinov law; Voinov Reference Voinov1976; Cox Reference Cox1986) on the outer region.

In the present context, we can interpret the necking dynamics as a kind of forced droplet-like disturbance that spreads into an effective precursor film left by the extending Taylor bubble. As noted above, a quasi-static theory requires the necking disturbance to grow much thicker than the precursor film. In our context, we recall that the dimensionless precursor-film thickness is unity and the pinch-off criterion is $h_m(t_*) = 1/B$ . Therefore, if $B \ll 1$ , the necking disturbance will necessarily grow much thicker than the precursor film prior to pinch-off, with $1 \ll h_m(t) \leqslant 1/B$ . A quasi-static dynamics can therefore be anticipated to arise for $B \ll 1$ over the ‘large’ temporal range $1 \ll t \leqslant t_*(H,B)$ . Within this dominant time interval, the necking disturbance, in accordance with the proposed quasi-static theory, forms a two-zone structure: an outer quasi-static zone of approximately parabolic form; and an inner zone localised near an apparent contact line $ x _N( t )$ that advances in accordance with a Cox–Voinov law. We proceed to develop the quasi-static theory and compare its predictions with the full time-dependent theory.

We begin by solving for the outer quasi-static region. Within the quasi-static framework, the flow adjusts rapidly to a near-static state with $q \approx 0$ and hence, from the expression for flux $q = ({1}/{3}) h^3 h_{\textit{xxx}}$ given below (2.23),

(3.1) \begin{equation} h_{\textit{xxx}} \approx 0. \end{equation}

The above represents the equation describing the ultimate shape formed by relaxation under linearised surface tension, $\kappa _x =0$ ; in other words, it is the linearised Young–Laplace equation describing a static two-dimensional meniscus. On the scales of the quasi-static outer region ( $h \gg 1$ ), the precursor-film thickness is effectively vanishing to leading order,

(3.2) \begin{equation} h( x _N( t ), t ) = 0. \end{equation}

The condition at the input nozzle (2.24) gives the further boundary condition

(3.3) \begin{equation} h(0, t ) = H/B. \end{equation}

Integration of (3.1) subject to (3.2) and (3.3) gives us the parabola

(3.4) \begin{equation} h(x,t) = \left ( A( t ) x +\frac {H}{Bx_N(t)} \right ) (x_N(t) - x), \end{equation}

where $A( t )$ is a constant of integration related to the height of the parabola. Since the volume must equal $ t$ in accordance with the dimensionless input flux $q(0, t )=1$ , the interface must satisfy the volume constraint

(3.5) \begin{equation} \int _0^{ x _N} h \; \textrm {d} { x } = \frac {1}{6}A(t) x _N^3+\frac {H}{2B}x_N = t , \end{equation}

where we have substituted (3.4) and evaluated the integral. Hence,

(3.6) \begin{equation} A( t ) = \frac {3}{x_N(t)^{3} } \left ( 2 t - \frac {Hx_N( t )}{B} \right )\!. \end{equation}

The only remaining unknown in the quasi-static solution (3.4) is now the position of the contact line $ x _N(t)$ . Thus, we require a further condition for closure.

The additional condition is given by a Cox–Voinov relation that matches the quasi-static solution to the precursor film via an inner travelling-wave state in the vicinity of $x \approx x_N(t)$ . Within the inner region, the solution transitions from its approximately linear form predicted by the outer solution (3.4) as $ x \to x _N( t )$ to the downstream precursor-film thickness. A step-by-step derivation of the associated matching condition

(3.7) \begin{equation} h_x^3 = -3\dot {x}_N\ln \left(\dot {x}_Nx_N^3\right) \qquad \big(x = x_N(t) \big), \end{equation}

reviewing its development from (2.23), is provided in Appendix A. The result relates the rate of advancement of the contact line $\dot x_N(t)$ to the instantaneous interfacial gradient of the outer solution at the contact line $h_x(x_N,t)$ . The length scale of the inner zone is of order $\sim (\dot x_N)^{-1/3}$ , and hence it is necessary that $(\dot x_N)^{-1/3} \ll x_N$ in order for the required length-scale separation defining the quasi-static regime to apply. In other words, the argument of the natural logarithm in (3.7) is intrinsically large as part of the consistency of the theory. As $x_N(t)$ grows, the condition $(\dot x_N)^{-1/3} \ll x_N$ will become ever more strongly satisfied with time ( $t \gg 1$ ), concurrently with the thickness becoming much larger than that of the precursor film, $h_m(t) \gg 1$ . The control of $x_N(t)$ implied by the matching condition (3.7) is solely responsible for introducing dependences of the quasi-static evolution on the viscous and capillary parameters.

With the system now fully closed, we seek an explicit evolution equation for the contact-line position $x_N(t)$ . First, we isolate the rate of change $\dot x_N(t)$ in (3.7) by writing that equation in the equivalent form

(3.8) \begin{equation} x_N^3 \dot {x}_N = \exp \left (W\left (- \frac {1}{3} {x_N^3h_x^3}\right )\right ) \qquad (x = x_N(t)), \end{equation}

where $W(z)$ is the Lambert $W$ -function (which, for our purposes, is sufficient to define as the unique positive root of ${\textit{We}}^W = z$ ). The criterion for quasi-static scale separation $(\dot x_N)^{-1/3} \ll x_N$ noted below (3.7) requires both sides of the equation above to be large, and hence the argument of the $W$ function in (3.8) is intrinsically large within the theory. The large-argument form of the Lambert- $W$ function (Abramowitz & Stegun Reference Abramowitz and Stegun1965) has the first three terms

(3.9) \begin{equation} W(z) \sim \ln z - \ln (\ln z) +\frac {\ln \left (\ln z\right )}{\ln z} \qquad (z \to \infty ), \end{equation}

and hence the first two terms in its exponentiation are

(3.10) \begin{equation} \exp (W(z) ) \sim \frac {z}{\ln z} \left (1 + \frac {\ln (\ln z )}{\ln z} \right ) \qquad (z \to \infty ). \end{equation}

With use of the first term in the expansion above, we determine that (3.8) has the maximally reduced leading-order asymptotic form

(3.11) \begin{equation} \dot x_N \sim \frac { - h_x^3 }{ 3 \ln \left (- \frac {1}{3} x_N^3 h_x^3 \right ) } \qquad (x = x_N(t)), \end{equation}

giving the rate of change of the contact line in terms of the frontal interfacial slope of the quasi-static region $h_x(x_N,t)$ . The frontal slope of the quasi-static parabolic solution (3.4) can be evaluated as

(3.12) \begin{equation} h_x(x_N,t) = -\frac {2}{x_N}\left (\frac {3 t}{x_N} - \frac {H}{B}\right )\!. \end{equation}

Substituting (3.12) into (3.11), we obtain the ordinary differential equation

(3.13) \begin{equation} \dot {x}_N = \dfrac {\dfrac {8}{3}\left (\dfrac {3 t}{x_N}-\dfrac {H}{B}\right )^3}{x_N^3\ln \left (\dfrac {8}{3}\left (\dfrac {3 t}{x_N}-\dfrac {H}{B}\right )^3\right )}, \end{equation}

yielding the desired explicit evolution equation for $x_N(t)$ .

The ordinary differential equation (3.13) describes the contact-line evolution of the quasi-static regime, analogously to corresponding differential equations developed in the context of droplets, as considered in axisymmetric (Hocking Reference Hocking1982) and two-dimensional (King & Bowen Reference King and Bowen2001; Kiradjiev et al. Reference Kiradjiev, Breward and Griffiths2019) cases. A distinctive element in the present case is the pinning of one side of the parabolic interface by the injection nozzle, which introduces terms involving the parameter $H/B$ , and results in an advancing centre of the parabola (as opposed to symmetrical spreading about a line source). We have thus drawn a link between the asymptotic structures that can underlie injection-driven pinch-off in capillaries and those encountered in the analysis of droplet spreading.

While (3.13) is the most reduced asymptotic form of the quasi-static ordinary differential equation (to be used in the asymptotic analysis of the next sub-section), for the numerical evaluation of the quasi-static theory, we opt here to integrate the unsimplified form (3.8) owing to the large $\ln (\ln z)/\ln z$ residual in the expansion (3.10). We solve (3.8) with (3.12) numerically subject to an initial condition on $x_N(t_0)$ using Runge–Kutta–Fehlberg integration. In order for the argument to be positive, as is required of the Lambert $W$ -function here, we require the initial condition to satisfy $x_N(t_0) \lt x_N^{(0)}\equiv 3Bt_0/H$ . As is common in intermediate asymptotics (Barenblatt Reference Barenblatt1996), the quasi-static states have an attractor to which the solutions converge that is independent of the initial condition (demonstrated here in Appendix B). Thus, to give space for the quasi-static attractor to establish, we initialise with $x_N(t_0) = 0.8 \, x_N^{(0)}$ with $t_0=5$ , and show the solution for $t \geqslant 50$ , for which the solutions have closely approached the attractor.

With $x_N(t)$ determined in this way, the parabolic profile (3.4) can be evaluated to give a prediction for the evolution of the interface, with pinch-off occurring critically once $h_m(t_*) = 1/B$ , where

(3.14) \begin{equation} h_m(t) = \frac {\left (3t-\dfrac {Hx_N(t)}{B}\right )^2}{3x_N(t) \left (2t-\dfrac {Hx_N(t)}{B} \right )} \end{equation}

is the maximum of the parabola. The prediction of the quasi-static theory is shown for $H=0.05$ and $H=0.4$ in figure 4. The left-hand panels show the quasi-static prediction for the (parabolic) interface profile as a red dashed line at the time of pinch-off, showing good agreement with the numerical prediction of the time-dependent necking-zone model. The evolution of the maximal thickness (3.14) is overlaid as a red, dashed line in the right-hand panels ( $d$ ) and ( $f$ ) of figure 4 (representing the relevant cases of small $B$ ), showing excellent agreement with the numerical predictions of the full necking-zone theory (2.23)–(2.25) over the required time scale, $1 \ll t \leqslant t_*$ . The results confirm that the quasi-static regime is capturing the interface evolution for cases of $B \ll 1$ .

3.2. Asymptotic solution to the quasi-static theory for $H \ll 1$

For $H \ll 1$ , the spacing between the nozzle and the sidewall is small relative to the size of the capillary. In this limit, the thickness of the film near the input nozzle can, similarly to the contact-line position (3.2), be approximated as zero on the scales of the outer quasi-static region for all times up to pinch-off $(h(0,t) = H/B \ll 1/B)$ . As a result, the parabola is approximately symmetric (cf. figure 4 $e$ , where $H = 0.05$ ), with a line of symmetry residing close to $x_N(t)/2$ . The solution thereby retains a self-similar interfacial shape, a simplification which affords a yet further reduced analytical prediction within the quasi-static theory.

The implication of $H \ll 1$ within the quasi-static theory is that we can neglect the $H$ term in the governing differential equation (3.11), which reduces it to the parameterless equation

(3.15) \begin{equation} \dot {x}_N \sim \frac {72t^3}{x_N^6\ln \left({72t^3}/{x_N^3}\right)}. \end{equation}

We now seek a leading-order asymptotic solution to (3.15) during the quasi-static time-scale ( $1 \ll t \leqslant t_*$ ). To derive this, we try the ansatz of the form

(3.16) \begin{equation} x _N( t ) \sim f( t ) t ^{\alpha } \qquad (t \gg 1), \end{equation}

where $\alpha$ is an unknown constant to be determined and $f( t )$ is an unknown function with the property that $\ln (f( t )) \ll \ln t$ for $ t \gg 1$ (in other words, $f(t)$ is assumed to be at most of logarithmic order in $ t$ or a power thereof). Substitution of (3.16) into (3.15) yields, on neglect of higher-order terms in $t \gg 1$ ,

(3.17) \begin{equation} \frac {{d} }{\textrm {d} t } \left( x _N^7 \right) = \frac {504 t ^3}{ 3(1- \alpha ) \ln t }. \end{equation}

Integrating, and using a change of variable in the integral $z = 4 \ln t$ , we obtain

(3.18) \begin{equation} x _N^7= \frac {504}{3(1- \alpha ) } \, E(4 \ln t ) , \end{equation}

where $E(z) \equiv \int _\infty ^z z^{-1} e^{z} \; \textrm {d} z$ is the exponential integral function, and we have set the constant of integration to zero in order to satisfy the early-time asymptotic condition $ x _N \to 0$ as $ t \to 0$ . With use of the large-argument asymptote of the exponential integral function $E(z) \sim e^z/z$ as $z \to \infty$ (Abramowitz & Stegun Reference Abramowitz and Stegun1965), (3.18) reduces to

(3.19) \begin{equation} x _N^7 \, = \frac {42 t ^4}{(1-\alpha )\ln t } \, \end{equation}

in the relevant limit $ t \gg 1$ . Comparing the equation above with the ansatz (3.16), we see that $\alpha = 4/7$ is necessary for consistency, and $f( t ) = (98/\ln t )^{1/7}$ . Then, since $\ln (f( t )) = O( \ln (\ln t ))$ , the derived function $f(t)$ satisfies the required asymptotic property stipulated in our ansatz (3.16) that $\ln (f( t )) \ll \ln t$ , confirming the asymptotic consistency of the derivation. We conclude that the position of the contact line evolves as

(3.20) \begin{equation} x _N \sim \left ( \frac {98 t ^4 }{\ln t } \right )^{1/7} \qquad ( 1 \ll t \leqslant t_*). \end{equation}

The result of (3.20) indicates a dominant $x_N \sim t^{4/7}/(\ln t)^{1/7}$ growth of the front of the necking disturbance for small $B$ . The corresponding prediction of the maximum of the parabolic profile (3.4), occurring at $ x = x _N( t )/2$ with (3.20), is

(3.21) \begin{equation} h_m( t ) \sim \frac {3}{2} \left (\dfrac { t ^3 \ln t }{98} \right )^{1/7}. \end{equation}

The above is overlaid in figure 4( $f$ ) as blue crosses in the example with $B=0.05$ and $H=0.05$ , showing excellent agreement with both the numerical prediction of the full necking-zone theory (2.23)–(2.25), and the prediction of the quasi-static theory derived from solving the differential equation (3.8) numerically.

A $t^{4/7}$ power component was likewise obtained for the contact-line position in the late-time regime of injection-driven spreading of a two-dimensional droplet along a precursor film, based on treating the logarithmic dependences that appear in (3.15) as constant (Kiradjiev et al. Reference Kiradjiev, Breward and Griffiths2019). Our analysis indicates that, in the context of the pinch-off problem, the additional slowly varying $\ln t$ factor we derive in (3.20) is necessary for a consistent leading-order pinch-off solution that is sufficient to encompass the complete quasi-static time interval preceding pinch-off, $1 \ll t \leqslant t_*(B)$ .

To determine the pinch-off time, we substitute (3.21) into the pinch-off criterion $h_m( t _*) = 1/B$ , giving

(3.22) \begin{equation} t _*^3 \ln \big( t _*^3\big)= a B^{-7}, \end{equation}

where $ a = ({112}/{27} )^2$ is a numerical constant. Hence,

(3.23) \begin{equation} t _* = \exp \left (\frac {1}{3}W\big( a {B^{-7}}\big)\right ) . \end{equation}

Since $B \ll 1$ , we can simplify the above using (3.9) to give the final result

(3.24) \begin{equation} t _* \sim \left ( \frac { a B^{-7}}{\ln \big( a B^{-7}\big) } \right )^{1/3}\sim 1.350 \left ( \frac {B^{-7}}{\ln (1/B) } \right )^{1/3}, \end{equation}

providing an explicit asymptotic prediction for the pinch-off time $t_*$ for small capillary number. The prediction is overlaid as a dashed red curve in figure 5, showing agreement with the pinch-off time predicted by the numerical solution to the unsimplified necking-zone theory of (2.23)–(2.25) in the relevant limit of $B \to 0$ . Despite being a theory based on $H \ll 1$ (since, if $H = O(1)$ , the parabolic interface is asymmetric during the dominant interval $1 \ll t \leqslant t_*(B)$ and the contact line cannot conform to a simple asymptote of the form (3.16)), the prediction nonetheless appears to capture the $B \to 0$ trend for cases of both $H=0.05$ and $0.4$ . The result of (3.24) implies that $ t _* \gg 1$ as $B \to 0$ and thus self-consistently predicts the existence of the large asymptotic time interval wherein the quasi-static regime occurs ( $1 \ll t \leqslant t_*(B)$ ).

3.3. Conditions for consistency of the necking-zone model

The original theory of the necking zone (2.23)–(2.25) was based on two underlying asymptotic modelling assumptions: first, that lubrication theory applies to leading order in the necking zone; and second, that a long (Taylor) bubble forms. We now consider the parametric conditions under which these conditions hold self-consistently.

Beginning with the lubrication approximation (2.10), we note that, in our non-dimensional variables, the condition of small interfacial slopes can be expressed as

(3.25) \begin{equation} |h_x| \lesssim \varepsilon \left ( \frac {B}{QC} \right )^{1/3}, \end{equation}

where $\varepsilon \ll 1$ is a small dimensionless parameter characterising the size of the aspect ratio (representing the tolerance of the approximation), and $Q=q_F/q_B$ is the ratio of the film flux to the bubble flux. Thus, the aspect ratio of the necking disturbance at the time of pinch-off can be characterised by

(3.26) \begin{equation} h_x \sim \frac {h_m(t_*)}{x_N(t_*)} \sim 0.50 \, \left (B \ln \left (1/B \right ) \right )^{1/3}, \end{equation}

upon using $h_m(t_*) = 1/B$ and

(3.27) \begin{equation} x_N(t_*)\sim 2.0\,\left (\frac {1}{B^{4}\ln (1/B)}\right )^{1/3}, \end{equation}

obtained by substituting the pinch-off time (3.24) into (3.20). Equating (3.25) and (3.26), and simplifying, we obtain the constraint on $Q$ and $C$ given by

(3.28) \begin{equation} QC \ln (1/C) \lesssim 12 \, \varepsilon ^3. \end{equation}

If the criterion above holds then the model prediction self-consistently maintains lubrication theory. The product $QC \equiv \mu q_F / \gamma d$ can be interpreted as an alternative capillary number derived from the injection flux of the viscous phase (as opposed to the flux of the injected inviscid bubble phase used to define $C$ ). The result indicates that the film capillary number $QC$ is the primary control in determining the self-consistency of lubrication theory. This finding is consistent with the fact that, similarly to a droplet driven by injection (e.g. Kiradjiev et al. Reference Kiradjiev, Breward and Griffiths2019), larger injection fluxes feeding the necking disturbance will cause it to thicken relatively faster than it spreads longitudinally, yielding steeper aspect ratios. There is a weak logarithmic dependence on $C$ , stemming from its role in controlling the precursor-film thickness and, in turn, the spreading rate via the Cox–Voinov law. The constraint above restricts the region of validity of the model predictions in the parameter space $(C,Q)$ below the locus plotted as a dashed curve in figure 6 for the illustrative value $\varepsilon = 0.2$ . Lubrication theory is thus more strongly satisfied for smaller film capillary numbers $QC \ll 1$ (verification that lubrication theory arises is confirmed by direct comparisons with full-Stokes simulations in Appendix C).

Figure 6.

The $C$ – $Q$ parameter space partitioned by the characteristic regions in which the two asymptotic assumptions of the model (3.28) and (3.30) hold self-consistently (coloured shading). The criteria, for illustrative tolerances of $\varepsilon =\delta = 0.2$ , are shown as dashed curves, with the Taylor-bubble condition lying entirely inside the criterion for lubrication theory. Within the region of model self-consistency, the sub-region where quasi-static theory applies to good approximation is shown in green.

Second, we examine the self-consistency of the Taylor-bubble regime, defined by the property that the bubble cap extends characteristically further downstream than the scale of the necking disturbance (2.3), forming a long bubble,

(3.29) \begin{equation} x_N(t_*) \ll x_F(t_*). \end{equation}

Using the dimensionless form of the prediction for the characteristic length of the Taylor bubble given below (2.9), namely $x_F \sim Bt/Q$ , and (3.27), we find that, after some simplification, the ratio reduces to

(3.30) \begin{equation} \frac {x_N(t_*)}{x_F(t_*)} \sim 1.5\, Q \lesssim \delta , \end{equation}

where $\delta \ll 1$ is a parameter representing the tolerance. Thus, it is indicated that the characteristic length of the bubble to the length of the necking disturbance is controlled independently by the flux ratio $Q$ . For $\delta = 0.2$ , the above yields $Q \lesssim 0.1$ , indicating that well defined Taylor bubbles are formed if the flux of the viscous fluid is less than one tenth of the bubble flux. The condition is indicated by a horizontal dashed line on figure 6.

The region of $(C,Q)$ space in which the Taylor bubble assumption (3.30) applies is entirely contained within the region where lubrication theory (3.28) applies. Since the self-consistency of the model requires both of these conditions to be satisfied, we conclude that, for tolerances of $\varepsilon = \delta = 0.2$ , model self-consistency (represented by the coloured regions) occurs sufficiently on the basis of the Taylor-bubble condition alone. A decrease in tolerances moves the dashed lines representing validity of each condition closer together, with the conclusion that they remain separate holding for tolerances of $\varepsilon \sim \delta \sim 0.1$ . For yet smaller tolerances, some overlap between the constraints of (3.28) and (3.30) is possible, with the lubrication constraint introducing some restriction on the flux ratio at moderate capillary numbers $\gtrsim 10^{-1}$ . The subregion highlighted in green on the parameter space (figure 6) defines where the small- $B$ asymptotic limit of the pinch-off time (3.24) is within 20 % of the pinch-off time predicted by the full unsimplified necking-zone model (2.23)–(2.25) for $H=0.05$ . Thus, for characteristic values of $C \lesssim 0.04$ , the analytical quasi-static theory (§ 3.1) begins to apply to good approximation.

4.1. Theory limitations

The theory provided in this paper describes an asymptotically self-consistent lubrication regime. As shown by the parametric conditions for asymptotic consistency (§ 3.3), the analysis here is limited to predict the time up to pinch-off of a long (i.e. Taylor) and inviscid bubble. First, we note that lubrication theory breaks down immediately after pinch-off, owing to the development of infinitely steep interfaces on either side of the pinch-off point. The newly detached interfaces at both the rear of the newly formed bubble, and at the nose of the remnant of bubble fluid attached to the input nozzle, will both exhibit infinitely steep gradients. Hence, solutions (like the one shown in figure 3) are valid up to the point of pinch-off, but lubrication theory cannot necessarily capture the dynamics immediately after pinch-off. Full (nonlinear) curvature and full-Stokes flow may be required to resolve the newly formed cap and bubble rear. Moreover, if we were to consider larger flux ratios for which short or approximately circular-cross-sectioned bubbles are generated at flux ratios of order unity (situations above the horizontal dashed line shown in figure 6), similar full-Stokes considerations may be required. In this case, the flow would no longer form a Taylor bubble to which we can match the necking zone via a precursor-film condition of the form (2.19). The necking dynamics would instead interact more directly with the developing bubble, likely necessitating both nonlinear curvature and full-Stokes resolution.

We acknowledge that the analysis presented in this work is restricted to the planar coflow geometry. Extension of this analysis to the axisymmetric geometry would be of interest, and the demonstration of the quasi-static regime in the two-dimensional case presented here provides a first step towards this. In the axisymmetric case, the analogue of the parabolic quasi-static states would be surfaces of minimum curvature connecting two concentric annuli. Analysis of similar states in the static context, described by the axial Young–Laplace equations (e.g. Slobozhanin, Alexander & Fedoseyev Reference Slobozhanin, Alexander and Fedoseyev1999; De Gennes, Brochard-Wyart & Quéré Reference De Gennes, Brochard-Wyart and Quéré2003; Collicott, Lindsley & Frazer Reference Collicott, Lindsley and Frazer2006; Lv & Hardt Reference Lv and Hardt2021), shows that the static states exhibit a rich variety of forms (relative to the universally parabolic states arising here), as well as the possibility for minimising surfaces not to exist, or to bifurcate.

Finally, we note that the underlying asymptotic structure of the quasi-static theory assumes an inviscid bubble, as is typical in classical studies of Taylor bubbles (e.g. Fairbrother & Stubbs Reference Fairbrother and Stubbs1935; Bretherton Reference Bretherton1961). If the bubble viscosity is appreciable, a zero pressure gradient no longer applies in regions where the interface is flat, and hence the parallel flow solution is no longer applicable. Consequently, the asymptotic structure of a Taylor bubble, containing a broad interior region with a near-uniform film thickness, no longer exists. An interesting theoretical question would be to explore the implications of finite bubble viscosity.

4.2. Comparisons

The generation of Taylor bubbles has received significant experimental attention (e.g. Cubaud et al. Reference Cubaud, Tatineni, Zhong and Ho2005; Xiong et al. Reference Xiong, Bai and Chung2007; Fu et al. Reference Fu, Funfschilling, Ma and Li2010; Lu et al. Reference Lu, Fu, Zhu, Ma and Li2016; Li, Wu & Chen Reference Li, Wu and Chen2021; Huang & Yao Reference Huang and Yao2022). Many experimental investigations have aimed to correlate bubble characteristics (such as pinch-off time or, equivalently, bubble length or volume) and dimensionless parameters of the system, such as capillary number. Cubaud et al. (Reference Cubaud, Tatineni, Zhong and Ho2005) determined a linear relationship between the length of a bubble formed in a cross-flow junction of square cross-section, and the inverse of the liquid volumetric flux fraction. The scaling is also validated in the coflow geometry by Xiong et al. (Reference Xiong, Bai and Chung2007), who use an experimental particle image velocimetry system to provide a visual depiction of the interface with time, with the interface appearing approximately parabolic near to the nozzle (figure 7d of Xiong et al. Reference Xiong, Bai and Chung2007), exhibiting qualitative agreement with the predictions of the necking-zone theory in § 2.2.

Although there are no identical experimental or numerical studies on Taylor bubble formation in a planar coflow geometry, we can draw qualitative comparisons with existing empirical scaling laws. Assuming the Taylor bubble takes an approximately cylindrical shape, the scaling law proposed by Cubaud et al. (Reference Cubaud, Tatineni, Zhong and Ho2005) is analogous to a prediction of pinch-off time given by $t_* \propto 1/q_F$ . This inverse scaling for pinch-off time is also evidenced by the experimental investigation of Salman et al. (Reference Salman, Gavriilidis and Angeli2006) into the formation of Taylor bubbles by coaxial injection of air into a small cylindrical channel filled with water. The progression of the air–water interface is recorded with a high-speed camera and, from photographic observations, the frequency of bubble production is concluded to increase with the flux of either fluid phase. They also conclude that the period of bubble formation exhibits only a weak dependence on the surface tension of the outer liquid, and is determined predominantly by the nozzle size and flow rates. The prediction for pinch-off time of Taylor bubbles as $t_* \sim q_F^{-1}$ , derived implicitly from the scaling law for bubble length given by Cubaud et al. (Reference Cubaud, Tatineni, Zhong and Ho2005), is qualitatively consistent with the findings of the present paper that $t_* \sim q_F^{-4/3}$ , and we speculate that the different power is due to the differing geometries. As noted above in § 4.1, our analysis could be extended to address the more complex geometries associated with a cylindrical capillary or microchannel of square cross-section, which could in principle be addressed by appropriately adapting or generalising the necking theory and quasi-static frameworks. Similarly, we anticipate some of the ideas presented herein could potentially be applied to other related bubble systems including Hele-Shaw cells (e.g. Hazel & Heil Reference Hazel and Heil2002; Gaillard et al. Reference Gaillard, Keeler, Le Lay, Lemoult, Thompson, Hazel and Juel2021; Lawless et al. Reference Lawless, Keeler, Hazel and Juel2024).