2.1. Problem statement and governing equations

We consider a rigid sphere with radius $R$ translating vertically with a prescribed constant velocity $-U\boldsymbol{e}_{x}$ through a linearly stratified fluid with a reference density $\rho _0$ and a vertical density gradient $\partial \tilde {\rho }/\partial x \lt 0$ , where $\boldsymbol{e}_{x}$ denotes the upward unit vector (see figure 1 a) and $\tilde {\rho }(x)$ is the undisturbed background density profile. The constant-velocity setting is adopted primarily to enable direct comparison with available experiments by Hanazaki et al. (Reference Hanazaki, Konishi and Okamura2009a ), Okino et al. (Reference Okino, Akiyama and Hanazaki2017) and Akiyama et al. (Reference Akiyama, Waki, Okino and Hanazaki2019); the same configuration was also employed in previous numerical studies (Torres et al. Reference Torres, Hanazaki, Ochoa, Castillo and Van Woert2000; Hanazaki et al. Reference Hanazaki, Nakamura and Yoshikawa2015; Zhang et al. Reference Zhang, Mercier and Magnaudet2019). The fluid has uniform dynamic viscosity $\mu$ and molecular diffusivity $\kappa$ . We normalise lengths, velocities, times, pressure and density variations using the characteristic scales $R$ , $U$ , $R/U$ , $\rho _0 U^2$ and $-R\,\partial \tilde {\rho }/\partial x$ , respectively. In the body-fixed reference frame, the fluid velocity is denoted by $\boldsymbol{u}(\boldsymbol{x}, t)$ at position $\boldsymbol{x}=(x,y,z)$ and time $t$ . The local fluid density is expressed as

(2.1) \begin{align} \tilde {\rho }(x) + \rho (\boldsymbol{x},t)= \rho _{0} + t - x + \rho (\boldsymbol{x},t)\,, \end{align}

where $\rho (\boldsymbol{x}, t)$ denotes the density disturbance resulting from the body motion.

Figure 1.

$(a)$ Sketch of the flow configuration and definition of some quantities; $(b)$ light fluid dragged down by the sphere at $ (\textit{Fr}, \textit{Re}, \textit{Pr}) = (1,100, 70)$ , with the iso-contours and colours highlighting isopycnals and density variations; $(c)$ stable jet observed in a 2-D axisymmetric simulation at $ (\textit{Fr}, \textit{Re}, \textit{Pr}) = (0.1,100, 700)$ ; $(d)$ unstable jet observed with the same set of parameters in a 3-D simulation. In $(c{,}d)$ , colours refer to the absolute vertical velocity $\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{e}_x - 1$ ; the inset provides an enlarged view of the upper part of the jet.

Figure 1. Long description

Panel A: A sketch of the flow configuration and definition of some quantities. It shows a sphere settling through a stratified fluid, with arrows indicating the direction of flow and density variations. The sphere is labeled with various parameters and vectors indicating directions. Panel B: A 2-D axisymmetric simulation showing a stable jet. The color gradient represents the absolute vertical velocity, with iso-contours and colors highlighting isopycnals and density variations. The inset provides an enlarged view of the upper part of the jet. Panel C: A 3-D simulation showing an unstable jet with the same set of parameters as Panel B. The color gradient represents the absolute vertical velocity, with iso-contours and colors highlighting isopycnals and density variations. Panel D: Another 3-D simulation showing an unstable jet with the same set of parameters as Panel B and C. The color gradient represents the absolute vertical velocity, with iso-contours and colors highlighting isopycnals and density variations.

The flow is described under the Boussinesq approximation, whose validity, even in the high-Prandtl-number regime ( $ \textit{Pr} \approx 700$ ), is supported by previous numerical studies through comparisons with either non-Boussinesq simulations (Wang et al. Reference Wang, Wang and Deng2023; Abdal et al. Reference Abdal, Kahouadji, Shin, Chergui, Juric, Caulfield and Matar2025) and experimental data such as the velocity and density distributions in the wake region (Okino et al. Reference Okino, Akiyama and Hanazaki2017, Reference Okino, Akiyama, Takagi and Hanazaki2021; Mercier et al. Reference Mercier, Wang, Péméja, Ern and Ardekani2020). Under the Boussinesq approximation, the non-dimensional governing equations for the velocity $\boldsymbol{u}$ and density disturbance $\rho$ take the form

(2.2) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0 , \\[-28pt] \nonumber \end{align}

(2.3) \begin{align} \partial _t\boldsymbol{u} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} = - \boldsymbol{\nabla }p - (\textit{Fr})^{-2}\rho {{\boldsymbol{e}}_x} + \textit{Re}^{-1}{{\nabla} ^2}{\boldsymbol{u}}, \\[-28pt] \nonumber \end{align}

(2.4) \begin{align} \partial _t\rho + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\rho = \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{e}_{x} - 1 + \textit{Pe}^{-1}{{\nabla} ^2}\rho. \\[-12pt] \nonumber \end{align}

Equations (2.3) and (2.4) involve the Reynolds and Péclet numbers, characterising the relative importance of advective and diffusive effects in the momentum and density balances, respectively. These two numbers are defined as $ \textit{Re} = \textit{UR}/\nu$ and $\textit{Pe} = \textit{UR}/\kappa$ , where $\nu = \mu /\rho _0$ is the kinematic viscosity. Their ratio yields the Prandtl number, $ \textit{Pr} = \nu /\kappa$ . Additionally, the ratio of inertial to buoyancy forces is quantified by the Froude number, $ \textit{Fr} = U/(\textit{NR})$ , with $N = [-(g/\rho _0)\partial _x\tilde {\rho }]^{1/2}$ being the Brunt–Väisälä frequency and $g$ denoting gravity. Several studies (Hanazaki et al. Reference Hanazaki, Nakamura and Yoshikawa2015; Okino et al. Reference Okino, Akiyama and Hanazaki2017; Akiyama et al. Reference Akiyama, Waki, Okino and Hanazaki2019) defined the Reynolds number based on the sphere diameter $2R$ , whereas the characteristic radius $R$ is selected here to maintain consistency with our previous work (Zhang et al. Reference Zhang, Mercier and Magnaudet2019) and align with the definition of $ \textit{Fr}$ . The hydrostatic pressure component, $-(gR/U^2)x + (\textit{Fr})^{-2}(x^2/2 - xt)$ , is incorporated into the modified pressure $p(\boldsymbol{x}, t)$ . At the sphere surface $\mathcal{S}$ , with unit normal $\boldsymbol{n}$ directed into the fluid, the no-flux and no-slip boundary conditions hold, implying

(2.5) \begin{align} \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\rho = \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_{x} \quad \text{and }\quad \boldsymbol{u}=0 \quad \text{on }(\mathcal{S})\, . \end{align}

