Puff turbulence in the limit of strong buoyancy

We provide a numerical validation of a recently proposed phenomenological theory to characterize the space–time statistical properties of a turbulent puff, both in terms of bulk properties, such as the mean velocity, temperature and size, and scaling laws for velocity and temperature differences both in the viscous and in the inertial range of scales. In particular, apart from the more classical shear-dominated puff turbulence, our main focus is on the recently discovered new regime where turbulent fluctuations are dominated by buoyancy. The theory is based on an adiabaticity hypothesis which assumes that small-scale turbulent fluctuations rapidly relax to the slower large-scale dynamics, leading to a generalization of the classical Kolmogorov and Kolmogorov–Obukhov–Corrsin theories for a turbulent puff hosting a scalar field. We validate our theory by means of massive direct numerical simulations finding excellent agreement. This article is part of the theme issue ‘Scaling the turbulence edifice (part 2)’.


Introduction
Turbulence in a fluid puff is a problem recently attracting a great deal of attention in relation to the global emergency caused by the COVID-19 pandemic. Its relevance stems from the key role of airborne virus transmission by viral particles released by an 2022 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Problem setting and ruling equations
We consider a fluid impulsively ejected in an undisturbed environment from a small circular opening of radius R, as shown in figure 1. The fluid motion is governed by the incompressible Navier-Stokes equations for the velocity u and pressure p fields, coupled to the advection-diffusion equation for the temperature field θ, i.e. by the Oberbeck-Boussinesq (OB) equations [12]:  Here, we assumed the fluid density ρ to be linearly dependent on the temperature θ : ρ(θ) = ρ a [1 − β(θ − θ a )], θ a and ρ a being the (constant) ambient temperature and density, respectively, and β the thermal expansion coefficient. In the previous equations, g = (0, 0, −g) is the gravitational acceleration, κ the thermal diffusion coefficient and ν (∼ κ) the kinematic fluid viscosity. In the following, we denote T as the difference between the fluid temperature and the ambient one: T = θ − θ a . Note that the momentum equation for the velocity field and the temperature advection-diffusion equation are fully coupled; indeed, the temperature field does not behave as a simple passive scalar field advected by the velocity field, but can alter the flow field itself by the back-reaction term βg(θ − θ a ) in the momentum equation.

