3.1. Qualitative plume structure and dynamics

The qualitative plume structure and dynamics are indicated in figure 3, which shows time series of the density field over a time span of roughly $\tau ^{*}$ for two different values of ${Ri}_0$, as well as three different values of ${Re}_0$ for each ${Ri}_0$. Note that $\tau ^{*}$ generally decreases with increasing ${Ri}_0$ and increasing ${Re}_0$, and therefore the time series in figure 3 span different physical times.

Figure 3.

Time series of the density field over a time span of roughly $\tau ^{*}$ for six different cases: (a) ${Ri}_0=2$, ${Re}_0=100$; (b) ${Ri}_0=2$, ${Re}_0=316$; (c) ${Ri}_0=2$, ${Re}_0=1000$; (d) ${Ri}_0=20$, ${Re}_0=100$; (e) ${Ri}_0=20$, ${Re}_0=316$; (f) ${Ri}_0=20$, ${Re}_0=1000$. Time progresses in increments of $\tau ^{*}/4$ from left to right. Panels (a i–f i) show where the AMR grid is refined, with grey indicating $\ell =2$, and black indicating $\ell =3$. Panels (a v–f v) include the in-plane direction of the velocity vector (not scaled with velocity magnitude).

In each of the simulations, figure 3 shows that large-scale vortical structures form at the base of the plume and then rise vertically due to buoyancy, entraining fluid as they propagate (this entrainment is indicated by the streamlines provided in the rightmost column of the figure). This process repeats for the duration of the simulations at a time scale of approximately $\tau ^{*}$, corresponding to the puffing instability discussed in § 1; measurements of $\tau ^{*}$ and the puffing frequency are given in § 3.5. Figure 3 further shows that the vortical structures are more diffuse for the ${Ri}_0=2$, ${Re}_0=100$ case, as compared to larger values of ${Re}_0$ at the same ${Ri}_0$. This is due to the slower propagation speed of the vortices, as compared to the rate of diffusion.

From a dynamical perspective, figure 3 shows that at ${Ri}_0=2$, ${Re}_0=1000$, the puffing oscillations undergo a transition where the vortices become more turbulent, as indicated by the development of motions over multiple scales within the vortices, consistent with experimental observations made by Subbarao & Cantwell (Reference Subbarao and Cantwell1992). The plume also appears to lean to one side as a result of a flapping mode, which has been observed previously in both non-reacting plumes (Cetegen Reference Cetegen1997b) and diffusion flames Cetegen & Dong (Reference Cetegen and Dong2000); this mode will be discussed in more detail in § 3.2.

For the ${Ri}_0=20$, ${Re}_0=100$ case shown in figure 3, there is still relatively little small-scale motion due to the low ${Re}_0$, but the formation of vortices is more sinuous (i.e. anti-symmetric) in nature, as compared to the varicose (i.e. symmetric) behaviour observed at the same ${Re}_0$ for ${Ri}_0=2$. At larger ${Re}_0$ for ${Ri}_0=20$, the flow is generally turbulent, even along the centreline near the bottom of the domain. This early transition to turbulence is correlated with the penetration of heavy air into the core of the plume, as can be seen in a number of snapshots. We call these Rayleigh–Taylor (RT) spikes due to their resemblance to classical RT instabilities. Overall, this leads to a puffing instability that is more turbulent than the puffing instability at lower ${Ri}_0$ and ${Re}_0$.

3.2. Temporal variability

Figure 4 shows time series of the vertical velocity $u_z/U_0$ and $Y_{He}$ at height $z/R_0=0.5$ along the centreline (i.e. $r/R_0=0$) over time span $10\tau ^{*}$. Note that time in figure 4 has been normalized by the expected time scale $\tau ^{*}$, based on (1.1) rather than on the computed puffing frequency, which is presented and discussed in § 3.5. In general, the time series of $u_z/U_0$ become increasingly complex with increasing ${Ri}_0$ and ${Re}_0$. For ${Ri}_0=2$ in figure 4(a), the time series are smooth and periodic for all ${Re}_0$, with a characteristic time scale close to $\tau ^{*}$, corresponding to the puffing motion. As either ${Ri}_0$ or ${Re}_0$ is increased, the local minima and maxima begin to vary more substantially from cycle to cycle. For sufficiently high ${Re}_0$, the signals become erratic and chaotic, indicating the transition to a fully turbulent plume as shown, for instance, in figures 4(c,d).

