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 infected person via coughing, sneezing, speaking or simply breathing [1]. Focusing on a cough for the sake of example, its typical duration is 200–500 ms, the typical mouth opening of male subjects is approximately , and the resulting Reynolds number is about [2,3] and even larger for a sneeze. The resulting flow field is thus fully turbulent. Very recently [4-6], the key role of turbulence in dictating the fate of virus-containing droplets in violent human expulsions has been elucidated. More specifically, it turns out that turbulence plays a crucial role in determining droplets evaporation time, resulting in errors up to 100% when the turbulent fluctuations are filtered or completely averaged out [4]. Because of the fact that the issue on the droplet evaporation time is crucial to establish whether viruses lingering on dry nuclei upon droplet evaporation retain their full potential of infection [5], turbulence in a puff needs to be fully understood and characterized in a statistical sense.

As a matter of fact, a consistent statistical theory is currently not available for puff turbulence. The reason is probably due to the fact that the system at hand does not possess all desirable symmetries characterizing the classical ideal turbulence setting. Current knowledge of puff turbulence is confined to the pioneering work by Kovasznay et al. [7] only dealing with bulk quantities involved in the puff dynamics: the puff bulk translational velocity and the average puff radius.

A first attempt to characterize from a statistical point of view the small-scale structure of turbulent fluctuations of both velocity and temperature, but neglecting the two-phase nature of the problem, is reported in [8] where an adiabatic generalization of the Kolmogorov–Obukhov picture of steady Navier–Stokes turbulence [9,10] has been presented in two relevant cases: one in which the buoyancy of the puff is negligible and the opposite situation where buoyancy dominates the puff dynamics. It is shown that the phenomenological theory gives sense to the concept of inertial range, energy flux and scaling behaviours, in space and time, of proper statistical observables. A similar approach have been successfully used also in describing inertial range scaling laws in decaying homogeneous anisotropic turbulence [11].

The proposed theory is still waiting for verification (in terms of accurate direct numerical simulation and/or experiments) in the second limit dominated by buoyancy. In [8], only the case of negligible buoyancy and the intermediate regime where buoyancy is just starting to play a role have been analysed. Here, we fill the gap by analysing, in terms of state-of-the-art direct numerical simulation, the buoyancy-dominated regime with the aim of supporting the theoretical prediction of [8].

The paper is organized as follow. In §2, we describe the governing equations governing a turbulent puff. In §3, we propose a phenomenological model to describe the spatial and temporal evolution of the turbulent puff, both in terms of bulk properties (mean velocity, temperature and size) and velocity differences in the inertial and viscous range of scales; in order to test the validity of the theory, we perform direct numerical simulations, whose details are provided in §4. The comparison of theory and simulations is reported in §5, and the main conclusions are summarized in §6.

Problem setting and ruling equations

We consider a fluid impulsively ejected in an undisturbed environment from a small circular opening of radius , as shown in figure 1. The fluid motion is governed by the incompressible Navier–Stokes equations for the velocity and pressure 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 , and being the (constant) ambient temperature and density, respectively, and the thermal expansion coefficient. In the previous equations, is the gravitational acceleration, the thermal diffusion coefficient and () the kinematic fluid viscosity. In the following, we denote as the difference between the fluid temperature and the ambient one: . 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 in the momentum equation.

Figure 1

Side view and cross section of the turbulent puff studied in the present work, together with the reference coordinate axis. (Online version in colour.)

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-of-the-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 , 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 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 , animated by a typical velocity , follows from the momentum conservation law of the puff, ; 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: and thus . Here, and 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 . The scaling laws of the puff bulk properties thus follow from simple power counting. Once the scaling behaviour of is known, one can easily predict the power-law decay for the large-scale temperature difference, , between the puff and the cooler ambient. is indeed proportional to the reciprocal of the cloud volume , from which: . 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 where is the air thermal expansion coefficient and is the gravitational acceleration. For , one has . Because the eddy viscosity term goes as , soon or later buoyancy dominates (i.e. its decay is slower) over the eddy viscosity. This happens provided that . For , the acceleration on the left-hand side of the Navier–Stokes equations now balances the buoyancy term: 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 large-scale 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 . Exploiting the scaling behaviours (3.1) one immediately gets and where is the energy flux a . 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-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, . Evaluating at in equation (3.2) the expressions for follows: Below velocity fluctuations are smooth, i.e. , and the resulting scaling behaviour for is