The first of (2.5) is readily obtained from the standard no-flux condition $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }(\tilde {\rho }+\rho )=0$ by noting that $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\tilde {\rho } = -\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_{x}$ . In addition, all disturbances vanish in the far field, so that $\rho$ and $\boldsymbol{u}$ satisfy

(2.6) \begin{align} \rho \xrightarrow {} 0 \quad \text{and}\quad \boldsymbol{u} \xrightarrow {} \boldsymbol{e}_{x} \quad \text{for}\quad ||\boldsymbol{x}||\xrightarrow {}\infty \,. \end{align}

We are also interested in examining how the forces acting on the sphere vary with the governing parameters. Since the sphere is constrained to move in the $\boldsymbol{e}_x$ -direction, the drag force is always parallel to $\boldsymbol{e}_x$ , while a transverse (or lift) force due to the flow instability may arise in the horizontal $(\boldsymbol{e}_y,\boldsymbol{e}_z)$ plane. Therefore, the force components are defined as

(2.7) \begin{eqnarray} F_D = \boldsymbol{e}_x\boldsymbol{\cdot }\int _{\mathcal{S}} \mathbb{T} \boldsymbol{\cdot }\boldsymbol{n}\,{\rm d}\mathcal{S}\,, \quad F_{L,y} = \boldsymbol{e}_{y}\boldsymbol{\cdot }\int _{\mathcal{S}} \mathbb{T} \boldsymbol{\cdot }\boldsymbol{n}\,{\rm d}\mathcal{S}\,,\quad F_{L,z} = \boldsymbol{e}_{z}\boldsymbol{\cdot }\int _{\mathcal{S}} \mathbb{T} \boldsymbol{\cdot }\boldsymbol{n}\,{\rm d}\mathcal{S}\,, \end{eqnarray}

where $\mathbb{T} = -p\mathbb{I} + Re^{-1}(\boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{u}^{{T}})$ is the stress tensor, $\mathbb{I}$ denoting the unit tensor. Normalising these force components by $({\pi }/{2} )R^2 \rho _0 U^2$ yields the drag and lift coefficients $C_D$ and $C_{L,y}$ (or $C_{L,z}$ ), respectively. As the sphere settles through the stratified fluid, the isopycnals are deflected around it, as illustrated in figure 1 $(b)$ , which leads to an increase in the drag coefficient $C_D$ . In the absence of flow instability, a steady-state density field is eventually established, in which advection and diffusion are in balance. When the stratification is sufficiently strong, a buoyant jet forms downstream of the sphere, as shown in the example of figure 1 $(c)$ . In this figure, resulting from a 2-D axisymmetric simulation, contours of the absolute vertical velocity, $\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{e}_x - 1$ , reveal that the fluid moves significantly faster than the background flow in the wake region. Figure 1 $(d)$ depicts a qualitatively different scenario obtained from fully 3-D simulations with the same parameter set. The tail region of the jet is now unstable and meanders weakly, as emphasised in the enlarged view. The present study thus aims to characterise the onset, evolution, and parametric sensitivity of this 3-D jet instability, and clarify the underlying physical mechanisms at stake.

2.2. Numerical methods and grid design

The JADIM code employs a finite-volume discretisation of the governing equations (2.2)–(2.4) on a staggered grid, where velocity components are defined at cell faces, while pressure and density are stored at cell centres (Magnaudet et al. Reference Magnaudet, Rivero and Fabre1995). Time integration combines a third-order explicit Runge–Kutta scheme for advective and source terms with a semi-implicit Crank–Nicolson method for diffusive terms, ensuring second-order accuracy in both time and space (Calmet & Magnaudet Reference Calmet and Magnaudet1997). Incompressibility is enforced at each time step via a projection method: an intermediate velocity field is corrected by solving a Poisson equation for the pressure increment, ensuring that the final velocity field is divergence-free to machine precision. At high Péclet number, i.e. for large $ \textit{Pr}$ , the central differencing of the advective density flux can generate non-physical oscillations. To avoid such oscillations and preserve the monotonicity of the density disturbance, the advective flux $\boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u})$ is discretised using a total-variation-diminishing (TVD) ‘monotonic-centred’ scheme based on the Van Leer limiter (Van Leer Reference Van Leer1977). We also tested several weighted essentially non-oscillatory (WENO) schemes (Balsara et al. Reference Balsara, Garain and Shu2016) for the discretisation of the density advective flux and found that they yield almost identical results (not shown). Since WENO schemes are substantially more expensive computationally, owing to their larger interpolation stencils, the TVD Van Leer scheme represents an appropriate choice for the present 3-D simulations. The JADIM solver summarised above was extensively validated by Zhang et al. (Reference Zhang, Mercier and Magnaudet2019). In particular, the corresponding 2-D axisymmetric simulations past a settling sphere were shown to accurately capture both the detailed structure of the thin jet that emerges at low Froude number and the associated drag enhancement.

In the present study, the same numerical strategy is extended to a 3-D framework. Computations are performed in the spherical coordinate system sketched in figure 1 $(a)$ , with non-uniform grid distributions in the radial ( $\xi$ ) and polar ( $\theta$ ) directions, and a uniform grid distribution in the azimuthal ( $\phi$ ) direction. Near the sphere surface, the minimum grid spacing in the radial direction, $\varDelta _\xi$ , and that in the polar direction near the upper pole, $\varDelta _\theta$ , are set to $\varDelta _\xi = \varDelta _\theta = 5\times 10^{-4}$ . For the most demanding case considered, $(\textit{Fr}, \textit{Re}, \textit{Pr}) = (0.02,150, 700)$ , the density boundary layer is significantly thinner than the momentum one, with a characteristic thickness $\delta _\rho \sim Re^{-1/2}Pr^{-1/3}=\mathcal{O}(10^{-2})$ , while the radius of the jet scales as $\delta _j \sim (Fr/RePr)^{1/2} =\mathcal{O}(10^{-3})$ (Hanazaki et al. Reference Hanazaki, Nakamura and Yoshikawa2015; Okino Reference Okino2023). Therefore, the above resolution ensures that both the density boundary layer and the inner jet structure are properly resolved. A grid independence study presented in Appendix A confirms that doubling $\varDelta _\xi$ and $\varDelta _\theta$ does not alter the results. Influence of the azimuthal resolution is also considered in that appendix. It is established that selecting a uniform distribution with $\varDelta _\phi =\pi /32$ ensures grid convergence. The spherical computational domain has an outer radius $\xi _{\textit{max}} = 40$ , a size that was shown by Zhang et al. (Reference Zhang, Mercier and Magnaudet2019) to be sufficiently large to avoid artificial confinement effects, especially in highly stratified cases (very low $ \textit{Fr}$ ) where the jet is very short. The total number of cells in the domain is $(N_\xi , N_\theta , N_\phi ) = (420, 200, 64)$ , with the $\xi$ - and $\theta$ - directions consistent with the previously validated 2-D axisymmetric setup. The overall grid arrangement is illustrated in figure 20 $(a)$ .

