Capillary breakup of a liquid bridge: identifying regimes and transitions

Abstract

Computations of the breakup of a liquid bridge are used to establish the limits of applicability of similarity solutions derived for different breakup regimes. These regimes are based on particular viscous–inertial balances, that is, different limits of the Ohnesorge number $Oh$. To accurately establish the transitions between regimes, the minimum bridge radius is resolved through four orders of magnitude using a purpose-built multiscale finite element method. This allows us to construct a quantitative phase diagram for the breakup phenomenon which includes the appearance of a recently discovered low-$Oh$ viscous regime. The method used to quantify the accuracy of the similarity solutions allows us to identify a number of previously unobserved features of the breakup, most notably an oscillatory convergence towards the viscous–inertial similarity solution. Finally, we discuss how the new findings open up a number of challenges for both theoretical and experimental analysis.

1 Introduction

The breakup of liquid volumes is a process that is ubiquitous throughout industry and nature. Recently, this process has attracted significant attention due to its importance for the functioning of a range of microfluidic technologies where one would like to be able to control the generation of uniform-sized droplets which then become building blocks, as in three-dimensional (3-D) printers (Derby 2010), or modes of transport for reagents, as in lab-on-a-chip devices (Stone, Stroock & Adjari 2004). In order to optimise the functioning of these technologies it is key to be able to understand the physical mechanisms controlling the breakup process and, ideally, to have a tool capable of predicting the response of the system to alterations in operating conditions.

From a theoretical perspective, the breakup phenomenon has attracted interest due to both its technological importance and its status as a free-surface flow which exhibits finite-time singularities (Eggers 1997). In the former case, motivated by a desire to describe engineering-scale systems, the focus has been on capturing the global dynamics of the process using multiphysics computational fluid dynamics (CFD) codes, with the actual breakup usually under-resolved due to the inherently multiscale nature of the problem. In contrast, in the latter case most interest has revolved around using asymptotic methods to study the micromechanics of the breakup process, with the larger-scale flow often not considered. Potentially, the advances made in understanding the micromechanics could be fed into large-scale CFD codes in order to exploit the advantages of both approaches, but this is yet to be achieved.

The failure of CFD packages to reliably capture breakup phenomena has been highlighted in Fawehinmi et al. (2005) for drop formation, where the appearance of satellite drops depended on grid resolution, and Hysing et al. (2009) for a two-dimensional (2-D) rising bubble in a regime where breakup was anticipated. In Hysing et al. (2009), six different codes were tested for the same problem and all gave different results ranging from no observed breakup through to multiple bubble detachments. In other words, the codes were unreliable for this process. The conclusion of the authors provides motivation for the current work: ‘Although the obtained benchmark quantities were in the same ranges, they did not agree on the point of break up or even what the bubble should look like afterwards, rendering these results rather inconclusive. To establish reference benchmark solutions including break up and coalescence will clearly require much more intensive efforts by the research community’.

Computational approaches cannot, in principle, resolve the thinning of a liquid volume down to the point at which breakup occurs, when the thread is infinitesimally small. Instead, at some point the pre-breakup process must be terminated, by ‘cutting’ the liquid thread joining two volumes, and a post-breakup state must be initiated based on the final solution in the pre-breakup phase. Whilst some computational codes do this process ‘automatically’, notably those using Eulerian meshes such as volume-of-fluid, there is no guarantee that this approach accurately represents the physical reality. Theoretical approaches to this problem have been considered in different regimes and are described in Eggers (2014). To implement an appropriate post-breakup solution, one must identify which pre-breakup regime the process is in, based on the scale for the cutoff. Therefore, knowing what regime a given breakup process is in and determining the transitions between these regimes are critical for the topological change to be correctly handled. The complexity of this procedure has been revealed in Castrejón-Pita et al. (2015); see § 2.2.3 for details, where multiple regimes were discovered for each breakup event.

Here, by computing the pre-breakup solution to the smallest possible scales that our simulation will allow, we will measure the accuracy of various similarity solutions proposed in the literature and, importantly, establish their limits of applicability. To do so, rather than just identifying power-law behaviour, methods will be developed to determine the regions of phase space where similarity solutions proposed for the different ‘regimes’ accurately approximate the full solution. The procedure is intended to leave no ambiguity as to when a certain regime is encountered. This attempt at a more rigorous approach will lead us to discover previously unobserved features of the breakup process that open up several new avenues of enquiry.

2 Regimes and transitions in the breakup phenomenon

This article focuses on the breakup of incompressible Newtonian liquids with constant viscosity ${\it\mu}$ and density ${\it\rho}$ , which have a constant surface tension ${\it\sigma}$ at their liquid–gas interface and are surrounded by a dynamically passive gas. The breakup is considered axisymmetric with the axis of the thread along the $z$ -axis of a cylindrical coordinate system ( $r,z$ ), as shown in figure 1 for a liquid bridge geometry, with characteristic length scale $R$ . In general, the flow is governed by the Navier–Stokes equations, but when inertial forces are weak the dominant balance of viscous and capillary forces gives a characteristic speed of breakup of $U={\it\sigma}/{\it\mu}$ so that the capillary number $Ca={\it\mu}U/{\it\sigma}=1$ . Then the appropriate dimensionless number characterising breakup is the Ohnesorge number $Oh$ , which is obtained by substituting $U$ into the Reynolds number $Re={\it\rho}UR/{\it\mu}={\it\rho}{\it\sigma}R/{\it\mu}^{2}=Oh^{-2}$ . The Ohnesorge number is the square root of the ratio of a viscous length scale $\ell _{{\it\mu}}={\it\mu}^{2}/({\it\rho}{\it\sigma})$ to the characteristic scale of the system $R$ . Then $Oh^{-1}=0$ is Stokes flow and $Oh=0$ gives inviscid flow.

Figure 1. Illustration of breakup in the liquid bridge geometry using the (dimensionless) coordinate system $(r,z)$ , with the black lines showing the computational domain. The minimum radius of the thread is $r_{min}$ , which can change vertical position as the breakup proceeds if satellite drops are formed (b), and the inward normal to the free surface is $\boldsymbol{n}$ .

The scales used so far are based only on global quantities. To analyse the breakup process locally it is useful to introduce a time-dependent local Reynolds number $Re_{local}(t)={\it\rho}U_{z}L_{z}/{\it\mu}$ based on instantaneous axial scales $U_{z}(t)$ and $L_{z}(t)$ (which will be larger than, or comparable to, radial scales). The practicalities of calculating this quantity are discussed in § 5.1. This local parameter will indicate the dominant forces during a breakup event and hence should identify the flow regime.

2.1 Regimes of breakup

Near breakup much research has focused on finding similarity solutions for the different regimes, all of which are categorised in § 4 of Eggers & Villermaux (2008). A summary of the key results is provided here, with particular attention paid to the dependence of the minimum thread radius $r_{min}$ , the maximum axial velocity $w_{max}$ and the local Reynolds number $Re_{local}$ on the time from breakup ${\it\tau}=t_{b}-t$ , where $t_{b}$ is the breakup time. Quantities have been made dimensionless with characteristic scales for lengths, velocities and time of $R$ , ${\it\sigma}/{\it\mu}$ and ${\it\mu}R/{\it\sigma}$ .

2.1.1 The inertial regime

In the ‘inertial regime’, henceforth the I-regime, where there is a balance between inertial and capillary forces, dimensional analysis (Keller & Miksis 1983; Brenner et al. 1997) of inviscid flow ( $Oh=0$ ) gives as ${\it\tau}\rightarrow 0$ that

(2.1a-c ) $$\begin{eqnarray}\displaystyle r_{min}=A_{I}(Oh\,{\it\tau})^{2/3},\quad Re_{local}\sim ({\it\tau}/Oh^{2})^{1/3},\quad w_{max}\sim (Oh^{2}/{\it\tau})^{1/3}, & & \displaystyle\end{eqnarray}$$

where $A_{I}$ is a constant of proportionality, given in Eggers & Villermaux (2008) as $A_{I}\approx 0.7$ . Notably, computations performed in the inertial regime have never listed $A_{I}$ and so it will be of interest to determine this value so that (2.1) is precisely determined and becomes a predictive tool. It will be useful to note, in terms of $r_{min}$ , that $Re_{local}\sim Oh^{-1}\,r_{min}^{1/2}$ and $w_{max}\sim Oh\,r_{min}^{-1/2}$ .

A feature of the inertial regime is the ‘overturning’ of the free surface, observed in both experiments (Chen, Notz & Basaran 2002) and simulations in this regime (Schulkes 1994), which invalidates attempts at a slender description of the thread. In Day, Hinch & Lister (1998) it was predicted that at breakup the free surface forms a double-cone shape with angles of $18.1^{\circ }$ and $112.8^{\circ }$ with the $z$ -axis. This phenomenon has been observed experimentally in Castrejón-Pita et al. (2012) and computationally in Wilkes, Phillips & Basaran (1999).

2.1.2 The viscous regime

The solutions derived for both the viscous and viscous–inertial regimes, henceforth referred to as the V-regime and the VI-regime, were derived in a one-dimensional approximation of the Navier–Stokes equations which relies on the thread remaining ‘slender’ as the breakup is approached. For a free surface represented by $r=r(z,t)$ this means that the gradient of the free surface must remain small, $\partial r/\partial z\ll 1$ .

For the V-regime, the similarity solution for Stokes flow ( $Oh^{-1}=0$ ) was derived in Papageorgiou (1995a , b ) and is given by

(2.2a-c ) $$\begin{eqnarray}\displaystyle r_{min}=0.0709{\it\tau},\quad Re_{local}\sim {\it\tau}^{2{\it\beta}-1}/Oh^{2},\quad w_{max}\sim {\it\tau}^{{\it\beta}-1}, & & \displaystyle\end{eqnarray}$$