Phenomenological theory for bulk and two-point puff statistics
In this section, we provide a quick review of known results for both bulk puff properties (e.g. typical puff velocity and typical puff size) and two-point observables characterizing the mixing zone of the puff. The results we are going to review for the bulk properties in the limit of vanishing buoyancy are known since the seminal paper [7]. The remaining results have been presented only very recently by the present authors [8] and are still waiting for a conclusive confirmation from numerical simulations and/or experiments. Providing such a confirmation in terms of state-ofthe-art direct numerical simulation is the main concern of the present paper.
Let us consider a fluid cloud impulsively ejected from a spatially localized source (e.g. at the origin of a three-dimensional reference system) in an otherwise undisturbed environment (e.g. air). The injection process is supposed to last until time t 0 , after which the cloud freely evolves into the ambient under the combined action of gravity and shear. In the freely evolving regime, the fluid cloud is named puff, to distinguish it from the initial jet phase [13]. Let us denote by T 0 > 0, the difference between the cloud temperature and the ambient temperature at the end of the jet phase.
When the buoyancy effect is negligible, the temporal behaviour of the typical size of the puff, denoted by L(t), animated by a typical velocity u L ∼ L/t, follows from the momentum conservation law of the puff, L 3 u L = const; note that, the typical translational velocity of the puff is determined in a Lagrangian frame of reference, following the point associated with the maximum velocity of the puff: L(t) ∼ L 0 (t/t 0 ) 1/4 and thus u L (t) ∼ u 0 (t/t 0 ) −3/4 . Here, L 0 and u 0 = L 0 /t 0 are the typical cloud size and cloud velocity at the end of the jet phase. These scaling behaviours have been obtained in [7] from the large-scale Navier-Stokes equations. These latter indeed reduce to a diffusive equation with a time-dependent eddy viscosity ν T (t) ∼ (t/t 0 ) −1/2 L 2 0 /t 0 . The scaling laws of the puff bulk properties thus follow from simple power counting. Once the scaling behaviour of L(t) is known, one can easily predict the power-law decay for the large-scale temperature difference, T L , between the puff and the cooler ambient. T L is indeed proportional to the reciprocal of the cloud volume L 3 , from which: T L ∼ T 0 (L/L 0 ) −3 = T 0 (t/t 0 ) −3/4 . Note that the theory discussed above is based on a time-dependent eddy viscosity closure and predicts the scaling behaviour for the (sole) bulk properties. As a results, we found the same scaling laws for the typical radius of the puff and for the distance from the origin travelled by the puff. It then follows that, apart prefactors, lengths in a puff obey the same scaling behaviours and the same applies also for the typical fluctuating velocities inside the puff and the mean velocity of the puff. They are indeed dimensionally related to the typical radius of the puff and to the distance it travels from the origin.
The eddy viscosity description obtained in [7] makes it possible to generalize the above scaling behaviour to the buoyancy-dominated case. Indeed, the order of magnitude of the buoyancy contribution to the puff evolution is b = βgT L where β is the air thermal expansion coefficient and g is the gravitational acceleration. For T L ∼ T 0 (t/t 0 ) −3/4 , one has b ∼ βgT 0 (t/t 0 ) −3/4 . Because the eddy viscosity term goes as t −7/4 , soon or later buoyancy dominates (i.e. its decay is slower) over the eddy viscosity. This happens provided that t t b ≡ u 0 /(βgT 0 ). For t t b , the acceleration on the left-hand side of the Navier-Stokes equations now balances the buoyancy term: u L /t ∼ u 2 L /L ∼ βgT 0 (L/L 0 ) −3 originating the following new scaling laws [8]: Let us now move to analyse the two-point statistics and assume that the bulk quantities, we have identified serve as large-scale properties of the flow, below which a cascade process may originate according to a generalized Kolmogorov-Obukhov picture. The idea is that, despite the fact that the problem at hand is not stationary, and thus very far from the classical playground of the Kolmogorov theory, small-scale turbulent fluctuations can rapidly relax to the slower largescale dynamics. This is the essence of the adiabaticity hypothesis leading to a generalization of the classical Kolmogorov theory for a turbulent puff. We refer to [8] for the phenomenological theory when the buoyancy is negligible. Here, we focus our attention on the sole regime where buoyancy dominates and the bulk puff properties are ruled by the scaling laws (3.1). The starting point is to assume the existence of an inertial range of scales characterized by a constant flux of energy given by (t) ∼ u 3 L /L ∼ (δ r u) 3 /r. Exploiting the scaling behaviours (3.1) one immediately It is worth observing that to get equation (3.2) we have supposed that buoyancy is only important to dictate large-scale balance, being instead negligible within the inertial range. A similar scenario has been found to hold in other convective systems as, e.g. Rayleigh-Taylor turbulence [14][15][16][17][18].
The energy flux is supposed to persist up to the generalized Kolmogorov viscous scale, η: the scale at which the inertial self-advection matches the viscous term in the Navier-Stokes equations. Namely, δ η u η ∼ ν. Evaluating δ r u at r = η in equation (3.2) the expressions for η follows: Below η velocity fluctuations are smooth, i.e. δ r u ∼ (r/η)δ η u, and the resulting scaling behaviour for r η is A similar way of reasoning for the two-point temperature fluctuations, δ r T, leads to the adiabatic generalization of Obukhov-Corrsin theory (OC51) of passive scalar advection [19,20]: where ε 0 = u 0 T 2 0 /L 0 is the flux of scalar variance at t = t 0 . The above scaling law is accompanied by the scaling relation within the temperature dissipative range (i.e. for scales r r d , r d being the dissipative scale coinciding with η because here ν ∼ κ): The scaling behaviours (3.2) and (3.5) fix, via simple power-counting, the mean-field predictions for the equal-time velocity and temperature inertial range structure functions. To account for intermittency corrections of both velocity [21] and temperature fluctuations [22][23][24][25], a further assumption is needed. Accordingly, in [8], it has been postulated that a turbulent puff possesses the same spatial scaling exponents as those of the stationary, homogeneous and isotropic turbulent system hosting a passively behaving scalar. This was expected to hold coherently with the adiabaticity hypothesis. The mean-field predictions for the structure functions of order p are thus corrected by a multiplicative factor [r/L(t)] −σ p , for the velocity, and [r/L(t)] −ξ p for the temperature. Note that because of the temporal dependence encoded in L(t) the intermittency corrections change both the spatial and the temporal structure-function scaling laws. The values of σ p and ξ p used are those of the stationary, homogeneous and isotropic turbulence advecting a scalar field [8] taken from state-of-the art numerical simulations of homogeneous isotropic turbulence [26]; in particular, we use σ 4 = 0.06, σ 6 = 0.27, ξ 4 = 0.24 and ξ 6 = 0.37. The resulting model for the structure functions of the (longitudinal) velocity and the temperature difference is thus and where A and B are unknown, non-universal (i.e. dependent on all details of the system) prefactors.The predictions for the structure functions within the viscous/dissipative range of scales follow from (3.4) and (3.6) by simple power counting without any intermittency correction.
In the above discussion, we assumed that the coupling between velocity and temperature is negligible everywhere in the inertial range of scales, except at the largest scales where buoyancy fixes the balance between inertia and buoyancy (i.e. it forces the velocity fluctuations).
Our assumption led to a Kolmogorov scenario which is consistent with the hypothesis that temperature is passive within the inertial range of scales. If one assumes that velocity and temperature are scale-by-scale balanced within the inertial range of scales, a new scaling behaviour would emerge (scaling a la Bolgiano), which is however not observed in our simulations. In conclusion, the fact that temperature behaves passively in the inertial range of scales is a work hypothesis which is successively verified via a consistency check.