Figure 4.

Time series of (a–d) $u_z/U_0$ and (e–h) $Y_{He}$ on the centreline ($r/R_0=0$) at height $z/R_0=0.5$ above the plume. Different line colours and widths correspond to different values of ${Re}_0$. Note that the time has been normalized by the expected time scale $\tau ^{*}$, based on (1.1) rather than the computed puffing frequency presented later in the text. Plots are for: (a,e) ${Ri}_0 = 2$, (b,f) ${Ri}_0 = 6.3$, (c,g) ${Ri}_0 = 20$, and (d,h) ${Ri}_0 = 63$.

Significantly, the transition to turbulent behaviour does not occur monotonically with increasing ${Re}_0$. For ${Ri}_0=2$ in figure 4(a), the local maxima in $u_z/U_0$ are much more consistent between each puffing cycle for ${Re}_0=100$ and ${Re}_0=1000$, while at ${Re}_0=316$, the local maxima between each puffing cycle vary more substantially. This increasing complexity is also evident in the time series of $Y_{He}$ for ${Ri}_0=6.3$ shown in figure 4(f); the only deviation from a uniform time series occurs momentarily for the intermediate Reynolds number ${Re}_0=316$. This behaviour is more consistent with nonlinear dynamical systems that have complex paths to chaos (e.g. period doubling phenomena, windows of less chaotic dynamics) rather than a simple broadening of scales typical of a fully turbulent flow.

The effects of the RT spikes noted in figure 3 are also evident in figure 4. For the two smallest values of ${Ri}_0$, in figures 4(e,f), $Y_{He}\approx 1$ at all times, implying that the more dense ambient air does not penetrate the centre of the plume at $z/R_0 = 0.5$. However, for ${Ri}_0=20$ and $63$ in figures 4(g,h), $Y_{He}$ frequently decreases below $1$, with the frequency, duration and strength of the decreases growing as ${Ri}_0$ and ${Re}_0$ increase. As indicated in figure 3, instances of $Y_{He} < 1$ are not a result of shear layer roll-up causing the ambient air to reach the centreline horizontally. Rather, these decreases are a result of a downward spike of heavier air. This phenomenon is particularly important when considering extensions to reacting plumes, since this would indicate the presence of an additional mechanism for supplying oxidizer to the flame beyond the mixing along the shear layer.

The flapping motion mentioned in § 3.1 and shown for the ${Ri}_0=2$, ${Re}_0=1000$ case in figure 3 is indicated quantitatively in figure 5 through time series of the radial velocity $u_r/U_0$ on the plume centreline at $z/R_0=0.5$. The radial velocity, in particular, corresponds to horizontal motions of the plume in the $x$–$y$ plane and is thus associated with the horizontal flapping behaviour shown in figure 3. For ${Ri}_0=2$ and ${Re}_0=1000$ in figure 5(a), $u_r/U_0$ oscillates at a frequency that is half that of the oscillations in $u_z/U_0$ shown in figure 4(a), although both oscillations are synchronized. That is, despite the difference in frequency, the sinuous mode associated with the plume flapping is synchronized with the varicose mode associated with the puffing. The flapping oscillation does not occur for the smaller values of ${Re}_0$ at ${Ri}_0=2$, although figure 5(b) shows that the oscillation occurs for all ${Re}_0$ at ${Ri}_0=6.3$. Variations in the radial velocity in the more turbulent simulations for higher ${Ri}_0$ are due primarily to the turbulent nature of the flow, rather than to a distinguishable oscillatory mode, hence we do not plot them here.

Figure 5.

Time series of $u_r/U_0$ on the centreline ($r/R_0=0$) at height $z/R_0=0.5$ above the plume for (a) ${Ri}_0=2$ and (b) ${Ri}_0=6.3$. Different line colours and widths correspond to different values of ${Re}_0$.

