Multi-scale analysis of respiratory droplets transport within the breathing cloud

Infectious diseases are transmitted primarily through pathogen laden droplets commonly exhaled during respiratory events such as breathing, coughing, or sneezing. The transport and evaporation of droplets are governed by the fundamental laws of fluid mechanics and convection-diffusion. From these, analytical models can be created helping in better understanding pathogens transmission from a mechanical perspective. The droplet transport within the humid air breath cloud and the local ventilation are crucial for accurately predicting the droplet fate. Different levels of complexity are possible for modelling the breath cloud, from simple 1D models to complex unsteady 3D simulations including discrete phase models for the droplets simulation. The former are too simple to capture the fluid dynamics of intermittent jets in detail, while the latter are too computationally expensive. The current work presents a novel multi-scale approach where an analytical model of the droplet transport and evaporation is coupled to unsteady CFD simulations of warm humid puffs of exhaled air. The proposed model has the advantage of the accuracy of an analytical model and the computational cost of a relatively standard unsteady CFD simulation, and can be used to predict the trajectory of the exhaled droplets for a variety of respiratory events.


Introduction
The outbreak of SARS-CoV-2 has resulted in significant changes to our daily lifes, including the use of face masks, regular hand and surface disinfections, and the need to maintain social distance.World Health Organization (WHO) recommended a social distance of at least 1 meter, but many researchers have shown how in most cases this distance is not necessarily enough to prevent infection.Respiratory droplets can range from approximately a millimetre down to a fraction of micrometre in size and can either deposit on surfaces or remain airborne, giving rise to different transmission routes, either by contact with an infected surface, or by droplet inhalation on behalf of a susceptible host.An important role is played by the ambient conditions in terms of relative velocity between the droplet and the surrounding air, and water vapour partial pressure.This is particularly true for smaller droplets that evaporate rapidly and whose inertia is negligible so that they can be carried arbitrarily far by the local air flow, dangerously contributing to the chance of airborne transmission, and easily breaking the one metre social distancing rule.In this framework, it is apparent how the droplet transport within the humid air breathing cloud and/or within the ventilation air flow in a confined environment becomes essential for the correct prediction of the droplet fate.While for the latter case results cannot be generalised easily as much depends on the room and on the specific ventilation system layouts, for the former case general models of the exhaled air cloud as a sequence of jets and puffs are possible for different reference respiratory events and could be of great aid in predicting the risks associated to short-distance direct inhalation airborne transmission route.
In the last few years several studies investigated the effects of indoor ventilation on the droplets trajectory.Ji et al. [1] showed how mixing ventilation -due to the augmented turbulent heat and mass transfer -induce a faster droplet evaporation with respect to displacement ventilation.Motamedi et al. [2] performed CFD simulations to evaluate the risk of infection with different types of indoor ventilation.The study shows how the risk, depends on both the position of the infected individual and the type of ventilation, whereas no ventilation strategy can be clearly deemed better than others.In all the cases considered, it has been demonstrated that social distancing alone is not enough to contain indoor virus spread.Other works focused on the breath cloud.Xie et al. [3] used an empirical model of non-isothermal jet to model the breath cloud and estimate the falling distance for droplets of different size in a quiescent environment.It was shown that droplets can be carried even more than six metres away in case of sneezing.Bourouiba et al. [4] after a series of experiments on human coughs and sneezes, highlight how respiratory events should not be modelled as steady jets but rather as sequence of turbulent jets and puffs.An analytical model based on momentum conservation is proposed.Other important aspects that influence the disease transmission in indoor enviroments are the ambient thermohygrometric conditions: Bahramian et al. [5] show how a higher indoor temperature results in an increase of aerosolized droplet lifetime, and hence an increased airborne transmission.Busco et al. [6] validated a series of experimental results on sneezes with a CFD model accounting for the angular head motion and the dynamic pressure response during the sneeze.Their model is compared to simpler sneeze models having costant head angle and ejection velocity.Their results show a much larger exhaled droplets spatial spread.
In the literature a rather net distinction is found between analytical and CFD models.The latter have superior breath cloud modelling capabilities, while the former, being characterised by a continuous approach, provide more accurate trajectory and evaporation predictions.This calls for a multi-scale approach that could exploit the advantages of the two methods.
In a previous work, the authors developed an analytical model of droplet transport and evaporation [7], also addressing the role of its non-volatile components, such as salts and enzymes, on the droplet evaporation.An elaborated, although simplified, numerical model of the buoyant breath cloud as a sequence of jets and puffs was included and provided the local velocity field and thermo-hygrometric boundary conditions required by the analytical model.
In this work, a multi-scale approach is proposed and preliminarily tested where the aforementioned breath cloud model is replaced by unsteady CFD simulations of intermittent warm humid air jets.In this way, the intrinsic limits of the 1D breath cloud model can be overcome, thus making the model prediction more accurate, yet without the need of computationally expensive Discrete Phase Model (DPM) CFD simulations.With the model proposed, the trajectory of exhaled droplets can be predicted for any kind of respiratory event and ambient conditions.As such, the model could be used to assess the risks associated to different pathogen transmission routes under different real life scenarios.Preliminary results for the case of periodic mouth breathing in a confined environment representing a room at a given temperature and relative humidity are shown and discussed.

