1. Introduction
Single bubbles, droplets and particles can experience complex paths when moving freely under the influence of gravity in an otherwise quiescent fluid (Magnaudet & Eames Reference Magnaudet and Eames2000; Clift, Grace & Weber Reference Clift, Grace and Weber2005; Ern et al. Reference Ern, Risso, Fabre and Magnaudet2012; Mathai, Lohse & Sun Reference Mathai, Lohse and Sun2020; Bonnefis et al. Reference Bonnefis, Sierra-Ausin, Fabre and Magnaudet2024). Understanding the origin and nature of these irregular paths has been a longstanding concern in multiple disciplines, including mechanical and chemical engineering, aerodynamics and meteorology. To a large extent, the onset of the first non-vertical path is closely related to the primary wake instability that occurs beyond a critical Reynolds number even if the body moves at a constant speed and orientation (Ern et al. Reference Ern, Risso, Fabre and Magnaudet2012). The first step in understanding path instability is to examine the conditions under which the axisymmetric wake of a fixed body with different boundary conditions (e.g. no-slip for particles and free-slip for bubbles) first becomes unstable (Dandy & Leal Reference Dandy and Leal1989; Natarajan & Acrivos Reference Natarajan and Acrivos1993; Johnson & Patel Reference Johnson and Patel1999; Ghidersa & Dušek Reference Ghidersa and Dušek2000; Magnaudet & Mougin Reference Magnaudet and Mougin2007; Yang & Prosperetti Reference Yang and Prosperetti2007; Tchoufag, Magnaudet & Fabre Reference Tchoufag, Magnaudet and Fabre2013; Chiarini, Gauthier & Boujo Reference Chiarini, Gauthier and Boujo2025). In this context, Magnaudet & Mougin (Reference Magnaudet and Mougin2007) demonstrated that, regardless of the boundary conditions at the surface of the body, the axisymmetric wake becomes unstable when the maximum vorticity generated on the external side of the body surface exceeds a critical Reynolds-number-dependent threshold. Beyond this threshold, the steady wake exhibits a pair of counter-rotating trailing vortices generating a non-zero transverse lift force. Once the body is free to move, this force induces an oblique motion in the symmetry plane of the wake, ultimately leading to a non-vertical path (Mougin & Magnaudet Reference Mougin and Magnaudet2002, Reference Mougin and Magnaudet2006; Horowitz & Williamson Reference Horowitz and Williamson2010).
The close connection between the primary wake instability behind a fixed body and the onset of the first non-vertical path when such a body is free to move has been well established for bubbles and particles, as well as droplets with a dynamic viscosity significantly higher than that of the surrounding fluid. For example, in the case of a high-Reynolds-number bubble whose surface is free of surfactants, the threshold for wake instability corresponds to a critical aspect ratio
$\unicode{x03C7}\approx 2.1$
(where
$\unicode{x03C7}$
is the ratio of major and minor axes of the body) (Yang & Prosperetti Reference Yang and Prosperetti2007), which is close to the threshold for path instability (
$\approx 2.0$
) of a bubble that freely rises (Duineveld Reference Duineveld1995; Zenit & Magnaudet Reference Zenit and Magnaudet2008; Bonnefis et al. Reference Bonnefis, Sierra-Ausin, Fabre and Magnaudet2024; Shi et al. Reference Shi, Zhang and Magnaudet2025c
). Similarly, for a solid sphere, where the no-slip boundary condition applies at the surface, the threshold for wake instability occurs at a critical Reynolds number of approximately 210, in good agreement with the range of the first path instability (
$\in [210, 260]$
) of a freely rising light sphere (Jenny, Dušek & Bouchet Reference Jenny, Dušek and Bouchet2004; Horowitz & Williamson Reference Horowitz and Williamson2010; Auguste & Magnaudet Reference Auguste and Magnaudet2018). Lastly, for highly viscous droplets (i.e. those with
$\mu ^\ast \gg 1$
, where
$\mu ^\ast$
denotes the drop-to-fluid viscosity ratio) that behave similarly to solid particles, the relationship between wake and path instabilities appears evident. Specifically, Albert et al. (Reference Albert, Kromer, Robertson and Bothe2015) used direct numerical simulations (DNS) to investigate the path of corn oil droplets rising in pure water (for which
$\mu ^\ast \approx 46$
) and found that the path transition (from vertical to steady oblique) occurs at a critical Reynolds number equal to approximately 198, close to the primary wake instability threshold for a solid sphere (
$\approx 210$
). The slightly lower critical Reynolds number observed is expected, as the droplet undergoes slight deformation (
$\unicode{x03C7} \approx 1.05$
at the threshold), leading to an increase of the external surface vorticity (Magnaudet & Mougin Reference Magnaudet and Mougin2007) and, consequently, a reduction of the critical Reynolds number for wake instability.
The relationship between wake and path instabilities differs for droplets with low-to-moderate viscosity ratios (i.e. those with
$\mu ^\ast ={\mathcal{O}}(0.1{-}1)$
). For such droplets, the onset of the first path instability occurs at a significantly lower external surface vorticity (and thus a lower Reynolds number at a fixed aspect ratio, or vice versa) than that predicted for the primary wake instability using the criterion proposed by Magnaudet & Mougin (Reference Magnaudet and Mougin2007). A representative example of this phenomenon is the experiment by Wegener (Reference Wegener2009) (see also Wegener, Kraume & Paschedag Reference Wegener, Kraume and Paschedag2010), which investigated the motion of single toluene droplets of various sizes rising in water. The corresponding drop-to-fluid viscosity ratio was
$\mu ^\ast = 0.62$
(see table 2 of Wegener et al. (Reference Wegener, Kraume and Paschedag2010) for detailed physical properties). For droplets with an equivalent radius exceeding
$R\approx 1.1\,\text{mm}$
, the rising speed initially increased but then experienced a sudden decrease of approximately
$30\,\%$
, followed by pronounced oscillations around this reduced mean value. Furthermore, after several cycles of rising-speed oscillations, the path evolved from rectilinear to oblique (see figure 5.4 of Wegener (Reference Wegener2009)). Notably, at this threshold radius, the surface vorticity generated at the external side of the droplet in a fixed-droplet configuration (based on present DNS results, to be outlined in § 5.2) is only one-third of that associated with the primary wake instability (Magnaudet & Mougin Reference Magnaudet and Mougin2007; Shi et al. Reference Shi, Climent and Legendre2025b
). Similar observations for droplets with
$\mu ^\ast ={\mathcal{O}}(0.1{-}1)$
– particularly the presence of a critical droplet size beyond which the (mean) terminal rising velocity undergoes a sudden reduction – have been reported in earlier experiments (Klee & Treybal Reference Klee and Treybal1956; Thorsen, Stordalen & Terjesen Reference Thorsen, Stordalen and Terjesen1968). A comprehensive review of related studies can be found in Abdelouahab & Gatignol (Reference Abdelouahab and Gatignol2011) and, more recently, in Zhang et al. (Reference Zhang, Wang, Stevens and Fei2019) and Godé et al. (Reference Godé, Charton, Rachih, Lamadie, Elyakime, Legendre and Climent2025).
Since the detailed experiments of Wegener (Reference Wegener2009), several attempts have been made using DNS to replicate the first path instability of toluene droplets rising freely in water. Early studies in this direction (Bäumler et al. Reference Bäumler, Wegener, Paschedag and Bänsch2011; Engberg & Kenig Reference Engberg and Kenig2014; Wegener Reference Wegener2014) carried out simulations in a two-dimensional axisymmetric configuration, assuming that the flow remains axisymmetric during the initial stage of the first path instability (i.e. before the path transitions to an oblique trajectory). In this constrained flow set-up, the only possible cause of oscillations in the rising speed is the onset of axisymmetric deformations about the droplet’s minor axis. Indeed, once such shape oscillation modes become unstable, the rising speed oscillates as well, since both the vertical drag and the vertical added mass depend on the droplet cross-section and, consequently, on its shape (Magnaudet Reference Magnaudet2011; Lalanne, Tanguy & Risso Reference Lalanne, Tanguy and Risso2013; Shi et al. Reference Shi, Zhang and Magnaudet2025c
). However, predictions from these axisymmetric simulations indicated that the first unstable shape oscillation mode occurs beyond a critical size of
$R\approx 2.2\,\text{mm}$
(Bäumler et al. Reference Bäumler, Wegener, Paschedag and Bänsch2011), which is twice the size of the first path instability reported by Wegener (Reference Wegener2009). This discrepancy motivated subsequent studies to perform fully resolved three-dimensional simulations to better capture the first path instability (Bertakis et al. Reference Bertakis, Groß, Grande, Fortmeier, Reusken and Pfennig2010; Eiswirth et al. Reference Eiswirth, Bart, Atmakidis and Kenig2011). However, due to the slow development of axisymmetry-breaking processes (and thus the long physical time required), it was not possible until the recent works of Engberg & Kenig (Reference Engberg and Kenig2015) and Charin et al. (Reference Charin, Lage, Silva, Tuković and Jasak2019) for the first path instability at the threshold droplet size (
$R\approx 1.1\,\text{mm}$
) to be reasonably replicated. Specifically, DNS results in the fully developed regime from Charin et al. (Reference Charin, Lage, Silva, Tuković and Jasak2019; figure 11 therein) revealed that for droplets with
$R\geqslant 1\,\text{mm}$
, two pairs of counter-rotating streamwise vortices form in the wake, as inferred from the three-dimensional velocity field at the rear of the droplet. This indicates that the axisymmetric wake has already broken down for
$R\approx 1\,\text{mm}$
. However, it is important to note that beyond this threshold, the rising path may still remain rectilinear, since the transverse force remains zero because the wake, consisting of two pairs of vortex threads, retains its biplanar symmetry, an observation confirmed by the DNS results of Charin et al. (Reference Charin, Lage, Silva, Tuković and Jasak2019).
The biplanar-symmetric wake structure revealed by the DNS of Charin et al. (Reference Charin, Lage, Silva, Tuković and Jasak2019) intuitively suggests that the instability responsible for the symmetry breaking of the flow around a droplet with a low-to-moderate viscosity ratio is associated with the azimuthal wavenumber
$m=2$
(Ghidersa & Dušek Reference Ghidersa and Dušek2000; Yang & Prosperetti Reference Yang and Prosperetti2007). The mathematical nature of this symmetry breaking differs fundamentally from that observed in cases involving rising bubbles and settling or rising particles. In these cases, the first non-straight path is typically triggered by a mode with azimuthal wavenumber
$m=1$
, leading to a transition from an axisymmetric to a uniplanar-symmetric flow state (Jenny et al. Reference Jenny, Dušek and Bouchet2004; Yang & Prosperetti Reference Yang and Prosperetti2007; Tchoufag et al. Reference Tchoufag, Magnaudet and Fabre2013; Bonnefis et al. Reference Bonnefis, Sierra-Ausin, Fabre and Magnaudet2024). Given this distinction, it is not surprising that the criterion for the primary wake instability established for the latter case by Magnaudet & Mougin (Reference Magnaudet and Mougin2007) fails to predict the first path instability of toluene droplets rising in water.
For a better understanding of the underlying physical mechanisms driving the first path instability of droplets with low-to-moderate viscosity ratios, it is first necessary to examine more systematically a simplified configuration – the wake instability of the flow past a fixed droplet over a wide range of viscosity ratio and Reynolds number. The first attempt in this direction appears to be the study by Edelmann, Le Clercq & Noll (Reference Edelmann, Le Clercq and Noll2017), which reported on three-dimensional simulations of uniform flow past a spherical droplet at Reynolds numbers of
${\mathcal{O}}(100)$
. Interestingly, those authors reported that at a viscosity ratio of 0.5, the wake exhibited a biplanar-symmetric structure, a feature that was later identified in the wake of freely rising droplets (Charin et al. Reference Charin, Lage, Silva, Tuković and Jasak2019). Subsequent and more systematic numerical investigations have been reported (Rachih Reference Rachih2019; Godé Reference Godé2024; Shi et al. Reference Shi, Climent and Legendre2024a
). Specifically, Godé (Reference Godé2024) highlighted the significant role of the internal flow (i.e. the flow inside the droplet) in triggering the primary wake instability. In that work, the author observed an internal flow bifurcation for
$\mu ^\ast$
up to two under the constraint that the external flow remained axisymmetric. In fact, when the external flow was allowed to respond to this internal flow bifurcation, as shown in Shi et al. (Reference Shi, Climent and Legendre2024a
), the axisymmetric wake broke down, and the entire flow transitioned into a biplanar-symmetric flow similar to that reported in Edelmann et al. (Reference Edelmann, Le Clercq and Noll2017). However, the critical viscosity ratio below which this internal bifurcation occurs is not fixed; rather, it varies within the range of 1–10 depending on the external and internal Reynolds numbers (definitions provided in the next section) (Shi et al. Reference Shi, Climent and Legendre2024a
).
Based on the previous studies mentioned above, two key issues remain open regarding the primary wake and path instabilities of droplets with low-to-moderate viscosity ratios:
-
(i) Fixed droplet in a uniform flow. The physical mechanisms driving the internal flow bifurcation (Godé Reference Godé2024; Shi et al. Reference Shi, Climent and Legendre2024a ) remain unclear. Moreover, a criterion is still lacking for determining whether the axisymmetric flow is stable for a given set of parameters (in terms of Reynolds numbers and viscosity ratio). This criterion cannot be based solely on a fixed viscosity ratio, as demonstrated in Shi et al. (Reference Shi, Climent and Legendre2024a ).
-
(ii) Freely rising/settling case. The direct relationship between the wake instability caused by the internal flow bifurcation in the fixed-droplet case and the first path instability of a freely rising/settling droplet has yet to be established. Specifically, for the well-documented first path instability of toluene droplets freely rising in water, can we reasonably predict the critical droplet size using a criterion for the primary wake instability derived from studying fixed-droplet cases?
In this work, our aim is to address these two open issues, with a particular focus on the first – namely the primary wake instability induced by internal flow bifurcation. To this end, we revisit the problem of flow instability past a spherical droplet and conduct DNS over a wide range of dimensionless numbers using the JADIM code developed at IMFT (Legendre et al. Reference Legendre, Rachih, Souilliez, Charton and Climent2019; Rachih Reference Rachih2019; Rachih et al. Reference Rachih, Legendre, Climent and Charton2020; Godé et al. Reference Godé, Charton, Climent and Legendre2023; Godé Reference Godé2024). The paper is organised as follows. In § 2, we formulate the problem and outline the numerical approach. Section 3 provides an overview of the results, highlighting the connection between the internal flow bifurcation and the vorticity generated on the internal side of the droplet surface. A typical transition sequence with increasing internal Reynolds number for fixed external Reynolds number and viscosity ratio is discussed in § 4. In § 5.1, we present a physical explanation for the mechanism driving the internal flow bifurcation. The relationship between the primary wake instability of a fixed droplet and the first path instability when the droplet is free to move is explored in § 5.2, with a particular focus on the case of toluene droplets rising in water. Based on the confirmed relationship, § 5.3 presents the threshold droplet size for the internal bifurcation of a nearly spherical droplet moving freely in water, using the criterion for internal bifurcation proposed in the present work. Conclusions, along with the perspectives arising from this study, are presented in § 6.
2. Problem statement and numerical approach
We consider a spherical droplet of radius
$R$
, density
$\rho ^i$
, and dynamic viscosity
$\mu ^i$
that is fixed in a Newtonian fluid of density
$\rho ^e$
and dynamic viscosity
$\mu ^e$
. Far from the droplet interface, the external flow is a uniform stream along
$\boldsymbol{e}_x$
, described by
$\boldsymbol{u}^\infty = u_{\textit{rel}}\boldsymbol{e}_x$
, where
$u_{\textit{rel}}$
represents the slip velocity of the external fluid relative to the droplet. The entire flow field is governed by the incompressible Navier–Stokes equations:
\begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^k=0, \quad \rho ^k\left (\frac {\partial \boldsymbol{u}^k}{\partial t}+ \boldsymbol{u}^k\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}^k\right )= - \boldsymbol{\nabla }p^k +\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\unicode{x03C4}}^k, \end{equation}
where
$\boldsymbol{\unicode{x03C4}}^k=\mu ^k ( \boldsymbol{\nabla }\boldsymbol{u}^k+{}^{T}\boldsymbol{\nabla }\boldsymbol{u}^k )$
is the viscous part of the stress tensor
$\boldsymbol{\varSigma }^k=-p^k\boldsymbol{I}+\boldsymbol{\unicode{x03C4}}^k$
and
$\boldsymbol{u}^k$
and
$p^k$
denote the disturbed velocity and pressure, respectively. Here,
$k=i$
(likewise,
$k=e$
) refers to the fluid inside (outside) the droplet.
The boundary conditions are outlined below. At the surface of the droplet, the normal velocity must vanish due to the non-penetration condition, whereas the tangential velocity and shear stress must be continuous. These constraints yield the following boundary conditions at the droplet surface
$r=R$
:
\begin{align} \boldsymbol{u}^i\boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{u}^e\boldsymbol{\cdot }\boldsymbol{n}=0\,, \quad \boldsymbol{n}\times \boldsymbol{u}^i= \boldsymbol{n}\times \boldsymbol{u}^e\,, \quad \boldsymbol{n}\times (\boldsymbol{\unicode{x03C4}}^i\boldsymbol{\cdot }\boldsymbol{n})= \boldsymbol{n}\times (\boldsymbol{\unicode{x03C4}}^e\boldsymbol{\cdot }\boldsymbol{n}), \end{align}
where
$r=(x^2+y^2+z^2)^{1/2}$
is the distance from the droplet centre and
$\boldsymbol{n}$
is the outward unit normal to the droplet surface. In the far field, we assume that the disturbance induced by the droplet vanishes, implying that
$\boldsymbol{u}^k=\boldsymbol{u}^\infty$
as
$r\to \infty$
. Unless otherwise stated, the initial velocity field throughout the domain corresponds to the undisturbed state; that is, we set
$\boldsymbol{u}^k = \boldsymbol{u}^\infty$
everywhere (i.e. both inside and outside the droplet) at
$t = 0$
. The boundary conditions on the droplet interface defined by (2.2a–c
) are satisfied starting from the first time step of the simulation.
The steady-state solution of the problem is characterised by three dimensionless numbers: the viscosity ratio
$\mu ^\ast =\mu ^i/\mu ^e$
, the external Reynolds number
${\textit{Re}}^e$
and the internal Reynolds number
${\textit{Re}}^i$
. The latter two are defined as
\begin{equation} {\textit{Re}}^e=\frac {\rho ^e u_{\textit{rel}} (2R)}{\mu ^e}, \quad {\textit{Re}}^i=\frac {\rho ^i u_{\textit{rel}} (2R)}{\mu ^i}. \end{equation}
The drop-to-fluid density ratio can be expressed in terms of these three parameters as
$\rho ^\ast =\mu ^\ast \,{Re}^i/{Re}^e$
. In what follows, the viscosity ratio is varied from 0.01 to 100, allowing us to examine the evolution of the flow structure from the clean-bubble limit (
$\mu ^\ast \to 0$
) to the solid-sphere limit (
$\mu ^\ast \to \infty$
). At a given
$\mu ^\ast$
, the two Reynolds numbers are varied independently, whereas in a real drop–liquid system, only one of them is independent once the drop-to-liquid density ratio
$\rho ^\ast$
is specified. Thus, arbitrarily varying
${\textit{Re}}^i$
and
${\textit{Re}}^e$
implies artificially changing the density ratio while keeping the viscosity ratio
$\mu ^\ast$
fixed (see § 5 to relate this analysis to experiments).
The three-dimensional time-dependent simulations were performed using the JADIM code developed at IMFT. This code has been previously applied to simulate the three-dimensional flow around spherical bubbles and particles, as well as the associated hydrodynamic forces (Legendre & Magnaudet Reference Legendre and Magnaudet1998; Adoua, Legendre & Magnaudet Reference Adoua, Legendre and Magnaudet2009; Shi et al. Reference Shi, Rzehak, Lucas and Magnaudet2020, Reference Shi, Rzehak, Lucas and Magnaudet2021) and has been extended to compute three-dimensional flows around and inside spherical droplets (Legendre et al. Reference Legendre, Rachih, Souilliez, Charton and Climent2019; Rachih et al. Reference Rachih, Legendre, Climent and Charton2020; Godé et al. Reference Godé, Charton, Climent and Legendre2023; Shi et al. Reference Shi, Climent and Legendre2024a , Reference Shi, Climent and Legendre2025a ). The reader is referred to Shi et al. (Reference Shi, Climent and Legendre2024a ) for details regarding the numerical implementation, including the mesh grid, boundary conditions and validation tests confirming the reliability of the numerical approach. Of particular note is the numerical implementation of the boundary conditions on the droplet interface defined by (2.2a–c ). Using finite-difference discretisation, the tangential condition that simultaneously enforces the continuity of viscous shear stress and tangential velocity is implemented with second-order accuracy (Legendre et al. Reference Legendre, Rachih, Souilliez, Charton and Climent2019). The pressure fields inside and outside the droplet are solved separately via a pseudo-Poisson equation, with reference pressures imposed in both domains. These reference values are in fact linked to the Laplace pressure jump because the droplet is spherical. The only difference between the present study and that considered in Shi et al. (Reference Shi, Climent and Legendre2024a ) is the type of undisturbed flow (the linear shear flow now being a uniform flow). It is worth noting that the mesh grid used in Shi et al. (Reference Shi, Climent and Legendre2024a ), which is also used in the present work, features highly refined grid cells near the droplet interface, with at least five nodes located within both the internal and external boundary layers for Reynolds numbers of up to 1000. This refinement ensures that both the internal and external flows are accurately resolved at high Reynolds numbers, particularly within the boundary layers on both sides of the droplet interface.
3. Overview of the results
3.1. Identification of internal and external flow bifurcations
The breaking of axisymmetry can be tracked by examining the perturbation energy, for which a convenient measure is the mean kinetic energy of the azimuthal velocity component (hereafter referred to as the azimuthal energy) (Thompson, Leweke & Provansal Reference Thompson, Leweke and Provansal2001; Magnaudet & Mougin Reference Magnaudet and Mougin2007):
\begin{equation} E^k=\frac {1}{\rho ^e V_s u_{\textit{rel}}^2} \int _{V^k}{\rho ^k ||\boldsymbol{u}_\varphi ^k||^2{\textrm{d}}V^k}, \end{equation}
where
$V_s=4 \pi\! R^3/3$
is the volume of the droplet and
$\boldsymbol{u}_\varphi ^k$
is the azimuthal component of the local velocity. Here,
$E^k$
(respectively
$V^k$
) with
$k=i$
or
$e$
denotes the azimuthal energy (respectively the domain) inside or outside the droplet. The azimuthal energy of the entire flow field is then given by
$E= E^i +E^e$
, which becomes positive as soon as the axisymmetry of the base flow breaks down.
Two different types of bifurcation can be identified based on the behaviour of
$E^i$
and
$E^e$
. To illustrate this, we consider a series of cases with
$({Re}^e,{Re}^i) = (300, 1000)$
but with varying
$\mu ^\ast$
. Figure 1(a) presents the time evolution of the three azimuthal energy components for a low-viscosity-ratio droplet (
$\mu ^\ast = 0.5$
). The axisymmetry of the base flow breaks down at
$t\approx 60\,R/u_{\textit{rel}}$
, as indicated by the onset of growth in the total azimuthal energy
$E$
, which peaks at
$t\approx 80\,R/u_{\textit{rel}}$
before stabilising at a slightly lower value beyond
$t\approx 140\,R/u_{\textit{rel}}$
. During this transition,
$E^i$
maintains a substantial magnitude relative to
$E$
, while
$E^e$
remains negligibly small in the initial stages (approximately for
${\textit{tu}}_{\textit{rel}}/R$
increasing from 60 to about 70). Figure 1(b) displays isosurfaces of the streamwise vorticity
$\omega _x$
at selected times during the transient. This vortical component becomes non-zero as soon as the bifurcation occurs. In the early stages of the bifurcation (
$t=66\,R/u_{\textit{rel}}$
),
$\omega _x$
is significant only inside the droplet, consistent with the initially negligible
$E^e$
observed in figure 1(a). As time progresses, the disturbance, represented by
$\omega _x$
isocontours, grows and spreads outside the droplet at
$t\approx 70\,R/u_{\textit{rel}}$
, leading to the formation of four vortex threads in the droplet wake. For the case under consideration, this wake structure remains stable in the fully three-dimensional developed state.
Characteristics of an internal flow bifurcation for
$(\mu ^\ast ,{Re}^e,{Re}^i) = (0.5, 300, 1000)$
. (a) Total, internal and external azimuthal energy as a function of time. (b) Isosurfaces of the streamwise vorticity,
$\omega _x R/u_{\textit{rel}} = \pm 0.2$
, at three selected time instants (indicated by numbers in b and marked as circles in a). Grey and black threads correspond to positive and negative
$\omega _x$
, respectively.
Same as figure 1, but for an external flow bifurcation in the case
$(\mu ^\ast ,{Re}^e,{Re}^i) = (20, 300, 1000)$
. In (b), the isosurfaces correspond to
$\omega _x R/u_{\textit{rel}} = \pm 0.1$
. The left part displays the vortical structure only inside the droplet and in the downstream half-space where the sign of
$\omega _x$
in the wake is positive.
For comparison, figure 2 presents the corresponding evolution for a droplet with a viscosity ratio of
$\mu ^\ast = 20$
, which is 40 times larger than in the previous case. Compared with the low-
$\mu ^\ast$
case, the key difference in the evolution lies in the internal energy
$E^i$
, which now remains vanishingly small throughout the transition (figure 2
a). Additionally, the wake structure during the transition differs between the two cases. As shown in figure 2(b), the wake now comprises only one pair of streamwise vortices, in contrast to the two pairs observed in the low-
$\mu ^\ast$
case. From the perspective of linear dynamical systems theory (Ghidersa & Dušek Reference Ghidersa and Dušek2000; Yang & Prosperetti Reference Yang and Prosperetti2007), this wake structure indicates a symmetry breaking driven by a mode with azimuthal wavenumber
$m=1$
. In contrast, in the low-
$\mu ^\ast$
case, the presence of four vortex threads in the wake suggests that the symmetry breaking is caused instead by a mode with an azimuthal wavenumber
$m=2$
.
Hereafter, we denote by an internal bifurcation the type of bifurcation that occurs in the low-
$\mu ^\ast$
case, where the internal azimuthal energy
$E^i$
remains significant throughout the transition. Conversely, we refer to as an external bifurcation the instability occurring in the high-
$\mu ^\ast$
case, where
$E^i$
remains vanishingly small throughout the transition. Figure 3(a) summarises the results for
$E^i$
and
$E^e$
in the fully developed state obtained at
$({Re}^e,{Re}^i) = (300, 1000)$
over a wide range of
$\mu ^\ast$
. Based on this classification, the internal bifurcation occurs for
$\mu ^\ast$
smaller than approximately 12, while the external bifurcation occurs for
$\mu ^\ast$
larger than about 15. (Note that, due to the distinct physical origins of the internal and external flow bifurcations, there is no reason to expect the critical viscosity ratio below which internal bifurcation occurs and that beyond which wake instability arises to coincide. The fact that, at
$({Re}^e, {\textit{Re}}^i) = (300, 1000)$
, both transitions appear to take place within the interval
$12 \leqslant \mu ^\ast \leqslant 15$
is a coincidence rather than a general result.)
Results for (a) the azimuthal energies
$E^i$
and
$E^e$
in the fully developed state and (b) the maximum surface vorticity as a function of viscosity ratio
$\mu ^\ast$
obtained at steady state at
$({Re}^e,{Re}^i)= (300, 1000)$
. In both panels, filled symbols in red (green) denote the onset of an internal (external) flow bifurcation.
For a uniform flow past a solid particle or a clean bubble, the axisymmetry of the flow breaks down via the external bifurcation when the maximum vorticity generated on the external side of the body surface exceeds a critical
${\textit{Re}}^e$
-dependent value (Magnaudet & Mougin Reference Magnaudet and Mougin2007). To determine whether this empirical criterion also applies to a droplet, we examine the evolution of the maximum surface vorticity with the viscosity ratio. Following Magnaudet & Mougin (Reference Magnaudet and Mougin2007), we decompose the vorticity
$\boldsymbol{\omega }^k = \boldsymbol{\nabla }\times \boldsymbol{u}^k$
into a normal component
$(\boldsymbol{\omega }^k \boldsymbol{\cdot }\boldsymbol{n}) \boldsymbol{n}$
and an azimuthal component
$\boldsymbol{\omega }^k - (\boldsymbol{\omega }^k \boldsymbol{\cdot }\boldsymbol{n}) \boldsymbol{n}$
. Since the base flow is axisymmetric, the primary vorticity field contains only an azimuthal component. Moreover, due to the viscosity contrast (except when
$\mu ^\ast = 1$
), the surface vorticity (and in particular, its azimuthal component) is discontinuous across the interface. We therefore define two distinct azimuthal surface vorticities:
\begin{equation} \boldsymbol{\omega }_S^i = \lim _{r \rightarrow R^-} \boldsymbol{\omega }^i - (\boldsymbol{\omega }^i \boldsymbol{\cdot }\boldsymbol{n}) \boldsymbol{n}, \quad \boldsymbol{\omega }_S^e = \lim _{r \rightarrow R^+} \boldsymbol{\omega }^e - (\boldsymbol{\omega }^e \boldsymbol{\cdot }\boldsymbol{n}) \boldsymbol{n}, \end{equation}
where
$\boldsymbol{\omega }_S^i$
(
$\boldsymbol{\omega }_S^e$
) denotes the azimuthal surface vorticity on the internal (external) side of the droplet surface. The maximum values of these quantities are denoted as
$\omega _s^i$
and
$\omega _s^e$
, respectively, where
$\omega _s^k = \max { (||\boldsymbol \omega _S^k|| )}$
. Unless stated otherwise, the results for
$\omega _s^i$
and
$\omega _s^e$
hereafter correspond to values obtained in the fully developed state of an imposed axisymmetric configuration, which are representative of the flow just prior to the onset of bifurcation. Besides, we refer to
$\omega _s^e$
(
$\omega _s^i$
) as the maximum external (internal) surface vorticity since
$\boldsymbol{\omega }^k = \boldsymbol{\omega }_S^k$
in an axisymmetric configuration.
Figure 3(b) (green symbols) presents the maximum external surface vorticity
$\omega _s^e$
as a function of
$\mu ^\ast$
. Clearly,
$\omega _s^e$
increases with increasing
$\mu ^\ast$
and exceeds approximately
$13.5 u_{\textit{rel}}/R$
(corresponding to the value at
$\mu ^\ast = 15$
), beyond which the external bifurcation occurs. This threshold is close to the critical value (
$13.8 u_{\textit{rel}}/R$
for
${\textit{Re}}^e = 300$
) predicted by Magnaudet & Mougin (Reference Magnaudet and Mougin2007; see (4.1) therein). Also shown in figure 3(b) (red symbols) is the corresponding maximum internal surface vorticity
$\omega _s^i$
. In contrast to
$\omega _s^e$
,
$\omega _s^i$
decreases with increasing
$\mu ^\ast$
. The absence of an internal bifurcation for
$\mu ^\ast \gtrsim 12$
can be attributed to
$\omega _s^i$
falling below a critical threshold. We discuss this threshold in more detail in the next section.
3.2. The internal flow bifurcation
In the previous section, we observed that for the series of cases with
$({Re}^e, {\textit{Re}}^i) = (300, 1000)$
, an internal bifurcation sets in when the viscosity ratio
$\mu ^\ast$
is smaller than approximately 12. However, the threshold
$\mu ^\ast$
can vary significantly with
${\textit{Re}}^e$
and
${\textit{Re}}^i$
(Edelmann et al. Reference Edelmann, Le Clercq and Noll2017; Godé Reference Godé2024; Shi et al. Reference Shi, Climent and Legendre2024a
), making it difficult to determine the internal bifurcation regime based solely on
$\mu ^\ast$
. In this section, we demonstrate that the internal bifurcation is closely related to the maximum internal surface vorticity
$\omega _s^i$
in the base flow and occurs when
$\omega _s^i$
exceeds a critical value, which can be satisfactorily fitted using only the internal Reynolds number
${\textit{Re}}^i$
.
(a) Internal azimuthal energy,
$E^i$
, in the fully developed state as a function of the viscosity ratio,
$\mu ^\ast$
, for various
${\textit{Re}}^i$
at
${\textit{Re}}^e = 200$
. (b) Maximum internal surface vorticity as a function of
${\textit{Re}}^i$
. In (b), in addition to the data at selected
${\textit{Re}}^i$
values shown in (a), an additional data series with increasing
${\textit{Re}}^i$
for
$(\mu ^\ast , {\textit{Re}}^e) = (0.5, 200)$
is also included (denoted by a thin dashed line and star symbols). In both panels, filled symbols indicate the onset of internal flow bifurcation. In (b), for each iso-
${\textit{Re}}^i$
data series,
$\mu ^\ast$
increases from top to bottom, and the thick black line represents the prediction from (3.3).
We begin by examining the regime of internal bifurcation in the parameter space
$(\mu ^\ast , {\textit{Re}}^i)$
. To this end, we fix
${\textit{Re}}^e$
at 200, ensuring that no external bifurcation occurs even in the solid-sphere limit
$\mu ^\ast \to \infty$
(Johnson & Patel Reference Johnson and Patel1999; Citro et al. Reference Citro, Tchoufag, Fabre, Giannetti and Luchini2016). Figure 4(a) presents the internal azimuthal energy
$E^i$
in the fully developed state as a function of the viscosity ratio
$\mu ^\ast$
for various values of
${\textit{Re}}^i$
. Under the selected
${\textit{Re}}^e$
, no internal bifurcation occurs for
${\textit{Re}}^i \leqslant 300$
. However, at higher
${\textit{Re}}^i$
, the threshold
$\mu ^\ast$
– below which the bifurcation sets in – increases from 3 at
${\textit{Re}}^i = 400$
to 15 at
${\textit{Re}}^i = 1500$
. Figure 4(b) shows the corresponding maximum internal surface vorticity,
$\omega _s^i$
. Unlike figure 4(a), the results are plotted against
${\textit{Re}}^i$
to highlight the dependence of the critical
$\omega _s^i$
on
${\textit{Re}}^i$
. For each iso-
${\textit{Re}}^i$
data series,
$\omega _s^i$
decreases with increasing
$\mu ^\ast$
, following a trend similar to that observed in figure 3(b). These results indicate that the critical
$\omega _s^i$
decreases as
${\textit{Re}}^i$
increases. To evaluate the consistency of this trend, we carried out an additional series of simulations with increasing
${\textit{Re}}^i$
while keeping
$(\mu ^\ast , {\textit{Re}}^e)$
fixed at
$(0.5, 200)$
. The corresponding results for
$\omega _s^i$
are shown in figure 4(b) (dashed line and star symbols). Under this condition, the internal bifurcation takes place as
${\textit{Re}}^i$
exceeds approximately 326, corresponding to a critical
$\omega _s^i$
of about
$6.56\,u_{\textit{rel}}/R$
.
We collect the critical values of
$\omega _s^i$
, denoted as
$\omega _c^i$
, for the four data series with
${\textit{Re}}^i \gt 300$
and fit them to a power-law relation in
${\textit{Re}}^i$
. The resulting empirical relation is
\begin{equation} \omega _c^i\, R/u_{\textit{rel}} \approx 1+0.33\big({\textit{Re}}^i/1000\big)^{-2.5}. \end{equation}
Although (3.3) is obtained for
${\textit{Re}}^e = 200$
, the predicted
$\omega _c^i$
is generally applicable to different values of
${\textit{Re}}^e$
. To verify this, we formed two additional series of simulations: one at
${\textit{Re}}^i=500$
and the other at
${\textit{Re}}^i=1000$
, varying
${\textit{Re}}^e$
from 50 to 500 in both cases. Figure 5 presents the resulting
$\omega _s^i$
as a function of the viscosity ratio, with cases involving the onset of internal bifurcation marked by filled symbols. According to (3.3), the critical
$\omega _s^i$
for the internal bifurcation is approximately 2.9 at
${\textit{Re}}^i=500$
and 1.3 at
${\textit{Re}}^i=1000$
. These two predictions are represented by horizontal dashed lines in figure 5(a,b), satisfactorily distinguishing the cases with internal bifurcation from those without.
Maximum internal surface vorticity as a function of viscosity ratio
$\mu ^\ast$
for various
${\textit{Re}}^e$
(distinguished by coloured symbols) at (a)
${\textit{Re}}^i=500$
and (b)
${\textit{Re}}^i=1000$
. In both panels, filled symbols denote cases where internal bifurcation occurs, and the horizontal dashed line represents the corresponding
$\omega _c^i({Re}^i)$
according to (3.3).
Based on the discussion above, we may state that correlation (3.3) can serve as an empirical criterion to determine whether the internal flow corresponding to a given set
$(\mu ^\ast , {\textit{Re}}^e, {\textit{Re}}^i)$
is stable or not. Specifically, given the maximum internal surface vorticity
$\omega _s^i$
at the internal Reynolds number
${\textit{Re}}^i$
under consideration, the internal flow is unstable (stable) if
$\omega _s^i(\mu ^\ast , {\textit{Re}}^e, {\textit{Re}}^i)$
is greater (smaller) than
$\omega _c^i({Re}^i)$
. This correlation helps to explain the significant variation in the threshold viscosity ratio,
$\mu _c^\ast$
, for the internal bifurcation. In particular, since
$\omega _c^i$
decreases with increasing
${\textit{Re}}^i$
, internal bifurcation is more likely to occur for droplets with larger
${\textit{Re}}^i$
. This explains why, in all previous studies (Edelmann et al. Reference Edelmann, Le Clercq and Noll2017; Rachih Reference Rachih2019; Godé Reference Godé2024; Shi et al. Reference Shi2024a
), internal bifurcation has generally been observed at relatively large
${\textit{Re}}^i$
(typically,
${\textit{Re}}^i \geqslant 300$
). On the other hand, since
$\omega _c^i$
according to (3.3) is independent of
${\textit{Re}}^e$
, while
$\omega _s^i$
increases with
${\textit{Re}}^e$
, the threshold viscosity ratio
$\mu _c^\ast$
at a given
${\textit{Re}}^i$
increases with increasing
${\textit{Re}}^e$
, as shown in figure 5.
4. Transition sequence
In this section, we focus on the series of cases with
$(\mu ^\ast , {\textit{Re}}^e) = (0.5, 200)$
and examine how the flow structure evolves with increasing
${\textit{Re}}^i$
. This corresponds to varying the density ratio. In § 5.2, we show that the phenomena described below can indeed be observed for realistic physical properties of liquid–liquid systems. Since the whole problem depends on
$(\mu ^\ast , {\textit{Re}}^e, {\textit{Re}}^i)$
, similar asymmetric flow structures can also arise with increasing
${\textit{Re}}^e$
at a fixed
$(\mu ^\ast , {\textit{Re}}^i)$
or with decreasing
$\mu ^\ast$
at fixed
$({Re}^e, {\textit{Re}}^i)$
. We did not explore in detail the bifurcation sequence under the latter two conditions, as doing so would require significant computational resources. However, it should be noted that while the critical
${\textit{Re}}^e$
or
$\mu ^\ast$
for the first internal flow bifurcation can be reasonably predicted using the criterion (3.3), the sequence of higher-order bifurcations with increasing
${\textit{Re}}^e$
or decreasing
$\mu ^\ast$
may differ from that observed with increasing
${\textit{Re}}^i$
.
4.1. Axisymmetric flow regime
For the series of cases with
$(\mu ^\ast , {\textit{Re}}^e) = (0.5, 200)$
, our three-dimensional simulations indicate that the axisymmetry of the flow breaks down through an internal bifurcation as
${\textit{Re}}^i$
exceeds a critical value approximately equal to 326 (see the star symbols in the inset of figure 4
b).
(a) Streamlines and (b) isocontours of the azimuthal vorticity
$\omega _\phi$
around the droplet for
${\textit{Re}}^i = 50$
(top panels) and
${\textit{Re}}^i = 325$
(bottom panels). For both cases,
$(\mu ^\ast , {\textit{Re}}^e) = (0.5, 200)$
. In (a), the vertical red line denotes
$x = 0$
. In (b), coloured lines represent
$-\omega _\phi R/u_{\textit{rel}} = 1$
(red), 2 (green), 3 (blue), 4 (magenta) and 5 (navy).
Figure 6(a) illustrates the streamlines around the droplet for
${\textit{Re}}^i = 50$
(top panel) and
${\textit{Re}}^i = 325$
(bottom panel), with the latter corresponding to the case just before the onset of bifurcation. In both cases, the external streamlines at the rear remain attached to the droplet, indicating the absence of a standing eddy in the wake prior to the onset of internal bifurcation. The internal flow structure resembles a Hill (spherical) vortex (Hill Reference Hill1894), although a fore–aft asymmetry with respect to
$x = 0$
(marked by the red vertical line in figure 6
a) can be inferred from the results at
${\textit{Re}}^i = 50$
. This fore–aft asymmetry becomes more evident when examining the isocontours of the azimuthal vorticity
$\omega _\phi$
, as shown in figure 6(b). Near the front of the droplet, the internal isocontours tilt towards the stagnation point, instead of aligning horizontally along the symmetry axis as in a Hill vortex. This tilting is more pronounced at higher
${\textit{Re}}^i$
: close to the stagnation point, the edges of the internal isocontours align almost parallel to the droplet surface. This strong tilting of azimuthal vorticity has also been observed in the wake of an oblate spheroidal bubble (Magnaudet & Mougin Reference Magnaudet and Mougin2007; Yang & Prosperetti Reference Yang and Prosperetti2007), where the external
$\omega _\phi$
isocontours tilt so that they align nearly perpendicular to the symmetry axis at the threshold of the external bifurcation.
4.2. Biplanar-symmetric flow regime
The internal flow bifurcation sets in beyond
${\textit{Re}}^i \approx 326$
for
$(\mu ^\ast , {\textit{Re}}^e) = (0.5, 200)$
. Following this bifurcation, the flow transits from an axisymmetric to a biplanar-symmetric structure (figure 1
b). Once the axisymmetry breaking has saturated, the flow in all cases within this regime remains steady, indicating that the bifurcation is regular. Figure 7(a) presents the total azimuthal energy in the final state as a function of
${\textit{Re}}^i$
. The results clearly indicate that the bifurcation is supercritical. Figure 7(b) shows the time evolution of the total azimuthal energy at
${\textit{Re}}^i = 345$
. After linear transient growth at a normalised growth rate of approximately 0.025 (indicated by a dashed straight line), the initial deviation from linearity levels off with a decreasing growth rate, further confirming the supercritical nature of the bifurcation (Strogatz Reference Strogatz2018, p. 82).
(a) Variation of the total azimuthal energy of the steady state,
$E$
, with the internal Reynolds number
${\textit{Re}}^i$
close to the threshold. (b) Plot of
$E$
as a function of time for
${\textit{Re}}^i = 345$
. In both panels, the straight dashed line highlights the linear scaling. In (b), the dashed line represents a normalised energy growth rate of 0.025.
Figure 8 presents the streamwise vorticity structure in the fully developed state for
${\textit{Re}}^i = 345$
. The resulting configuration, consisting of four vortex threads of equal intensity, closely resembles that observed in figure 1(b). Due to the entrainment of fluid elements by these vortex threads, two distinct symmetry planes exist. The first, denoted as the
$y=0$
plane (coloured blue), is characterised by the inward motion of fluid elements towards the symmetry axis. The second plane, denoted as the
$z=0$
plane (coloured green), is associated with an outward motion away from the symmetry axis. Note that the positions of the two symmetry planes are determined by the initial disturbance. In our numerical set-up, a weak streamwise linear shear flow of
$10^{-4}\, y (u_{\textit{rel}}/R)\, \boldsymbol{e}_x$
is imposed to trigger the bifurcation, thereby prescribing the locations of the two symmetry planes. The velocity variation across the droplet scale is
$10^{-4}u_{\textit{rel}}$
, negligibly small compared with the ambient flow.
Isosurfaces of the streamwise vorticity,
$\omega _x R/u_{\textit{rel}} = \pm 0.05$
, past a droplet at
${\textit{Re}}^i=345$
(grey and black threads correspond to positive and negative values, respectively). The flat surface in green (blue) highlights the symmetry plane in which the flow diverges (converges).
To further illustrate the biplanar symmetry of the flow structure, figure 9 presents the two-dimensional streamlines of the disturbance in selected cross-stream planes. The disturbance is obtained by subtracting the streamwise velocity component from the full velocity field, i.e. the plotted streamlines correspond to
$\boldsymbol{u}^k - u_x^k \boldsymbol{e_x}$
. At a distance of
$2R$
upstream of the droplet (figure 9
a), all streamlines radiate outward from the symmetry axis connected to the front stagnation point, indicating that the flow remains axisymmetric at this location. At the cross-stream plane passing through the droplet centre (figure 9
b), the axisymmetry of the inside base flow effectively breaks into four energetic vortices. The internal flow structure closely resembles that associated with the azimuthal wavenumber
$m=2$
in the wake of a solid sphere (see figure 13c in Ghidersa & Dušek (Reference Ghidersa and Dušek2000)). As this asymmetric internal disturbance influences downstream flows (figure 9
c), vortex pairs at the same
$y$
level gradually deviate from the symmetry plane
$y=0$
while simultaneously converging towards each other.
Two-dimensional streamlines of the disturbance
$\boldsymbol{u}^k - u_x^k \boldsymbol{e_x}$
in selected cross-stream planes for
${\textit{Re}}^i=345$
with
$x/R =$
(a) −2, (b) 0 and (c) 5. In each panel, the red circle represents the boundary of a cylindrical surface
$(y^2+z^2)^{1/2} = R$
. The thick horizontal blue line (vertical green line) denotes the symmetry plane
$y=0$
(
$z=0$
), as shown in figure 8.
4.3. Uniplanar-symmetric flow regime
A secondary bifurcation occurs as
${\textit{Re}}^i$
exceeds
$\approx 370$
for
$(\mu ^\ast , {\textit{Re}}^e) = (0.5, 200)$
. This new bifurcation disrupts the biplanar-symmetric flow structure, resulting in a flow with a single plane of symmetry in the fully developed state.
Taking the case at
${\textit{Re}}^i=375$
as an example, figure 10(a) shows the time evolution of the total energy
$E$
during this transition. The primary bifurcation sets in at
${\textit{tu}}_{\textit{rel}}/R \approx 300$
, from which
$E$
initially exhibits linear growth before decreasing and temporarily reaching a first plateau value as
${\textit{tu}}_{\textit{rel}}/R$
exceeds
$\approx 400$
. Figure 11(a) illustrates the vortical structure at this time. The flow remains biplanar-symmetric, similar to that observed in the fully developed state for
${\textit{Re}}^i = 345$
(see figure 8). Supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.10565 provides more detail on the time evolution of the vortical structure for this case. The evolution clearly shows that, during the primary bifurcation, the vortical perturbations are propagated outward, i.e. from the droplet interior to the exterior. Also, the vortical structures on the internal and external sides of the interface do not perfectly align, owing to the viscosity contrast, which influences the transfer of vorticity across the interface. Until
${\textit{tu}}_{\textit{rel}}/R \approx 400$
, the evolution of energy and flow structure closely resembles that of the biplanar flow regime, confirming that the primary bifurcation remains supercritical. However, at
${\textit{Re}}^i=375$
, the biplanar-symmetric flow structure is unstable. As shown in figure 10(a), shortly after stabilising at the first plateau,
$E$
again increases and when
${\textit{tu}}_{\textit{rel}}/R$
exceeds 650 it re-stabilises at a second plateau value approximately twice as high as the first one. Figure 11(b,c) depicts the vortical structures at two instants: one during the secondary transition (
${\textit{tu}}_{\textit{rel}}/R=530$
; figure 11
b) and the other in the fully developed state (
${\textit{tu}}_{\textit{rel}}/R=2000$
; figure 11
c). These results reveal that during the transition, the flow symmetry with respect to
$y=0$
breaks down. Specifically, one of the vortex pairs at the same
$y$
level (the pair with
$y\gt 0$
in figure 11) shrinks while the other grows over time. Throughout this evolution, the flow symmetry with respect to
$z=0$
persists. More details on the evolution of the vortical structure during the secondary bifurcation can be found in supplementary movie 1.
(a) The azimuthal energy,
$E$
, and the lift coefficient,
$C_L$
, as functions of time for
${\textit{Re}}^i = 375$
. (b) Variation of
$E$
and
$C_L$
in the fully developed state with the internal Reynolds number
${\textit{Re}}^i$
for
${\textit{Re}}^i$
ranging from 350 to 400. In (a), the two circles mark
$t u_{\textit{rel}}/R = 400$
and 530, respectively. The dashed black line near
$E(t)$
shows that the initial deviation from linearity levels off with a decreasing growth rate, indicating the supercritical nature of the primary bifurcation. The dashed black line near
$C_L(t)$
indicates that the deviation levels off with an increasing growth rate, highlighting the subcritical nature of the secondary bifurcation.
Isosurfaces of the streamwise vorticity,
$\omega _x R/u_{\textit{rel}} = \pm 0.05$
, past a droplet at selected time instants for
${\textit{Re}}^i=375$
with
${\textit{tu}}_{\textit{rel}}/R = $
(a) 400, (b) 530 and (c) 2000. More details about the time evolution of the vortical structure can be found in supplementary movie 1.
The flow asymmetry with respect to
$y=0$
during the second transition results in a non-zero lift force
$F_L$
along the
$y$
axis, directed from the primary vortex pair towards the shrinking one. We define this direction as the positive
$y$
axis. The resulting lift force can be quantified using a lift coefficient, defined as
$F_L = C_L \pi R^2 \rho ^e u_{\textit{rel}}^2 / 2$
. Figure 10(a) (red line) shows the time evolution of the lift coefficient. The growth of the lift force begins at
${\textit{tu}}_{\textit{rel}}/R \approx 375$
(not shown), well before
$E$
reaches its first plateau. This indicates that the secondary bifurcation sets in before the primary bifurcation fully saturates. However, the initial growth rate of
$C_L$
remains small. Beyond
${\textit{tu}}_{\textit{rel}}/R \approx 450$
,
$C_L$
starts to increase progressively, reaching approximately 0.065 in the fully developed state. Note that this value is close to that induced by the external flow bifurcation in the case of a solid sphere moving at
${\textit{Re}}^e = 250$
(Johnson & Patel Reference Johnson and Patel1999; Shi et al. Reference Shi, Rzehak, Lucas and Magnaudet2021). The evolution of the lift coefficient also provides insight into the nature of the secondary bifurcation (Citro et al. Reference Citro, Tchoufag, Fabre, Giannetti and Luchini2016). Specifically, in figure 10(a), the dashed line near
$C_L(t)$
clearly shows that the initial deviation from linearity is followed by an increasing growth rate of
$C_L(t)$
. Hence, unlike the primary bifurcation, which is supercritical, the secondary bifurcation is subcritical – similar to, for example, the secondary bifurcation observed in the wake of a circular cylinder (Henderson & Barkley Reference Henderson and Barkley1996).
Figure 10(b) shows the values of the total azimuthal energy and the lift coefficient in the fully developed state for
${\textit{Re}}^i$
up to 400. Notably, both quantities exhibit a sharp increase as
${\textit{Re}}^i$
exceeds approximately 370. This further supports the subcritical nature of the secondary bifurcation, which allows for the presence of jumps and hysteresis as the control parameter (here,
${\textit{Re}}^i$
) is varied (Strogatz Reference Strogatz2018, p. 61). This latter aspect is discussed in the next section, where we examine the stability of the dynamical system under finite-amplitude initial disturbances.
4.4. Bistable flow regime
As
${\textit{Re}}^i$
increases further while keeping
$(\mu ^\ast , {\textit{Re}}^e) = (0.5, 200)$
, the primary bifurcation occurs earlier in time. For instance, at
${\textit{Re}}^i = 450$
, the azimuthal energy
$E$
starts to increase at approximately
$150\,R/u_{\textit{rel}}$
(not shown), compared with about
$300\,R/u_{\textit{rel}}$
for
${\textit{Re}}^i = 375$
(see figure 10
a). In contrast, the secondary bifurcation becomes progressively slower and eventually ceases to occur for
${\textit{Re}}^i \gt 468$
. This latter behaviour is highlighted in figure 12(a), which shows the time evolution of
$C_L$
for cases near this transition. For
${\textit{Re}}^i \leqslant 468$
,
$C_L$
initially rises slowly to a plateau before gradually increasing to its final value. However, for larger
${\textit{Re}}^i$
cases,
$C_L$
settles at a plateau that decreases with increasing
${\textit{Re}}^i$
and remains at this level as
${\textit{tu}}_{\textit{rel}}/R \to \infty$
. The
$C_L(t)$
evolutions resemble systems close to a saddle-node bottleneck (Strogatz Reference Strogatz2018, § 4.3), where a saddle-node remnant or ghost induces slow passage. To better illustrate this, figure 12(b) presents the rate of change of the lift coefficient,
${\textrm{d}}C_L(t)/{\textrm{d}} (t u_{\textit{rel}}/R )$
, as a function of
$C_L(t)$
during the interval where
$C_L$
varies slowly over time. A bottleneck causing slow passage is clearly visible at
${\textit{Re}}^i=466$
, shrinking as
${\textit{Re}}^i$
increases. For
${\textit{Re}}^i \gt 468$
, the ‘path’ originating from
$C_L=0$
terminates at a stable fixed point with a small but finite lift coefficient (indicated by filled symbols in figure 12
b).
(a) Time evolution of the lift coefficient for internal Reynolds numbers near the transition to the bistable regime. (b) Variation of
${\textrm{d}}C_L(t)/{\textrm{d}}(t u_{\textit{rel}}/R )$
as a function of
$C_L(t)$
over the time interval where
$C_L(t)$
evolves slowly (results for
$t u_{\textit{rel}}/R \leqslant 300$
are omitted). For all considered
${\textit{Re}}^i$
, the simulation starts from an initially unperturbed state. For
${\textit{Re}}^i=469$
, an additional simulation (labelled as ‘perturbed’) was performed, starting from an initially asymmetric state based on a result from
${\textit{Re}}^i=467$
at
$t u_{\textit{rel}}/R = 675$
(denoted by an open blue circle in both panels). In (b), the two cases with
${\textit{Re}}^i \gt 468$
approach stable fixed points with small but finite
$C_L$
(denoted by filled symbols), whereas in the remaining cases,
$C_L$
after escaping the bottleneck continues to increase with time (as indicated by the arrows).
Figure 13 illustrates the streamwise vorticity structure in the fully developed state for
${\textit{Re}}^i=469$
. Given the small
$C_L$
in this case, the flow structure remains nearly biplanar-symmetric, albeit with a slight up–down asymmetry. Figure 14 presents the azimuthal energy and lift coefficient in the fully developed state for
${\textit{Re}}^i$
increasing from 450 to 500. Due to the bottleneck effects mentioned above, both
$E$
and
$C_L$
undergo an abrupt decrease as
${\textit{Re}}^i$
exceeds approximately 468.
Isosurfaces of the streamwise vorticity,
$\omega _x R/u_{\textit{rel}} = \pm 0.05$
, in the fully developed state for
${\textit{Re}}^i=469$
corresponding to different initial conditions. (a) Simulation starting from an initially axisymmetric state. (b) Simulation starting from a slightly asymmetric state derived from the transient result for
${\textit{Re}}^i=467$
, where
$C_L=0.0102$
(denoted by an open circle in figure 12).
Variation of azimuthal energy
$E$
and lift coefficient
$C_L$
in the fully developed state for
${\textit{Re}}^i$
increasing from 450 to 500. All simulations started with an initially axisymmetric flow.
We are now in a position to examine the response of the flow in the presence of a finite-amplitude initial disturbance. This is necessary because, as mentioned in the last part of § 4.3, the secondary bifurcation is not supercritical, and multiple stable states may coexist depending on the initial conditions. To begin with, we revisit the cases near
${\textit{Re}}^i = 468$
. All cases were initialised from an axisymmetric state where
$C_L = 0$
. To investigate the possible existence of an additional stable state, we performed the simulation with
${\textit{Re}}^i=469$
again starting from a slightly asymmetric initial state with
$C_L=0.0102$
. This ‘initial’ state was derived from the transient result for
${\textit{Re}}^i=467$
at the instant when the system had just passed the bottleneck (denoted by an open symbol in figure 12
b). As seen in figure 12(b) (thick line), even with such a small initial disturbance, the system was able to ‘escape’ the bottleneck. Following this escape,
$C_L$
increased sharply to its final level, nearly identical to the values reached at slightly smaller
${\textit{Re}}^i$
(thick line in figure 12
a). Figure 13(b) shows the structure of the streamwise vorticity in the fully developed state. The flow exhibits a strong up–down asymmetry, significantly different from the flow structure obtained when starting from an initially axisymmetric state (figure13
a), but similar to that shown in figure 11(c).
The test above indicates that for
${\textit{Re}}^i \gt 468$
, the flow may evolve towards at least two distinct asymmetric branches. One, in which the flow remains nearly biplanar-symmetric, appears to be stable only to small disturbances. The other, in which the flow exhibits a uniplanar-symmetric structure, may be triggered by larger amplitude of ‘up–down’ disturbances (hence
$C_L \gt 0$
). Hereafter, these two asymmetric branches are referred to as the biplanar and uniplanar branches, respectively. So far, all cases discussed (except one) started from an initially axisymmetric state, corresponding to a vanishingly small initial disturbance. To examine the response of the flow system to ‘large’ disturbances, we carried out these cases again starting from a uniplanar-symmetric state reached in the fully developed stage of
${\textit{Re}}^i=450$
. For this new ‘initial’ state, the up–down asymmetry of the flow leads to a lift coefficient equal to
$C_L = 0.088$
, providing an initial disturbance that we believe is large enough to trigger the transition to the uniplanar branch.
Azimuthal energy
$E$
(green symbols), lift coefficient
$C_L$
(red symbols) and drag coefficient
$C_D$
(blue symbols) in the fully developed state for
${\textit{Re}}^i$
increasing from 300 to 550 with
$(\mu ^\ast , {\textit{Re}}^e) = (0.5, 200)$
. Circles: simulations starting with an initially axisymmetric flow; crosses: simulations starting with a uniplanar-symmetric flow corresponding to the fully developed state at
${\textit{Re}}^i = 450$
. Vertical dashed lines highlight the critical
${\textit{Re}}^i$
values marking regime transitions. In (a), the two shaded grey regions correspond to the two bistable regimes. In (b), the horizontal dashed line denotes the drag coefficient obtained by enforcing axisymmetry regardless of
${\textit{Re}}^i$
.
Figure 15(a) (cross symbols) presents the total energy and lift coefficient in the fully developed state for
${\textit{Re}}^i$
increasing from 300 to 550. To obtain these results, the simulations were initiated from a uniplanar-symmetric flow corresponding to the fully developed state at
${\textit{Re}}^i = 450$
. The corresponding results obtained from an initially axisymmetric state (open circles) are shown as well. Based on these results, the following picture of the transition sequence emerges. For
${\textit{Re}}^i \lt 326$
, the flow remains axisymmetric. For
${\textit{Re}}^i \in [326, 366]$
and
${\textit{Re}}^i \in [370, 468]$
, axisymmetry breaks down, but the system evolves towards a unique asymmetric branch depending on
${\textit{Re}}^i$
. Specifically, regardless of the amplitude of initial disturbance, the stable asymmetric branch is always biplanar within the first interval and uniplanar within the second. In the narrow range
${\textit{Re}}^i \in (366, 370)$
and for
${\textit{Re}}^i \gt 468$
, both asymmetric branches coexist, and the selected branch depends on the initial disturbance amplitude. In particular, the biplanar branch is stable only to small disturbances, while the uniplanar branch is more likely to emerge when larger disturbances are present. A detailed discussion on the evolution of the lift coefficient in the first bistable regime, where
${\textit{Re}}^i \in (366, 370)$
, is provided in Appendix A. Determining a quantitative threshold for the disturbance required to promote a transition between the two branches is not straightforward. Nevertheless, we believe that in the regime where
${\textit{Re}}^i \gt 468$
, this threshold increases with increasing
${\textit{Re}}^i$
. Indeed, in our restarted simulations, the uniplanar branch is no longer stable when
${\textit{Re}}^i$
exceeds approximately 1500, whereas it remains stable for
${\textit{Re}}^i$
up to 3000 if the simulations start from the fully developed uniplanar state at
${\textit{Re}}^i = 1000$
, where the initial disturbance is larger (
$C_L = 0.16$
at
${\textit{Re}}^i = 1000$
, compared with
$C_L = 0.088$
at
${\textit{Re}}^i = 450$
). Of particular note is that for the uniplanar branch, the wake becomes unsteady when
${\textit{Re}}^i$
exceeds about 1000. In these high-
${\textit{Re}}^i$
cases, the wake oscillates in time – leading to temporal fluctuations in the drag and lift coefficients – while the flow retains its uniplanar symmetry.
Figure 15(b) presents the evolution of the drag coefficient
$C_D$
, defined as
$F_D = C_D \pi R^2 \rho ^e u_{\textit{rel}}^2 / 2$
, in the considered
${\textit{Re}}^i$
range. The drag remains virtually unchanged for
${\textit{Re}}^i$
up to 326, where the flow remains axisymmetric. This behaviour is consistent with previous findings (Feng & Michaelides Reference Feng and Michaelides2001; Edelmann et al. Reference Edelmann, Le Clercq and Noll2017; Shi et al. Reference Shi, Climent and Legendre2024a
), which demonstrated that the drag coefficient
$C_D$
depends only weakly on the internal Reynolds number when the flow is constrained to be axisymmetric. Thus, the horizontal dashed line in figure 15(b), representing
$C_D$
at
${\textit{Re}}^i = 325$
, can be regarded as the reference drag coefficient for the corresponding axisymmetric configuration (hereafter referred to as
$C_D^{axi}$
) at higher
${\textit{Re}}^i$
for
$(\mu ^\ast , {\textit{Re}}^e) = (0.5, 200)$
. Using this reference value, it becomes clear that the transition from the axisymmetric to the asymmetric branch (whether biplanar or uniplanar) is always accompanied by an increase in drag. The magnitude of this increase, denoted as
$\Delta C_D = C_D - C_D^{axi}$
, initially rises smoothly as
${\textit{Re}}^i$
surpasses 326, where the bifurcation is supercritical. As
${\textit{Re}}^i$
increases further,
$\Delta C_D$
shows a jump for the three critical
${\textit{Re}}^i$
values (366, 370 and 468), where the flow system transits between asymmetric branches. Specifically, at
${\textit{Re}}^i = 366$
and 370, the switch from biplanar to uniplanar symmetry leads to an increase in
$\Delta C_D$
of approximately 0.02, whereas at
${\textit{Re}}^i = 468$
, the transition results in a decrease of approximately 0.03. This distinction highlights the fundamental difference between the flow regime transition at the last critical
${\textit{Re}}^i$
and those at the first two.
5. Discussion
5.1. Mechanism of the primary wake instability
We have shown in § 3.2 and further corroborated in § 4.1 that the rotational symmetry of the base flow breaks down when the maximum internal surface vorticity
$\omega _s^i$
exceeds a
${\textit{Re}}^i$
-dependent threshold
$\omega _c^i$
. Since the whole problem is governed by three dimensionless parameters, namely the viscosity ratio
$\mu ^\ast$
, the external Reynolds number
${\textit{Re}}^e$
and the internal Reynolds number
${\textit{Re}}^i$
, the threshold
$\omega _c^i$
may also be interpreted as a critical internal Reynolds number when considering a series of cases with fixed
$\mu ^\ast$
and
${\textit{Re}}^e$
. This scenario has been explored in detail in § 4, where we set
$(\mu ^\ast , {\textit{Re}}^e)=(0.5, 200)$
and determined the response of the flow as
${\textit{Re}}^i$
(hence
$\omega _s^i$
) increased. The axisymmetric base flow was found to become unstable at
${\textit{Re}}^i \approx 326$
, beyond which it first transits to a steady, biplanar-symmetric flow through a supercritical bifurcation. The biplanar flow structure post-bifurcation suggests that the steady mode with azimuthal wavenumber
$m=2$
is responsible for breaking the axisymmetry. In the following, we elaborate in more detail on the mechanism by which, once produced into the base axisymmetric flow, the azimuthal vorticity may lead to its destabilisation.
As internal vorticity plays a key role in the onset of instability, it is relevant to examine its distribution within the droplet in the base flow, particularly near the instability threshold. As seen in figure 6(b), for cases with fixed
$\mu ^\ast$
and
${\textit{Re}}^e$
, the isocontours of the azimuthal vorticity
$\omega _\phi$
tilt towards the front stagnation point as
${\textit{Re}}^i$
increases. By analysing additional results, we observed the same trend when increasing
${\textit{Re}}^e$
(for fixed
$\mu ^\ast$
and
${\textit{Re}}^i$
) or decreasing
$\mu ^\ast$
(for fixed
${\textit{Re}}^e$
and
${\textit{Re}}^i$
) (not shown). These findings, combined with the fact that the maximum internal vorticity
$\omega _s^i$
increases with increasing
${\textit{Re}}^e$
and
${\textit{Re}}^i$
and decreasing
$\mu ^\ast$
, indicate that a higher
$\omega _s^i$
is associated with a more pronounced tilting of the
$\omega _\phi$
isocontours inside the droplet.
Isocontours of the normalised streamwise gradient of the azimuthal vorticity,
$\partial \omega _\phi /\partial x \,(R^2/u_{\textit{rel}})$
, inside the droplet. The three numbers in parentheses at the bottom of each panel correspond to
$(\mu ^\ast , {\textit{Re}}^e, {\textit{Re}}^i)$
. Specifically,
${\textit{Re}}^i$
increases from 50 to 325 from (a) to (c),
${\textit{Re}}^e$
increases from 5 to 25 from (d) to (f) and
$\mu ^\ast$
decreases from 10 to 4 from (g) to (i). In each row, the last panel corresponds to the case closest to the onset of instability. In all panels, the ambient flow is directed from left to right.
Maximum values of (a) the normalised azimuthal vorticity
$\omega _\phi \,(R/u_{\textit{rel}})$
and (b) its streamwise gradient
$\partial \omega _\phi /\partial x \,(R^2/u_{\textit{rel}})$
inside the droplet as a function of the internal Reynolds number for
$(\mu ^\ast , \,{Re}^e)=(0.5, \,200)$
. In both panels, coloured dashed lines denote results from an axisymmetric flow configuration. Circles: simulations starting from an initially axisymmetric flow; crosses: simulations starting from a uniplanar-symmetric flow corresponding to the fully developed state at
${\textit{Re}}^i = 450$
.
The tilting of the internal
$\omega _\phi$
isocontours towards the front stagnation point as
$\omega _s^i$
increases is an insightful observation, as a similar topological change has been reported for the external
$\omega _\phi$
isocontours in the near wake of oblate spheroidal bubbles (Magnaudet & Mougin Reference Magnaudet and Mougin2007; Yang & Prosperetti Reference Yang and Prosperetti2007). In that configuration, the external
$\omega _\phi$
isocontours tend to align nearly perpendicular to the symmetry axis as
${\textit{Re}}^e$
or the aspect ratio
$\unicode{x03C7}$
(the ratio of the major to minor axes) increases. As this trend progresses, the streamwise gradient of
$\omega _\phi$
must become increasingly pronounced for the viscous term
$\nu ^e \partial ^2\omega _\phi /\partial x^2$
to counterbalance the inertial terms. Such a scenario is inherently unstable, leading to the breakdown of axisymmetry in the base flow when
${\textit{Re}}^e$
exceeds a critical value for a given
$\unicode{x03C7}$
(or vice versa).
In analogy to the argument above, the tilting of the internal
$\omega _\phi$
isocontours observed in figure 6(b) should correspond to a gradual increase in the streamwise gradient of
$\omega _\phi$
as
${\textit{Re}}^i$
increases. To confirm this, we present in figure 16(a–c) the isocontours of
$\partial \omega _\phi /\partial x$
for increasing
${\textit{Re}}^i$
. Clearly, the maximum
$\partial \omega _\phi /\partial x$
(normalised by
$u_{\textit{rel}}/R^2$
) near the front of the droplet increases from 4.6 to 35.4 as
${\textit{Re}}^i$
rises from 50 to 325. A similar trend is observed with increasing
${\textit{Re}}^e$
(see figure 16
d–f, where
$(\mu ^\ast , {\textit{Re}}^i) = (0.5, 500)$
) and with decreasing
$\mu ^\ast$
(see figure 16
g–i, where
$({Re}^e, {\textit{Re}}^i) = (200, 500)$
). These results confirm that, as
$\omega _s^i$
increases, both the tilting of the
$\omega _\phi$
isocontours and the streamwise gradient of the azimuthal vorticity within the droplet become more pronounced. Following the argument by Magnaudet & Mougin (Reference Magnaudet and Mougin2007), the internal flow can no longer remain stable if this gradient becomes sufficiently large. Indeed, although not explicitly shown, an inspection of the flow field at the onset of internal bifurcation reveals that the region where disturbances initially grow is closely aligned with the axial and meridional coordinates of the maximum
$\partial \omega _\phi /\partial x$
in the corresponding axisymmetric configuration (see, for instance, the vortical threads near the front part of the droplet surface in the first inset of figure 1
b). This serves as indirect evidence supporting the proposed relation between internal flow instability and the maximum streamwise gradient of the azimuthal vorticity.
While the relationship between internal flow instability and the presence of a sufficiently large
$\partial \omega _\phi /\partial x$
proposed here is specific to axisymmetric flows (relevant to the first instability), we suspect that
$\partial \omega _\phi /\partial x$
may still play a role in subsequent flow bifurcations, where the flow is already three-dimensional. Figure 17 shows the evolution of the maxima of both
$\omega _\phi$
and
$\partial \omega _\phi /\partial x$
inside the droplet as a function of
${\textit{Re}}^i$
for the series of cases with
$(\mu ^\ast ,\,{Re}^e)=(0.5,\,200)$
. Clearly,
$\partial \omega _\phi /\partial x$
continues to increase significantly with
${\textit{Re}}^i$
throughout the considered range, even when the flow is imposed as axisymmetric. In contrast, the maximum azimuthal vorticity remains nearly constant once the flow first transits from the biplanar to the uniplanar branch. These different behaviours are likely related to variations in
$C_D$
as the flow transits at the two larger critical
${\textit{Re}}^i$
values (see discussion in the last paragraph of the previous section) and highlight the consistent role of
$\partial \omega _\phi /\partial x$
in triggering subsequent flow instabilities. Of course, once the flow becomes three-dimensional, both the streamwise and polar vorticity components become non-zero, and for a comprehensive understanding of flow instability, it may not be sufficient to examine only the azimuthal vorticity component. A rigorous stability analysis of the internal flow field, such as those conducted by Yang & Prosperetti (Reference Yang and Prosperetti2007) and Tchoufag et al. (Reference Tchoufag, Magnaudet and Fabre2013), would be required to fully elucidate the mechanisms at play. Such an analysis would also help verify whether, at the first bifurcation, it is indeed the steady azimuthal mode with wavenumber
$m=2$
that is amplified, leading to a non-axisymmetric but still steady flow. Conducting such an analysis, though beyond the scope of the present study, would be a valuable endeavour for future research.
5.2. Relation between wake and path instabilities of a freely moving droplet of low-to-moderate
$\mu ^\ast$
In the previous sections, the droplet was considered to be spherical, with its centroid assumed to remain fixed. However, under physically realistic conditions, droplets can deform while moving freely under the effect of buoyancy and gravity. A key question is whether the internal flow bifurcation persists in this more general configuration and whether it plays a role in triggering the first path instability.
To provide insights into this question, we carry out three-dimensional time-dependent simulations of a buoyant, deformable droplet rising freely in an immiscible liquid that is otherwise at rest. The physical properties of the liquid–liquid system are selected to match those of a toluene droplet rising in water, as investigated experimentally by Wegener et al. (Reference Wegener, Kraume and Paschedag2010). The droplet radius is set to
$R=1.2\,\text{mm}$
, which is slightly larger than the threshold reported in the experiment. Unlike all simulations discussed so far, which were performed using the JADIM code, the simulations discussed below are carried out using the open-source flow solver Basilisk (Popinet Reference Popinet2009, Reference Popinet2015). This choice allows us to take advantage of the adaptive mesh refinement technique built into the solver (van Hooft et al. Reference van Hooft, Popinet, van Heerwaarden, van der Linden, de Roode and van de Wiel2018) and to maintain a resolution at the interface comparable to that provided by the boundary-fitted mesh used earlier.
The one-fluid approach, combined with the geometric volume-of-fluid method, is used to track and evolve the liquid–liquid interface. At the interface, the local density (dynamic viscosity) of the fluid medium is approximated using the arithmetic (harmonic) averaging rule. Note that the harmonic model for the dynamic viscosity does not guarantee a continuous variation of the shear stresses across the interface (Kothe Reference Kothe1998; Magnaudet et al. Reference Magnaudet, Bruhier, Mer and Bonometti2025), and a sufficiently high resolution is therefore required to eliminate this numerical artefact. To model surface tension, the balanced-force surface-tension formulation is used (Francois et al. Reference Francois, Cummins, Dendy, Kothe, Sicilian and Williams2006), which is based on the continuum surface force approach originally proposed by Brackbill, Kothe & Zemach (Reference Brackbill, Kothe and Zemach1992). In addition, a second-order-accurate calculation of curvature is performed using the height-function technique developed by Popinet (Reference Popinet2009). Details of the numerical schemes can be found in Popinet (Reference Popinet2018). The computational domain is a cubic box with an edge length of
$240R$
. A free-slip condition is imposed on all four lateral boundaries, while a periodic condition is applied to the top and bottom boundaries (Zhang, Ni & Magnaudet Reference Zhang, Ni and Magnaudet2021). The droplet is initially spherical and is released from rest at a position midway between the four lateral boundaries and
$15R$
above the bottom of the simulation domain. The spatial resolution is refined to approximately
$1/68 R$
near the interface and to about
$1/17R$
in the far wake (starting approximately
$10R$
downstream of the droplet). Illustration of the grid structure in the vicinity of the droplet as well as in the far wake may be found in figure 18 of Shi et al. (Reference Shi, Zhang and Magnaudet2024b
). The accuracy provided by the above grid resolution is confirmed through a grid-independence study detailed in Appendix A of Shi et al. (Reference Shi, Zhang and Magnaudet2025c
), where the path oscillations of a bubble rising near a wall at Reynolds numbers (both
${\textit{Re}}^e$
and
${\textit{Re}}^i$
) up to approximately 1000 have been considered. To ensure that the resolution remains sufficient for the present case of a freely rising droplet, a corresponding grid study in a two-dimensional axisymmetric configuration has been conducted. There, reducing the minimum grid size from approximately
$1/68\,R$
to
$1/136\,R$
produces negligible changes in both the time evolution and the terminal value of the droplet rising velocity, for toluene droplets with radii ranging from
$0.5$
to
$2.1\,\text{mm}$
. This confirms the reliability of the numerical approach for the current problem.
Time evolution of (a) the vertical velocity (
$V_v$
, red line, left axis) and droplet aspect ratio (
$\unicode{x03C7}$
, green line, right axis) as well as (b) the horizontal velocity (
$V_h$
) for a single toluene droplet of radius
$R=1.2\,\text{mm}$
rising in quiescent water (for detailed physical parameters, see table 2 of Wegener et al. (Reference Wegener, Kraume and Paschedag2010)). Both
$V_v$
and
$V_h$
are normalised by
$(gR)^{1/2}$
. In both panels, solid lines represent the present simulation results. In (a), red open symbols denote experimental data of
$V_v$
from Wegener et al. (Reference Wegener, Kraume and Paschedag2010) for
$t/(R/g)^{1/2}$
up to 300, beyond which wall effects in the experiment significantly influenced the rising speed. The insets display the isosurfaces of the vertical component of the vorticity,
$\omega _v (R/g)^{1/2}=\pm 0.5$
(grey and black threads denote positive and negative values, respectively), at selected time instants (indicated at the top of each panel; values normalised by
$(R/g)^{1/2}$
). In all insets, the gravitational acceleration points vertically downwards, such that the droplet initially rises vertically upwards until approximately
$t/(R/g)^{1/2} = 200$
.
Figure 18(a) compares the time evolution of the vertical velocity of the droplet,
$V_v$
, obtained from our three-dimensional simulation (solid line in red) with the corresponding experimental data (red symbols) from Wegener et al. (Reference Wegener, Kraume and Paschedag2010). In both cases, the rising speed initially increases to a maximum at
$t/(R/g)^{1/2} \approx 75$
before decreasing and oscillating around a mean value of approximately
$1.0(gR)^{1/2}$
for
$t/(R/g)^{1/2} \leqslant 250$
. The first maximum of the rising speed and the reduced frequency (or Strouhal number,
$St = 2fR/V_m$
, where
$V_m$
is the mean rising speed over the first two oscillation cycles and
$f$
is the corresponding mean oscillation frequency) are
$V_v^{max } = 1.34(gR)^{1/2}$
and
$St = 0.043$
in the simulation, closely consistent with the experimental values of
$V_v^{max } = 1.39(gR)^{1/2}$
and
$St = 0.037$
. The slight difference in the oscillation frequency is likely due to confinement effects in the experiment. There, the droplet initially rises along the axis of a cylindrical domain of radius approximately
$21R$
(Wegener Reference Wegener2009, p. 89), compared with
$120R$
in the simulation. As a result of this difference, the subsequent rising behaviour of the droplet in the experiment differs significantly from that observed in the simulation. In the experiment, the droplet begins to migrate laterally at
$t/(R/g)^{1/2} \approx 225$
and collides with the wall at
$t/(R/g)^{1/2} \approx 360$
(Wegener Reference Wegener2009; see figure 5.4 therein) (experimental data near the collision are omitted from figure 18
a). As expected, the wall imposes a retarding effect on the rise speed, which causes a damping of velocity oscillations (Magnaudet, Takagi & Legendre Reference Magnaudet, Takagi and Legendre2003; Zeng, Balachandar & Fischer Reference Zeng, Balachandar and Fischer2005; Shi Reference Shi2024). In contrast, the lateral migration in the simulation begins earlier, at
$t/(R/g)^{1/2} \approx 200$
, as seen in figure 18(b). Beyond this point, the horizontal velocity,
$V_h$
, grows in time while the rising velocity undergoes damped oscillations. Both velocity components stabilise beyond
$t/(R/g)^{1/2} \approx 600$
, with the droplet rising steadily along an oblique path at an angle of approximately
$\tan ^{-1} (V_h/V_v) \approx 15^\circ$
. Notably, in this terminal state, the external Reynolds number is approximately
${\textit{Re}}^e = 256$
, and the force balance among drag, lift, buoyancy and gravity yields drag and lift coefficients of
$C_D = 0.374$
and
$C_L = 0.10$
, respectively. These coefficients agree well with their counterparts in the corresponding spherical fixed-drop configuration, where we obtain
$C_D = 0.372$
and
$C_L = 0.093$
.
From the fully resolved three-dimensional simulation results, it is also possible to examine the time evolution of the droplet shape. Assuming that the droplet remains almost oblate spheroidal during its rise, its shape can be characterised using the aspect ratio
$\unicode{x03C7}$
, which represents the ratio of major and minor axes. The evolution of
$\unicode{x03C7}(t)$
according to our numerical results is shown in figure 18(a) (green line and right vertical axis). Throughout the evolution, the level of deformation remains marginal. Specifically, the maximum aspect ratio is approximately 1.16 and is reached shortly after the onset of the internal bifurcation. Thereafter,
$\unicode{x03C7}$
oscillates around 1.1 before stabilising at 1.08. Prior to stabilisation, the radian frequency of the shape oscillations is about
$0.132(g/R)^{1/2}$
. On the other hand, for the same liquid–liquid system, the radian frequency associated with the fundamental shape mode (i.e. the mode associated with the polar wavenumber
$l = 2$
) of a nearly spherical droplet is (Lamb Reference Lamb1932; Lalanne et al. Reference Lalanne, Tanguy and Risso2013)
\begin{equation} \omega _s(l) = \sqrt {\frac {(l-1)l(l+1)(l+2)\gamma }{[\rho ^\ast (l+1) + l]\rho ^e R^3}} \approx 6.6(g/R)^{1/2}, \end{equation}
which is about 50 times higher than the oscillation frequency observed in the aspect ratio. Hence, deformation instability is unlikely to be responsible for the observed oscillations. Rather, we suspect that the shape oscillations observed during the transient stage are primarily driven by variations in the rise velocity. This is consistent with the small amplitude of droplet deformation, which scales linearly with the Weber number and thus primarily with the square of the rising speed. As a result,
$\unicode{x03C7}$
reaches a local maximum shortly after
$V_v$
peaks, and the same correspondence holds for the local minima of
$\unicode{x03C7}$
and
$V_v$
.
We now elaborate, based on the numerical results from Basilisk, the relationship between the internal flow bifurcation and the evolution of the droplet motion. By examining the time evolution of the internal and external azimuthal energies (not shown), we found that the axisymmetry of the flow breaks down due to an internal flow bifurcation at
$t/(R/g)^{1/2} \approx 60$
, i.e. before the abrupt decrease in
$V_v(t)$
. This feature is highlighted by the insets in figure 18(a), which illustrate the evolution of the streamwise vorticity during the interval
$60 \leqslant t/(R/g)^{1/2} \leqslant 70$
. Note that since the initial configuration is axisymmetric, the transition to a non-axisymmetric state is triggered by the amplification of random numerical disturbances (Zhang et al. Reference Zhang, Ni and Magnaudet2021). Similar to what was observed in figure 1(b), the disturbance associated with the asymmetry initially develops and grows only inside the droplet. Moreover, this disturbance is constrained within four vortex threads of equal intensity, maintaining a biplanar-symmetric structure during the transition. The insets on the left-hand side of figure 18(b) show the vortical structure during the time oscillations of
$V_v$
. These results indicate that the velocity oscillations are closely related to the unsteady development of the wake. Specifically,
$V_v$
reaches a local maximum at, for example,
$t/(R/g)^{1/2} = 125$
and
$172$
, where the vortices shrink to their minimum extent, and the reverse occurs when
$V_v$
reaches a local minimum, such as at
$t/(R/g)^{1/2} = 96$
and
$143$
. The flow transits from biplanar-symmetric to uniplanar-symmetric at
$t/(R/g)^{1/2} \approx 200$
. As seen in the three rightmost insets of figure 18(b), during this secondary transition, the left–right asymmetry grows over time, similar to the up–down asymmetry observed for a fixed droplet in figures 11 and 13(b). This asymmetry leads to a lift force that drives the lateral migration of the droplet.
Characteristic parameters obtained in an axisymmetric configuration for toluene droplets of sizes close to the threshold of the first path instability. Variation of the (a) internal and (b) external Reynolds number (horizontal axis) as a function of droplet radius
$R$
(vertical axis). Maximum (c) internal and (d) external) surface vorticity (red line) as a function of
${\textit{Re}}^i$
and
${\textit{Re}}^e$
, respectively. In (a), the tick labels at the bottom match those in (c), and the same correspondence holds between (b) and (d). In (c), the black solid line represents the criterion for internal flow bifurcation (3.3), while in (d), the black solid line corresponds to the criterion for external flow bifurcation from Magnaudet & Mougin (Reference Magnaudet and Mougin2007).
Based on the confirmed relationship between the internal flow bifurcation and the first path instability, we now assess the applicability of the empirical criterion proposed in § 3.2 (i.e. (3.3)) for predicting the threshold droplet size at which path instability occurs first. To apply criterion (3.3), data on the maximum internal vorticity,
$\omega _s^i$
, and the internal Reynolds number,
${\textit{Re}}^i$
– both of which depend on the droplet radius,
$R$
– are required. To obtain these data, we conducted an additional series of simulations using Basilisk, considering the same liquid–liquid system while imposing axisymmetry on the flow field. The droplet radius was varied from
$0.5$
to
$1.5\,\text{mm}$
in increments of
$0.05\,\text{mm}$
. Figure 19(a) shows the variation of
${\textit{Re}}^i$
(horizontal axis) as a function of the droplet radius (vertical axis) for
$R$
near the path instability threshold. The corresponding results for
$\omega _s^i$
are shown in figure 19(c) (red line) as a function of
${\textit{Re}}^i$
. The intersection of the curve
$\omega _s^i({Re}^i)$
with the criterion (3.3) (black solid line in figure 19
c) corresponds to a critical internal Reynolds number of 330, which translates to a critical radius of
$R=0.97\,\text{mm}$
, as indicated in figure 19(a). This closely matches the experimentally observed threshold of
$R\approx 1.1\,\text{mm}$
(Wegener et al. Reference Wegener, Kraume and Paschedag2010). Furthermore, we recall that the critical radius for the onset of path instability was found to be
$R\approx 1.0\,\text{mm}$
in recent DNS by Charin et al. (Reference Charin, Lage, Silva, Tuković and Jasak2019).
The numerical results for
$\omega _s^e(R)$
and
${\textit{Re}}^e(R)$
, also obtained from these axisymmetric simulations, allow us to examine whether the external flow remains stable. Figure 19(b) presents the external Reynolds number as a function of
$R$
near the threshold of the internal flow bifurcation, with the critical radius corresponding to
${\textit{Re}}^e=236$
. Figure 19(d) compares the maximum external surface vorticity (red line) with the threshold predicted by Magnaudet & Mougin (Reference Magnaudet and Mougin2007) for external flow bifurcation (black solid line). At
${\textit{Re}}^e=236$
, the normalised maximum external surface vorticity,
$\omega _s^e/(u_{\textit{rel}}/R)$
, is approximately
$4.7$
, which is only one-third of the threshold value (
$13.5$
). Hence, the external flow remains stable at the critical droplet size of path instability.
5.3. Threshold droplet radius for internal bifurcation of a nearly spherical droplet moving in water
Based on the key findings above, reference values for the critical droplet size required for the onset of internal flow bifurcation in a fluid–fluid system can be estimated. Specifically, for an immiscible liquid droplet freely moving in water, the critical radius, denoted as
$R_c$
, beyond which internal bifurcation may occur, can be determined by evaluating the maximum internal surface vorticity at steady state, where buoyancy, gravity and drag forces are in equilibrium. A relatively coarse estimate of
$R_c$
can also be obtained by considering the internal and external Reynolds numbers at steady state. Specifically, our fixed-droplet simulations indicate that internal bifurcation occurs typically for
${\textit{Re}}^i \gtrsim 350$
, provided that
${\textit{Re}}^e = {\mathcal{O}}(100)$
. Using these observations, we determined
$R_c$
for various viscosity and density ratios (for details, see Appendix B), with both
$\mu ^\ast$
and
$\rho ^\ast$
ranging from 0.1 to 10, a parameter range commonly encountered for real liquid–liquid systems (Balla et al. Reference Balla, Kavuri, Tripathi, Sahu and Govindarajan2020). These results are summarised in figure 20(a), which applies to droplets with small deformation, say, for
$\unicode{x03C7}\leqslant 1.1$
. Together with the prerequisite
${\textit{Re}}^e = {\mathcal{O}}(100)$
, the range of validity may be further interpreted as an upper bound on the Morton number,
$Mo$
, up to approximately
$3\times 10^{-9}$
, based on the empirical correlation by Myint, Hosokawa & Tomiyama (Reference Myint, Hosokawa and Tomiyama2007; see (14) therein), where
$Mo = g(\mu ^e)^4|1-\rho ^\ast |/(\rho ^e\gamma ^3)$
with
$\gamma$
denoting the interfacial tension. Notably, within the considered range of
$\mu ^\ast$
and
$\rho ^\ast$
, the critical
$R_c$
varies from approximately
$0.2$
to
$2\,\text{mm}$
, which is well within the typical size range of most practical systems (Clift et al. Reference Clift, Grace and Weber2005). Figure 20(b) provides a zoomed-in view for
$\mu ^\ast , \rho ^\ast \in [0.5,2]$
. Within this refined parameter range, the minimum
$R_c$
at a given
$\mu ^\ast$
is generally larger for a light droplet (
$\rho ^\ast \lt 1$
) than for a heavy one (
$\rho ^\ast \gt 1$
). For example, at
$\mu ^\ast = 1$
, the minimum
$R_c$
is approximately
$0.95\,\text{mm}$
for a light droplet (attained at
$\rho ^\ast = 0.5$
), while it is only about
$0.5\,\text{mm}$
for a heavy droplet (attained at
$\rho ^\ast = 2.0$
). Furthermore, for a toluene droplet rising in water, a case with
$(\mu ^\ast , \rho ^\ast ) = (0.62, 0.86)$
, figure 20(b) indicates a critical radius
$R_c$
of approximately
$1.05\,\text{mm}$
, which closely agrees with the experimental threshold reported by Wegener et al. (Reference Wegener, Kraume and Paschedag2010).
Threshold droplet radius
$R_c$
(in
$\text{mm}$
) in the
$(\mu ^\ast ,\,\rho ^\ast )$
phase plane for internal bifurcation in the case of a nearly spherical droplet rising (dashed lines) or settling (solid lines) in water. (a) Results for
$\mu ^\ast$
and
$\rho ^\ast$
varying from 0.1 to 10. (b) Same as (a) but for
$\mu ^\ast , \rho ^\ast \in [0.5,2]$
. In both panels, coloured lines denote iso-
$R_c$
contours, with values (in
$\text{mm}$
) indicated in the figure. The filled circle, located at
$(\mu ^\ast , \rho ^\ast ) = (0.62, 0.86)$
, corresponds to the case of toluene droplets in water, for which the critical
$R_c$
is approximately
$1.05\,\text{mm}$
, as determined from the present regime map. In (a), the dotted line in black corresponds to
${\textit{Re}}^e = 100$
(or equivalently
${\textit{Re}}^i = 350$
), indicating that to the right (left) of this line, a sufficiently large
${\textit{Re}}^i$
(
${\textit{Re}}^e$
) is required for internal flow bifurcation to occur. For details on the constraint of internal bifurcation in terms of
${\textit{Re}}^e$
and
${\textit{Re}}^i$
, see (B1) in Appendix B.
6. Summary
We carried out three-dimensional numerical simulations of a uniform flow past a fixed spherical droplet over a wide range of governing parameters, namely the viscosity ratio
$\mu ^\ast$
, the external Reynolds number
${\textit{Re}}^e$
and the internal Reynolds number
${\textit{Re}}^i$
. Our results show that for droplets with low-to-moderate viscosity ratios, the axisymmetric wake becomes unstable because of an internal flow bifurcation. This behaviour is absent for bubbles, particles and droplets with large viscosity ratios, where the internal flow does not influence wake instability. The internal flow bifurcation is linked with the surface vorticity produced at the external side of the droplet interface. By varying
$\mu ^\ast$
,
${\textit{Re}}^e$
and
${\textit{Re}}^i$
independently, we found that the critical condition for the onset of internal bifurcation can be characterised in terms of the maximum vorticity on the internal side of the droplet surface,
$\omega _s^i$
. This leads to an empirical criterion based on a threshold
$\omega _s^i$
, denoted as
$\omega _c^i$
, which was found to depend solely on
${\textit{Re}}^i$
, to determine whether the axisymmetric flow remains stable.
Then, we selected a particular series of cases where
$(\mu ^\ast , {\textit{Re}}^e)=(0.5, 200)$
to study the flow evolution with an increasing internal Reynolds number. Starting from an initially axisymmetric velocity field, the flow first undergoes a supercritical bifurcation at
${\textit{Re}}^i=326$
, yielding a steady non-axisymmetric flow that retains biplanar symmetry. In this configuration, the wake consists of two pairs of counter-rotating vortex threads of equal intensity, suggesting that symmetry breaking is associated with an azimuthal mode with wavenumber
$m=2$
(Ghidersa & Dušek Reference Ghidersa and Dušek2000; Yang & Prosperetti Reference Yang and Prosperetti2007). With an additional increase in
${\textit{Re}}^i$
beyond 370, a secondary bifurcation occurs, which is found to be subcritical. Following this transition, the flow loses its biplanar symmetry and becomes uniplanar-symmetric, characterised by a single pair of counter-rotating vortices in the wake. Consequently, the droplet experiences a sizeable lift force, with
$C_L$
showing an abrupt increase from a vanishingly small value to approximately 0.06 as
${\textit{Re}}^i$
reaches 370. The secondary bifurcation persists up to
${\textit{Re}}^i=468$
, where
$C_L$
further increases from 0.06 to about 0.09. For
${\textit{Re}}^i \gt 468$
, the flow reverts to a biplanar-symmetric configuration. However, due to the subcritical nature of the secondary bifurcation, the final-state flow also depends on the initial disturbance amplitude. Specifically, by restarting all simulations from an asymmetric initial condition using the fully developed flow at
${\textit{Re}}^i=450$
(where
$C_L=0.088$
), we identified two bistable regimes: a narrow range where
${\textit{Re}}^i\in (366,370)$
and another for
${\textit{Re}}^i\gt 468$
. In both bistable regimes, the final-state flow transits from biplanar- to uniplanar-symmetric when the initial disturbance exceeds a certain threshold.
Based on these findings, we proposed a physical explanation for the mechanism driving the primary wake instability. Examination of the azimuthal vorticity field in the base flow close to the threshold revealed that the isocontours inside the droplet tilt significantly towards the front, particularly when
$\omega _s^i$
is large. This tilting is accompanied by a marked increase in the streamwise gradient of the internal azimuthal vorticity. Drawing an analogy with the argument proposed by Magnaudet & Mougin (Reference Magnaudet and Mougin2007) for an external flow bifurcation, we suggested that if this streamwise gradient becomes sufficiently large, the internal flow cannot remain stable, leading to axisymmetry breakdown. Although this criterion probably provides only a sufficient condition for the primary wake instability, it aligns quantitatively with our numerical observations. A detailed stability analysis of the base flow in this regime is of course required to confirm the above scenario and obtain a more accurate criterion.
Finally, we examine the relationship between the primary wake instability observed for a fixed droplet and the first path instability when the droplet is free to move. To this end, we conducted additional simulations of freely rising droplets, selecting physical parameters corresponding to those of a toluene droplet rising in quiescent water (Wegener et al. Reference Wegener, Kraume and Paschedag2010). Results from the three-dimensional simulations of a droplet with a radius of
$R=1.2\,\text{mm}$
confirmed the onset of an internal flow bifurcation prior to the emergence of the first path instability, thereby establishing a direct connection between wake and path instabilities. Using the data for
$(\omega _s^i, {\textit{Re}}^i)$
obtained from constrained axisymmetric simulations over a wide range of droplet radii, along with the empirical criterion
$\omega _c^i({Re}^i)$
derived from the fixed-droplet simulations, we found that the predicted threshold droplet size for the primary wake instability closely matches the experimental and numerical thresholds for the onset of the first path instability (Wegener et al. Reference Wegener, Kraume and Paschedag2010; Charin et al. Reference Charin, Lage, Silva, Tuković and Jasak2019). This further validates the proposed criterion for the internal flow bifurcation. Building on this confirmed relationship, we also estimated the threshold droplet size for the internal bifurcation of a nearly spherical droplet moving freely in water, using the criterion proposed in the present work.
One key aspect not addressed in this study is the mathematical nature of the secondary bifurcation that drives the transition from biplanar to uniplanar symmetry. Understanding this bifurcation is particularly important, as it leads to a lift force that, for a freely moving droplet, causes the transition from a vertical to an oblique path. Addressing this issue requires the development of a suitable global linear stability approach. Recent numerical techniques have made it possible to determine the threshold and nature of bifurcations in freely rising bubbles (Bonnefis et al. Reference Bonnefis, Sierra-Ausin, Fabre and Magnaudet2024). Extending this approach to systems where the internal and external flow fields are coupled through kinematic and dynamic boundary conditions appears to be a promising next step to gain deeper insight into this fundamental problem. In this context, a systematic azimuthal mode decomposition of the DNS results, which we have deliberately omitted here, would also help to rigorously characterise the primary and secondary bifurcation modes, and will be pursued in future work.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2025.10565.
Acknowledgements
P.S. gratefully acknowledges many fruitful discussions with J. Magnaudet, whose valuable insights have greatly contributed to understanding the possible mechanism of the primary wake instability. We also extend our gratitude to J. Zhang at Xi’an Jiaotong University for his generous support in implementing the periodic boundary condition in Basilisk. The computations were carried out on the HPC cluster hemera at HZDR.
Funding
This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) (P.S.; grant number 501298479).
Declaration of interests
The authors report no conflict of interest.
Appendix A. The first bistable regime encountered with increasing
${\boldsymbol{Re}}^{\boldsymbol{i}}$
for
$\boldsymbol{(\mu} ^{\boldsymbol{\ast}} , {\boldsymbol{Re}}^{\boldsymbol{e}}\boldsymbol{)=(0.5,200)}$
For the transition sequence discussed in § 4, a bistable regime exists within the narrow interval
${\textit{Re}}^i\in (366,370)$
, where the final-state flow structure can be either biplanar-symmetric or uniplanar-symmetric, depending on the amplitude of the initial disturbance.
Lift coefficient obtained from different initial conditions for
${\textit{Re}}^i$
increasing from 360 to 375. Circles correspond to cases initialised from an axisymmetric flow, while crosses denote cases starting from an initially asymmetric flow based on the final-state result for
${\textit{Re}}^i=450$
(where
$C_L=0.088$
).
Same as figure 12 but for
${\textit{Re}}^i$
increasing from 367 to 375. (a) The time evolution of the lift coefficient
$C_L$
for cases starting from an initially axisymmetric velocity field. (b) The variation of
${\textrm{d}}C_L(t)/{\textrm{d}}(t u_{\textit{rel}}/R )$
as a function of
$C_L(t)$
, illustrating the emergence of a local fixed point near
$C_L \approx 0.01$
for
${\textit{Re}}^i = 369$
.
Figure 21 shows the lift coefficient obtained from different initial states for
${\textit{Re}}^i$
increasing from 360 to 375. For cases starting from an axisymmetric flow, the final-state flow remains weakly biplanar-symmetric up to
${\textit{Re}}^i=369$
, where
$C_L$
is merely equal to 0.01. As
${\textit{Re}}^i$
increases further,
$C_L$
undergoes an abrupt rise to approximately 0.06 at
${\textit{Re}}^i=370$
, beyond which it shows a weak increase. Now, the above cases were carried out again but starting from an initially asymmetric velocity field. Using the final-state result from
${\textit{Re}}^i=450$
(where
$C_L=0.088$
), we find that the abrupt increase starts at
${\textit{Re}}^i\approx 367$
, slightly lower than in the previous scenario. However, for both
${\textit{Re}}^i\leqslant 366$
and
${\textit{Re}}^i\geqslant 370$
, the resulting
$C_L$
in the final state (and hence the corresponding flow structure) is independent of the initial conditions.
Figures 22(a) and 22(b) show the time evolution of
$C_L$
and its rate of change (as a function of
$C_L$
), respectively, obtained from the series of simulations starting from an axisymmetric flow. The evolution of
$C_L(t)$
highlights a bottleneck effect as
${\textit{Re}}^i$
decreases to 370, similar to that observed in § 4.4 for
${\textit{Re}}^i$
increasing beyond 465 (see figure 12). As
${\textit{Re}}^i$
decreases slightly to 369, a local fixed point with small
$C_L$
(approximately 0.01) emerges (figure 22
b). According to figure 22(b), this local fixed point shifts rapidly towards
$C_L=0$
as
${\textit{Re}}^i$
is decreased further. For
${\textit{Re}}^i\leqslant 366$
, this local fixed point becomes globally stable, meaning that it cannot be eliminated even if the simulations are initialised from a highly asymmetric flow corresponding to the final stage of
${\textit{Re}}^i=450$
(where
$C_L=0.088$
).
Appendix B. Regime map of internal bifurcation of a nearly spherical freely rising or settling droplet in water
The discussion in § 5.2 confirmed the close relationship between internal bifurcation and the first path instability of a freely rising droplet with a low-to-moderate viscosity ratio. Moreover, the criterion (3.3), based on the maximum internal surface vorticity, was found to predict reasonably well the threshold droplet size reported in experiments for a toluene droplet rising in water. This agreement motivates us to construct a regime map for internal bifurcation in a general liquid–liquid system. To narrow the scope, we focus on the case of a freely rising or settling droplet in water under gravity. The key question we seek to address is: given the viscosity and density of the droplet (hence
$\mu ^\ast$
and
$\rho ^\ast$
) a priori, what is the threshold droplet radius corresponding to the onset of internal bifurcation?
To answer this question, one could first establish an empirical correlation to estimate the maximum internal surface vorticity,
$\omega _s^i$
, in the parameter space
$(\mu ^\ast , {\textit{Re}}^e, {\textit{Re}}^i)$
. The corresponding external and internal Reynolds numbers,
$({Re}^e, {\textit{Re}}^i)$
, could then be determined for a droplet moving in water with radius
$R$
with given viscosity and density ratios
$(\mu ^\ast , \rho ^\ast )$
. By equating
$\omega _s^i$
, now expressed as a function of
$(R, \mu ^\ast , \rho ^\ast )$
, with the critical vorticity threshold
$\omega _c^i$
– which itself depends on
$(R, \mu ^\ast , \rho ^\ast )$
through
${\textit{Re}}^e$
and
${\textit{Re}}^i$
– one would obtain the critical droplet radius,
$R_c$
. While this approach is essential for practical applications, it requires substantial effort to derive a reliable empirical correlation for
$\omega _s^i$
over the three-parameter space
$(\mu ^\ast , {\textit{Re}}^e, {\textit{Re}}^i)$
. Given this complexity, we opt instead for a simplified estimate of
$R_c$
based on the fixed-droplet DNS results available from the present study.
We begin by examining the dependence of the critical internal surface vorticity,
$\omega _c^i$
, on the internal Reynolds number,
${\textit{Re}}^i$
. Inspection of (3.3) reveals that the most rapid variation occurs for
${\textit{Re}}^i \lesssim 300$
. At
${\textit{Re}}^i = 350$
,
$\omega _c^i$
decreases to approximately
$5.5\,u_{\textit{rel}}/R$
, which is generally lower than the resulting
$\omega _s^i$
for droplets with
$\mu ^\ast = {\mathcal{O}}(0.1{-}1)$
moving at
${\textit{Re}}^e = {\mathcal{O}}(100)$
. Thus, as a rule of thumb, internal bifurcation is likely to occur when
${\textit{Re}}^i$
exceeds approximately 350. Given this, the next step is to determine the critical droplet radius,
$R_c$
, required for
${\textit{Re}}^i$
to exceed this threshold. More specifically, since
${\textit{Re}}^i = {\textit{Re}}^e \rho ^\ast /\mu ^\ast$
, the problem reduces to finding
$R_c$
such that
\begin{equation} \left \{ \begin{aligned} &{Re}^e \rho ^\ast /\mu ^\ast \approx 350 & & \text{for} \quad {\textit{Re}}^e \gtrsim 100, \\ & & &\text{or} \\ &{Re}^e \rho ^\ast /\mu ^\ast \gtrsim 350 & & \text{for} \quad {\textit{Re}}^e \approx 100. \end{aligned} \right . \end{equation}
We now examine the relationship between
${\textit{Re}}^e$
and the droplet radius
$R$
. Assuming that the flow remains axisymmetric in the terminal state and that no shape oscillations occur, the balance among the drag, gravity and buoyancy forces reads
\begin{equation} \rho ^e |\rho ^\ast - 1| \frac {4}{3} \pi R^3 g = \frac {1}{2} \pi R^2 \rho ^e C_D\, u_{\textit{rel}}^2. \end{equation}
Rearranging this expression, we obtain
\begin{equation} R = \sqrt [3]{\frac {3}{32} \, \frac {\left ( \mu ^e \right )^2}{\left ( \rho ^e \right )^2 |\rho ^\ast - 1| g} \, C_D \left ( {\textit{Re}}^e \right )^2}. \end{equation}
To employ (B3), an appropriate correlation for the drag coefficient,
$C_D$
, is required, as it generally depends on
$(\mu ^\ast , \rho ^\ast , {\textit{Re}}^e)$
. Assuming weak deformation, a reliable correlation for
$C_D$
can be found in our recent work (Shi et al. Reference Shi, Climent and Legendre2024a
), where the hydrodynamic force on spherical droplets was examined. It was shown that, in the absence of internal bifurcation,
$C_D$
depends only weakly on
$\rho ^\ast$
and can be approximated as
\begin{equation} C_D = C_D^B + (R_\mu ^m - 1) \frac {C_D^S - C_D^B}{(3/2)^m - 1}, \end{equation}
where
$R_\mu = (2 + 3\mu ^\ast )/(2 + 2\mu ^\ast )$
represents the intensity of the Stokeslet and
$C_D^B$
and
$C_D^S$
denote the drag coefficients in the clean-bubble (
$\mu ^\ast \to 0$
) and solid-sphere (
$\mu ^\ast \to \infty$
) limits, respectively. The exponent
$m$
is a fitted function of
${\textit{Re}}^e$
and, along with
$C_D^B$
and
$C_D^S$
, takes the following expressions:
\begin{align} C_D^B &= \frac {16}{{Re}^e} \left \{ 1 + \left [ \frac {8}{{Re}^e} + \frac{1}{2} \left(1 + 3.315\left ({Re}^e\right )^{-1/2}\right) \right ]^{-1} \right \}, \end{align}
\begin{align} C_D^S &= \frac {24}{{Re}^e} \left[ 1 + 0.15\left ({Re}^e\right )^{0.687}\right], \end{align}
\begin{align} m &= 1 + 0.01\left ({Re}^e\right )^{1.1}. \end{align}
Now, substituting the expression for the critical
${\textit{Re}}^e$
(given by (B1)) and the drag coefficient correlation ((B4) together with (B5)) into (B3), and noting that the viscosity and density of water under standard conditions are
$\mu ^e \approx 10^{-3} \, \text{Pa s}$
and
$\rho ^e \approx 1000 \, \text{kg} \, \text{m}^{-3}$
, respectively, we obtain solutions for
$R_c$
over a broad range of viscosity and density ratios,
$(\mu ^\ast , \rho ^\ast )$
, with the results summarised in figure 20.
Open access