where ${\it\beta}\approx 0.175$ . Recently, Eggers (2012) showed that this is the only similarity solution of infinitely many possible ones which remains stable to small perturbations. The free-surface profile associated with this similarity solution remains symmetric about the pinch point and the local profile depends on an external length scale (i.e. it is not ‘universal’), so that it is a self-similar problem of the second kind (Barenblatt 1996). The persistence of the V-regime for high-viscosity liquids has been observed experimentally in McKinley & Tripathi (2000).

2.1.3 The viscous–inertial regime

As ${\it\tau}\rightarrow 0$ , from (2.1) and (2.2) we find in the I-regime that $Re_{local}\rightarrow 0$ and in the V-regime that $Re_{local}\rightarrow \infty$ . This contradicts the initial assumptions behind their derivations (Lister & Stone 1998). Therefore, neither the I- nor the V-regime is valid right up to breakup, and a regime where viscous and inertial effects are in balance, so that $Re_{local}\sim 1$ , must be considered. In this VI-regime, the Navier–Stokes equations are required. It was shown in Eggers (1993) that for the VI-regime a ‘universal’ set of exponents exist for which

(2.3a-c ) $$\begin{eqnarray}\displaystyle r_{min}=0.0304{\it\tau},\quad Re_{local}\sim 1,\quad w_{max}=1.8Oh\,{\it\tau}^{-1/2}, & & \displaystyle\end{eqnarray}$$

with a highly asymmetric free-surface shape joining a thin thread to a much steeper ‘drop-like’ profile. It was later demonstrated in Brenner, Lister & Stone (1996) that this solution is the most favourable of a countably infinite set of similarity solutions. Notably $w_{max}=\max (w)$ is the maximum velocity out of the thin thread ( $w>0$ in what follows), rather than $\text{max}|w|$ , which scales in the same way but has a pre-factor $3.1$ (Eggers 1993). In terms of $r_{min}$ this gives $w_{max}=0.3Oh\,r_{min}^{-1/2}$ .

2.2 Transitions between regimes

It has been established that whichever regime the breakup process starts in, it will eventually end up in the VI-regime. However, it is possible that the transition from either the V- or the I-regime will occur at such small scales that this regime is irrelevant, at least from a practical perspective. For example, experiments in Burton, Rutledge & Taborek (2004) for millimetre-sized mercury drops show that the similarity solution in the I-regime adequately describes the breakup down to the nanoscale. Therefore, determining the transitions between the regimes is an important part of understanding the breakup process.

Phase diagrams which identify the regions of $(Oh,r_{min})$ space associated with different regimes have been constructed for the related problem of drop coalescence in Paulsen et al. (2012), Paulsen (2013) and Sprittles & Shikhmurzaev (2014b ) – though, intriguingly, these publications disagree on the number of regimes – whereas for the breakup phenomenon an equivalent diagram is currently lacking. Most progress in this direction has been made in Castrejón-Pita et al. (2015), where sketches of two characteristic phase space trajectories were presented (in their figure 1E). Here, we provide a systematic exploration of $(Oh,r_{min})$ space in order to construct a phase diagram for breakup which is equivalent to those obtained for coalescence. To do so, the scalings proposed for the transitions between the various regimes as $r_{min},{\it\tau}\rightarrow 0$ are considered.

2.2.1 Viscous to viscous–inertial transition (V $\rightarrow$ VI)

When $Oh\gg 1$ , the V-regime adequately describes the initial stages of breakup. However, the inertial term in the Navier–Stokes equations does not remain negligible as $r_{min}\rightarrow 0$ and so a transition to the VI-regime occurs at a bridge radius $r_{min}^{V\rightarrow VI}$ . In Basaran (2002) and Eggers (2005), by balancing the solutions in the V- and VI-regimes, i.e. by taking $Re_{local}\sim 1$ in (2.2), it was shown that this should occur when

(2.4) $$\begin{eqnarray}\displaystyle r_{min}^{V\rightarrow VI}\sim Oh^{2/(2{\it\beta}-1)}\sim Oh^{-3.1}, & & \displaystyle\end{eqnarray}$$

However, experimental evidence in Rothert, Richter & Rehberg (2003) appears to contradict this result, finding instead that $r_{min}^{V\rightarrow VI}$ is constant, i.e. independent of $Oh$ . Investigating these regimes computationally should provide new insight into the transition.

2.2.2 Inertial to viscous–inertial transition (I $\rightarrow$ VI)

When $Oh\ll 1$ the Euler equations accurately capture the initial stages of breakup until the local Reynolds number in (2.1) drops to $Re_{local}\sim 1$ , which occurs when ${\it\tau}\sim Oh^{2}$ so that

(2.5) $$\begin{eqnarray}\displaystyle r_{min}^{I\rightarrow VI}\sim Oh^{2}. & & \displaystyle\end{eqnarray}$$

This cross-over was first confirmed computationally in Notz, Chen & Basaran (2001) and experimentally in Chen et al. (2002) for the case of a water–glycerol mixture with $Oh=0.16$ , but has never been studied systematically across a range of $Oh$ .

2.2.3 The discovery of multiple regime transitions in Castrejón-Pita et al. (2015)

From the predicted transitions of (2.4) and (2.5) the phase diagram may be expected to look qualitatively like figure 2(a). In this figure, to see which regimes are encountered during breakup at a given $Oh$ , one follows a vertical line (at the given $Oh$ ) from the top axis downwards (through decreasing $r_{min}$ ), e.g. see paths 1 and 2. For small $Oh$ (path 1) one has an I-regime crossing into a VI-regime when $r_{min}\sim Oh^{2}$ , whilst for large $Oh$ (path 2) one has a V-regime crossing into the VI-regime at $r_{min}\sim Oh^{-3.1}$ .

However, a recent publication by Castrejón-Pita et al. (2015) shows that the picture can be more complex: the breakup can pass transiently through multiple different regimes. For example, Castrejón-Pita et al. (2015) show that, counter-intuitively, at $Oh=0.23$ there is no I $\rightarrow$ VI transition but rather an I $\rightarrow$ V $\rightarrow$ VI transition, so that a new unexpected low- $Oh$ V-regime is encountered. Similarly, a high- $Oh$ I-regime is observed. This behaviour is sketched out qualitatively in a phase diagram (their figure 1E) showing how $Re_{local}$ varies as breakup is approached for different values of $Oh$ .

In this work, we will use computations to build the first fully computed and quantitatively constructed phase diagram for the breakup so that the nature of the transitions can be established. To orientate the reader over the forthcoming sections, figure 2(b) gives a preview of the computed phase diagram, in a form that is only asymptotically accurate as $r_{min}\rightarrow 0$ but serves to neatly illustrate the computed flow transitions. The result shows clearly how the new phase diagram differs from the expected behaviour (figure 2 a) due to the appearance of the low- $Oh$ V-regime sandwiched between the I- and VI-regimes, so that path 1 now first encounters an I $\rightarrow$ V transition. In contrast to Castrejón-Pita et al. (2015), no evidence for the high- $Oh$ I-regime can be seen. These features will be considered in further detail over the forthcoming sections.

In §§ 7–9 we will show how the phase diagram in figure 2(b) was actually constructed, but before doing so we must specify the problem considered.

Figure 2. Phase diagrams of Ohnesorge number $Oh$ against minimum bridge radius $r_{min}$ : (a) expected phase diagram; (b) computed phase diagram. The paths of typical breakup events are given by arrowed lines 1 and 2, which show how different regimes are encountered as $r_{min}$ decreases. At higher $Oh$ (path 2) a V $\rightarrow$ VI transition is both expected and computed, whilst at lower $Oh$ (path 1), a single I $\rightarrow$ VI transition was expected but a more complex behaviour was computed, with a low- $Oh$ V-regime encountered before the VI-regime, so that I $\rightarrow$ V $\rightarrow$ VI transitions are found, as discovered in Castrejón-Pita et al. (2015). The dashed lines show the scaling of the transitions between the different regimes as $r_{min}\rightarrow 0$ , with I $\rightarrow$ V and V $\rightarrow$ VI transitions scaling as $r_{min}\sim Oh^{2}$ and V $\rightarrow$ VI occurring when $r_{min}\sim Oh^{-3.1}$ .

3 Problem formulation

A liquid bridge geometry (figure 1) is used with the liquid trapped between two stationary solid discs of (dimensional) radius $R$ a distance $2R$ apart and surrounded by a dynamically passive gas. This set-up allows us to isolate the breakup dynamics, limit the elongation of the domain (e.g. in contrast to dripping phenomena) and retain an experimentally realisable set-up. A future work will consider the effects of geometry through elongation by allowing the plates to move apart. The contact line where the liquid–gas free surface meets the solid remains pinned at the disc’s edge throughout. Assuming gravitational effects are negligible allows us to consider a plane of symmetry at $z=0$ of a cylindrical polar coordinate system $(r,z,{\it\theta})$ so that, using $R$ as a characteristic scale for lengths, the solid is located at dimensionless position $z=1$ and the free surface is pinned at $(r,z)=(1,1)$ . The extension to include gravity is not a difficult one, but it does not add significant value to a study focused on the small-scale breakup.

Using $U={\it\sigma}/{\it\mu}$ as a scale for velocity, $T=R/U$ as a scale for time and ${\it\mu}U/R$ for pressure, the (dimensionless) Navier–Stokes equations are

(3.1a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\quad \frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}=Oh^{2}\,\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}, & & \displaystyle\end{eqnarray}$$

where the stress tensor is

(3.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64B}=-p\unicode[STIX]{x1D644}+[\boldsymbol{{\rm\nabla}}\boldsymbol{u}+(\boldsymbol{{\rm\nabla}}\boldsymbol{u})^{\text{T}}]. & & \displaystyle\end{eqnarray}$$