The Dirichlet boundary condition (2.6) is imposed to the velocity field over the part of outer boundary $\xi = \xi _{\textit{max}}$ extending from the lower pole to the edge of the wake region. Within this wake region (defined arbitrarily as the cone originating from the sphere centre and making a $60^\circ$ semi-angle with the upper part of the $\boldsymbol{e}_{x}$ -axis), the non-reflecting outlet condition proposed by Magnaudet et al. (Reference Magnaudet, Rivero and Fabre1995) is imposed. To avoid spurious reflections of internal waves from the outer boundary, a linear Rayleigh damping strategy is applied to density disturbances within a sponge layer extending over the last five cells adjacent to the outer boundary in the $\xi$ -direction (Slinn & Riley Reference Slinn and Riley1998; Chongsiripinyo et al. Reference Chongsiripinyo, Pal and Sarkar2017). Specifically, a damping term $-\psi (\xi )\rho$ is added to the right-hand side of (2.4), where the weighting function $\psi (\xi )$ increases quadratically from zero at the inner edge of the sponge layer to unity at $\xi = \xi _{\textit{max}}$ . This method was already shown to be effective by Zhang et al. (Reference Zhang, Mercier and Magnaudet2019).

Some numerical aspects specific to 3-D simulations deserve to be mentioned. In the 2-D axisymmetric study of Zhang et al. (Reference Zhang, Mercier and Magnaudet2019), the time history of the drag coefficient was monitored, and the flow field was considered converged when $C_D$ -variations dropped below $0.2\,\%$ over the final 2000 time steps. In contrast, 3-D simulations are inherently more complex due to the possible emergence of the jet instability, which introduces temporal and spatial fluctuations. To address this difficulty, we established distinct convergence criteria based on the flow stability. For stable flows, e.g. at $ \textit{Fr} \geqslant 0.3$ for $(Re, Pr) = (100, 700)$ , the same $C_D$ -based strategy is still applied. Conversely, for strongly unstable flows, e.g. at $ \textit{Fr} \leqslant 0.05$ , we rely on a criterion based on the lift coefficient, $C_L$ . The flow is considered statistically stationary once $C_L$ exhibits periodic saturation over the last ten oscillation cycles (see panels $(c)$ and $(f)$ in figure 3). In contrast, in intermediate regimes, e.g. $0.05 \lt \textit{Fr} \lt 0.3$ , the instability may evolve very slowly, and the flow does not reach convergence even beyond $t \gt 80$ . These cases are therefore excluded from the quantitative analysis presented later. Moreover, we define the onset of jet instability based on a threshold on the azimuthal velocity component: the jet is said to be unstable when $|u_\phi ^{\textit{max}}| \gt 10^{-2}$ close to its centreline. Cases with lower values of $|u_\phi ^{\textit{max}}|$ are considered stable.

It is important to mention that no direct quantitative validation of the present 3-D simulations of the jet instability is currently possible because neither experimental data nor previous 3-D simulations under similar conditions are available for this problem, as already discussed in § 1. Existing experiments (Hanazaki et al. Reference Hanazaki, Kashimoto and Okamura2009a ; Akiyama et al. Reference Akiyama, Waki, Okino and Hanazaki2019; Okino et al. Reference Okino, Akiyama, Takagi and Hanazaki2021) clearly demonstrate the occurrence of jet meandering in strongly stratified regimes, but do not provide detailed 3-D data. Nevertheless, some indirect evidence can be used to support the reliability of the present 3-D results. In particular, figure 2 shows that the predicted $(\textit{Fr},\textit{Re})$ -domain where the jet instability manifests itself agrees well with the experimental observations of Akiyama et al. (Reference Akiyama, Waki, Okino and Hanazaki2019), while figure 14 shows that 2-D axisymmetric and full 3-D simulations provide vertical distributions of the axial velocity along the thin upward jet in close agreement with each other before the instability sets in.

In most simulations, the jet instability arises naturally due to the amplification of ambient numerical disturbances. Nevertheless, in regimes close to the onset of instability, e.g. $ \textit{Fr} \approx 0.1$ for $(Re, Pr) = (100, 700)$ , the natural development of disturbances is so slow that it is desirable to accelerate the growth of the instability with the help of an external disturbance. Similarly, in the kinetic energy budget analysis presented in § 4.2, artificial disturbances are required to properly estimate the growth rate. In these cases, the disturbance is introduced only during the first 5000–20 000 time steps beyond the time required to reach saturation in the equivalent 2-D axisymmetric case (depending on the time step, it is applied over a time period $ \tau \sim 1-2$ ); then the computation proceeds naturally. The disturbance is imposed on the vertical velocity field $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{e}_{x}$ in the form $u_x' = u_0' \exp \{-[(y - y_0)^2 + (x - x_0)^2] / \epsilon ^2\}$ , with amplitude $u_0' = 10^{-4}$ , central location $x_0 = 2, y_0 = 0.005$ and decay width $\epsilon = 0.05$ . Appendix B establishes that this artificial disturbance does not alter the intrinsic nature of the flow, even when it is given a much larger amplitude, i.e. $u_0' = 10^{-2}$ : the jet initially exhibits transient oscillations but subsequently relaxes back to its initial axisymmetric structure after the disturbance is removed.

Last, we employed a post-processing technique based on the tracking of trajectories of massless Lagrangian particles to investigate the flow behaviour in the sphere’s wake. Each particle evolves according to the kinematic relation ${d}_t\boldsymbol{x}_p=\boldsymbol{u}(\boldsymbol{x}_p,t)$ , where $\boldsymbol{x}_p$ denotes the instantaneous particle position. The velocity $\boldsymbol{u}(\boldsymbol{x}_p,t)$ at the particle location is obtained by interpolating the flow field over the surrounding $2 \times 2 \times 2$ grid cells. Then, starting from the particle position at time $n\Delta t$ (with $\Delta t$ denoting the time step), the position at time $(n+1)\Delta t$ is obtained explicitly as $\boldsymbol{x}_p^{\,n+1}=\boldsymbol{x}_p^{\,n}+\boldsymbol{u}(\boldsymbol{x}_p,t)\,\Delta t.$ The results reported below were obtained with $\Delta t = 1\times 10^{-2}$ . We verified that reducing the time step to $1\times 10^{-3}$ does not alter the particle trajectories, thereby confirming temporal convergence of the particle-tracking procedure.