Analytical model
A brief overview of the mathematical model, that describes the behaviour of the exhaled droplets, is now given.In the following, it is assumed that the droplets have spherical shape and uniform 40th UIT International Heat Transfer Conference (UIT 2023) Journal of Physics: Conference Series 2685 (2024) 012066 thermophysical properties.

Transport
The droplet trajectory is computed from the balance of the forces acting on the droplet itself as it moves through the air.These are its weight (including buoyancy) and drag, while other forces are assumed negligible [8].From Newton's second law the droplet acceleration can be written as where where x and z are horizontal and vertical directions respectively, u r =u p −u m is the relative velocity between the droplet and the air stream, τ d the characteristic time associated to drag calculated as the ratio between the vertical momentum of the droplet and the vertical component of the drag force, and the subscript 't' denotes the terminal condition.By terminal the zero acceleration condition the particle will attain in the limit is meant.If a spherical droplet with initial relative velocity u r,0 falling in a moist air stream having velocity u m is considered, by integrating equation ( 1) For the case of Stokes' flow, i.e. when the Reynolds number is low, drag force is a linear function of the relative velocity, as a consequence τ d is constant and equal to τ d,t .In case of non Stokes' flow, instead, τ d changes during the droplet flight and tends to τ d,t .By further integrating equation ( 2) with initial position s p,0 , the droplet trajectory results From equation (2) it can be noted that as time passes, the influence of the initial condition on the droplet velocity vanishes.In addition, in case of small terminal velocity the particle is entirely transported by the local air flow, and therefore can travel over arbitrarily long distances.

Evaporation
Considering convection-diffusion equation in stationary form without source or sink terms and neglecting cross-diffusion of air in water, the mass transfer from the droplet surface to the surrounding moist air due to evaporation can be written in terms of the droplet diameter reduction rate in time as where the subscripts 's' and '∞' indicate that the related quantity is to be evaluated by the droplet surface or far from it in the surrounding air, and the subscripts 'a' and 'v' denote dry air and water vapour respectively.D p is the droplet diameter, Sh the Sherwood number (that can be evaluated using Ranz-Marshall empirical correlation), p the pressure, ρ p the droplet density, R v the water vapour ideal gas constant, and T the temperature.Finally, D v is the mass diffusivity of water in air that can be evaluated with the correlation proposed in [9] D v = 21.By introducing the evaporation rate K defined as separating the variables, and integrating equation ( 6) in time from an initial droplet diameter D p,0 , is it possible to write the droplet diameter as a function of time In general K is not constant during the droplet evaporation, yet it can be shown that for small diameters K tends to a constant value, function of the particle density and the ambient thermohygrometric conditions alone.

Energy
When a droplet is released in the air, its evaporation starts immediately.At the beginning the droplet is gradually cooled (or warmed) from its initial temperature to the ambient wet bulb temperature.At a later stage, the droplets keeps shrinking at constant temperature.In case the droplet contains a non-volatile fraction, the evaporation slows down as it progresses according to Raoult's law, and the droplet temperature gradually grows up to when thermal equilibrium with the ambient is met and evaporation terminates.
In terms of thermal energy balance, at any time the sensible thermal power stored by the particle is balanced by evaporative cooling and convective heating Qsn = Qev + Qcv .
Expanding these terms and recasting equation (8) in terms of the droplet temperature rate of change is found, where Nu is the Nusselt number, λ the thermal conductivity of the moist air, c p the droplet specific heat, and h lt the water latent heat of vaporisation.To be noted that in case the droplet contains some non-volatile compounds its density and specific heat will change over time and should be computed as weighted average between the corresponding values of the solid and the liquid fractions.
The fact that the droplet temperature variation is small and limited in time allows the dependency of K upon the droplet temperature to be neglected, where T p =T s .equation ( 9) can thus be integrated assuming an initial droplet temperature T s,0 to give where τ ev is the characteristic time associated to evaporation and T s,t is the droplet terminal temperature, i.e. the ambient wet bulb temperature when the droplet is evaporating.To be noted how the solutions in equation ( 2) and in equation ( 10) are of the same form.