Here, $\boldsymbol{u}$ and $p$ are, respectively, the velocity and pressure in the liquid, and the Ohnesorge number is $Oh={\it\mu}/\sqrt{{\it\rho}{\it\sigma}R}$ .

On the free surface, whose location $f(r,z,t)=0$ must be obtained as part of the solution, the kinematic equation

(3.3) $$\begin{eqnarray}\displaystyle \frac{\partial f}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}f=0 & & \displaystyle\end{eqnarray}$$

is applied alongside the usual balance of fluid stresses with capillarity in the directions tangential and normal to the free surface:

(3.4a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}\boldsymbol{\cdot }(\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})=\mathbf{0},\quad \boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}\boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{n}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{n}$ is the inward normal, $\unicode[STIX]{x1D644}$ is the metric tensor of the coordinate system and $p$ is taken relative to the (constant) gas pressure.

At the liquid–solid boundary, no-slip and impermeability conditions are applied, $\boldsymbol{u}=0$ and the free surface is pinned at the contact line $f(1,1,t)=0$ .

3.1 Initial conditions and initiation of breakup

The dynamics of breakup is studied by generating a thread shape that is taken to the edge of its Rayleigh–Plateau stability limit (Slobozhanin & Perales 1993; Lowry & Steen 1995), so that small perturbations will trigger breakup. The movie ‘FullBreakup’ in the supplementary material available at http://dx.doi.org/10.1017/jfm.2016.276 shows a typical breakup event resulting from this approach. In this way, the breakup dynamics can be computed independently of a driving force such as the flux from a nozzle, as in jetting, or the pull of gravity, as in dripping (Rubio-Rubio, Sevilla & Gordillo 2013). This reduces the dimensionality of the system to just one controlling parameter $Oh$ so that constructing a phase diagram becomes a tractable task.

The required initial thread shape can be generated by either solving the Young–Laplace equation for the free-surface position (Fordham 1948; Brown & Scriven 1980; Thoroddsen, Takehara & Etoh 2005) or using the finite element code, described below, in a quasi-static mode, to generate this shape by gradually removing fluid from the thread in a manner similar to the experimental set-up described in Meseguer & Sanz (1985). The latter method is used, as this also enables the shape on the stable–unstable boundary to be established. The result is the first profile in figure 3(a). The instability can be triggered either by continuing to gradually remove fluid from the thread or by very slowly moving the plates apart, and both approaches were shown to produce indistinguishable results, thus confirming that the only control parameter is $Oh$ .

4 Multiscale finite element computations

The mathematical model formulated in § 3 is for an unsteady free-surface flow influenced by the forces of inertia, viscosity and capillarity. To study this system over the entire range of $Oh$ requires computational methods.

4.1 Previous computational studies

The breakup of liquid volumes has been studied computationally using the Stokes (Gaudet, McKinley & Stone 1996; Pozrikidis 1999), Euler (Schulkes 1994; Day et al. 1998) and Navier–Stokes equations for both Newtonian (Ashgriz & Mashayek 1995; Popinet 2009) and non-Newtonian (Li & Fontelos 2003; Bhat et al. 2010) liquids. Considerable progress has been made by the group at Purdue University using finite element methods for axisymmetric flows such as dripping (Ambravaneswaran, Phillips & Basaran 2000), jetting (Ambravaneswaran et al. 2004), drop-on-demand (Chen & Basaran 2002) and liquid bridge breakup (Suryo & Basaran 2006) with an array of different physical effects such as non-Newtonian fluids (Yildirim & Basaran 2001), surfactant dynamics (McGough & Basaran 2006) and electric field effects (Collins et al. 2008). These algorithms were the first to demonstrate a number of experimental features such as overturning of the free surface for a viscous fluid (Wilkes et al. 1999), transitions between scaling regimes (Chen et al. 2002) and the identification of multiple regime transitions (Castrejón-Pita et al. 2015).

Notably, it has been shown that with sufficient care, the results of sharp interface methods, such as those implemented at Purdue, can sometimes be recovered by the diffuse interface approach, in which topological changes are easier to handle computationally. In particular, the diffuse interface method developed in Yue et al. (2004) was shown in Zhou, Yue & Feng (2006) to accurately recover the sharp interface results of Wilkes et al. (1999). Furthermore, in Zhou et al. (2006), the flexibility of the diffuse interface method was demonstrated by computing compound drop formation, where additional complexity arises from the need to track two free surfaces. Despite these advances, we will implement a sharp interface approach (see § 4.2) whose reliability and accuracy has been repeatedly confirmed.

Despite the numerous successes of computational schemes in accurately predicting many of the global features of breakup, there have been fewer investigations of the local dynamics and its comparison with the similarity solutions. One of the most impressive works in this direction is Suryo & Basaran (2006), where scales comparable to those resolved in the current work ( $r_{min}<10^{-4}$ ) were captured and the scaling behaviour of both Newtonian and power-law liquids was investigated. Other progress includes results in Notz et al. (2001) and Chen et al. (2002), where I $\rightarrow$ VI transitions were shown, and Castrejón-Pita et al. (2015), in which multiple transitions were observed (see § 2.2.3). Our work will build on these previous investigations by systematically quantifying the accuracy of the similarity solutions across the entire parameter space (of Newtonian liquids) in order to build the first quantitatively constructed phase diagram for breakup.

In principle, computations are the ideal tool with which to map a phase diagram for the breakup process and to investigate in which regions each of the similarity solutions is accurate; a similar procedure was used for the coalescence phenomenon in Sprittles & Shikhmurzaev (2014b ). The difficulty of this procedure is that many decades of $r_{min}$ are required in order to reliably determine the scaling behaviour of different regime boundaries (if they exist). Although one may hope to compute $r_{min}$ until its behaviour falls into the VI-regime, in practice this may occur at scales below the possible computational resolution, e.g. $r_{min}<10^{-5}$ (Burton et al. 2004). Therefore, establishing the range of applicability of the similarity solutions using computations and then using these solutions to carry the dynamics to scales below realisable computational resolution seems to be a promising strategy for capturing breakup.

4.2 Computational approach

A focus of this work is to resolve the spatial and temporal dynamics of the pinch-off process to the smallest scales possible in order to compare with the similarity solutions proposed in the literature for the final stages of breakup. The approach used is based on the finite element framework originally developed in Sprittles & Shikhmurzaev (2012b , 2013) to capture dynamic wetting problems and subsequently used to study the coalescence of liquid drops (Sprittles & Shikhmurzaev 2012a , 2014a ), including two-phase calculations (Sprittles & Shikhmurzaev 2014b ); drop impact phenomena (Sprittles & Shikhmurzaev 2012a ); the detachment of bubbles from an orifice (Simmons, Sprittles & Shikhmurzaev 2015); and dynamic wetting in a Knudsen gas (Sprittles 2015). These flow configurations are all inherently multiscale, either due to the disparity of length scales in the problem formulation or because of the dynamics of the process itself, which generates small scales during its evolution, as is the case for breakup phenomena. The framework allows the global flow in these problems to be resolved alongside localised smaller scales, without relying on particular limits of $Oh$ .

As a step-by-step user-friendly guide to the implementation of this computational framework has been provided in Sprittles & Shikhmurzaev (2012b ), only the main details will be recapitulated here alongside some aspects which are specific to the current work. The code uses the arbitrary Lagrangian Eulerian scheme, based on the method of spines (Ruschak 1980; Kistler & Scriven 1983), to capture the evolution of the free surface in 2-D or 3-D axisymmetric flows. In order to keep the problem computationally tractable, the mesh is graded so that small elements can be used near the pinch-off region whilst larger elements are used where scales associated with the global flow are present.

The mesh starts with ${\approx}4000$ triangular elements ( ${\approx}500$ surface nodes) and then adaptively refines as the breakup progresses by adding spines whenever elements become too deformed. At the end of the computation, the number of elements is in the range $7000$ – $50\,000$ ( $800$ – $7000$ surface nodes), depending on the particular breakup requirements. For example, for a computation at $Oh=10^{-3}$ , the mesh is refined due to the formation of a ‘corner’ in the free surface (see figure 15), resulting in 869 surface nodes with a minimum surface element length of $3\times 10^{-5}$ . In contrast, at $Oh=0.16$ , refinement is required due to the formation of a thin thread of liquid (figure 9) which requires 5263 surface nodes with surface elements in this region all having length ${\approx}10^{-4}$ . Overturning of the free surface is permitted by using angled spines (i.e. not only horizontal) whose slope varies along the thread.

Refinement of the mesh is restricted by controlling the smallest permitted element size $h_{min}$ , which affects the final number of elements. Reducing $h_{min}$ allows the code to converge to smaller $r_{min}$ so that more of the breakup is recovered. By reducing $h_{min}$ from $10^{-2}$ down to $10^{-5}$ , convergence of the scheme under spatial refinement has been established, with each decrease in $h_{min}$ revealing more of the solution (smaller $r_{min}$ ) but producing curves which are graphically indistinguishable from those obtained on cruder meshes, where both solutions exist.

The result of the spatial discretisation is a system of nonlinear differential algebraic equations of index two (Lötstedt & Petzold 1986), which are solved using the second-order backward differentiation formula, whose application to the Navier–Stokes equations is described in detail in Gresho & Sani (1999). The time step automatically adapts during a simulation to capture the appropriate temporal scale at each instant.

In the problem considered, the smallest scale that needs resolving is the minimum neck radius $r_{min}$ , and eventually this becomes so small that accurate converged solutions can no longer be obtained. For the entire range of $Oh$ considered, computations remain accurate down to $r_{min}=10^{-4}$ and in certain cases they can go as far as $r_{min}=10^{-5}$ .

5 Quantitatively identifying regimes