The dynamics in the shear layer are indicated by the time series of $u_z/U_0$ and $Y_{He}$ shown in figure 6. These time series are taken at the same height (i.e. $z/R_0=0.5$) as in figure 4, but at radial location $r=\delta _Y/2$, where $\delta _Y/2$ is the half-width of the flow based on where the temporal and azimuthal average of $Y_{He}$ is 0.5. By probing this radial location, we ensure that there are substantial and coherent fluctuations in the density field for all simulations, unlike at the centreline.

Figure 6.

Same as figure 4, but at location $r=\delta _Y/2$, $z/R_0=0.5$ within the shear layer.

The time series of $u_z/U_0$ in figures 6(a–d) are generally consistent with the centreline results in figure 4. Namely, the complexity of the time series increases with increasing ${Ri}_0$ and ${Re}_0$, with consistently repeating signals observed only for the smallest value of ${Ri}_0$ shown in figure 6(a). By contrast to the centreline results, however, the increase in complexity occurs monotonically for each ${Ri}_0$ with increasing ${Re}_0$; these trends will be discussed in more detail in § 3.3.

The shear layer time series of $Y_{He}$ in figures 6(e–h) vary far more substantially at each value of ${Ri}_0$ and ${Re}_0$ than the corresponding centreline time series in figure 4. Most notably, for ${Ri}_0=2$ as ${Re}_0$ increases, figure 6(e) shows a transition from a smooth, almost sinusoidal, variation in $Y_{He}$ to a time series with more distinct plateaus at $Y_{He}\approx 1$. This is the result of the slower time scale associated with convection for the lower ${Re}_0$, allowing for increased broadening of the helium–air interface as a result of diffusion. At larger ${Re}_0$, the flow is accelerated more rapidly, and as a result, can produce gradients that are even stronger than those initiated at the inlet. This transition can also be seen at ${Ri}_0=6.3$ in figure 6(f) for lower ${Re}_0$. In the simulations with ${Ri}_0\geq 20$, while the signals are not as regular, we also see sharper changes in $Y_{He}$ within the shear layer, as compared with the centreline results in figure 4, consistent with previous observations of scalar mixing in turbulent flows (Sandoval Reference Sandoval1995; Buch & Dahm Reference Buch and Dahm1998).

Finally, at the location in the shear layer, there are clear and consistent phase relationships between the different variables. Because density (related directly to $Y_{He}$) is an active scalar for vertical momentum transport, it is intuitive that local changes in $Y_{He}$ would also be reflected in $u_z$. In particular, comparison of figures 6(a–d) with figures 6(e–h) reveals that large peaks in $u_z/U_0$ are preceded by increases in $Y_{He}$. This makes sense intuitively because larger density discrepancies lead to larger buoyant forces that accelerate the flow. This is true for both the laminar and turbulent plumes, and these correlations are examined in more detail in the next subsection.

3.3. Transition from laminar to turbulent flow

A transition from laminar to chaotic dynamics in the near field of the plumes can be observed in the time series in § 3.2 as ${Ri}_0$ and ${Re}_0$ increase. In this subsection, we examine this transition more closely by using the state space dynamics of the vertical velocity $u_z/U_0$. We first focus on this transition in the near-field region where the puffing instability is most evident. It is important to note that for buoyant plumes, the flow will transition from laminar to turbulent sufficiently far downstream (except for very small ${Re}_0$), since potential energy is continually converted to kinetic energy, resulting in an increase in the local Reynolds number (Hunt & Van den Bremer Reference Hunt and Van den Bremer2011). Therefore, we first focus on this transition at a fixed spatial location within the near-field region, then show the transition from laminar to chaotic dynamics as downstream distance is varied.

To begin, figure 7 shows the state space of $u_z(t+\tau ^{*}/4)/U_0$ versus $u_z(t)/U_0$ along the centreline at $z/R_0=0.5$ for each simulation summarized in table 3. Each trajectory is coloured according to the simulation time such that lighter yellow indicates earlier times and darker blue indicates later times. White space indicates a state that the trajectory does not reach for the given spatial location during the analysis period. For the two smallest values of ${Re}_0$ in the ${Ri}_0=2$ and $6.3$ cases, the state-space trajectories form a stable limit cycle corresponding to the puffing behaviour. However, as ${Re}_0$ increases to $316$ for ${Ri}_0=2$, the trajectory is not as consistent compared to other trajectories for this value of ${Ri}_0$, even compared to higher ${Re}_0$, indicating increased flow complexity. This ‘island’ of increased variability, as indicated by the larger amount of state space occupied by the trajectory, is consistent with the non-monotonicity with ${Re}_0$ noted in figure 4(a), and is likely due to the complex coupling between the puffing oscillations and shear layer roll up. Overall, however, for ${Ri}_0=2$, figure 7 does not reveal a transition to fully chaotic (i.e. turbulent) dynamics for the range of ${Re}_0$ examined here.