Methods
From the analytical model outlined above, it is apparent how the local air velocity and thermohygrometric conditions in terms of temperature and relative humidity play a fundamental role on the droplet transport, evaporation, and thermal equilibrium.These quantities, for instance, are found in the definitions of the characteristic times and terminal quantities in equation ( 2) and in equation ( 10), and largely affect the partial pressures in equation ( 6).If the local air properties were constant and uniform things would be easier, but this is not the case for respiratory droplets.These in fact are exhaled within the breath cloud, a warm and humid buoyant air intermittent puff that gradually mixes with the colder ambient air, so that the ambient conditions seen by the droplet change continuously, at least as long as the droplet remains trapped inside the breath cloud.A mean is needed to retrieve the required information on the ambient air and feed it to the analytical model.In the following a preliminary coupling between the analytical model and an unsteady CFD simulation is introduced and discussed.

Case study
The case study for the CFD analysis consists of a simplified 2D open room geometry, 6 m long and 4 m high, where a breathing mannequin is placed.The breathing process is modelled as a time-periodic inlet/outlet boundary having period of 4 s and composed of a 2.4 s long exhalation followed by a 1.6 s inhalation [10].The flow is normal to the mouth surface, a circle of 2 cm in diameter, and its velocity varies sinusoidally whereas the period of the positive and negative halves of the sinusoid is not the same as mentioned (see figure 1b).Maximum positive velocity during exhalation is 1.042 m/s while minimum negative velocity during inhalation is −1.563 m/s, so that a tidal volume of 0.5 l per breath is obtained.

CFD model
The CFD model is built in OpenFOAM.A multiphase solver is adopted to account for moisture in the air, together with k-ω SST turbulence model.For the preliminary results presented in this work a 2D structured mesh was chosen.To minimise the impact of the boundaries on the results, the outlet patch (named open ceil in figure 1a) is located sufficiently far from the area of interest, at a height of 8 m from the floor.Slip condition is imposed at all solid boundaries but the floor and the mannequin.At the mouth inlet the velocity function in figure 1 is imposed, inlet temperature and relative humidity being 307.15K and 95% respectively.A fixed temperature is imposed on the mannequin while the remaining solid patches are assumed adiabatic.The mannequin surface temperature was estimated using the formula T clo ≈T skin +q ′′ skin /(0.155I clo ) where the insulation coefficient (I clo ) for coverall clothes and the body heat flux for an individual standing at rest (q ′′ skin ) 40th  [11].Assuming I clo =0.72 clo, T skin =36 • C, q ′′ skin =70 W/m 2 the clothes surface temperature is T clo ≈28 • C. Ambient initial conditions are still air at 293.15 K and 50% relative humidity.Since the breath cloud and the body surface are warmer than the surroundings, both cloud buoyancy and thermal plume by the body occur.
The total simulated time is 50 s, and data are saved every 0.1 s.Table 1 displays the results of the mesh sensitivity test.The first column reports the mesh elements number, while the second gives the minimum and maximum diagonal cell size in the area of interest, where it must be considered that the mesh is finer close to the mouth and gets coarser with the distance.The remaining columns show the mean values, for a coupled simulation of a 10 µm droplet, of the distance travelled by the droplet (d), and the local temperature (T ) and relative humidity (RH) encountered at three instants of time taken at the end of the first three exhalations (i.e. at 2.4, 6.4, and 10.4 s).The intermediate mesh is chosen as a reasonable trade-off between accuracy and computational time.

Coupling
The CFD simulation is run once before the analytical model.When all the results are saved, the analytical model is launched to simulate a series of droplets of different size.The boundary conditions for the analytical model is limited to the droplet initial position, velocity, temperature, size, and solid fraction properties, whereas all the information related to the ambient and the breath cloud is taken from CFD.
At each time step the analytical model provides the simulated time and the droplet coordinates to the coupling script.Since the current time and position will not match one of the CFD saves and mesh cell centres, interpolation will be needed both in space and time: the CFD results of the two closest saves in time are read and queried using the probe filter of the VTK module [12].This filter locates the closest cell centres to the query point and takes care of the interpolation returning the desired ambient air values.This information is then linearly interpolated by the coupling script before being fed back to the analytical model where all the remaining simulation parameters are updated.These include the droplet position, velocity, diameter and temperature with their rate of change, chemical composition, characteristic times, and forces, together with the thermophysical properties of both droplet and ambient air.
Figure 2 shows the logical flow of the multi-scale method proposed.Of course, the accuracy of the method is proportional to the number of CFD saves, and the size of the grid.