3.1. Evolution of the transverse force stemming from the unstable jet

We first fix the Reynolds and Prandtl numbers at $(Re, Pr) = (100, 700)$ and progressively reduce the Froude number, starting from the critical curve shown in figure 2. Figure 3 $(a{-}d)$ presents colour maps of the absolute streamwise velocity $u_x = \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{e}_x-1$ projected onto the vertical cross-sectional plane $(x,y)$ for $ \textit{Fr} = 0.3$ , $0.1$ , $0.05$ and $0.02$ , respectively. In all cases, a high-speed jet forms in the wake, consistent with previous numerical findings (Torres et al. Reference Torres, Hanazaki, Ochoa, Castillo and Van Woert2000; Hanazaki et al. Reference Hanazaki, Nakamura and Yoshikawa2015; Zhang et al. Reference Zhang, Mercier and Magnaudet2019) and experimental observations (Hanazaki et al. Reference Hanazaki, Kashimoto and Okamura2009a ). A bell-shaped structure (dark blue region) is also visible along the jet axis, as previously described by Hanazaki et al. (Reference Hanazaki, Nakamura and Yoshikawa2015). These authors showed that this peculiar structure is associated with internal waves generated in the wake. As stratification intensifies, i.e. $ \textit{Fr}$ decreases, the jet becomes progressively shorter and thinner, in agreement with earlier observations.

Figure 3.

Examples of the jet evolution and instability characteristics observed by varying $ \textit{Fr}$ , with $ \textit{Re}=100$ and $ \textit{Pr}=700$ in all cases. $(a-d)$ Contours of the vertical velocity projected onto the vertical cross-sectional $(x,y)$ plane, with $(a)$ $ \textit{Fr} = 0.3$ (stable); $(b)$ $ \textit{Fr} = 0.1$ (chaotic instability); $(c)$ $ \textit{Fr} = 0.05$ (spiral instability); $(d)$ $ \textit{Fr} = 0.02$ (standing-wave instability). $(e-g)$ Time histories of the $y$ - and $z$ -components of $C_L$ during the saturated oscillation stage, with $(e)$ $ \textit{Fr} = 0.1$ ; $(f)$ $ \textit{Fr} = 0.05$ ; $(g)$ $ \textit{Fr} = 0.02$ .

As shown in figure 3 $(a{-}d)$ , while the jet and bell-shaped structure remain axisymmetric at $ \textit{Fr} = 0.3$ , oscillations emerge for $ \textit{Fr} \leqslant 0.1$ , indicating the onset of the 3-D instability, consistent with experimental observations (Hanazaki et al. Reference Hanazaki, Kashimoto and Okamura2009a ; Akiyama et al. Reference Akiyama, Waki, Okino and Hanazaki2019). At $ \textit{Fr} = 0.1$ , a weak meandering is visible near the jet tip at the vertical position $x \approx 4.5$ , while the bell-shaped structure remains axisymmetric. As stratification increases further, i.e. $ \textit{Fr} = 0.05$ and $0.02$ , the jet exhibits stronger oscillations that take place closer to the sphere, at $x \approx 2.5$ and $1.6$ , respectively. The meandering behaviour may be characterised by the lift coefficient, $C_L$ . Figure 3 $(e{-}g)$ shows the time histories of the $y$ - and $z$ -components of $C_L$ during the saturated oscillatory regime for the three unstable cases. These plots reveal sinusoidal oscillations in both $C_{L,y}$ and $C_{L,z}$ , with amplitudes increasing from $\mathcal{O}(10^{-4})$ at $ \textit{Fr} = 0.1$ – $\mathcal{O}(10^{-2})$ at $ \textit{Fr} = 0.05$ and $0.02$ . This increase likely arises from two factors, namely the shorter distance between the meandering jet and the sphere and the increase in the oscillation amplitude as $ \textit{Fr}$ decreases.

Figure 4.

Wake evolution tracked in the $(C_{L,y},C_{L,z})$ phase plane in three distinct unstable regimes, with $ \textit{Re}=100$ and $ \textit{Pr}=700$ in all cases: $(a)$ $ \textit{Fr} = 0.1$ (chaotic regime); $(b)$ $ \textit{Fr} = 0.05$ (spiral regime); $(c)$ $ \textit{Fr} = 0.02$ (standing-wave regime); $(d)$ fast Fourier transforms (FFT) of the lift coefficient $C_{L,y}$ for different stratification levels. In $(a-c)$ , the colour along the path transitions from grey to black following time progression.

Figure 4. Long description

The image contains four panels showing different aspects of particle motion in stratified fluids. Panel A is a scatter plot with axes labeled C_L,z and C_L,y, showing chaotic regime data. Panel B is a circular plot with axes labeled C_L,z and C_L,y, showing spiral regime data. Panel C is another scatter plot with the same axes as Panel A, showing standing-wave regime data. Panel D is a line graph with axes labeled FFT and 2πFrf, showing Fourier transforms of the lift coefficient for different stratification levels. The color transitions from grey to black along the path in Panel A, indicating time progression.

To better quantify the jet instability, we examine the time evolution of the lift coefficients in the $(C_{L, y},C_{L, z})$ phase space, as shown in figure 4 $(a{-}c)$ . At $ \textit{Fr} = 0.1$ , the path exhibits chaotic excursions throughout the phase plane. Although the lift coefficient remains small ( $\mathcal{O}(10^{-4})$ ), this behaviour is not due to numerical noise, since simulations with weaker stratifications (higher $ \textit{Fr}$ ) indicate that $C_L$ -variations remain of $\mathcal{O}(10^{-6})$ . A similar chaotic dynamics has been reported in unstratified flows past bluff bodies, such as circular disks (Shenoy & Kleinstreuer Reference Shenoy and Kleinstreuer2008; Auguste et al. Reference Auguste, Fabre and Magnaudet2009) and rigid spheroids (Chrust et al. Reference Chrust, Bouchet and Dušek2010). Such transitions are often associated with irregular vortex shedding patterns. We did not examine this possibility in detail, nor did we attempt to confirm more rigorously the chaotic nature of the observed motion.