Figure 7.

State-space trajectories of the vertical velocity at $r/R_0=0$, $z/R_0=0.5$, generated by plotting the same signal on different axes with an offset by $\tau ^{*}/4$. The lines are coloured by time elapsed during the simulation. Plots are for: (a–e) $Ri_0 = 2$, (f–j) $Ri_0 = 6.3$, (k–o) $Ri_0 = 20$, (p–r) $Ri_0 = 63$; and (a,f,k,p) $Re_0=100$, (b,g,l,q) $Re_0=178$, (c,h,m,r) $Re_0=316$, (d,i,n) $Re_0=562$, (e,j,o) $Re_0=1000$.

For ${Ri}_0=6.3$, by contrast, figure 7(f–j) shows that there is a substantial increase in complexity between ${Re}_0=178$ and ${Re}_0=316$. Because ${Re}_0$ controls the prevalence of small-scale motions and the onset of nonlinear chaotic dynamics, after a critical ${Re}_0$, the puffing instability becomes sufficiently chaotic such that the state-space trajectories are not as consistent. When this occurs, the trajectories occupy more of the state-space locations within a bounded region, eventually exploring most points within the region. Additional data and analysis are required to claim definitively whether the trajectories occupy all of the state space within the bounded region. For the two highest values of ${Ri}_0$ examined here, figure 7 shows that the state-space dynamics are fully chaotic even for the lowest ${Re}_0$, indicating that the transition to chaotic dynamics in these cases occurs for ${Re}_0$ less than 100.

Based on these results, the critical Reynolds number, denoted ${Re}_{crit}$, at which the transition to chaotic dynamics occurs depends on ${Ri}_0$, with smaller values of ${Re}_{crit}$ for increasing ${Ri}_0$. Among the present simulations, ${Re}_{crit}>1000$ for ${Ri}_0=2$, $178<{Re}_{crit}<316$ for ${Ri}_0=6.3$, and ${Re}_{crit}<100$ for both ${Ri}_0=20$ and $63$. In general, for a given ${Ri}_0$, as ${Re}_0$ increases, the flow transitions from steady laminar to unsteady laminar, to unsteady transitional, to unsteady turbulent behaviours.

The state-space trajectories in figure 7 are all calculated at a single centreline location, but the transition from laminar to turbulent behaviours also occurs with increasing vertical location due to the conversion of potential (i.e. buoyant) energy into kinetic energy. Figure 8 shows this transition for the ${Ri}_0=2$, ${Re}_0=1000$ case. At $z/R_0=0.25$ and $1$, the state-space trajectories approximately follow stable limit cycles, but at $z/R_0=2$, the region of state space occupied by the trajectory has grown substantially and the trajectory explores essentially all of the states within this region.

Figure 8.

State-space trajectories of the vertical velocity for ${Ri}_0=2$, ${Re}_0=1000$ along the centreline at different vertical locations. The lines are coloured the same as in figure 7: (a) $z/R_0=0.25$, (b) $z/R_0=1$, (c) $z/R_0=2$.

This dynamical transition with downstream location is expected in most plumes as the flow convects vertically (Hunt & Van den Bremer Reference Hunt and Van den Bremer2011) and has been noted before by Subbarao & Cantwell (Reference Subbarao and Cantwell1992). It is possible that a robust criterion could be derived to show more precisely the spatial location where this transition occurs. For example, estimates of the correlation dimension can be computed using standard algorithms in nonlinear dynamics (Hegger, Kantz & Schreiber Reference Hegger, Kantz and Schreiber1999), where most laminar cases would be essentially one-dimensional but turbulent ones would be larger; this analysis has been done for the closely related flame flickering instability (Gotoda, Kawaguchi & Saso Reference Gotoda, Kawaguchi and Saso2008). Application of these criteria to the present simulations is left as a direction for future research.