A similar way of reasoning for the two-point temperature fluctuations, , leads to the adiabatic generalization of Obukhov–Corrsin theory (OC51) of passive scalar advection [19,20]: where is the flux of scalar variance at . The above scaling law is accompanied by the scaling relation within the temperature dissipative range (i.e. for scales , 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-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 are thus corrected by a multiplicative factor , for the velocity, and for the temperature. Note that because of the temporal dependence encoded in the intermittency corrections change both the spatial and the temporal structure-function scaling laws. The values of and 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 , , and . The resulting model for the structure functions of the (longitudinal) velocity and the temperature difference is thus and where and 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 , 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 . 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 the flow solver Fujin, an in-house code, extensively validated in a variety of problems [27-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 neighbouring ones (40 preceding and 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 () and one where its effect is secondary (red symbols) with (). As shown in the figure, while in the former the puff starts raising due to buoyancy soon after the ejection phase ends at , the latter instead initially maintains an almost straight trajectory, and only at later times when starts rising. We can thus safely consider the two cases as good representative of the buoyancy-dominated and shear-dominated regimes.

Figure 2

Time history of the vertical displacement of the puff centre of mass. The red circles and blue squares represent the data from the case with low and high buoyancy. (Online version in colour.)

For the two configuration above, we start reporting in figure 3 the bulk property of the puff evolution, namely its bulk size , the streamwise velocity and the temperature difference . 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 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 , while those for the buoyancy-dominated regime for . 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 from which we obtain in our case . Thus, we conclude that up to the viscous contribution is negligible at large scales and the momentum is conserved with good accuracy. Note that, after , we do not make any assumption on momentum conservation.

Figure 3

Scaling laws for the bulk properties of the puff: (filled symbols), (empty symbols) and (half-filled symbols) (divided by a factor two for graphical reasons) for (a) shear-induced fluctuations and (b) buoyancy-driven fluctuations. The data reported in the left panel are taken from [8]. In the figures, the solid lines represent the proposed scaling laws, while the symbols the results of our simulations. (Online version in colour.)

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 and the Rayleigh number . 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 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 4

Scaling laws for the non-dimensional numbers of the puff: (a) and (b) for shear-induced fluctuations (red circles) and buoyancy-driven fluctuations (blue squares). In the figures, the solid lines represent the proposed scaling laws, while the symbols the results of our simulations. (Online version in colour.)

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. . 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).

Figure 5

Second-order structure function for the (a,b) velocity and (c,d) temperature difference: (a,c) raw data and (b,d) data scaled by the predicted inertial scaling. In the figures, the solid lines represent the proposed scaling laws, while the symbols the results of our simulations. The data refer to the time interval , with colours indicating different time instants separated by an interval of in the following order: dark-blue, blue, grey, light-green, green and brown. (Online version in colour.)

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 in each of them ( for the viscous range and 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. , , and , , , 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.

Figure 6

Time histories of (a,b) (filled symbols), (empty symbols) and (half-filled symbols) and of (c,d) (filled symbols), (empty symbols) and (half-filled symbols) for two separations taken in the inertial (a,c) and viscous (b,d) range of scales, for the case with buoyancy-driven fluctuations. In the figures, the solid lines represent the proposed scaling laws, while the symbols the results of our simulations. (Online version in colour.)

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-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. ), 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 prediction; the clear departure of these lines from our data indicates both the importance of including intermittency corrections in the model, and the validity of our approach to do so [41].

Figure 7

(a,b) (left) and (right) as a function of and (c,d) (left) and (right) as a function of . All the structure functions are divided by their expected temporal scalings. The symbols represent the data from the case with high buoyancy, while the solid and dotted lines represent the ratio of the expected spatial scaling laws with and without the intermittency corrections. (Online version in colour.)

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 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 or . 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 buoyancy-dominated 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 spatio-temporal 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 expulsions [4]. Another issue left for the future is the application and extension of the present theory to a multiphase problem.