At $ \textit{Fr} = 0.05$ , the $(C_{L, y},C_{L, z})$ trajectory exhibits a spiral structure, preceded by planar oscillations. We refer to this behaviour as the ‘spiral’ mode. Interestingly, similar spiral paths of transverse force coefficients have been reported in counter-flows past heated spheres (Kotouč et al. Reference Kotouč, Bouchet and Dušek2009), where local buoyancy forces induce a similar symmetry-breaking behaviour. In that case, two longitudinal vortex threads were found to rotate slowly around each other, corresponding to a gradual rotation of the wake’s symmetry plane. Notably, the peak amplitude observed here in the lift coefficients ( $C_L^{\textit{max}} \approx 0.014$ ) is comparable to that reported by the previous authors, namely $C_L^{\textit{max}} = 0.019$ (see their figure 13). As the stratification becomes stronger ( $ \textit{Fr} = 0.02$ ), the $(C_{L, y},C_{L, z})$ oscillations transition to a periodic regime with zero-mean value, confined to a fixed symmetry plane. We refer to this regime as the ‘standing-wave’ mode. Similar patterns have been observed in unstratified flows past thin disks and spheres (Shenoy & Kleinstreuer Reference Shenoy and Kleinstreuer2008; Chrust et al. Reference Chrust, Bouchet and Dušek2010), and are distinct from the classical planar ‘zig-zag’ mode, which also preserves planar symmetry but exhibits a non-zero-mean lift force (Fabre et al. Reference Fabre, Auguste and Magnaudet2008; Meliga et al. Reference Meliga, Chomaz and Sipp2009). Nevertheless, the peak amplitude in the lift coefficients is significantly smaller in the present case, with $C_L^{\textit{max}} \approx 0.02$ , to be compared with $C_L^{\textit{max}} = 0.042$ in the uniform flow past a circular cylinder at $ \textit{Re} = 90$ (Shenoy & Kleinstreuer Reference Shenoy and Kleinstreuer2008) or $C_L^{\textit{max}} = 0.069$ for a sphere at $ \textit{Re} = 150$ (Johnson & Patel Reference Johnson and Patel1999). Figure 4 $(d)$ presents the Fourier spectra of the lift coefficients for different levels of stratification. In all cases, the dominant peak closely matches the Brunt–Väisälä frequency, such that the normalised oscillation frequency $2\pi Fr(fR/U)$ is close to unity. This correspondence suggests that the fluctuations in the transverse force, hence the jet meandering, are intimately linked to the internal waves generated by the body as it descends through the fluid.

Figure 5.

Variation with the Froude number of the maximum amplitude of the transverse force, $C_L^{\textit{max}}$ (blue line, right axis), and the axial location $h_x^{\textit{max}}$ of the peak vertical velocity on the jet axis in the axisymmetric configuration (red line, left axis).

Additional simulations for stratification levels in between $ \textit{Fr} = 0.05$ and $ \textit{Fr} = 0.02$ reveal that the $(C_{L, y},C_{L, z})$ trajectory exhibits a flattened spiral shape at $ \textit{Fr} = 0.03$ (not shown). The successive transitions observed from $ \textit{Fr}=0.1$ – $ \textit{Fr}=0.02$ indicate that, as stratification intensifies, the meandering jet becomes increasingly confined to a 2-D vertical plane. Following the interpretation of Kotouč et al. (Reference Kotouč, Bouchet and Dušek2009) in the context of heated spheres, the emergence of the spiral mode is attributed to an instability in the far-wake region, where the vorticity threads are less constrained and allowed to twist due to a slow rotation of the wake’s symmetry plane. This motivates an investigation of the relationship between the observed instability ‘style’ and the axial location $h_x^{\textit{max}}$ at which the velocity reaches its maximum on the jet axis in the corresponding axisymmetric configuration. Figure 5 quantifies this connection across a range of Froude numbers. As $ \textit{Fr}$ decreases, so does $h_x^{\textit{max}}$ . The maximum amplitude of the transverse force, $C_L^{\textit{max}}$ , also decreases monotonically with $ \textit{Fr}$ , reinforcing the view that stronger stratification leads to earlier jet destabilisation and greater intensity of jet oscillations. The $h_x^{\textit{max}}(\textit{Fr})$ and $C_L^{\textit{max}}(\textit{Fr})$ curves exhibit a change of slope in the range $0.05\lesssim \textit{Fr}\lesssim 0.06$ . This change coincides with a transition in the instability behaviour - from a spiral mode originating in the far wake to a standing-wave mode initiated closer to the sphere – suggesting that the spatial location at which the jet destabilises plays a critical role in shaping the ensuing dynamics.

3.2. Wake structure

Figure 6 presents a typical example of the instantaneous wake structure in the spiral and standing-wave modes. Visualisations in panels $(a{-}b)$ reveal the helical structure of the spiral mode, with two longitudinal vortex threads wrapped around each other. A similar structure was identified by Kotouč et al. (Reference Kotouč, Bouchet and Dušek2009) in the wake of a heated sphere. According to their interpretation, advective effects in the near wake act to stabilise the flow and suppress vortex shedding, whereas further downstream the wake becomes more prone to instability.

Figure 6.

Vortical structure in the sphere’s wake. $(a{-}b)$ Spiral mode at $ \textit{Fr} = 0.05$ shown at two successive time instants, $t = 30$ and $t = 30.1$ ; $(c{-}d)$ standing-wave mode at $ \textit{Fr} = 0.02$ shown at $t = 17.09$ and $t = 17.11$ . In each row, panel $(\textrm{i})$ presents contours of the axial vorticity $\omega _x$ in a colour scale ranging from dark blue ( $\omega _x=-3$ ) to dark red ( $\omega _x=+3$ ), and streamlines in the cross-sectional plane $x = 1.8$ in $(a{-}b)$ and $x = 1.6$ in $(c{-}d)$ ; panel $(\textrm{ii})$ shows the 3-D vortical structure visualised by iso-surfaces $\omega _x = \pm 3$ in $(a{-}b)$ and $\omega _x = \pm 20$ in $(c{-}d)$ , with the vertical velocity iso-contour $(u_x - 1) = 1.2$ highlighted in red; panel $(\textrm{iii})$ is a close-up view of the wake structure shown in $(\textrm{ii})$ .

Here, at $ \textit{Fr}=0.05$ , stratification effects remain moderate enough for the jet to start meandering only some distance downstream of the sphere, thereby providing sufficient spatial extent for the wake symmetry plane to rotate slowly. This slow rotation is no longer observed at $ \textit{Fr}=0.02$ (panels $(c{-}d)$ ), the two counter-rotating longitudinal vortices remaining in a fixed plane in that case. The vortex structure exhibits a clear periodic shedding, in agreement with the oscillatory behaviour of the transverse force observed in figure 4 $(c)$ . The symmetry plane is seen to align with an azimuthal angle of $\theta = 45^\circ$ in the $(y,z)$ plane, consistent with the dominant direction of the lift oscillations observed in figure 4 $(c)$ . This further supports the classification of the observed instability as a standing-wave mode. The wake structure closely resembles that found in some regimes in unstratified flows past spheres and disks (Shenoy & Kleinstreuer Reference Shenoy and Kleinstreuer2008; Kotouč et al. Reference Kotouč, Bouchet and Dušek2009; Chrust et al. Reference Chrust, Bouchet and Dušek2010). These findings establish the one-to-one connection between the various oscillatory modes observed in the transverse force and the corresponding vortex structures in the wake. In the present study, the sphere is forced to translate vertically with a constant speed. If it were free to move under buoyancy effects, its path would be influenced by this unsteady wake dynamics. Therefore, it is expected that the distinct oscillatory wake modes identified above would translate into different styles of rise or fall of the body.