Computations are unable to resolve down to $r_{min}=0$ , i.e. the point at which the topological change occurs, so that the precise time $t=t_{b}$ of pinch-off is unknown. This can be problematic when comparing to scaling laws of the form $r_{min}=A{\it\tau}^{B}+C$ , where ${\it\tau}=t_{b}-t$ requires the time of breakup, as small changes in $t_{b}$ can have a large effect on a curve of $r_{min}$ versus ${\it\tau}$ , which is then used to determine which regime a breakup process is in. Similar conclusions were reached for coalescence in Thoroddsen et al. (2005). In the final regime one has $C=0$ , otherwise for ‘transitional’ or ‘transient’ regimes $C$ must be determined as well. In general, there is no easy method to overcome these issues without assuming a particular functional form for the breakup.

However, in the V- and VI-regimes it is possible to exploit $r_{min}$ being a linear function of ${\it\tau}$ ( $B=1$ ) so that the speed at which the minimum thread radius evolves is constant, i.e. independent of both $t_{b}$ and $C$ . Therefore, to identify these regimes the ‘speed of breakup’ given by ${\dot{r}}_{min}=\text{d}r_{min}/\text{d}t$ can be compared to the predictions of (2.2) that ${\dot{r}}_{min}=-0.071$ and (2.3) that ${\dot{r}}_{min}=-0.030$ as a function of $r_{min}$ without needing to know $t_{b}$ . This has the further advantage that $C$ does not have to be fitted when a regime is observed transiently, as in Castrejón-Pita et al. (2015).

The same approach cannot be used to identify the inertial regime where $B=2/3$ , so that the speed of breakup is no longer constant. Therefore, we try a similar idea to that above and determine an $Oh$ -independent quantity which should be linear in the inertial regime. This relies on $C\approx 0$ , so that this regime must be close to breakup. Then, given that $r_{min}=A_{I}(Oh\,{\it\tau})^{2/3}$ , we introduce $l_{min}=Oh^{-1}r_{min}^{3/2}$ , which will satisfy $l_{min}=A_{I}^{3/2}{\it\tau}$ if the expected scaling holds. Differentiating this quantity with respect to time gives $\dot{l}_{min}=\text{d}l_{min}/\text{d}t=-A_{I}^{3/2}$ , i.e. a time-independent constant.

A key aspect of this work is to develop techniques to quantify when the breakup dynamics is in a certain regime (Sprittles & Shikhmurzaev 2014b ). To do so, we must define what it means to be ‘in a regime’. We choose to define a breakup process to be in the V- or VI-regime when the speed of breakup ${\dot{r}}_{min}$ is within $0.015$ of that predicted by (2.2) or (2.3), respectively. In other words, speeds in the range ${\dot{r}}_{min}=-0.071\pm 0.015$ indicate V-regime dynamics whilst ${\dot{r}}_{min}=-0.030\pm 0.015$ for the VI-regime. In a similar way, the inertial regime is defined by points at which $\dot{l}_{min}$ is within 0.15 of $A_{I}^{3/2}$ , where we will later see that $A_{I}^{3/2}=0.5$ . Choosing smaller (larger) margins to define each regime will shrink (expand) the area of a given regime in the phase diagram but should not alter the qualitative picture. The particular values chosen avoided regimes overlapping and enabled us to satisfactorily map the phase diagram. The method for identifying regimes is illustrated in the Appendix.

Rather than assuming that the regimes fill all of phase space, we will be careful to split the entrance to and exit from a particular regime. For example, the exit from the V-regime will be labelled as $r_{min}^{V\rightarrow }$ and the entrance to the VI-regime as $r_{min}^{\rightarrow VI}$ , rather than assuming a priori that these coincide at $r_{min}^{V\rightarrow VI}$ .

5.1 Analysis of the regimes

To understand the dynamics in the different regimes and determine when the assumptions behind the similarity solutions of § 2.1 are valid, local to the breakup point we analyse the assumptions made about (i) the relative magnitudes of viscosity and inertia and (ii) the slenderness of the thread. In practice, this means calculating (i) the local Reynolds number $Re_{local}$ and (ii) a measure for the slenderness of the free-surface profile.

To measure these characteristics, one needs to define an appropriate axial length scale $L_{z}$ and velocity scale $U_{z}$ near the pinch point. The obvious choice for $U_{z}$ is the maximum axial velocity in the thread, $U_{z}=w_{max}$ , which will occur near, but not at, the position where $r=r_{min}$ . Therefore, a characteristic axial scale can be defined as the vertical distance between the positions of minimum radius and maximum velocity, so that $L_{z}=|z(r=r_{min})-z(w=w_{max})|$ . This defines (i) the local Reynolds number as $Re_{local}=Oh^{-2}\,U_{z}L_{z}$ and (ii) the slenderness as the ratio of the characteristic radial length scale $L_{r}=r_{min}$ to $L_{z}$ , so that the geometry is slender when $L_{r}/L_{z}\ll 1$ . Importantly, $w_{max}=\max (w)$ is used rather than $\max |w|$ (see figure 9 c) to avoid discontinuities in $L_{z}$ , and hence $Re_{local}$ , which occur when the position of $\max |w|$ jumps from being above the pinch point to below it (as in the V $\rightarrow$ VI transition).

Other possible definitions of $U_{z}$ and $L_{z}$ exist (see Castrejón-Pita et al. 2015), so it is not a priori clear we have made the correct choices. However, computations will show that our definitions allow us to reproduce the expected scalings for $Re_{local}$ in both the I-regime (2.1) and the V-regime (2.2), which serves to validate our approach.

6 Overview

In §§ 7–9, respectively, the dynamics of breakup for large ( $Oh>1$ ), intermediate ( $0.1<Oh<1$ ) and small ( $Oh<0.1$ ) Ohnesorge numbers are analysed separately and used to establish the limits of applicability of the different similarity solutions described in § 2.1. This will allow us to define the regimes of breakup on a phase diagram.

In each section, the initial focus will be on a breakup event which highlights the features of the prominent regime, namely the V-regime ( $Oh=10$ ), VI-regime ( $Oh=0.16$ ) and I-regime ( $Oh=10^{-3}$ ), with a movie of each breakup in the supplementary material. In each case, snapshots of the free-surface shape, axial velocity and pressure distributions in the breakup region will be presented. Evidence that the breakup is in a given regime will be provided by comparing the computational results with the similarity solution’s predictions for the evolution of the free-surface shape, minimum bridge radius $r_{min}$ , maximum axial velocity $w_{max}$ and scaling of the local Reynolds number $Re_{local}$ .

Figure 3. For $Oh=10$ , characteristic of the V-regime, plots show (a) evolution of the entire free surface, (b) a close-up of the breakup region, (c) the axial velocity at the free surface, $w_{fs}$ , and (d) the pressure at the free surface, $p_{fs}$ , at $r_{min}=10^{-1}$ (1, blue), $r_{min}=10^{-2}$ (2, red), $r_{min}=10^{-3}$ (3, green) and $r_{min}=10^{-4}$ (4, brown).

Having confirmed the existence of a regime, we then identify its boundaries on the phase diagram to determine the regime transitions. As explained in § 5, these boundaries are specified by calculating the values of ${\dot{r}}_{min}$ and $\dot{l}_{min}$ as a function of $r_{min}$ . Having discovered the transitional behaviour, the local measures $Re_{local}$ and slenderness $L_{r}/L_{z}$ (introduced in § 5.1) are studied in an attempt to rationalise these findings.

The analysis in §§ 7–9 will not only allow us to construct a phase diagram (cf. figure 18) for the process but also open up new questions about the breakup process; to avoid distracting from the main focus in §§ 7–9, these aspects will be considered in more detail in § 11.

7 Viscous-dominated flow (large $Oh$ )

For $Oh=10$ free-surface profiles in figure 3(a,b) show the development of a thin thread of liquid whose smallest minimum radius remains at $z=0$ throughout, so that no satellite drops form. The axial velocity (figure 3 c) at the free surface $w_{fs}$ , which is approximately $r$ -independent, shows a rapid drainage from the thread driven by a pressure gradient from $z=0$ (figure 3 d).

Evidence that this breakup is in the V-regime can be seen in figure 4(b), where the free-surface profile obtained at $r_{min}=10^{-4}$ is almost indistinguishable from the similarity solution for the V-regime, and figure 5(a), in which the computed axial velocity (curve 3) follows the predicted (singular) form $w_{max}\sim r_{min}^{-0.825}$ as $r_{min}\rightarrow 0$ .

Figure 4. A comparison of the computed free-surface profiles with the similarity solutions for the (a) VI-regime and (b) V-regime. Computed profiles are at the minimum radius $r_{min}=10^{-4}$ for (a) $Oh=0.16$ and (b) $Oh=10$ , typical of breakup in the VI- and V-regimes, respectively. These profiles compare well to the similarity solutions (dashed lines) in (a) with ( ${\it\xi}_{VI},{\it\phi}_{VI}$ ) provided in Eggers (1993), giving $(r,z)=({\it\tau}{\it\phi}_{VI},z_{0}+Oh\,{\it\tau}^{1/2}{\it\xi}_{VI})$ with no free parameters, and in (b) with ( ${\it\xi}_{V},{\it\phi}_{V}$ ) from Papageorgiou (1995b ) and Eggers (1997), so that $(r,z)=({\it\tau}{\it\phi}_{V},Oh^{2-2{\it\beta}}\,{\it\tau}^{{\it\beta}}\ell _{z}{\it\xi}_{V})$ with $\ell _{z}=1.85\times 10^{-3}$ . In (a), the height of the pinch point $z_{0}=0.268$ for the VI-regime is calculated using the expression in Eggers (1993) for the drift of this point $z(r=r_{min})=z_{0}-1.6Oh\,{\it\tau}^{-1/2}$ (shown in figure 11 as the dashed line).