It should be noted that within the two distinct regimes of puffing indicated by figures 7 and 8 (i.e. laminar and turbulent), the state-space trajectories continue to change with increasing ${Re}_0$. For example, in the ${Ri}_0=2$ simulations at the two lowest values of ${Re}_0$, we see ‘period doubling’ behaviour where the periodic trajectory changes from requiring one orbit at ${Re}_0=100$ to two orbits at ${Re}_0=178$ to return to the same state-space location. For the two highest values of ${Ri}_0$, the flow is turbulent for all ${Re}_0$, but the bounded region within which the state-space trajectories reside grows with increasing ${Re}_0$. This is consistent with many turbulent flows where increases in the Reynolds number correspond to an increasing prevalence of intermittent rare events that expand the occupied state-space region. At a certain value of ${Re}_0$, however, results for the two largest values of ${Ri}_0$ shown in figure 8 indicate that this region approaches a fixed size, corresponding to fully turbulent behaviour.

3.4. Temporal correlations

The results in § 3.2 suggest that the flow variables are correlated in potentially complicated ways. These correlations and the temporal correspondence between variables can be examined in more detail using state-space trajectories such as those shown in figure 9, where trajectories are calculated in the $p$–$u_z$ state space along the centreline at $z/R_0=0.5$ for three different values of ${Re}_0$ and for each ${Ri}_0$. For ${Ri}_0=2$ and ${Re}_0=100$, $p$ and $u_z$ are highly correlated, with $p$ lagging behind $u_z$ (i.e. the trajectory in figure 9(a) goes counterclockwise as time evolves). In the context of these state-space trajectories, correlations are indicated by consistent orbital motions through state space, where $p$ and $u_z$ vary together in a regular way.

Figure 9.

State-space trajectories of the vertical velocity $u_z/U_0$ and mechanical pressure $p/\rho _\infty U_0^{2}$ at $r/R_0=0$, $z/R_0=0.5$. The lines are coloured by time elapsed during the simulation. Plots are for: (a) $Ri_0=2$, $Re_0=100$; (b) $Ri_0=6.3$, $Re_0=100$; (c) $Ri_0=20$, $Re_0=100$; (d) $Ri_0=63$, $Re_0=100$; (e) $Ri_0=2$, $Re_0=316$; (f) $Ri_0=6.3$, $Re_0=316$; (g) $Ri_0=20$, $Re_0=316$; (h) $Ri_0=63$, $Re_0=316$; (i) $Ri_0=2$, $Re_0=1000$; (j) $Ri_0=6.3$, $Re_0=1000$; (k) $Ri_0=20$, $Re_0=1000$.

As ${Ri}_0$ increases from 2 to 6.3, figure 9(b) shows that the state-space trajectories become more complex, with period doubling and slight broadening of the limit cycles, although $p$ still lags behind $u_z$. Further increases in ${Re}_0$ or ${Ri}_0$ result in a chaotic flow for which the correlation between $p$ and $u_z$ is much less pronounced. In particular, results for the larger values of ${Re}_0$ and ${Ri}_0$ in figure 9 show that some orbital trajectories do still exist, but the trajectories explore an increasingly large region of state space. State-space trajectories including $u_r$ and $Y_{He}$ are less interesting at this spatial location and are thus not shown here.

Figure 10 shows corresponding state-space trajectories within the shear layer at $z/R_0=0.5$, $r=\delta _Y/2$. Given the more substantial variations in $u_r$ and $Y_{He}$ at this location (as is also indicated by the variations in figure 6), here we show state-space trajectories for all possible combinations of $Y_{He}$, $u_r$ and $p$. Once again, for small ${Ri}_0$ and ${Re}_0$, these trajectories form closed and compact orbits, with $u_r$ and $Y_{He}$ lagging behind $p$. Large values of $u_r$ and $Y_{He}$ are generally highly correlated since, when $u_r$ is large and positive (i.e. the plume is not entraining strongly), the flow is generally accelerated vertically by buoyancy, corresponding to higher values of $Y_{He}$.

Figure 10.