Direct numerical simulations
In order to verify the above phenomenological theory, we rely on fully resolved direct numerical simulations. The equations of motion (2.1)-(2.3) are solved numerically within a cubic domain box of size L = 84R, R being the radius of the circular opening from which air is injected according to the time-varying profile representative of human cough proposed by Gupta et al. [3], with a peak Reynolds number of Re = u max 2R/ν ≈ 20 000. The flow is assumed to be periodic at the four sides and we prescribe a convective outlet boundary condition at the outflow boundary. We use   [27][28][29][30][31][32][33][34], based on the (second order) finite-difference method for the spatial discretization and the (second order) Adams-Bashfort scheme for the temporal discretization. See also: https://groups.oist.jp/ cffu/code.
In the performed simulations, the domain is discretized with 1260 grid points per side with a uniform spacing in all directions, resulting in a total number of 2 billion grid points. We verified the convergence of the results by comparing them with those obtained with different grid resolutions. The structure function presented below from the results of the simulations are obtained as follows: first, we identify the distance from the emission point where the maximum velocity fluctuation is present. In that plane, and in 80 neighbouring ones (40 preceding and 40 following it), we compute the velocity difference at the desired power and average the results. Note that only the points inside the puff are considered; these are selected by choosing those in which the temperature is 1% greater than the ambient one.