4.1. Stage $II$ : varicose instability

Figure 8.

Stage II: varicose instability. $(a)$ : trajectory of the transverse force components in the $(C_{L,y},C_{L,z})$ phase plane for $2 \lt t \lt 3$ , with circles indicating the time instants corresponding to the snapshots shown in $(b)$ ; $(b)$ horizontal cross-sections of the iso-surface $(u_x - 1) = 1.2$ at $x = 2.1$ at selected times instants; $(c{-}h)$ distribution of the vertical velocity in the cross-sectional $(x,y)$ plane at the same instants of time. In $(c{-}h)$ , the solid red and blue dashed lines denote the iso-contours $(u_x - 1) = 1.2$ and $(u_x - 1) = -1.4$ , respectively.

Figure 8 $(a)$ shows that, in the varicose stage, both components of the transverse force remain below $10^{-6}$ up to $t =2.6$ , indicating a nearly axisymmetric wake. This is corroborated by the horizontal cuts of the iso-surfaces of the absolute vertical velocity displayed in panel $(b)$ , which retain circular shapes undergoing an expanding-contracting-expanding sequence. The velocity contours in panels $(c{-}h)$ capture this process: the initial sword-like tip region (panels $(c{-}d)$ ) transitions to a tip region with a pronounced neck (panels $(e{-}f)$ ) before exhibiting clear signs of pinch-off (panels $(g{-}h)$ ). Such an axisymmetric bulging–necking–bulging dynamics, typical of varicose modes, is widely documented for liquid jets (Sato Reference Sato1960; Hussain & Thompson Reference Hussain and Thompson1980; Huang & Hsiao Reference Huang and Hsiao1999) and flickering diffusion flames (Cetegen & Dong Reference Cetegen and Dong2000; Zhang et al. Reference Zhang, Xia and Gao2021, Reference Zhang, Yang, Li, Lin, Qi and Xia2024). The behaviour observed here closely resembles that of such flames, in which the buoyancy-induced outer vortex ring plays a key role in triggering the varicose instability (Zhang et al. Reference Zhang, Yang, Li, Lin, Qi and Xia2024).

Figure 9.

Formation of vortex rings and their interaction with the jet during the varicose stage (stage II). $(a$ – $c)$ Streamlines in the laboratory frame outside the jet at $t=2.1$ , $2.2$ and $2.3$ , respectively; $(d)$ axial positions of the jet neck (symbols) and of the vortex rings (lines) as functions of time; solid and dashed lines denote the lower (A) and upper (B) rings in each pair, respectively; $(e)$ FFT of the drag coefficient during this stage. In panels $(a-b)$ , ‘VR’ stands for vortex ring, ‘A’ and ‘B’ denote the lower and upper rings within a given pair, and ‘I’–‘IV’ indicate the successive vortex-ring pairs arising during this stage.

Figure 9 $(a{-}c)$ shows the emergence of vortex rings around the jet. At $t=2.1$ , the velocity streamlines remain open. By $t=2.2$ , two vortex rings form. The first of them, VR-IA (notations used for the successive vortex rings are defined in the caption of figure 9), centred at $x = 1.92$ , results from the pinching of the neighbouring streamlines, whereas VR-IB, centred at $x = 2.0$ , is likely induced by the squeezing of fluid elements in between VR-IA and the jet. At $t=2.3$ , the previous vortex pair and the next one (VR-II) have been advected downstream and are no longer visible in the spatial window displayed in panel $(c)$ . In contrast, the next two pairs are, and the jet is seen to neck in between the two rings of the VR-III pair. The evolution of the axial position of the jet’s neck and those of the centres of the successive vortex rings are displayed in figure 9 $(d)$ , revealing that the former closely follows the latter. This in turn suggests that these vortex structures drive the necking and pinch-off of the jet. Figure 9 $(e)$ shows the FFT of the drag force acting on the sphere. It reveals two dominant peaks, one at $f_1=1.67$ corresponding to the shedding frequency of the vortex rings, and a second one at $f_2=8.0$ matching the Brunt–Väisälä frequency (see figure 4 d). The drag fluctuations induced by the shedding of vortex rings highlight the role of these specific flow structures in the stress distribution at the sphere surface.

Figure 10.

Flow field in the wake region. $(a)$ Temporal evolution of the radial velocity along the vertical line $y=0.05$ , $z=0$ ; $(b)$ time history of the density disturbance at two axial locations, $x=1.60$ and $1.88$ , along the same vertical line; $(c)$ isopycnal lines $\rho (x,r,t)-x=\mathrm{const.}$ in the axisymmetric base flow at two instants of time indicated by the dashed lines in $(b)$ . In $(c)$ , the isopycnals enclosed in the dashed ellipses are seen to be more widely spaced at the upper location than at the lower one.

Figure 10. Long description

Panel A: A line graph shows the temporal evolution of the radial velocity along the vertical line. The x-axis represents the position x in centimeters, and the y-axis represents the radial velocity u_y in centimeters per second. Multiple lines represent different time points, with a legend indicating specific times. Panel B: A line graph displays the time history of the density disturbance at two axial locations along the same vertical line. The x-axis represents time t in seconds, and the y-axis represents the density disturbance -ρ. Two lines represent different axial positions, with a legend indicating specific positions. Panel C: A line graph shows isopycnal lines in the axisymmetric base flow at two instants of time. The x-axis represents the position x in centimeters, and the y-axis represents the position y in centimeters. Multiple lines represent different times, with a legend indicating specific times. The isopycnals enclosed in the dashed ellipses are more widely spaced at the upper location than at the lower one.