Figure 5. Curves are for 1: $Oh=10^{-3}$ , 2: $Oh=0.16$ and 3: $Oh=10$ . (a) Evolution of the maximum axial velocity $w_{max}$ as the minimum radius $r_{min}$ decreases, showing excellent agreement with the similarity solutions (dashed lines) from the I-regime, $w_{max}\sim r_{min}^{-0.5}$ , the VI-regime, $w_{max}=0.3Oh\,r_{min}^{-0.5}$ , and the V-regime, $w_{max}\sim r_{min}^{-0.825}$ . (b) The minimum bridge radius $r_{min}$ against time from breakup ${\it\tau}$ . Computations converge to the similarity solutions for the I-, VI- and V-regimes, respectively, with the lower dashed line being the similarity solution for the I-regime ((2.1) with $A_{I}=0.63$ ), the middle one for the VI-regime (2.3) and the highest one for the V-regime (2.2).

Assuming that the breakup remains in the V-regime, the data can be extrapolated to obtain the breakup time $t_{b}$ and, as a result, plot $r_{min}$ against ${\it\tau}$ in figure 5(b) (curve 3), where there is excellent agreement with the similarity solution of (2.2) for $r_{min}<10^{-2}$ .

7.1 Identification of the V-regime

When $Oh\gg 1$ the breakup is expected to start in the V-regime before transitioning into the VI-regime. Figure 6(a) shows that when $Oh=10$ (curve 5), the breakup speed ${\dot{r}}_{min}$ reaches the value predicted in the V-regime, $-0.071$ , at approximately $r_{min}=10^{-2}$ and remains there until $r_{min}=10^{-4}$ , i.e. until the end of the computation. For higher $Oh$ , curves are graphically indistinguishable from the case of $Oh=10$ , so this is the Stokes flow solution. Therefore, for $r_{min}\geqslant 10^{-4}$ the VI-regime is not observed at large $Oh$ .

Figure 6. (a) Evolution of the breakup speed ${\dot{r}}_{min}$ for 1: $Oh=1$ , 2: $Oh=1.5$ , 3: $Oh=2.5$ , 4: $Oh=5$ and 5: $Oh=10$ . The similarity solution in the V-regime ( ${\dot{r}}_{min}=-0.071$ ) is shown as a dashed line. At higher $Oh$ (e.g. curve 5) the speed remains at $-0.071$ , whilst at lower values (e.g. curve 1) there is a transition towards the VI-regime where ${\dot{r}}_{min}=-0.030$ . (b) Variation of the local Reynolds number for 1: $Oh=0.25$ , 2:  $Oh=0.5$ , 3: $Oh=1$ , 4: $Oh=2.5$ and 5: $Oh=10$ . Curves show that in the V-regime $Re_{local}$ gradually increases until a critical value is achieved, at which point the V-regime is exited. The increase of $Re_{local}$ in the V-regime is shown to follow the predicted scaling from (2.2) of $Re_{local}\sim r_{min}^{-0.65}$ (lower dashed line). The transition out of the V-regime ( $r_{min}=r_{min}^{V\rightarrow }$ marked as circles) is found to occur when $Re_{local}\approx 0.85$ (horizontal dashed line).

Figure 7 shows the region of $Oh$ – $r_{min}$ phase space where V-regime dynamics is encountered, with the circles marking the computationally determined boundaries of this regime. The procedure for identifying this regime is shown in the Appendix and, as discussed in § 5, it is based on the speed of breakup ${\dot{r}}_{min}$ remaining between $-0.071\pm 0.015$ . Next, the changes in behaviour triggering the entrance and exit from this V-regime are considered.

7.1.1 Entrance into the V-regime

The V-regime is entered ( $r_{min}^{\rightarrow V}$ ) at a constant $Oh$ -independent value $r_{min}\approx 10^{-2}$ , which suggests that this transition may occur when the thread can be considered slender, so that the assumptions behind the similarity solution (2.2) become valid. To determine whether the geometry near the breakup point is ‘slender’ or not, we plot $L_{r}/L_{z}$ in figure 8. For $Oh=10$ (curve 5), and generally for larger $Oh$ , initially $L_{r}/L_{z}\sim 1$ (no scale separation), but as $r_{min}$ shrinks the breakup region becomes slender. Defining the region to be slender when $L_{r}/L_{z}<0.1$ , it is found that for all $Oh>1$ slenderness occurs at $r_{min}\approx 10^{-2}$ , in agreement with our result for $r_{min}^{\rightarrow V}$ . Therefore, the transition into the V-regime appears to be dictated by the geometry of the thread.

Figure 7. Phase diagram showing the existence of a V-regime where there is a balance between viscous and capillary forces. The transition $r_{min}^{\rightarrow V}$ is where the breakup enters this regime and $r_{min}^{V\rightarrow }$ is where it leaves. Computational results (circles) show that the entrance is at a constant $Oh$ -independent $r_{min}^{\rightarrow V}=0.014$ (horizontal dashed line) whilst the exit follows $r_{min}^{V\rightarrow }=5.5\times 10^{-4}Oh^{-3.1}$ (lower dashed line).

Figure 8. Evolution of the slenderness of the free surface in the breakup region for 1: $Oh=10^{-3}$ , 2: $Oh=10^{-2}$ , 3: $Oh=0.14$ , 4: $Oh=0.16$ , 5: $Oh=10$ . This quantity, which must be small for the thread to be slender (as is required by the V- and VI-regimes), is defined using $L_{r}=r_{min}$ and $L_{z}=|z(r=r_{min})-z(w=w_{max})|$ .

7.1.2 Exit from the V-regime

The V-regime is left when the speed of breakup ${\dot{r}}_{min}$ diverges from $-0.071\pm 0.015$ (e.g. curve 1 in figure 6 a). Figure 7 shows that this transition follows the proposed scaling for V $\rightarrow$ VI that $r_{min}=A\,Oh^{-3.1}$ , with computations finding $A=5.5\times 10^{-4}$ .

The transition out of the V-regime ( $r_{min}^{V\rightarrow }$ ) occurs when inertial effects become non-negligible, so that (2.2) is no longer valid. Figure 6(b) shows that the increase in $Re_{local}$ follows remarkably well the scaling predicted by the similarity solution for the V-regime that $Re_{local}\sim r_{min}^{-0.65}$ . This also serves to validate our definition of $Re_{local}$ . Circles placed at $r_{min}^{V\rightarrow }$ on curves in figure 6(b) show that the transition out of the V-regime occurs when $Re_{local}\approx 0.85$ . One may expect that at this point, where $Re_{local}\sim 1$ , the VI-regime will be encountered, and this will be the focus of the next section.

8 Viscous–inertial balance (intermediate $Oh$ )

At $Oh=0.16$ in figure 9(a,b) one can see the appearance of a satellite drop centred at $z=0$ . The close-up images show how this is driven by symmetry breaking about the pinch point after $r_{min}\approx 10^{-2}$ (curve 2), as observed experimentally in Rothert, Richter & Rehberg (2001). Consequently, two pinch points occur at a finite distance either side of the symmetry plane ( $z=0$ ) with a large pressure acting to push fluid out of the thinnest region (figure 9 d). In contrast to profiles in the I-regime (cf. § 9), the geometry near the pinch point remains slender; see curve 4 in figure 9(b) and the slenderness data in figure 8 (curve 4).

Features that indicate the breakup is occurring in the VI-regime are as follows. Figure 4(a) shows that the free-surface profile at $r_{min}=10^{-4}$ compares well with the similarity solution (dashed line) derived in Eggers (1993) with no free parameters. Below the pinch point ( $z<0.25$ ) in figure 4, deviations of the computed free-surface shape from the similarity solution occur as there is no ‘far field’ in our set-up (as the plane of symmetry is at $z=0$ ). This also affects the velocity profiles (see figure 9 c), which agree with the similarity profiles in Eggers (1993) above the pinch point but not below – a point that will be discussed further in § 11.2.1. However, curve 2 in figure 5(a) shows that the maximum axial velocity follows $w_{max}=0.3Oh\,r_{min}^{-0.5}$ from (2.3), again with no adjustable parameters.

Figure 9. For $Oh=0.16$ plots show (a) the entire free surface, (b) a close-up of the breakup region, (c) the axial velocity at the free surface $w_{fs}$ (with $w_{max}$ when $r_{min}=10^{-3},~10^{-4}$ shown) and (d) the pressure at the free surface $p_{fs}$ at 1: $r_{min}=10^{-1}$ (blue), 2: $r_{min}=10^{-2}$ (red), 3: $r_{min}=10^{-3}$ (green) and 4: $r_{min}=10^{-4}$ (brown).

In figure 10(a), the breakup speed is shown for $Oh>0.15$ , values whose significance will become apparent. In this range, all curves tend towards the value of $-0.030$ , indicating that breakup enters the VI-regime.

Figure 10. (a) Evolution of the speed of breakup when $Oh>Oh_{c}$ . Curves are for 1: $Oh=0.16$ , 2: $Oh=0.25$ , 3: $Oh=0.5$ . At $Oh>Oh_{c}$ curves converge towards the similarity solution in the VI-regime given by (2.3), i.e.  ${\dot{r}}_{min}=-0.030$ . (b) Different dynamics are observed for the breakup speed in the VI- and V-regimes with 1: $Oh=0.16$ , 2: $Oh=10$ . For these values of $Oh$ , the computed solutions converge in an oscillatory (1) and a monotonic (2) manner towards similarity solutions in the VI- and V-regimes, given by (2.3) and (2.2), respectively, which are the lower and upper dashed lines.