Results
We start our analysis by showing in figure 2 the time evolution of the centre of mass vertical position for two reference configurations: one where buoyancy is dominant (blue symbols) with b ≈ 0.5 (t b ≈ 2t 0 ) and one where its effect is secondary (red symbols) with b ≈ 0.05 (t b ≈ 14t 0 ). As shown in the figure, while in the former the puff starts raising due to buoyancy soon after the ejection phase ends at t = t b ≈ t 0 , the latter instead initially maintains an almost straight trajectory, and only at later times when t ≈ t b t 0 starts rising. We can thus safely consider the two cases as good representative of the buoyancy-dominated and shear-dominated regimes.
For the two configuration above, we start reporting in figure 3 the bulk property of the puff evolution, namely its bulk size L, the streamwise velocity u L and the temperature difference T L . In both configurations, after the initial transient phase of the ejection which is dependent on the emission profile, the bulk properties of the puff attain power law behaviours. The puff size L grows in size, and accordingly it slows down and adapt to the ambient temperature. In the figures, we report the theoretical predictions with solid lines and the simulations data with symbols. As can be appreciated, the agreement between the two is very good. Note that, the theory is not expected to hold in the early stage of the evolution, but only at later stages after the information of the initial condition is lost. In particular, from the figure, we note that the scaling laws for the shear-dominated regime are valid for t t 0 , while those for the buoyancy-dominated regime for t t b . Furthermore, to derive the scalings for the bulk properties, we assumed momentum conservation; in order to verify the validity of this assumption, we can define a viscous time as t ν ≈ L 2 /ν from which we obtain in our case t ν /t b ≈ 8000 1. Thus, we conclude that up to t b     the viscous contribution is negligible at large scales and the momentum is conserved with good accuracy. Note that, after t b , we do not make any assumption on momentum conservation.
A complete description of the bulk properties of the puff can be obtained by considering the time evolution of two non-dimensional numbers: the Reynolds number Re = u L L/ν and the Rayleigh number Ra = βgT L L 3 /(νκ). The time evolution of these quantities is shown in figure 4, together with their theoretical predictions. We observe that both theory and numeric confirm that in the shear-dominated regime the Reynolds number Re decays with time, while the Rayleigh number remains constant over time. On the other hand, in the buoyancy dominated regime both are constant and independent of time. These results further corroborates the validity of the phenomenological theory for the bulk properties of the flow. Based on these, we can now proceed to verify the predictions for the two-point statistics of the turbulent puff. Figure 5 shows the second-order structure function for the velocity and temperature differences for the time instants marked in figure 3 after the end of the jet phase, i.e. t > t 0 . As can be seen by the left plots, the magnitude of the structure function decreases as time passes, providing a scattering of the curves. Also, consistently with the bulk predictions for this regime from equation (3.1), the variation in the temperature data is larger than what observed for the velocity structure function. The time variation can be successfully compensated by scaling the data with the theoretical predictions of equations (3.7) and (3.8), as reported in the right plots. When doing so, all the data reasonably collapse onto a single curve in the inertial range of scales, while remaining scattered in the viscous one. An opposite effect could be obtained by scaling the data with the theoretical prediction for the viscous range of scales (not shown).  As an additional step towards the verification of our scaling predictions for the structure functions in the viscous and inertial range of scales, we choose a separation r/L 0 in each of them (r/L 0 ≈ 0.007 for the viscous range and r/L 0 ≈ 0.18 for the inertial range) and report the time history of the second, fourth-order and sixth-order structure functions for the velocity and temperature difference, i.e. S 2 (r), S 4 (r), S 6 (r) and S 2 (r), S 4 (r), S 6 (r), respectively. These are shown in figure 6 for the case with buoyancy-driven fluctuations. The continuous lines in the plots are the corresponding scaling laws from our theoretical predictions and we can observe in all cases an excellent agreement between the theoretical prediction and simulation data. As stated above, for this flow the intermittency correction affect not only the spatial scalings but also the temporal ones reported here, thus the good agreement of our results is also suggesting that our predictions are correctly modelling the intermittency corrections.
Finally, we verify the theoretical predictions for both the temporal and spatial scalings altogether by showing in figure 7 the data in the form of extended self-similarity [35,36]. Structure functions of any order in turbulence are expected to behave as power laws; Benzi and co-workers [35] noted that even when the Reynolds number is not really large, any two structure functions show a mutual power law behaviour, i.e. when plotted one against the other, a straight line appears with an exponent equal to the ratio between the two exponents., see e.g. [37][38][39][40]. Here, we show the fourth and sixth structure functions for velocity and temperature as a function of the second order one in the buoyancy dominated regime for several times following the end of the ejection (i.e. t > t 0 ), the same of figure 5. First, the good collapse of different time instants further corroborate the validity of our temporal scalings; next, the spatial one are verified by comparing the slope of the curve with the continuous lines with slopes deduced from our predictions, showing again a very good agreement. Note also that the dashed lines are the resulting spatial scalings that do not account for the intermittency corrections, thus representing a pure mean-field