Correlations between $Y_{He}$ and $u_r$ (green circles), $Y_{He}$ and $p$ (blue squares), and $u_r$ and $p$ (red triangles), at $r=\delta _Y/2$, $z/R_0=0.5$ within the shear layer for several of the simulations. Plots are for: (a) $Ri_0=2$, $Re_0=100$; (b) $Ri_0=6.3$, $Re_0=100$; (c) $Ri_0=20$, $Re_0=100$; (d) $Ri_0=63$, $Re_0=100$; (e) $Ri_0=2$, $Re_0=316$; (f) $Ri_0=6.3$, $Re_0=316$; (g) $Ri_0=20$, $Re_0=316$; (h) $Ri_0=63$, $Re_0=316$; (i) $Ri_0=2$, $Re_0=1000$; (j) $Ri_0=6.3$, $Re_0=1000$; (k) $Ri_0=20$, $Re_0=1000$.

3.5. Puffing frequency

The frequency at which large-scale vortices are shed during the buoyant plume evolution is commonly referred to as the puffing frequency. This frequency is the most widely used statistical measure of the puffing instability, and both experimental and computational approaches can easily measure and compare this simple metric. In our simulations, the periodic puffing motion is present qualitatively in essentially all the data presented in §§ 3.1–3.4.

The most typical and simplest method of measuring the puffing frequency is to extract a time series of a flow variable (typically the mechanical pressure or vertical velocity) and perform a fast Fourier transform (FFT) to obtain the power spectral density (PSD) of the signal. The peak of the PSD is generally classified as the puffing frequency. In figure 11, we show the PSD of the vertical velocity $u_z$ taken along the centreline at $z/R_0=0.5$ for six different simulations. Each of the PSDs has been normalized by the maximum magnitude, and the Strouhal number ${St}_0$ is normalized by the expected value ${St}_0^{*}$ based on Wimer et al. (Reference Wimer, Lapointe, Christopher, Nigam, Hayden, Upadhye, Strobel, Rieker and Hamlington2020) and (1.1).

Figure 11.

Normalized power spectral densities (PSDs) of $u_z$ along the centreline at $z/R_0 = 0.5$ for (a–c) ${Ri}_0=2$ and (d–f) ${Ri}_0=20$, at (a,d) ${Re}_0=100$, (b,e) ${Re}_0=316$, and (c,f) ${Re}_0=1000$. The observed Strouhal number ${St}_0$ is normalized by the predicted value ${St}^{*}_0$ obtained from (1.1).

For ${Ri}_0=2$ and all values of ${Re}_0$, figures 11(a–c) show that each PSD has a distinct peak at ${St}_0/{St}_0^{*} \approx 0.84$, indicating that the puffing frequency is slightly below that predicted by Wimer et al. (Reference Wimer, Lapointe, Christopher, Nigam, Hayden, Upadhye, Strobel, Rieker and Hamlington2020) (there are also smaller peaks at higher harmonics of the puffing frequency). For each of the ${Ri}_0=20$ cases, by contrast, we do not observe as clear a peak in the PSD. When ${Re}_0=100$ in figure 11(d), there are actually two peaks that could conceivably be identified as the puffing frequency. At ${Re}_0=316$ in figure 11(e), the peak in the signal is fairly clear but far from the expected value. At ${Re}_0=1000$ in figure 11(f), it is difficult to say with any certainty that a single distinct peak even exists.

Given the results in figure 11, measuring a single distinct puffing frequency for each simulation using a probing technique can be difficult. Subbarao & Cantwell (Reference Subbarao and Cantwell1992) also noted this difficulty for higher ${Re}$ flows when conducting experiments. Because puffing is a global instability, we should be able, in principle, to probe the flow anywhere within the plume and determine the puffing frequency. In previous studies, the probe location was typically placed where the largest fluctuations occur, roughly one radius vertically along the centreline (Cetegen & Kasper Reference Cetegen and Kasper1996) or slightly off-centre at the base of the plume (Bharadwaj & Das Reference Bharadwaj and Das2017). A more robust technique that better captures the global nature of the instability was used by Wimer et al. (Reference Wimer, Day, Lapointe, Meehan, Makowiecki, Glusman, Daily, Rieker and Hamlington2021), where a singular value decomposition (SVD) was performed on two-dimensional slices of vertical velocity.