Figure 11. Evolution of the vertical position of the pinch point, i.e. minimum thread radius $z(r=r_{min})$ , indicating the formation of satellite drops once the pinch point moves away from the centre line $z(r=r_{min})>0$ . Curves are for 1: $Oh=10^{-3}$ , 2: $Oh=0.16$ , 3: $Oh=0.5$ and 4: $Oh=1$ . The dashed line shows the expression from Eggers (1993) at $Oh=0.16$ for the drift of the pinch point, $z(r=r_{min})=z_{0}-1.6Oh\,{\it\tau}^{-1/2}$ with $z_{0}=0.268$ .

What is immediately striking about the curves in figure 10(a) is that they appear to oscillate around the value of  $-0.030$ , converging towards it as $r_{min}\rightarrow 0$ . To highlight this behaviour, figure 10(b) compares curves for $Oh=10$ , which remains in the V-regime, and $Oh=0.16$ , where the VI-regime is most prominent. Whilst the value of  $-0.071$ is approached monotonically from below by curve 2, the approach of curve 1 towards $-0.030$ is entirely different. The motion in the VI-regime is clearly oscillatory in the radial direction, and on closer inspection of figure 5(a) and figure 11 oscillations in axial quantities $w_{max}$ (and hence $Re_{local}$ in figure 6 b) and $z(r=r_{min})$ can also be seen. Having identified this behaviour, it can also be recognised in the plot of $r_{min}$ against ${\it\tau}$ in curve 2 of figure 5(b), with the bridge radius oscillating around and converging towards the dashed line predicted by the similarity solution (2.3). As far as we are aware, this is the first time this oscillatory behaviour has been observed and it will be discussed further in § 11.2.

Figure 12(a) focuses on the narrow range $Oh=0.14-0.16$ to highlight a sharp transition in the dynamics of breakup which occurs at a critical $Oh_{c}\approx 0.15$ . For $Oh=0.16>Oh_{c}$ (curve 2), the transition into the VI-regime occurs at $r_{min}\approx 10^{-2}$ , whilst for $Oh=0.14<Oh_{c}$ (curve 1) the transition into this regime is delayed until $r_{min}\approx 4\times 10^{-4}$ . Remarkably, we will show that this is due to the existence of a ‘low- $Oh$ V-regime’, first discovered in Castrejón-Pita et al. (2015), which precedes entry into the VI-regime for $Oh<Oh_{c}$ . In figure 12(a) for $Oh=0.14$ (curve 1), it can be seen that this recently discovered regime is entered when $r_{min}^{\rightarrow V}=3\times 10^{-3}$ and exited again when $r_{min}^{V\rightarrow }=6\times 10^{-4}$ , during which period the speed of breakup is in the V-regime range of $-0.071\pm 0.015$ .

Figure 12. Curves are for 1: $Oh=0.14$ and 2: $Oh=0.16$ . (a) Changes in the evolution of the breakup speed ${\dot{r}}_{min}$ , against $r_{min}$ , around the critical point $Oh_{c}=0.15$ . Curve 1 shows the appearance of the low- $Oh$ V-regime around $r_{min}\approx 10^{-3}$ where the breakup speed dips to ${\dot{r}}_{min}\approx -0.071$ , in contrast to curve 2 where the speed immediately tends towards the value from the VI-regime of $-0.030$ . (b) The local Reynolds number $Re_{local}$ for the same values of $Oh$ shows that the low- $Oh$ V-regime coincides with a drop in $Re_{local}$ . The dashed line $Re_{local}=0.85$ is the value below which V-regime dynamics is to be expected.

The transition in flow behaviour at $Oh_{c}$ is also seen from the free-surface shapes in the vicinity of the breakup region. In figure 13(b) significant differences in the breakup geometry are shown for Ohnesorge numbers just below ( $Oh=0.14$ ) and above ( $Oh=0.16$ ) the critical value. At $r_{min}=10^{-4}$ , for $Oh=0.16$ (curve 1a) a long slender thread of length ${\approx}0.22$ connects the hemispherical volume above ( $z>0.27$ ) to a small satellite drop below ( $z<0.05$ ), whilst at $Oh=0.14$ (curve 1) the thread is much shorter at ${\approx}0.07$ . Consequently, the satellite drops formed in each case differ considerably, with $Oh=0.16$ forming a much fatter drop connected to a longer thin thread.

The free-surface profile in the low- $Oh$ V-regime (curve 2) does not have the slender symmetric shape in the breakup region that is characteristic of the conventional V-regime, which suggests that this regime is encountered in a weaker sense than the high- $Oh$ V-regime which satisfies all expected characteristics (see § 7). Using more stringent criteria for how the regimes are defined could result in this becoming a region of phase space where none of the similarity solutions are accurate; however, our definition is based on the breakup speed alone and this identifies it as a low- $Oh$ V-regime.

Figure 13. (a) Evolution of the free surface for $Oh=0.14$ , showing the development of a corner-like geometry (curve 3) while the flow is in the low- $Oh$ V-regime followed by the development of a thin thread once the VI-regime is entered (curve 4). (b) Comparison of free-surface shapes below the critical point, at $Oh=0.14$ (curves 1, 2), and above it, at $Oh=0.16$ (curves 1a, 2a), for 1, 1a: $r_{min}=10^{-4}$ and 2, 2a: $r_{min}=10^{-3}$ .

8.1 Identification of regimes

In figure 14 the VI- and low- $Oh$ V-regimes are shown on a phase diagram. The extreme change in behaviour at $Oh_{c}=0.15$ is clearly visible.

Figure 14. Phase diagram showing the low- $Oh$ V-regime and the VI-regime. The transition $r_{min}^{\rightarrow VI}$ follows the scaling $r_{min}=2.3\times 10^{-4}Oh^{-3.1}$ (curve 1) for $Oh>Oh_{c}=0.15$ . For $Oh<Oh_{c}$ the low- $Oh$ V-regime appears and is bounded by curves 2: $r_{min}=0.2Oh^{2}$ and 3: $r_{min}=10Oh^{2}$ as $r_{min}\rightarrow 0$ .

8.1.1 Entrance into the VI-regime for $Oh>Oh_{c}$

For $Oh>Oh_{c}$ , the transition into the VI-regime follows the expected scaling for V $\rightarrow$ VI of $r_{min}=B\,Oh^{-3.1}$ , with $B=2.3\times 10^{-4}$ . Curves 1–3 in figure 6(b) show that as expected $Re_{local}$ increases until it reaches ${\approx}1$ , a value characteristic of the VI-regime. The sharp changes in the gradients of these curves at the point where $Re_{local}$ is a maximum marks the beginning of the formation of satellite drops, i.e. the axial shift of the minimum bridge radius to $z(r=r_{min})>0$ , as can be seen from figure 11. This signals the appearance of satellite drops and entry into the VI-regime. Although $L_{r}/L_{z}$ increases at this point (figure 8), the thread remains slender. Once in the VI-regime, we do not see any evidence of an exit from this regime.

8.1.2 Entrance into the low- $Oh$ V-regime for $Oh<Oh_{c}$

When $Oh<Oh_{c}$ , appearance of the VI-regime is delayed by the presence of a low- $Oh$ V-regime. The entrance to this new regime follows a $r_{min}\sim Oh^{2}$ scaling as $r_{min}\rightarrow 0$ , which is characteristic of the I $\rightarrow$ VI transition, although here, as shown in § 9, this is an I $\rightarrow$ V transition.

Figure 12(b) corroborates the argument for a low- $Oh$ V-regime by showing, counter-intuitively, that once $Oh<Oh_{c}$ lower values of $Re_{local}$ are recovered. For example, for $Oh=0.16$ there is a dip to $Re_{local}=0.59$ , whilst for $Oh=0.14$ it falls as low as $Re_{local}=0.18$ . For $Oh=0.14$ , there is then a substantial period ( $4\times 10^{-4}<r_{min}<7\times 10^{-3}$ ) during which $Re_{local}$ remains at a value characteristic of the V-regime ( $Re_{local}<0.85$ ). Identifying the V-regime from the breakup speed in figure 12(a) results in a slightly smaller period ( $6\times 10^{-4}<r_{min}<3\times 10^{-3}$ ), a mismatch which appears to be caused by an initial lack of slenderness in the thread (figure 8), with $L_{r}/L_{z}<0.1$ only once $r_{min}<9\times 10^{-4}$ . This may explain why although ${\dot{r}}_{min}\approx -0.071$ in the new regime, there is no asymptotic approach to this value but rather a transient passing through speeds associated with a V-regime.

Entrance into the low- $Oh$ V-regime is triggered by the minimum radius moving away from $z=0$ as a satellite drop is formed. This is accompanied by other features that are characteristic of an I-regime, namely a corner-like free-surface shape (curve 3 in figure 13 a) and a decrease in the axial length scale $L_{z}$ due to the position of maximum velocity approaching the pinch point. However, in the newly created breakup region the axial velocity $U_{z}$ is still relatively small, so that $U_{z}L_{z}$ (i.e.  $Re_{local}$ ) is not large enough to trigger I-regime dynamics and instead a V-regime is encountered, a feature discovered in Castrejón-Pita et al. (2015). What remains unclear is why this mechanism persists at smaller $Oh$ and prevents transitions from the I-regime directly into the VI-regime.

8.1.3 Entrance into the VI-regime for $Oh<Oh_{c}$

As can be seen from figure 14, the transition out of the low- $Oh$ V-regime is quickly followed by a transition into the VI-regime for $Oh<Oh_{c}$ . This boundary appears to follow a $r_{min}\sim Oh^{2}$ scaling as $r_{min}\rightarrow 0$ , although for $Oh>4\times 10^{-2}$ it is approximately constant. The ${\sim}Oh^{2}$ scaling was predicted for the I $\rightarrow$ VI transition and observed above for the entrance into the low- $Oh$ V-regime. However, why this scaling is followed for V $\rightarrow$ VI is unclear.