Results
Figure 3 compares the results of the previous analytical model developed by the authors with those of the current multi-scale approach in terms of droplet trajectory and evaporation progress.Droplets of three different sizes are simulated: namely, 10 µm, 50 µm, and 100 µm, so that the latter will be mostly governed by gravity, the first by drag, and the remaining, to some extent, by   both.Simulations either stop when the droplet reaches the ground, or exits from the domain, or remains airborne (i.e. the thermal equilibrium is reached out of the breath cloud area of influence while the droplet relative velocity becomes sufficiently small), or the end time is reached.Simulation times can be rather different depending on the droplet size.On average, the analytical model requires ≈1 minute of CPU time per simulated droplet on an ordinary workstation equipped with Intel Core i7-2600 CPU at 3.4 GHz, while the multi-scale model takes ≈5 minutes in addition to the time required by the preliminary CFD run.The particle trajectories in an equivalent DPM CFD analysis are highly dependent on the mesh size and the time step, thus implying much larger computational costs to achieve a comparable accuracy with respect to the multi-scale model.
While the agreement is remarkable for what concerns the droplet evaporation, the trajectories computed by the two methods are rather different.The more realistic trajectory of the multiscale approach is apparent for the droplet in figure 3a which due to its small size is largely carried around by the vortices generated in the breath cloud; vortices that cannot be modelled analytically.Even though, both methods predict that the droplet is brought upward by the buoyant breath cloud, the path and the distance travelled are much different, and so the risks of virus transmission.Similarly, the droplet in figure 3b remains lingering in front of the exhaler before being inspired again.This particular behaviour cannot be caught by the analytical approach where inspiration is not modelled.Last, the droplet in figure 3c is caught by entrainment vortices suddenly reversing its course after crossing the jet shear layer, it is then attracted towards the no-slip mannequin surface due to the Coandȃ effect, reaching the ground before complete evaporation occurs.Also in this case, the difference in the predicted trajectory is due to physical phenomena that a simple analytical model cannot handle.

Conclusions
Exhaled droplets pose a potential risk for infectious diseases spreading.Larger droplets, with a diameter of 100 µm or more follow a parabolic trajectory governed by gravity, while smaller droplets are mainly transported by drag by the local airflow.Several works in the literature have investigated the role of the breath cloud on the droplet transport, either developing analytical models or performing complex DPM CFD analysis.
The weak point of analytical models is clearly the difficulty in implementing a simple but reliable breath cloud model.In this work, a multi-scale approach has been proposed to overcome this difficulty by means of a one-way coupling between an analytical model previously developed by the authors, and a CFD simulation of an unsteady breath cloud in a large room carried out with OpenFOAM; the local air velocity and thermohygrometric conditions needed by the analytical model are thus provided by the CFD analysis.This approach combines the advantages of CFD modelling capabilities in capturing the jet dynamics more in detail while dropping the need for DPM, and that of continuous analytical model accuracy in predicting the droplet transport and evaporation.
Preliminary results are presented where the multi-scale approach is compared to the purely analytical model for the case of three droplets of different diameters.While the two approaches agree very well on the prediction of the droplet evaporation, things are rather different on transport prediction, where the role of turbulence within the breath cloud or in the jet shear layer may have a dramatic impact over the droplet fate.
To the author's knowledge, such a multi-scale approach has never been investigated in the literature, and in this preliminary analysis it has been proven to be an interesting feasible alternative to complex and computationally intensive DPM CFD analyses, or simpler analytical methods with limited modelling capabilities.
In future developments, the authors plan to implement the analytical model coupling to a 3D CFD model featuring a more realistic mannequin geometry.Investigation on possible means to further speed up the CFD query procedure will also be carried out.

Figure 1 .
Figure 1.CFD model geometry and boundary patches (a), and normal velocity component at the mouth (b).

Figure 2 .
Figure 2. Logical flow of the multi-scale method.

Figure 3 .
Figure 3. Droplets trajectory (top) and evaporation progress (bottom): multi-scale approach (solid red line) vs analytical model (dashed blue line).From left to right the images refer to the 10 µm, the 50 µm, and the 100 µm droplet, respectively.