For the present simulations, neither a single probe nor an SVD approach produced consistent results that reflected accurately the puffing frequency for all simulations. The complexity of measuring the puffing frequency is indicated in figure 12. At each spatial location, we computed the FFT of the time series of vertical velocity, identified which frequency was associated with the peak in the spectra, computed ${St}_0$ corresponding to that frequency, then normalized based on the expected ${St}_0^{*}$ from (1.1). This was done for each simulation, and one radial slice is shown in figure 12. Further, we used a binary logarithm, $\log _2({St}_0/{St}_0^{*})$, such that a value of 0 (white in figure 12) corresponds to the peak ${St}_0$ being identical to $St_0^{*}$, and integer values correspond to harmonics of ${St}_0^{*}$.

Figure 12.

Frequency corresponding to the peak in the PSD of the fluctuating streamwise velocity computed at each spatial location within a radial slice for all simulations performed in this study. Images are coloured using $\log _2({St}_0/{St}_0^{*})$ such that 0 (or white) corresponds to the computed ${St}_0$ being identical to ${St}_0^{*}$ from (1.1), with integer values corresponding to harmonics. Plots are for (a–e) $Ri_0 = 2$, (f–j) $Ri_0 = 6.3$, (k–o) $Ri_0 = 20$, (p,q,r) $Ri_0 = 63$, with (a,f,k,p) $Re_0=100$, (b,g,l,q) $Re_0=178$, (c,h,m,r) $Re_0=316$, (d,i,n) $Re_0=562$, (e,j,o) $Re_0=1000$.

At ${Ri}_0=2$ in figure 12, almost the entire $r$–$z$ plane shows a consistent puffing frequency that is slightly below the expected value ${St}_0^{*}$. This is perhaps not surprising given that in figure 3, the flow is laminar for $z<1.5R_0$. The regions of dark blue and red are harmonics of ${St}_0$. At ${Ri}_0=6.3$, the region ${St}_0 \lesssim {St}_0^{*}$ becomes more confined for increasing ${Re}_0$, with values of ${St}_0$ tending towards $St_0^{*}$ as ${Re}_0$ increases. A similar variation with ${Re}_0$ is observed for ${Ri}_0=20$, except that in this case ${St}_0$ increases above ${St}_0^{*}$, most notably for ${Re}_0=316$ and ${Re}_0=562$. For ${Re}_0=100$ and ${Re}_0=178$ at ${Ri}_0=62$, the region that is close to ${St}_0^{*}$ is very small, even though the flow fields in figure 3 indicate the presence of pulsatile flow.

It is interesting to note that figure 12 shows that a low-frequency region forms near the centreline at the base of the plume for ${Ri}_0=20$, ${Re}_0\geq 178$, and for all simulations at ${Ri}_0=63$. Connecting this back to figure 4, these are also the simulations where RT spikes form along the centreline. Because these spikes cause large fluctuations in the vertical velocity on a time scale much larger than the puffing instability, the computed ${St}_0$ in this area is much smaller than the ${St}_0$ associated with the puffing.

To determine quantitatively the puffing frequency dependence on ${Ri}_0$ and ${Re}_0$, we apply the probe technique to the 3-D volume $x\in [-2R_0,2R_0]$, $y\in [-2R_0,2R_0]$ and $z\in [0,5R_0]$, resulting in fields of peak frequencies similar to those shown in figure 12. If the frequency at a particular spatial location satisfies $|{\log _2({St}_0/{St}_0^{*})}| < \log _2(4/3)$, then this frequency is retained for further analysis. This bound was chosen because the puffing frequency must fall within this range based on many previous measurements (Wimer et al. Reference Wimer, Lapointe, Christopher, Nigam, Hayden, Upadhye, Strobel, Rieker and Hamlington2020), and if a computed frequency falls outside this range, then it is most likely not associated with the puffing instability. Among the spatial locations that fall within this range, the most common value of ${St}_0$ (i.e. the mode) is shown in figure 13 for all simulations.

Figure 13.

Puffing frequency measured by the most commonly computed ${St}_0$ value in each simulation. Symbols correspond to different ${Ri}_0$, and colours correspond to different ${Re}_0$. (a) Puffing frequency as a function of ${Ri}_0$. (b) Puffing frequency as a function of ${Re}_0$. The dashed lines are the empirical correlation predicted by Wimer et al. (Reference Wimer, Lapointe, Christopher, Nigam, Hayden, Upadhye, Strobel, Rieker and Hamlington2020) and indicated in (1.1).