9 Inertia-dominated flow (small $Oh$ )

In figure 15 flow profiles for $Oh=10^{-3}$ are shown. As seen in § 8, a satellite drop is formed, but now the free surface forms a corner at the pinch point (figure 15 b). In a narrow region nearby a rapid increase in $w_{fs}$ and $p_{fs}$ (figure 15 c,d) can be seen. Notably, $w$ is no longer $r$ -independent. Below the pinch point the free surface forms an angle close to $18.1^{\circ }$ (predicted for the I-regime in Day et al. (1998)) with the $z$ -axis (dashed line), but there is no agreement with the predicted angle above ( $112.8^{\circ }$ ), which is found to be $78^{\circ }$ . Therefore, no ‘overturning’ of the free surface is observed, as discussed further in § 11.3.

Figure 15. For $Oh=10^{-3}$ plots are for (a) the entire free surface, (b) a close-up of the breakup region (with dashed line showing the angle $18.1^{\circ }$ from the $z$ -axis predicted in Day et al. (1998)), (c) the axial velocity at the free surface, $w_{fs}$ , and (d) the pressure at the free surface, $p_{fs}$ , at 1: $r_{min}=10^{-1}$ (blue), 2: $r_{min}=10^{-2}$ (red), 3: $r_{min}=10^{-3}$ (green) and 4: $r_{min}=10^{-4}$ (brown).

The maximum axial velocity is shown in figure 5(a) to scale as ${\sim}r_{min}^{-1/2}$ as predicted by the similarity solution (2.1). In contrast to the other two regimes, the radial velocity $-{\dot{r}}_{min}$ scales in the same way and is close to $w_{max}$ when multiplied by a factor of 5, showing that in the I-regime there is no separation of scales in the $r$ and $z$ directions. This is supported by observations that the thread is not slender as pinch-off is approached, as can be seen from curve 1 in figure 8.

Recalling that regions where $\dot{l}_{min}$ is approximately constant result in a $r_{min}\sim {\it\tau}^{2/3}$ scaling indicative of the I-regime, curve 1 in figure 16(a) shows the appearance of this regime for $Oh=10^{-3}$ around $10^{-3}<r_{min}<10^{-2}$ . Computations show that $\dot{l}_{min}=-A_{I}^{3/2}\approx -0.5$ so that $A_{I}=0.63$ , which can be used to extract the breakup time and plot $r_{min}$ against ${\it\tau}$ in curve 1 of figure 5(b). The value of $A_{I}$ is close to that given in Eggers & Villermaux (2008) of $0.7$ .

Figure 16. Curves are for 1: $Oh=10^{-3}$ , 2: $Oh=4\times 10^{-3}$ and 3: $Oh=10^{-2}$ . (a) Identification of the I-regime, which is defined when $\dot{l}_{min}=-0.5\pm 0.15$ (dashed lines). (b) Evolution of the local Reynolds number: transition points out of the I-regime at $r_{min}=r_{min}^{I\rightarrow }$ are marked as circles and show that this transition occurs when $Re_{local}\approx 14$ (horizontal dashed line). The decrease of $Re_{local}$ in the I-regime follows the predicted scaling from the similarity solution (2.1) that $Re_{local}\sim r_{min}^{0.5}$ .

9.1 Identification of the I-regime

Using $\dot{l}_{min}$ to define the I-regime, and starting from the point at which the satellite drop is formed, i.e. $r_{min}<3\times 10^{-2}$ (figure 11), where the geometry of the I-regime starts to take shape, the boundaries of the I-regime are shown in figure 17.

Figure 17. Phase diagram showing the I-regime. The transition into this regime $r_{min}^{\rightarrow I}$ is constant, ${\approx}0.019$ , whilst the exit $r_{min}^{I\rightarrow }$ scales as $Oh^{2}$ with $r_{min}=45Oh^{2}$ .

9.1.1 Entrance into the I-regime

Similar to the V-regime, transitions into the I-regime $r_{min}^{\rightarrow I}$ are due to geometry, with a certain time required for the thread to form the corners which are characteristic of this regime. As one can see from curves 1 and 2 in figure 8, the breakup region is not slender and $L_{r}/L_{z}\approx 0.6$ remains approximately constant as $r_{min}\rightarrow 0$ .

9.1.2 Exit from the I-regime

The exit from the I-regime $r_{min}^{I\rightarrow }$ occurs when viscous effects become significant and the local Reynolds number drops below a critical value. Figure 16(b) shows that the reductions in $Re_{local}$ follow the ${\sim}r_{min}^{0.5}$ scaling predicted for the I-regime (2.1), which confirms that our definition of $Re_{local}$ is a good one. It is found that the I-regime is exited when $Re_{local}\approx 14$ . From § 8 we know that at this point the breakup will transition into the low- $Oh$ V-regime, and figure 17 shows that this entire boundary follows an ${\sim}Oh^{2}$ scaling.

10 A phase diagram for breakup

Having computed transitions into and out of the three regimes across parameter space, the results can be stitched together to produce the phase diagram shown in figure 18. Notably, once a regime is reached, which is around $r_{min}\approx 10^{-2}$ across $Oh$ , the transitions between regimes are relatively sharp, so that there are not vast regions of parameter space where none of the scaling laws are applicable. There are also no overlaps between the regimes, partially due to our choices of how to define being ‘in’ each regime. Consequently, rather than worry about transitions into and out of the different regimes, it is now sensible to talk about transitions between different regimes. This is how figure 18 was converted to produce the neat phase diagram figure 2(b) in § 2.2.3, where (i) transitions between the regimes are provided and (ii) only the asymptotic form of these transitions as $r_{min}\rightarrow 0$ are kept (i.e. using only the dashed lines in figure 18). It is of course possible that there is further complexity in the phase diagram below $r_{min}=10^{-4}$ . Here, analytical work is required as computational studies are inherently confined to a finite resolution below which the breakup behaviour can only be implied.

The most surprising results are the appearance of a low- $Oh$ V-regime which prevents I $\rightarrow$ VI transitions and the dramatic change in behaviour around a critical value $Oh_{c}=0.15$ . For $Oh>Oh_{c}$ the scaling for the V $\rightarrow$ VI transition follows that expected from theory ( $r_{min}^{V\rightarrow VI}\sim Oh^{-3.1}$ ), whilst both the transitions from the I-regime into the low- $Oh$ V-regime and the exit into the VI-regime follow the $r_{min}\sim Oh^{2}$ scaling predicted for the I $\rightarrow$ VI transition, but only as $r_{min}\rightarrow 0$ . Notably, it was shown that many of the transitions between regimes can be predicted from our measures $Re_{local}$ and $L_{r}/L_{z}$ .

Figure 18. Phase diagram for the breakup phenomenon which stitches together results from figure 7, figure 14 and figure 17, where specific expressions for the transitions between regimes can be found.

An experimentally verifiable prediction for the existence of a sharp transition at $Oh_{c}=0.15$ is that going from $Oh=0.16$ down to $Oh=0.14$ decreases the length of the thin liquid thread which connects the satellite drop to the main volume by a factor of three; see figure 13(b).

11 Discussion

Computations using a multiscale finite element method have allowed us to accurately construct a phase diagram for breakup and determine the transitions between different regimes. During the analysis, a number of features have been encountered which are worthy of further attention, and these are now considered.

11.1 Unexpected transitional behaviour

The possibility of multiple flow transitions was suggested in Eggers (2005) and first observed in Castrejón-Pita et al. (2015). The low- $Oh$ V-regime can be seen in figure 2 of Castrejón-Pita et al. (2015) at similar values of $Oh$ to those found here. However, in contrast to Castrejón-Pita et al. (2015), no evidence of a transient high- $Oh$ I-regime (their figure 4 for the case of $Oh=1.81$ ) has been observed. This could be because our method of defining the I-regime does not pick up transient regimes far from breakup; however, as can be seen from figure 18, our phase diagram is full for $r_{min}<10^{-2}$ so that any I-regime that has been missed will not appear in the later stages of breakup. Moreover, for all cases where a high- $Oh$ I-regime could exist, $Re_{local}$ remained at values well below that seen for the I-regime, where $Re_{local}>14$ , giving strong evidence for the absence of a high- $Oh$ I-regime.

To make unambiguous conclusions about the dominant forces as breakup is approached, one could directly extract the size of the inertial and viscous forces from the computed terms in the Navier–Stokes equations and compare their magnitudes in the breakup region. This would be a more advanced method for calculating $Re_{local}$ that may provide additional insight into the transitions between different regimes. However, we have found such measurements difficult to obtain in the breakup region due to the tiny elements found there and the singular nature of the dynamics. Such issues are amplified when determining second derivatives of the velocity, which itself is only approximated quadratically across each element, in order to calculate the viscous forces. This will be the focus of some of our future work in this area.

An intriguing possibility is that the low- $Oh$ V-regime is responsible for the transitional behaviour observed in Rothert et al. (2003), where water–glycerol mixtures in the approximate range $0.06<Oh<3$ were shown to transition V $\rightarrow$ VI at a constant $r_{min}$ rather than following the predicted scaling $r_{min}\sim Oh^{-3.1}$ . Given the relatively low $Oh$ considered in these experiments, it seems possible that it was the low- $Oh$ V-regime that was being encountered rather than the usual one. Furthermore, figure 18 shows that in the range $0.04<Oh<0.14$ the transition from the low- $Oh$ V-regime into the VI-regime ( $r_{min}^{V\rightarrow VI}$ ) is approximately constant. Unfortunately, no direct comparison to Rothert et al. (2003) is possible due to the different flow configurations used, so at present we can only speculate about the possible observation of a low- $Oh$ V-regime in their experiments.

11.2 Oscillatory convergence to the VI-regime