Conclusion
A phenomenological theory for a turbulent puff has been recently presented in [8]. The main focus of the theory is the two-point statistics of both velocity and temperature. This latter field is a scalar (turbulent) field carried by the puff turbulent velocity fluctuations. The theory, based on an adiabatic generalization of the classical Kolmogorov-Obukhov theory, identifies two distinct dynamical regimes. In the first, turbulence is mechanically driven and buoyancy only has a subleading role; the situation in reversed in the second regime where buoyancy drives turbulent fluctuations. Associated with these different dynamical regime are different scaling laws, both in space and in time, for the two-point velocity/temperature structure functions. While the numerical verification has been provided in [8] for the mechanically driven puff turbulence, the verification of the second regime was still lacking. Here, we have provided such validation in terms of high-resolution direct numerical simulation. As a result of our simulations, we first confirmed the existence of an inertial range of scales at the end of which a viscous (dissipative) region takes place. A similar conclusion has been drawn for the temperature field. The resulting phenomenology appears very similar to the classical Kolmogorov-Obukhov (and Kolmogorov-Obukhov-Corrsin for temperature fluctuations) scenario with turbulence fluctuations now relaxing almost instantaneously to the time-varying flow large-scale properties, with buoyancy effects only confined to the large-scale of motion. This observation (which was a working-hypothesis in [8]) led ourselves in [8] to postulate the emergence of intermittency corrections in the present system as those known (numerically and experimentally) in the classical  homogeneous, isotropic and stationary turbulence. The validity of such model has been verified here with results clearly confirming the guess. We have also verified the space-time scaling behaviours of both velocity and temperature fluctuations resulting from the quasi-adiabatic scenario finding good agreement. Our theory is limited to the case with Pr ≈ 1 and this choice is motivated by the fact that air thermal diffusivity and viscosity are very close to each other. From a theoretical point of view, the regime we have analysed is also the most rich: indeed, in this case, an inertial range of scales for both velocity and temperature fluctuations exists, while this is not the case when Pr 1 or Pr 1. In the former case, the scalar evolves in a smooth turbulent background without intermittency, the so-called Batchelor regime, for which the scaling behaviour are reported in this work; in the latter case, despite an inertial range of scales develops for the velocity fluctuations, this is not for the scalar which exhibits a trivial bare diffusive regime.
As an interesting issue left for future research, we mention the understanding of the buoyancydominated regime in two dimensions. This regime seems to be interesting because it might be characterized by a different dynamical regime with buoyancy now balancing, scale by scale, the nonlinear terms in the Navier-Stokes equations. A similar situation is known to happen in Rayleigh-Taylor turbulence and it causes the emergence of scaling laws à la Bolgiano (instead of à la Kolmogorov) [42]. Understanding whether or not this is the case also in puff turbulence is an interesting issue left for future.
The present results can help in future modelling efforts to predict the distance reached by viral load in human cough events. Indeed, it was found that droplet evaporation is mainly controlled by the combined effect of turbulence and droplet inertia [5], and thus, the knowledge of the spatiotemporal evolution of turbulent fluctuations presented in this work can be used to obtain more realistic coarse-grained model of the original Navier-Stokes problem including the evaporation processes of small liquid droplets carried by the flow as, e.g., those involved in violent human