Figure 13(a) shows that the strongest variations in ${St}_0$ occur with changes in ${Ri}_0$. In particular, ${St}_0$ increases with increasing ${Ri}_0$, and the dependence on ${Re}_0$ is relatively weak. As ${Ri}_0$ increases, there is an increase in ${St}_0$ relative to the scaling proposed in (1.1). Nevertheless, the present results do lie close to this relationship, which was proposed in Wimer et al. (Reference Wimer, Lapointe, Christopher, Nigam, Hayden, Upadhye, Strobel, Rieker and Hamlington2020) based on a wide range of data sources.

Although the dependence of ${St}_0$ on ${Re}_0$ is not as strong as the dependence on ${Ri}_0$, figure 13(b) does reveal that there are variations in ${St}_0$ with ${Re}_0$. For the two highest values of ${Ri}_0$, in particular, figure 13(b) shows that ${St}_0$ reaches a maximum value at an intermediate value of ${Re}_0$ (316 in the case ${Ri}_0=20$, and 178 for ${Ri}_0=63$), but decreases for both smaller and larger values of ${Re}_0$. The dependence of ${St}_0$ on ${Re}_0$ is less pronounced for the two smaller values of ${Ri}_0$,

Based on the spatially varying puffing frequency fields shown in figure 12, there are several points to note regarding the single values of the puffing frequency shown in figure 13. First, as ${Ri}_0$ increases, figure 12 shows that the spatial region where the probe can be used to measure the puffing frequency narrows substantially. Additionally, within this region, there is also some uncertainty in the precise value of the puffing frequency. For example, at ${Ri}_0=20$ and ${Re}_0=1000$, as well as at ${Ri}_0=63$ and ${Re}_0=316$, there are different peak frequencies all within the same spatial region; any of these frequencies could be associated reasonably with the puffing frequency. These variations indicate greater uncertainty in the measurement of a single global puffing frequency associated with the plume dynamics.

Taken together, figures 12 and 13 suggest the following dependence of the puffing frequency on ${Ri}_0$ and ${Re}_0$. At low ${Re}_0$, the puffing frequency is consistently predicted below the scaling proposed in (1.1). As ${Re}_0$ increases, there is more uncertainty in the computed ${St}_0$, but the estimated ${St}_0$ approaches the proposed scaling. This trend is most clear for ${Ri}_0=20$: for ${Re}_0=100\unicode{x2013}178$, ${St}_0$ is below the proposed scaling; for ${Re}_0=316\unicode{x2013}562$, ${St}_0$ is above the scaling; and at ${Re}_0=1000$, the most commonly computed ${St}_0$ is near the proposed scaling, although with some uncertainty (based on the variability shown in figure 12).

Finally, it should be noted that the slight dependence of the puffing frequency on ${Re}_0$ observed here does not contradict previous works (Subbarao & Cantwell Reference Subbarao and Cantwell1992; Bharadwaj & Das Reference Bharadwaj and Das2017), where a dependence on ${Re}_0$ has been neglected or was not observed. Based on our simulations, we start to see the increase in ${St}_0$ with ${Re}_0$ roughly where the flow transitions from laminar to turbulent, implying that this transitional regime induces a Reynolds number dependence. Additionally, it was shown experimentally by Cetegen & Kasper (Reference Cetegen and Kasper1996) that there is a change in the scaling of ${St}_0$ at high Richardson numbers. However, these high ${Ri}_0$ values were achieved by varying the inlet velocity $U_0$, thus ${Re}_0$ was varying simultaneously. From our analysis of the data presented in figure 13, if the critical ${Re}_0$ at which transition occurs continues to decrease with increasing ${Ri}_0$, then we expect that the data shown in Cetegen & Kasper (Reference Cetegen and Kasper1996) would be near this critical regime where the frequency was found here to be dependent upon ${Ri}_0$. Of course, in order to make the connection precise, many more simulations with different ${Ri}_0$ and ${Re}_0$ would need to be conducted, but this Reynolds number dependence offers a possible explanation for the change in scaling observed previously at large ${Ri}_0$.