Figure 10(b) clearly illustrated that convergence towards the V- and VI-regimes is fundamentally different: at $Oh=10$ the bridge speed increases monotonically towards the value predicted by the similarity solution in the V-regime of $-0.071$ , whilst for $Oh=0.16$ the bridge speed converges in an oscillatory manner towards the speed predicted by the similarity solution in the VI-regime of $-0.030$ (preliminary results for the case where the plates confining the liquid bridge move apart are similar except the VI-regime is entered earlier, i.e. at larger $r_{min}$ ). The results are reinforced by figure 19, which shows the same oscillations when ${\dot{r}}_{min}$ is plotted against the time from breakup ${\it\tau}$ . Notably, the wavelength of the oscillations is approximately constant in logarithmic time $\ln {\it\tau}$ , a natural choice of variable, used e.g. in Eggers (1997), in the derivation of similarity solutions for breakup. Our results suggest that the similarity solutions in both the V- and VI-regimes are stable to perturbations, but that in the former case the associated eigenvalues of the perturbation are negative and purely real whilst in the latter case they are complex with negative real part.

Figure 19. Speed of breakup ${\dot{r}}_{min}$ against time from breakup ${\it\tau}$ for the case of $Oh=0.16$ , where oscillations characteristic of the VI-regime can be observed most strongly.

For the V-regime our findings agree with those in Eggers (2012), where it is shown that out of an infinite hierarchy of similarity solutions, only (2.2) is stable with purely real eigenvalues. Our findings for the VI-regime are entirely new and have not been predicted computationally, analytically or experimentally. The likely reasons for this are as follows:

1. (1) Computations. Very few simulations have captured the small scales resolved in our work for the axisymmetric Navier–Stokes system. However, in cases where they have, e.g. in Castrejón-Pita et al. (2015), the focus has been on $r_{min}$ against ${\it\tau}$ , whereas the oscillations are most clear when plotting ${\dot{r}}_{min}$ against ${\it\tau}$ or $r_{min}$ . What is surprising is that this behaviour has not been observed in slender jet codes which recover the similarity solution in Eggers (1993). Our own preliminary simulations indicate that the same oscillations can be observed in this setting too.

2. (2) Analysis. The current situation is summarised in p. 39 of Eggers & Villermaux (2008): thus far no one has analytically calculated the eigenvalues of the Navier–Stokes system for breakup. Consequently, it is unknown (i) whether the similarity solution is stable and (ii) if the associated eigenvalues are purely real or complex. The framework for the stability analysis is provided in appendix B of Brenner et al. (1996) and progress on similar flows has been achieved in Bernoff, Bertozzi & Witelski (1998), but extending these results to the Navier–Stokes pinch-off remains an open problem.

3. (3) Experiments. As with computations, so far the focus has been on $r_{min}$ against ${\it\tau}$ , where the oscillations can be missed, so that converting some of this data into ${\dot{r}}_{min}$ against ${\it\tau}$ or $r_{min}$ would be useful. However, here one faces the additional complication of noisy data, which could make the accurate extraction of derivatives a tricky procedure. Despite this, the oscillations observed are quite significant; at $Oh=0.16$ , between $r_{min}=10^{-3}$ – $10^{-2}$ the bridge speed doubles, from $-0.02$ to $-0.04$ (figure 10 b). If the length scale is $R=1~\text{mm}$ this occurs when the dimensional bridge radius is around $1{-}10~{\rm\mu}\text{m}$ , i.e. near the limits of optical resolution. Therefore, there is hope that experimental data can identify these changes.

11.2.1 Bumps on the VI similarity solution

Interestingly, continuing some of our simulations in the VI-regime beyond our self-imposed minimum radius of $r_{min}=10^{-4}$ resulted in the growth of a number of interesting features which, as far as we can see, appear unrelated to the aforementioned oscillations. These should be treated with caution, as the accuracy of the numerical scheme is close to its limit at such small radii, but the features appear to be robust.

As can be seen from figure 20(a) for $Oh=0.16$ , at around $r_{min}=5\times 10^{-5}$ the free surface has developed waves close to the pinch point. These feature in plots of axial velocity in figure 20(b), where one can see that the maximum absolute velocity (which will have $w<0$ ) shifts from around $z=0.08$ (in curve 1) to $z=0.23$ (curve 4). The velocity profile in curve 4, with the maximum absolute velocity near the pinch point, is closer to that predicted by the similarity solution in Eggers (1993).

Figure 20. Curves are for $Oh=0.16$ at 1: $r_{min}=2\times 10^{-4}$ , 2: $r_{min}=1.5\times 10^{-4}$ , 3: $r_{min}=10^{-4}$ and 4: $r_{min}=5\times 10^{-5}$ : (a) bumps on the free-surface profiles; (b) axial velocity distribution.

It is possible these bumps could be related to the iterated instabilities observed experimentally in Brenner, Shi & Nagel (1994) and interpreted as successive instabilities of the similarity solution of Eggers (1993). In this work, the instability was attributed to thermal fluctuations which were added to the lubrication formulation. However, it was pointed out that because this formulation only approximates the full Navier–Stokes system, such external noise may be unnecessary. The iterated instabilities have been seen computationally in McGough & Basaran (2006) for a surfactant-covered jet, but our results suggest that the instabilities can be generated without requiring additional effects such as noise or surfactants.

There could also be a relation to the destabilisation of similarity solutions due to finite-amplitude perturbations observed in Brenner et al. (1996), with intriguing similarities between our figure 20(a) and their figure 5 for the evolution of an unstable similarity solution.

11.3 Overturning of the free surface in the I-regime

Figure 15 showed that at $Oh=10^{-3}$ no overturning of the free surface was observed and, as far as we are aware, this does not contradict any previous findings in the liquid bridge geometry. Although simulations in Suryo & Basaran (2006) show overturning in the same geometry, this phenomenon only occurs for non-Newtonian cases. In contrast, for drop formation from an orifice, experimental studies show cases that do (Castrejón-Pita et al. 2012) and do not (Brenner et al. 1997) exhibit overturning, whilst simulations using the inviscid equations do show this feature (Schulkes 1994) and those using slender jet theory are unable to (Brenner et al. 1997).

The most comprehensive computational study of overturning is given in Wilkes et al. (1999), where the finite element method is used to study the effect of $Oh$ and the Bond number $Bo$ (i.e. gravitational forces) on this phenomenon. It is shown that the critical $Oh$ at which overturning is observed and the angles which the free-surface shape makes with the $z$ -axis at the pinch point depend on $Bo$ , which alters the geometry of the breakup. Furthermore, the angles recovered vary, with a maximum angle of $95^{\circ }$ observed, in contrast to the value of $112.8^{\circ }$ predicted from the inviscid theory (Day et al. 1998).

In the drop formation geometry, threads of arbitrary large lengths can be formed, whereas in our set-up the liquid is confined between two stationary plates. This appears to be the reason why no overturning has been observed even in cases of small $Oh$ which fall below the overturning limit in Wilkes et al. (1999). In particular, rather than a thread connected to a free drop, our geometry has a short thread (the satellite drop) connected to a hemispherical volume pinned to a plate.

11.4 Implications for breakup in CFD

Computational approaches that ‘go through’ the topological change of breakup all rely on cutting the thread, either manually (as in finite element methods) or automatically (as in volume-of-fluid approaches – VoF) once it reaches a specific $r_{min}$ , which may be implicitly related to the grid size (as in VoF). Therefore, it is of interest to know how accurate this approach is and what could be lost from a premature truncation. To do so, for the sake of argument consider a scheme with a fixed resolution (mesh size) of ${\rm\Delta}h=5\times 10^{-3}$ , giving a (relatively large) grid of $200\times 200$ for our problem, which cuts the thread wherever $r_{min}<{\rm\Delta}h$ .

The most obvious feature that can be lost from the cutoff is the satellite drop, as seen in Fawehinmi et al. (2005). Looking at figure 11, satellite drop formation, indicated by $z(r=r_{min})$ becoming non-zero, can clearly be seen for both $Oh=0.16$ (curve 2) and $Oh=0.5$ (curve 3). This occurs when the symmetric V-regime transitions into the asymmetric VI one, as can be seen clearly from free-surface profiles at $Oh=0.16$ in figure 9. For $Oh=0.16$ the pinch point moves when $r_{min}\approx 10^{-2}$ , whilst for $Oh=0.5$ this does not happen until $r_{min}\approx 10^{-3}$ . Therefore, using our cutoff, at $Oh=0.16$ the thread would be cut near to the correct pinch point, whilst for $Oh=0.5$ the thread would be severed at $z=0$ so that the formation of a satellite drop would be missed. Most worryingly, in the latter case the post-breakup state would be entirely wrong, as instead of having two breakup points a distance ${\rm\Delta}z=0.35$ apart (giving three distinct volumes) there would be just one at the centre of the drops (and hence two volumes). Notably, for $R=1~\text{mm}$ the minimum radius $r_{min}=10^{-3}$ corresponds to a dimensional radius of $1~{\rm\mu}\text{m}$ , so that these are truly macroscopic quantities.

To capture the features of breakup, such as satellite drops, one has to either (i) develop codes with huge levels of resolution/refinement, which will be extremely costly (particularly in 3-D) due to the multiscale nature of this phenomenon, or (ii) intelligently utilise the similarity solutions to take the codes through the topological change and up to a suitable post-breakup state. Such a scheme would be complex, but could lead to increased accuracy at a lower computational cost. By identifying in which regions of parameter space the different similarity solutions are accurate and devising methods to analyse in which flow regime a given breakup is occurring we have taken just some of the first steps in the development of such a scheme. Forthcoming works will extend these results to cases in which there is a strong externally driven elongation of the thread and/or interface formation physics which creates a singularity-free breakup; see Shikhmurzaev (2007).