The formation of recirculating eddies results from stratification-induced differences in the density-relaxation process, as figure 10 helps reveal. Figure 10 $(a)$ displays the temporal evolution of the radial velocity along a vertical line lying some distance away from the jet axis ( $y=0.05$ ). At $t\leqslant 2.179$ , the radial velocity remains positive at all vertical positions, indicating an outward flow. Then, at $t = 2.181$ , it becomes negative over a short range of vertical positions, $1.8\lt x\lt 1.87$ , signalling a local flow reversal associated with the onset of a vortex ring. The origin of this reversal is clarified in figure 10 $(b)$ , which displays the evolution of the density disturbance at two different heights along the above vertical line. While $\rho$ quickly stabilises at the lower position ( $x=1.6$ ), it continues to evolve over a much longer time at the upper position ( $x=1.88$ ). The imbalance in the density balance (2.4) responsible for this longer relaxation arises from the weaker radial diffusion at the upper position. Indeed, as figure 10 $(c)$ indicates, the isopycnals remain further apart at the upper location than at the lower one (see the distribution within the dashed ellipses, especially at $t=1.0$ ), which implies weaker radial density gradients at $x=1.88$ , hence a longer time for the density disturbance to reach a stationary value. This delayed adjustment of the upper fluid layer promotes a vertical gradient in buoyancy effects: the fluid in the upper region experiences a larger buoyancy force and therefore reaches a higher relative velocity, causing the streamlines to separate. Ultimately, this leads to the formation of a recirculating eddy and thereby to the onset of the varicose mode.

The above discussion makes it clear that the vortex-ring formation mechanism requires strong stratification conditions (low Froude number) and weak scalar diffusivity (high Prandtl number) for the lower fluid layer to equilibrate significantly faster than the upper layer, thereby promoting the emergence of a recirculating eddy. This is why the varicose instability only occurs in a specific parameter range, namely $ \textit{Fr} \lesssim 0.08,\ \textit{Re} \gtrsim 50,\ \textit{Pr}\gtrsim 70$ , as the regime map displayed in figure 19 will establish.

4.2. Stage III: onset of sinuous instability

As time proceeds, the varicose instability ceases and is succeeded by the sinuous instability corresponding to stage III in figure 7 $(b)$ . Figure 11 summarises this transition with the help of different indicators. The trajectories of the lift coefficients are seen to form spiral patterns in the $(C_{L,y},C_{L,z})$ phase plane (panel $(a)$ ), with $C_L$ growing to $\mathcal{O}(10^{-3})$ , signalling the loss of axial symmetry. Similarly, panel $(b)$ shows that the centre of the iso-contours of the vertical velocity drifts away from the vertical axis $(y,z)=(0,0)$ beyond $t \approx 5$ , an information confirmed by panels $(c{-}h)$ in which the meandering of the jet’s tip is seen to gradually increase over time. The instability that sets in resembles the sinuous modes observed in diffusion flames (Boulanger Reference Boulanger2010; Zhang et al. Reference Zhang, Xia and Gao2021; Xiao et al. Reference Xiao, Gupta, Macfarlane, Kennedy, Dunn, Kourmatzis, Torero and Masri2023). In this case, the instability arises from interactions between buoyancy and inertial forces – each dominating in different regimes – that destabilise the vortex sheet enveloping the flame, ultimately leading to the onset of the sinuous pattern.

Figure 11.

Stage III: onset of the sinuous instability. $(a)$ Trajectories of the transverse force components in the $(C_{L,y},C_{L,z})$ phase plane; $(b)$ iso-surfaces of the absolute vertical velocity $(u_x -1) = 1.2$ at $x = 2.1$ ; $(c{-}h)$ same as figure 8 $(c{-}h)$ in the time interval $3 \lt t \lt 6.4$ ; the successive snapshots are taken at times instants spotted by open circles in $(a)$ .

Figure 12.

Time evolution of the various terms in the disturbance kinetic energy budget (4.3) at different $ \textit{Fr}$ for $(\textit{Re},\textit{Pr}) = (100,700)$ : $(a)$ $ \textit{Fr} = 0.02$ ; $(b)$ $ \textit{Fr} = 0.05$ ; $(c)$ $ \textit{Fr} = 0.5$ .

Figure 12. Long description

Three line graphs depict the time evolution of various terms in the disturbance kinetic energy budget at different Froude numbers. Panel A: The line graph shows the time evolution of various terms in the disturbance kinetic energy budget for a Froude number of 0.02. The x-axis represents time (t) ranging from 12 to 16, and the y-axis represents the values of the terms ranging from -10 to 5. The legend indicates different terms: Advection (Advec) in blue, Shear in red, Buoyancy (Buoy) in green, Dissipation (Diss) in purple, Transport (Trans) in black, and Total in dashed black. The graph shows distinct trends for each term over time. Panel B: The line graph shows the time evolution of various terms in the disturbance kinetic energy budget for a Froude number of 0.05. The x-axis represents time (t) ranging from 20 to 35, and the y-axis represents the values of the terms ranging from -5 to 2. The legend is the same as in Panel A. The graph shows different trends for each term over time. Panel C: The line graph shows the time evolution of various terms in the disturbance kinetic energy budget for a Froude number of 0.5. The x-axis represents time (t) ranging from 40 to 100, and the y-axis represents the values of the terms ranging from -5 to 5. The legend is the same as in Panel A. The graph shows distinct trends for each term over time.

To better understand the mechanisms at stake, we examine the kinetic energy budget of the non-axisymmetric disturbance, which we assume to be small in the initial stage. This budget reads (Lombardi et al. Reference Lombardi, Caulfield, Cossu, Pesci and Goldstein2011; Pal et al. Reference Pal, Sarkar, Posa and Balaras2017)

(4.3) \begin{eqnarray} \partial _t e' = && -\overline {u}_j\partial _j e' -u'_iu'_{\kern-1pt j}\partial _j \overline {u}_i - (Fr)^{-2}\rho 'u'_x -{(2Re)}^{-1}(\partial _ju'_i + \partial _i u'_{\kern-1pt j})^2 \nonumber \\&& - \partial _i \big[u'_ie' + u'_ip' - 2Re^{-1} u'_{\kern-1pt j}(\partial _j u'_i +\partial _i u'_{\kern-1pt j})\big], \end{eqnarray}

where $e' = u_i'u'_i/2$ is the kinetic energy associated with the disturbance flow $u'_i = u_i - \overline {u}_i$ (with $\overline {u}_i$ the base flow corresponding to the axisymmetric solution), and $\rho '$ and $p'$ representing the non-axisymmetric density and pressure disturbances, respectively. As mentioned in § 2, the various terms involved in this budget are evaluated numerically by superimposing an artificial disturbance with a magnitude of $1\times 10^{-4}$ onto the $x$ -component of the velocity field. The first two terms in the right-hand side of (4.3) represent advection by the mean flow and generation by the mean shear, respectively, while the third term (sometimes referred to as the buoyancy flux) corresponds to generation by density disturbances and the last two terms represent dissipation and transport by the disturbance, respectively.

Figure 12 shows the time evolution of the volume-integrated terms of (4.3), evaluated over the wake region defined by the cylindrical domain $0 \leqslant r \leqslant 8$ and $1 \leqslant x \leqslant 10$ . Panels $(a{-}b)$ make it clear that, at low enough Froude number and sufficiently early times (up to $t\approx 26$ at $ \textit{Fr}=0.05$ ), the buoyancy flux (green line) is responsible for the generation of $e'$ . Generation by mean shear (red line) then starts to become significant, and both sources reach a similar magnitude after some time. The buoyancy flux is linked to the baroclinic vorticity generated by the tilting of velocity and density iso-surfaces near the base of the jet. The lower the diffusivity of the scalar, the more efficient this mechanism due to the strong horizontal density gradients that can then be maintained. Thus, increasing $ \textit{Fr}$ or decreasing $ \textit{Pr}$ reduces this contribution. Similarly, the shear-induced generation term decreases with increasing $ \textit{Fr}$ or decreasing $ \textit{Pr}$ since the jet then weakens and becomes thicker. These variations are confirmed in figure 12. Indeed, according to panels $(a{-}b)$ , both generation terms remain significant at $ \textit{Fr} = 0.05$ , although their magnitude is smaller than at $ \textit{Fr} = 0.02$ . Then, at $ \textit{Fr} = 0.5$ (panel $(c)$ ), they are insufficient to support the growth of the disturbance beyond some initial transient, so that the flow remains axisymmetric. We shall examine in more detail the mechanisms responsible for the onset of the sinuous instability in § 5.

4.3. Stages IV and V: transition and saturation of the sinuous instability

During stage IV, i.e. $6 \leqslant t \lt 9$ in figure 7 $(b)$ , the kinetic energy associated with the azimuthal velocity continues to increase but nonlinear effects make its growth deviate gradually from the previous exponential trend. Simultaneously, the height of the jet drops from $H^{{max} }= 2.2$ to $\approx 1.7$ , as figure 13 $(c{-}h)$ reveals. Panel $(a)$ in the same figure shows that the trajectories of the transverse force components evolve from a flattened spiral to a planar zig-zag as the jet shortens, indicating that the sinuous motion becomes increasingly confined to a vertical plane. Meanwhile, the beatings at the top of the jet reach a large amplitude, yielding strongly non-isotropic distributions of the vertical velocity in the horizontal plane (panel $(b)$ ). These beatings are responsible for the loss of coherence of the top part of the jet, and thus for its shortening.

Figure 13.

Stage IV: gradual saturation of the sinuous instability. $(a)$ Trajectories of the transverse force components in the $(C_{L,y},C_{L,z})$ phase plane; $(b)$ iso-surfaces of the vertical velocity $(u_x -1) = 1.2$ at $x = 2.0$ ; $(c-h)$ same as figure 8 $(c{-}h)$ in the time interval $6 \lt t \lt 8.8$ .

Figure 14.

Evolution of the absolute vertical velocity, $u_x - 1$ , along the vertical axis: $(a)$ $ \textit{Fr} = 0.02$ ; $(b)$ $ \textit{Fr} = 0.2$ . The inset in $(a)$ shows the time history of the streamwise and azimuthal enstrophy components during stage IV.

Figure 14. Long description

Panel A: A line graph depicts the evolution of absolute vertical velocity, u_z - 1, along the vertical axis, x, for a Froude number of 0.02. The x-axis ranges from 1.0 to 3.0, and the y-axis ranges from 0 to 25. Multiple lines represent different time points: Axisym. (black), t = 5.7 (red), t = 6 (blue), t = 6.798 (green), t = 7.479 (orange), t = 8.02 (purple), and t = 8.448 (brown). The inset shows the time history of the streamwise and azimuthal enstrophy components during stage IV, with the x-axis ranging from 0 to 12 and the y-axis ranging from 0 to 0.15 for enstrophy and from 2.8 to 3.8 for energy. Panel B: A line graph depicts the evolution of absolute vertical velocity, u_z - 1, along the vertical axis, x, for a Froude number of 0.2. The x-axis ranges from 1 to 8, and the y-axis ranges from 0 to 10. Multiple lines represent different time points: Axisym. (black), t = 10 (red), t = 15 (blue), t = 20 (green), t = 25 (orange), t = 30 (purple), and t = 35 (brown).

The gradual saturation process at work also alters the evolution of the vertical velocity distribution along the vertical axis, as illustrated in figure 14 $(a)$ . Variations in this distribution become visible at the jet’s tail around $t = 5.7$ and then propagate upstream, leading to oscillations at positions $2\lesssim x \lesssim 2.3$ . These oscillations grow over time, leading to the breakdown of the upper portion of the jet, characterised by near-zero values of $u_x$ , beyond $t = 6.8$ . The time evolution of the enstrophy associated with the streamwise vorticity component, $E_{\omega _x}=\int _{\varOmega }\omega _x^2\,{\rm d}\varOmega$ , is plotted in the inset enclosed in the same subfigure. This quantity is a relevant metric of the sinuous instability, since $\omega _x$ remains zero as long as the jet remains axisymmetric. $E_{\omega _x}$ is seen to increase sharply from near-zero values from $t \approx 6.5$ , while its azimuthal counterpart, $E_{\omega _\phi }=\int _{\varOmega }\omega _\phi ^2\,{\rm d}\varOmega$ , decreases simultaneously. These concurrent evolutions suggest that part of the primary azimuthal vorticity, $\omega _\phi$ , is converted into streamwise vorticity, $\omega _x$ , through a vortex tilting mechanism, a hypothesis that will be confirmed in § 5.1. For comparison, figure 14 $(b)$ displays the counterpart of the evolution on the jet velocity for a ten times larger Froude number, $ \textit{Fr} = 0.2$ . In this case, the $u_x$ -perturbations remain localised in the tail region ( $x \geqslant 5$ ) over a much longer period of time ( $10 \leqslant t \leqslant 35$ ), without significantly affecting the overall structure of the jet.

Figure 15.

Stage V: saturation of the sinuous instability. $(a)$ Trajectories of the transverse force coefficients in the $(C_{L,y},C_{L,z})$ phase plane; $(b)$ iso-surfaces of the absolute vertical velocity $(u_x -1) = 1.2$ at $x = 1.5$ ; $(c{-}h)$ same as figure 8 $(c-h)$ in the time interval $15.95 \lt t \lt 16.04$ corresponding to one period of the jet oscillations.

The sinuous instability reaches a saturated state in stage V. The corresponding flow characteristics are illustrated in figure 15. Both the evolution of the transverse force components in the $(C_{L,y},C_{L,z})$ plane (panel $(a)$ ) and the iso-surfaces of the vertical absolute velocity (panel $(b)$ ) indicate a planar zig-zagging oscillation in a vertical plane close to the $y = z$ diagonal. Unlike in stage IV, the height of the jet remains nearly constant throughout the sequence displayed in panels $(c{-}h)$ , confirming that the sinuous instability has reached saturation, leading to a statistically stationary wake structure.