This paper presents and discusses the results of the “2022 International Computational Fluid Dynamics Challenge on violent expiratory events” aimed at assessing the ability of different computational codes and turbulence models to reproduce the flow generated by a rapid prototypical exhalation and the dispersion of the aerosol cloud it produces. Given a common flow configuration, a total of 7 research teams from different countries have performed a total of 11 numerical simulations of the flow dispersion by solving the Unsteady Reynolds Averaged Navier–Stokes (URANS) or using the Large-Eddy Simulations (LES) or hybrid (URANS-LES) techniques. The results of each team have been compared with each other and assessed against a Direct Numerical Simulation (DNS) of the exact same flow. The DNS results are used as reference solution to determine the deviation of each modeling approach. The dispersion of both evaporative and non-evaporative particle clouds has been considered in 12 simulations using URANS and LES. Most of the models predict reasonably well the shape and the horizontal and vertical ranges of the buoyant thermal cloud generated by the warm exhalation into an initially quiescent colder ambient. However, the vertical turbulent mixing is generally underpredicted, especially by the URANS-based simulations, independently of the specific turbulence model used (and only to a lesser extent by LES). In comparison to DNS, both approaches are found to overpredict the horizontal range covered by the small particle cloud that tends to remain afloat within the thermal cloud well after the flow injection has ceased.

## I. INTRODUCTION

The Coronavirus Disease 2019 (COVID-19) pandemic has revealed the need for a better understanding of the flow physics that govern the airborne transmission of infectious diseases via pathogen-laden aerosols, such as Severe Acute Respiratory Syndrome (SARS), human influenza H1N1, avian influenza (H5N1), and tuberculosis. In response to such a need, Computational Fluid Dynamics (CFD) tools have been massively used to predict the short-term short-range flow and particle dispersion produced by expiratory events to investigate the underlying physical mechanisms and to predict the risk of infection 1 or 2 s after exhalation. Given the complexity of the flow and the huge range of parameters to be explored, expensive Direct Numerical Simulations (DNS) are rarely affordable and Reynolds Averaged Navier–Stokes (RANS)- or Large-Eddy Simulations (LES)-based approaches are most often used. When this is the case, literature results show a significant scattering in the predictions, indicating that no general consensus has been achieved yet on the optimal (if any) modeling approach.

The outbreak of the COVID-19 pandemic in 2020 has dramatically boosted the number of research papers in which CFD is used to investigate aspects of the physical mechanisms of the virus transmission. As an example, a Scopus search conducted at the end of October 2022 for articles with the words “CFD” and “COVID” or “SARS” in the title, in the keyword list, or in the abstract yielded 306 papers with 26 entries in 2020, 99 in 2021, and 141 in 2022. Some of these articles have been collected in special issues, including “Flow and the virus” in Physics of Fluids (2021), “Numerical and experimental investigation of airborne pathogen transmission, and associated heat and mass transfer processes” in International Communications in Heat and Mass Transfer (2021), “COVID-19 and indoor environment” in Indoor Air (COVID-19 and indoor environment, 2022), “Prevention and control of COVID-19 transmission in the indoor environment” in Indoor and Built Environment (Xu and Yu, 2022), or “CFD and COVID-19” in International Journal of Computational Fluid Dynamics (Saad, 2021). It is not the objective of this introduction to exhaustively review the different studies published in the literature so the interested reader is referred directly to these special issues for more details. Also, recent reviews on the role of CFD for modeling aerosol pathogen transmission can be found in Peng (2020), Mohamadi and Fazeli (2022), Sheikhnejad (2022), and Rayegan (2022), among others.

CFD has been used to investigate the flow and the particle dispersion released during expiratory events under different situations and conditions. Simulations of the short-term, short-range dispersion have been conducted to determine the range and shapes of the flow ejection and of the particle cloud a few seconds after the end of an isolated cough or sneeze or after a sequence of violent expiratory events (Abkarian , 2020; Diwan , 2020; Pendar and Páscoa 2020; Chong , 2021; Fabregat , 2021a; 2021b; Liu , 2021; Trivedi , 2021; and Wang , 2021). These simulations usually considered the ejection of warm fluid into a quiescent colder ambient and results have been used, first, to establish specific short-term safety distances and, second, to impose initial conditions for long-term simulations of the cloud dispersion within forced or naturally ventilated spaces. Long-term, long-range dispersion in specific indoor ventilated scenarios has been also simulated to investigate the effect of the background air flow on the turbulent dispersion of the aerosol cloud in a wealth of different indoor scenarios, such as classrooms (Foster and Kinzel, 2021; Narayanan and Yang, 2021), hospitals (Arjmandi , 2022), offices (Yu , 2018; Bhat , 2022), halls (Shao and Li, 2020), general ventilated spaces (Zhang 2021), cars (Arpino , 2022; Mathai , 2022), buses (Duchaine , 2021; Ho and Binns, 2021; and Zhang , 2021), airliner cabins (Yang , 2018; Talaat , 2021), restaurants (Ho, 2021; Liu 2021; and Wu 2021), elevators (Dbouk and Drikakis, 2021), and parkings (Nazari , 2021).

To simulate aerosol dispersion, DNS, LES, the numerical solution of the Unsteady Reynolds-Averaged Navier–Stokes equations (URANS), and hybrid URANS-LES methods have been used. These three techniques can be understood as three different levels of modeling to simulate turbulent flows. While DNS solves directly for all spatial and temporal scales, LES models the effect of the turbulent scales smaller than the grid on the larger scales that are explicitly solved by the computational grid. URANS-based simulations are based on a higher level of turbulence modeling, since they numerically solve the Unsteady Reynolds Averaged equations with models for the closure terms (Pope, 2000). Although some DNS have been reported in the literature so far (Diwan , 2020; Chong , 2021; and Fabregat 2021a), generally, LES have been preferred for predictions of the short-term dispersion because of the relatively small computational costs it requires (Abkarian , 2020; Pendar and Páscoa, 2020; Calmet , 2021; Liu , 2021; and Wang 2021). On the other hand, URANS, and to a lesser extent LES coupled with hybrid methods, has been used to simulate the long-term dispersion, which usually involves larger computational domains and longer simulations times compared to the short-term dispersion.

In this paper, we focus our attention on the effect of the turbulence modeling and of the numerical methods used in the simulations for the prediction of the short-term dispersion of an isolated violent expiratory event. The warm and relatively humid flow ejection into a typically colder and drier ambient usually lasts less than 1 s (Gupta , 2009; Tang , 2013) and has a maximum Reynolds number of order 10^{4} (Bourouiba, 2021). During the flow injection, a new turbulent jet is formed, with a penetration length from the source that scales with time as t^{1/2} (Chaudhuri , 2020). When the flow injection ceases, the jet evolves into a turbulent buoyant puff, with a time scaling of t^{1/4} as the turbulence progressively decays (Scorer, 1997). This transient turbulence level is challenging for simulations because many turbulence models for URANS or some subgrid-scale models for LES exhibit a decrease in their performance when dealing with decaying unsteady turbulent flows with buoyancy effects. This is especially true when the flow progressively relaminarizes during the puff stage. Aerosol particles are expelled during the flow injection, and most simulations found in the literature track the particle dispersion following individual evaporating or non-evaporating solid spherical particles by means of Lagrangian methods assuming low particle volume fraction ( $\varphi \u223c 10 \u2212 5$ Duguid, 1946) and, consequently, neglecting the particle collisions and adopting the one-way coupling between the phases. Coalescence and breakup phenomena are also neglected in almost all the available studies. The rationale behind this choice is that none of these mechanisms is well understood yet for the specific problem of pathogen-laden droplets transported in an ambient air flow (Zhou and Zou, 2021). This implies uncertainty in the modeling and (expected) high variability in the simulated results.

The larger inertial particles expelled during the exhalation follow quasi-ballistic trajectories that are almost independent of the flow, whereas small and less inertial particles tend to remain afloat within the wake of the decaying new jet and the frontal puff, even when they are larger than the Kolmogorov scale (Pourfattah , 2021). The Lagrangian tracking of these small particles is also challenging when using turbulence or subgrid-scale models because their trajectories are more sensitive to the surrounding flow conditions and to the local velocity. In these cases, random-walk models, tuned according to the local level of turbulence, may be used to account for the effect of the modeled fluctuations on the particle trajectory (Mofakham and Ahmadi, 2020).

The usual grid distribution strategy adopted in the simulations that are reported in the literature concentrates the mesh elements toward the flow inlet and near the growing shear layers of the new jet flow. This strategy produces skewed mesh elements that can affect the convergence of the numerical solution and may reduce its accuracy. Also, the use of a relatively coarse mesh away from the source can influence the dynamics of the frontal puff and the capability of LES to capture the decaying turbulent flow due to excessive filtering. These complex flow features of the spatially and temporally evolving flow associated with a violent expiratory event make the construction of the grid also challenging since one needs to correctly capture the time evolutions of the turbulent intensities in different regions of the flow. The application of Adaptive Mesh Refinement (AMR) techniques can have advantages to capture with enough spatial resolution, flow structures that appear at different times and at different locations along the computational domain. At the same time, however, we remark that AMR may affect the comparison between different simulations, which is what we aim to do in this paper. For example, AMR can make it more difficult to directly identify the limitations of a turbulence model because AMR dynamically adjusts the mesh resolution in different parts of the domain based on the local solution features, thus masking some of the deficiencies of a specific turbulence model. The turbulence model may appear to perform better in some regions of the domain simply because the mesh has been refined in those areas and not because the model is more accurate.

The relevance that the prediction of the risk of transmission of infectious diseases via pathogen-laden aerosols has, together with the complexity of the simulations of the flow and particle dispersion generated in a violent expiratory event, motivated the organization of this international CFD challenge, which was launched in October 2021. The aim of the challenge is to assess the accuracy with which the different modeling approaches can reproduce the dynamics of a prototypical violent expiratory event. In this paper, we compare and analyze the different results submitted by the participants, who used different CFD codes and different turbulence modelization (LES, URANS, and Hybrid), and we summarize the outcomes. In Sec. II, we thus present the framework of the Challenge, and we describe the main characteristics of the codes, meshes, and numerical methods used by the different participants. The results are presented, compared, and discussed in Sec. III, and the final conclusions of the study and recommendations are outlined in Sec. IV.

## II. CHALLENGE FRAMEWORK AND PROBLEM DEFINITION

The objective of the challenge, organized by the members of Universitat Rovira i Virgili (Spain) and University of Udine (Italy), is to evaluate the ability of different CFD codes and different turbulence models to reproduce the short-term flow and the particle short-range dispersion that characterizes a prototypical violent expiratory event, which has a duration of 0.4 s (typical of a mild cough expelled in an indoor environment) and is dominated by the inertia of the air injection, such that the dispersion process is assumed to be essentially independent of the background air currents produced by a forced or natural ventilation system. Specifically, the Challenge aims at addressing two tasks: (1) assess the ability of each combination of numerical method and turbulence modeling used for the simulation to reproduce predictions of DNS of the unsteady jet flow of an idealized exhalation event, characterized by the rapid, but limited in time injection of warm air into an initially quiescent ambient (stage I of the Challenge); (2) estimate the impact of the particle size and evaporation on the dispersion of the aerosol cloud, generated in the exhalation event, and compare the results with the predictions obtained combining DNS of the flow with a one-way coupling Lagrangian tracking scheme (stage II of the Challenge). Each team was free to carry out just one or both tasks. However, since the transmission of airborne diseases depends on the dispersion of pathogen-laden particles, the teams were strongly encouraged to address both.

To prevent bias among participants, the challenge was planned as a blind test. The participants had no information about choices made by the other teams before the deadline for the data submission and the file exchange between the participants and the organizers was made through a private Google Drive folder for each team. The DNS data published in Fabregat (2021a; 2021b), which were already available prior to the challenge, were suggested as a reference. The data requested had to be sent in ASCII or binary vtk files to facilitate the post-processing and visualization steps, which were performed using the open-source multi-platform ParaView software (Paraview, 2022).

The physical model, the computational domain, and the boundary conditions proposed for the challenge were the same as those used by Fabregat (2021a; 2021b). The interested reader is referred to these two publications for information about the specific set of physical parameters selected for the simulation. Figure 1 shows the cylindrical computational domain and the frame of reference. The dimensions of the cylindrical computational domain were $H$ = 1.60 and $D$ = 1.00 m, and the flow inlet consisted of a cylindrical pipe (mimicking the oral cavity) with axial length $ H p$ = 0.04 m and diameter $d$ = 0.02 m. The Cartesian coordinate system uses ( $x$, $y$, $z$) as the spanwise, vertical, and streamwise directions with the origin located on the symmetry axis and placed at the downstream end of the inlet pipe. Gravity acts in the negative y direction (see Fig. 1).

To facilitate the transition to turbulence, a Gaussian bump of height $ H d$ = 0.001 m and width $\sigma $ = 0.002 m centered at $z$ = −0.01 m was placed inside the inlet pipe (see inset in Fig. 1). The exhalation event was modeled as an unsteady air injection at 34 °C into an initially quiescent environment at 15 °C. The air injection velocity ( $w$) was assumed to be uniform at the inlet circular section and followed a saw tooth profile equal to zero at $t$ = 0, then linearly increasing up to a maximum value of $w$ = 4.8 m/s at $t$ = 0.15 s, and finally linearly decreasing back to $w$ = 0 at $t$ = 0.4 s. The physical properties of the air were assumed to be constant except for the linear density variation with temperature, which only was taken into account in the vertical momentum equation, according to the Boussinesq approximation. The values for viscosity, density, thermal conductivity, heat capacity, and the thermal expansion coefficient were, respectively, $\mu $ = 1.95 × 10^{−5 }Pa s, $ \rho 0$ = 1.22 kg/m^{3}, $k$ = 0.0277 W/m K, $Cp$ = 1010 J/kg K, and $\beta $ = 3.36 × 10^{−3} K^{−1}. The constant pressure or outflow boundary conditions were prescribed at the far field exit ( $z=H$) and at the lateral cylindrical surface $ ( D / 2 ) 2= x 2+ y 2$ for 0 < $\u2009z$ < $H$. At the surface $z$ = 0, ( $D$/2)^{2 }> $x$^{2 }+ $\u2009y$^{2} > ( $d$/2)^{2} and at the cylindrical inlet pipe wall, non-slip boundary conditions had to be used. Uniform Dirichlet boundary conditions were imposed at the pipe inlet ( $z$ = − $ H d$) according to the time evolution of the air injection velocity described above. This yields $u$ = 0, $v$ = 0, and $w$ = $w$( $t$). All the boundary surfaces were assumed to be adiabatic except for the circular inlet where a constant uniform Dirichlet condition for the temperature, namely, $T$ = 34 °C, had to be prescribed. Each team was free to select the boundary values for the relevant turbulence quantities at the circular inlet. The physical time covered by the simulations was also prescribed and set equal to $t$ = 1.60 s.

For stage II, we considered seven different spherical particle diameters: 4, 8, 16, 32, 64, 128, and 256 *μ*m. Particles had to be released continuously at the inlet with the same velocity of the surrounding fluid over the entire duration of the exhalation ( $t$ ≤ 0.4 s). Each team was free to decide how to model the evaporation of the aqueous fraction of the particles. If included, the minimum allowable diameter due to evaporation had to be one-third of the initial value, as in Fabregat (2021b). This corresponds to evaporating droplets laden with a 3% volume fraction of nonvolatile matter (Wang , 2021). Exhaled and ambient air relative humidities were taken as 85% and 65%, respectively. The physical properties of the particles were assumed to be those of water at ambient temperature ( $ \rho p$ = 1000 kg/m^{3}, $ k p$= 0.606 W/m K, $ C p p$ = 4180 J/kg K), whereas a value $\Delta H v$ = 2.257 × 10^{6} J/kg was selected for the enthalpy of evaporation.

The deliverables for stage I were vtk files, which are readable with ParaView (2022), with 2D instantaneous scalar fields on the symmetry plane $x$ = 0 at three different times: at $t$ = 0.25 s, during the exhalation, at $t$ = 0.40 s, corresponding to the end of the flow injection and at $t$ = 1.50 s, about 1 s after the exhalation. Each file contained the temperature (°C), the axial (horizontal) velocity component (m/s), and the vertical velocity component (m/s).

An international call for participation in the Challenge was released in October 2021, and initially 16 teams from around the world showed interest and asked for information about the instructions to contribute. Six months later at the deadline for data submission (May 1, 2022), seven teams submitted the requested results. The results obtained by each team were presented and discussed during an online workshop held on June 27, 2022. Table I shows, in alphabetical order, the list of teams that submitted a complete set of data for stage I, and their country of origin.

Barcelona Supercomputing Center | Spain |

Westmont College-Federal University of Uberlândia | USA-Brazil |

Otto von Guericke Universität Magdeburg | Germany |

Peter the Great St. Petersburg Polytechnic University | Russia |

RMIT University-The University of Sydney | Australia |

Universitat Politècnica de Catalunya | Spain |

University of Maribor-University of Erlangen Nuremberg | Slovenia-Germany |

Barcelona Supercomputing Center | Spain |

Westmont College-Federal University of Uberlândia | USA-Brazil |

Otto von Guericke Universität Magdeburg | Germany |

Peter the Great St. Petersburg Polytechnic University | Russia |

RMIT University-The University of Sydney | Australia |

Universitat Politècnica de Catalunya | Spain |

University of Maribor-University of Erlangen Nuremberg | Slovenia-Germany |

In total, 11 different simulations were submitted for stage I. Table II summarizes the list of teams, in random order, and the details of their simulations. Each simulation is identified by an alphanumeric code, reported in the last column of Table II. This code classifies each team with a randomly assigned letter (A–G) followed by an integer representing the simulations performed by that team (the integer is increased in the case of multiple simulations by a single team) and a further letter: U for a URANS-based simulation, L, for LES or H, for hybrid URANS-LES method. Finally, a Roman numeral is used to indicate the stage the simulation refers to (I or II).

Team . | Code . | Turbulence . | Spatial discretization . | Temporal discretization . | Grid . | Boundary conditions . | Code . | |
---|---|---|---|---|---|---|---|---|

A | Ansys Fluent v19.3 (Ansys Inc., 2022) | k-ε | Finite volume, second order | Implicit, second order, Δt = 10^{–5} s | 0.6 × 10^{6} hexahedral | Turbulence intensity 5% Turbulent viscosity ratio 10 | A1-U-I | |

URANS | k-ε RNG | A2-U-I | ||||||

k-ω SST | A3-U-I | |||||||

LES | WMLES S-Ω (Pr_{t }= 0.85) | 35 × 10^{6} polyhedral | Vortex method | A4-L-I | ||||

B | Alya-HPC mechanics (Vázquez , 2016) | LES of low-Mach equations (Vreman, 2004) (Pr_{t} = 0.7) | Finite element, second order | Explicit. Δt = 4.4 × 10^{–6} s | 46 × 10^{6} tetrahedral/prism | Laminar inlet | B1-L-I | |

C | OpenFOAM v8 and v9 (OpenFoam, 2022) | URANS: k- ω SST | Finite volume, second order | Adaptative 3.5 × 10^{–4} ≤ Δt ≤2 10^{–3} s | 6.2 × 10^{6} hexahedral (792 × 243 × 36)a | k = 0.0864 m^{2}/s^{2} ω = 387 s^{−1} | C1-U-I | |

D | Ansys Fluent v2021R2 (Ansys Inc., 2022) | Hybrid: SBES, SST-LES with WALE | Finite volume, second order | Bounded second order implicit. | 20 × 10^{6} poly-hexcore | Turbulence intensity 5% | D1-H-I | |

Δt = 10^{–5} s for t < 0.8 s | ||||||||

Turbulent viscosity ratio 10 | ||||||||

10^{–5} ≤ Δt ≤ 10^{–4} s for t > 0.8 s | ||||||||

E | UNSCYFL3D (Fontes , 2019) | URANS two-layer k-ε | Finite volume, second order | Implicit, second order Δt = 5 × 10^{–4} s | 0.68 × 10^{6} hexahedral (120 × 110 × 52)a | k = 10^{–8} m^{2}/s^{2} ε = 10^{–8} m^{2}/s^{3} | E1-U-I | |

F | STAR-CCM+ 2020.1 and 2021.3 (Siemens, 2022) | URANS: SST k-ω model | Finite volume, second order | Second order implicit, Δt = 10^{–3} s | 1.6 × 10^{6} hexahedral (366 × 120 × 56)a | Turbulence intensity 1% Turbulent viscosity ratio 10 | F1-U-I | |

LES: WALE model (Pr_{t }= 0.9) | PISO unsteady. Δt = 10^{–5} s | 14.2 × 10^{6} hexahedral (602 × 274 × 136)a | Laminar inlet | F2-L-I | ||||

G | OpenFOAM v7 (OpenFoam, 2022) | Hybrid: k-ω-SST DES | Finite volume, second order | Implicit, 6 × 10^{–6} ≤ Δt ≤ 6 × 10^{–4} s | 30 × 10^{6} hexahedra, polyhedra | Variable turbulence intensity and k as a function of the inlet Re number | G1-H-I |

Team . | Code . | Turbulence . | Spatial discretization . | Temporal discretization . | Grid . | Boundary conditions . | Code . | |
---|---|---|---|---|---|---|---|---|

A | Ansys Fluent v19.3 (Ansys Inc., 2022) | k-ε | Finite volume, second order | Implicit, second order, Δt = 10^{–5} s | 0.6 × 10^{6} hexahedral | Turbulence intensity 5% Turbulent viscosity ratio 10 | A1-U-I | |

URANS | k-ε RNG | A2-U-I | ||||||

k-ω SST | A3-U-I | |||||||

LES | WMLES S-Ω (Pr_{t }= 0.85) | 35 × 10^{6} polyhedral | Vortex method | A4-L-I | ||||

B | Alya-HPC mechanics (Vázquez , 2016) | LES of low-Mach equations (Vreman, 2004) (Pr_{t} = 0.7) | Finite element, second order | Explicit. Δt = 4.4 × 10^{–6} s | 46 × 10^{6} tetrahedral/prism | Laminar inlet | B1-L-I | |

C | OpenFOAM v8 and v9 (OpenFoam, 2022) | URANS: k- ω SST | Finite volume, second order | Adaptative 3.5 × 10^{–4} ≤ Δt ≤2 10^{–3} s | 6.2 × 10^{6} hexahedral (792 × 243 × 36)a | k = 0.0864 m^{2}/s^{2} ω = 387 s^{−1} | C1-U-I | |

D | Ansys Fluent v2021R2 (Ansys Inc., 2022) | Hybrid: SBES, SST-LES with WALE | Finite volume, second order | Bounded second order implicit. | 20 × 10^{6} poly-hexcore | Turbulence intensity 5% | D1-H-I | |

Δt = 10^{–5} s for t < 0.8 s | ||||||||

Turbulent viscosity ratio 10 | ||||||||

10^{–5} ≤ Δt ≤ 10^{–4} s for t > 0.8 s | ||||||||

E | UNSCYFL3D (Fontes , 2019) | URANS two-layer k-ε | Finite volume, second order | Implicit, second order Δt = 5 × 10^{–4} s | 0.68 × 10^{6} hexahedral (120 × 110 × 52)a | k = 10^{–8} m^{2}/s^{2} ε = 10^{–8} m^{2}/s^{3} | E1-U-I | |

F | STAR-CCM+ 2020.1 and 2021.3 (Siemens, 2022) | URANS: SST k-ω model | Finite volume, second order | Second order implicit, Δt = 10^{–3} s | 1.6 × 10^{6} hexahedral (366 × 120 × 56)a | Turbulence intensity 1% Turbulent viscosity ratio 10 | F1-U-I | |

LES: WALE model (Pr_{t }= 0.9) | PISO unsteady. Δt = 10^{–5} s | 14.2 × 10^{6} hexahedral (602 × 274 × 136)a | Laminar inlet | F2-L-I | ||||

G | OpenFOAM v7 (OpenFoam, 2022) | Hybrid: k-ω-SST DES | Finite volume, second order | Implicit, 6 × 10^{–6} ≤ Δt ≤ 6 × 10^{–4} s | 30 × 10^{6} hexahedra, polyhedra | Variable turbulence intensity and k as a function of the inlet Re number | G1-H-I |

^{a}

Number of elements: axial × radial × angular.

Table II shows that three teams performed the simulations with commercial codes: two teams (A and D) with Ansys Fluent and one team (F) with STAR-CCM+. The open-source OpenFOAM solver was used by two teams (C and G), whereas the two remaining teams (B and E) used their own in-house codes. Note that team B considered the low-Mach version of the transport equations (Le Quéré , 2005) without the buoyancy term in the vertical momentum equations. All the simulations were carried out with second-order finite volume codes, except for those of team B, which used a code based on a second-order finite element technique for spatial discretization. Of the 11 simulations, 6 were based on the numerical solution of the URANS equations, 3 used the LES technique, and 2 used a URANS-LES hybrid methods. It should be noted that for the physical model considered in the Challenge, the role of the URANS in these hybrid methods is limited to the region near the non-slip wall of the cylindrical inlet with the Gaussian bump (see Fig. 1). Most of the teams performed the time-marching procedure with implicit schemes and the time-steps used for the URANS-based simulations ranged approximately from 10^{−3} s to 10^{−5} s, while the LES and hybrid URANS-LES simulations were carried out with time steps in the range 10^{−5}–10^{−6} s. Participants were asked to select the mesh resolution according to a test of grid independence. The resulting grids were preferably constituted by hexahedra and polyhedral volumes. URANS-based simulations were carried out with a number of grid elements ranging from 0.6 × 10^{6} to 6.2 × 10^{6}, while the number of elements for LES and hybrid methods ranged from 14.2 × 10^{6} to 46 × 10^{6} elements.

Five teams submitted the requested data for stage II and the details of the Lagrangian tracking methods used by each team are indicated in Table III, together with the identifying code assigned to each simulation. For stage II, the code's last integer indicates if evaporation is (1) or is not (0) accounted for. As can be seen from Table III, evaporation was considered in 3 simulations out of a total of 12 simulations of the particle dispersion submitted under different conditions. In seven of these simulations, the flow was computed by solving the URANS equations while there are three simulations in which LES was used and only one was run using a hybrid method. All the teams presented results for the seven selected particle diameters. Particle tracking was performed considering a force balance that included the drag force, the gravity and buoyancy forces, and a turbulent dispersion force or a subgrid-scale random walk for URANS or LES, respectively.

Team . | Turbulence . | Evaporation . | Forces . | Number of particles . | Numerical method . | Code . |
---|---|---|---|---|---|---|

A | URANS | No | Drag, gravity, discrete random walk model | 13 800 for each diameter | Lagrangian tracking Runge–Kutta method | A1-U-II-0 |

A2-U-II-0 | ||||||

A3-U-II-0 | ||||||

LES | Diameters: 4, 8, 16, 32, 64, 128, and 256 μm | A4-L-II-0 | ||||

C | URANS | No | Drag, gravity, turbulent dispersion (stochasticDispersionRAS) | 13 800 for each diameter | Lagrangian tracking. Euler implicit. | C1-U-II-0 |

Diameters: 4, 8, 16, 32, 64, 128, and 256 μm released in half of the inlet cross section | Interpolation method: Linear weighted interpolation using cell values. | |||||

E | URANS | No | Drag, gravity, turbulence dispersion force | 5600 for each diameter | Lagrangian tracking scheme: Analytical scheme, interpolation method for fluid velocity: Second order | E1-U-II-0 |

Diameters: 4, 8, 16, 32, 64, 128, and 256 μm | ||||||

F | URANS | No | Drag, gravity, turbulence dispersion force | 13 800 for each diameter | Lagrangian tracking. | F1-U-II-0 |

Yes | F1-U-II-1 | |||||

Diameters: 4, 8, 16, 32, 64, 128, and 256 μm released at x = 0, y = 0, and z = 0. | Second order tracking integration and interpolation. | |||||

LES | No | F2-L-II-0 | ||||

Yes | F2-L-II-1 | |||||

G | Hybrid URANS-LES | No | Drag, gravity, turbulent dispersion (stochasticDispersionRAS) | 27 738 for each diameter | Lagrangian tracking. Euler implicit. | G1-H-II-0 |

Yes | Diameters: 4, 8, 16, 32, 64, 128, and 256 μm | Interpolation method: Linear weighted interpolation using cell values. | G1-H-II-1 |

Team . | Turbulence . | Evaporation . | Forces . | Number of particles . | Numerical method . | Code . |
---|---|---|---|---|---|---|

A | URANS | No | Drag, gravity, discrete random walk model | 13 800 for each diameter | Lagrangian tracking Runge–Kutta method | A1-U-II-0 |

A2-U-II-0 | ||||||

A3-U-II-0 | ||||||

LES | Diameters: 4, 8, 16, 32, 64, 128, and 256 μm | A4-L-II-0 | ||||

C | URANS | No | Drag, gravity, turbulent dispersion (stochasticDispersionRAS) | 13 800 for each diameter | Lagrangian tracking. Euler implicit. | C1-U-II-0 |

Diameters: 4, 8, 16, 32, 64, 128, and 256 μm released in half of the inlet cross section | Interpolation method: Linear weighted interpolation using cell values. | |||||

E | URANS | No | Drag, gravity, turbulence dispersion force | 5600 for each diameter | Lagrangian tracking scheme: Analytical scheme, interpolation method for fluid velocity: Second order | E1-U-II-0 |

Diameters: 4, 8, 16, 32, 64, 128, and 256 μm | ||||||

F | URANS | No | Drag, gravity, turbulence dispersion force | 13 800 for each diameter | Lagrangian tracking. | F1-U-II-0 |

Yes | F1-U-II-1 | |||||

Diameters: 4, 8, 16, 32, 64, 128, and 256 μm released at x = 0, y = 0, and z = 0. | Second order tracking integration and interpolation. | |||||

LES | No | F2-L-II-0 | ||||

Yes | F2-L-II-1 | |||||

G | Hybrid URANS-LES | No | Drag, gravity, turbulent dispersion (stochasticDispersionRAS) | 27 738 for each diameter | Lagrangian tracking. Euler implicit. | G1-H-II-0 |

Yes | Diameters: 4, 8, 16, 32, 64, 128, and 256 μm | Interpolation method: Linear weighted interpolation using cell values. | G1-H-II-1 |

## III. RESULTS AND DISCUSSION

### A. Instantaneous thermal and velocity fields (stage I)

The teams submitted the instantaneous thermal and velocity fields in the symmetry plane $x$ = 0 at three different times ( $t$ = 0.25, $t$ = 0.4, and $t$ = 1.5 s). We recall that the prototypical exhalation prescribed as inlet flow reaches the maximum velocity at $t$ = 0.15 and ends at 0.4 s (see Sec. II). We grouped the results provided by the teams using URANS in Figs. 2 and 3, and the results provided by either LES or hybrid methods in Figs. 4 and 5. Figures 2 and 4 show the non-dimensional temperature contours, whereas Figs. 3 and 5 show the contours of the axial (z) component of the velocity vector in physical units. The non-dimensional temperature is defined as $\theta =(T\u2212 T o)/( T i\u2212 T o)$, where $ T i$ is the temperature of the flow at injection ( $ T i$ = 34 °C) and $ T o$ is the background temperature ( $ T o$ = 15 °C). This yields $\theta =0$ at the far-field and $\theta =1$ at the flow inlet. The code used to identify each simulation, indicated in the last column of Table II, is shown to the left of the panels of Figs. 2–4.

A DNS or LES of a given turbulent flow represents a single realization of the flow, and the variability between the instantaneous fields of different realizations can be determined by performing different simulations with slightly different boundary conditions (Trivedi , 2021). In the panels in the second row of Figs. 2–5, we included, for comparison purposes, the average temperature and axial velocity contours of 13 independent LES of the flow carried out by the organizers with different inlet velocity boundary conditions. In these 13 LES, only the flow dispersion was simulated. The details of these simulations are described in the Appendix. The average temperature and axial velocity contours can be compared directly with the corresponding fields predicted by the URANS simulations, shown in Figs. 2 and 3, because the URANS predictions can be considered as equivalent to the ensemble averaged fields of many different flow realizations. This allows to estimate the variability of the different realizations of the same turbulent flow. In all the figures, the instantaneous contours of the DNS reported by Fabregat (2021a) are included on the top row for comparison purposes.

The theoretical shapes of the envelope of the thermal cloud, estimated using the models of Scorer (1997) and Richards (1968), and applied to violent expiratory events in Pallares and Fabregat (2022), have been superimposed to the temperature contours in Figs. 2 and 4. The model developed by Scorer (1997) to obtain the time-evolution of the front of the new jet predicts an envelope with a quasi-conical shape with an elliptical cross section. At the plane $x$ = 0, the contour of the envelope corresponds to the white triangular shape shown in Figs. 2 and 4. The position and radius of the frontal spherical thermal puff (Richards, 1968) is indicated in these figures with a white circle. The theoretical triangular and circular envelopes frame well the instantaneous thermal puff of the DNS shown in the top three panels of Fig. 2, especially for $t$ = 0.4 and $t$ = 1.5 s. Note that the flow injection, that lasts up to 0.4 s, is still active at $t$ = 0.25 s, and the frontal quasi-spherical puff, completely formed only when the flow injection ceases, is not fully developed.

The six predictions of the thermal cloud by the URANS-based simulations at $t$ = 0.25, 0.5, and 1.5 s are shown in Fig. 2. Note that the vertical dispersion of the jet obtained with these simulations is smaller than the theoretical one, indicated by the triangular shapes, and the one predicted by the DNS. The positions of the frontal quasi-spherical thermal puff at $t$ = 0.4 s and $t$ = 1.5 s are very well captured by the simulations C1-U-I and F1-U-I, while the other simulations (A1-U-I, A2-U-I, A3-U-I, and E1-U-I) underpredict to some extent the range and the vertical displacement of the frontal thermal puff. The reason for these improved predictions of the thermal puff dynamics has to do mainly with the grid resolution: The C1-U-I and F1-U-I simulations were run on meshes of 6.2 × 10^{6} and 1.6 × 10^{6} elements, respectively, whereas the A1-U-I, A2-U-I, A3-U-I, and E1-U-I simulations used coarser meshes of about 0.6 × 10^{6} elements. Comparing the predictions of simulations A3-U-I and C1-U-I, both carried out with the $k$-ω-SST model against the predictions of simulations A1-U-I and A2-U-I, carried out with $k$-ε models (see Table II), it can be seen that the specific turbulence model has a very limited impact on the prediction of the position and extension of the frontal thermal puff. Similar conclusions can be obtained by comparing the predictions of the axial velocity component plotted in Fig. 3. In this case, simulations performed with finer meshes (C1-U-I and F1-U-I) better reproduce the range of DNS instantaneous velocity contours, which are characterized by higher vertical dispersion when compared to URANS simulations. It is also remarkable that simulation F1-U-I shows a frontal region of high velocity that is detached from the thermal cloud and is qualitatively similar to that found in the DNS (see Fig. 3).

We also remark here that the instantaneous axial velocity was not reported by team B (B1-L-I); therefore, Fig. 5 shows the modulus of the instantaneous velocity for this team. At $t$ = 0.25 and $t$ = 0.4 s, all the LES show very similar temperature and velocity contours. As for the URANS-based simulations (cases A1-U-I to F1-U-I in Figs. 2 and 3), it is evident that the vertical dispersion of the jet is slightly underestimated with respect to the DNS, but both the range and the extension of the frontal thermal puff are well reproduced in all the cases. At $t$ = 1.5 s, simulations F2-L-I and G1-H-I exhibit intense velocity activity in the frontal thermal puff, as observed in DNS too (see Fig. 5). This can be also appreciated in Fig. 4, where the range of the thermal cloud of these two simulations appears larger than the LES range. The comparison of the characteristics of these two simulations with those of the remaining simulations shown in Table II does not reveal a clear reason for this difference in terms of mesh resolution, inlet velocity conditions, or SGS model. It can be probably attributed to the stochasticity of the different realizations of the same turbulent flow. In fact, Trivedi (2021) reported a relatively large variability in the range of the thermal cloud in different realizations (see, for example, Fig. 2 of Trivedi , 2021). Also, as it will be shown in Subsection III B, the variability found in the horizontal range of the thermal cloud at $t$ = 1.5 s among the different LES shown in Fig. 4 is compatible with the values of the horizontal size of the thermal cloud of the 13 independent LES.

### B. Position and size of the thermal cloud (stage I)

As shown in Figs. 6(a) ( $t$ = 0.25) and 6(b) ( $t$ = 0.4 s), the center location and the size of the thermal cloud are mostly well predicted by the LES and URANS-based simulations. At larger times [Fig. 6(c)], the uncertainty of the predictions increases. For the three times considered, the vertical location of the centroid resulting from the URANS-based simulations, plotted with red symbols at the bottom of the panels, agrees well with those of the DNS and the 13 LES, but the widths of the thermal cloud are underpredicted especially at $t$ = 1.50 s. At this time, the variability of the vertical sizes of the thermal cloud, defined as $2 \sigma T y$, is relatively small among the URANS predictions, which ranges between 0.07 (A1-U-I) and 0.08 m (C1-U-I). These values are about 50% smaller than those corresponding to the DNS (0.15) and the 13 LES (0.12 m), indicating a reduced turbulent vertical mixing due to the use of turbulence models in the URANS simulations. In these cases, the effect of buoyancy in the transport equations of the turbulent kinetic energy and of the dissipation is to enhance/suppress turbulence in unstable/stable stratified regions of the flow. However, it is known that these models show some deficiencies in the predictions of the spread rates of vertical plumes (Kuma and Dewan, 2014) and buoyant horizontal jets (Alfaifi , 2019).

The position of the centroid of the thermal cloud predicted by LES, plotted at the top of the panels of Fig. 6 with blue symbols, falls within the variability of the 13 LES for $t$ = 0.25 s [Fig. 6(a)] and $t$ = 0.4 s [Fig. 6(b)]. At $t$ = 1.5 s, there are two simulations (B1-L-I and D1-H-I) that predict vertical positions of the centroid below $y$ = 0. The inspection of the instantaneous temperature distribution provided by the simulation B1-L-I at $t$ = 1.5 s (Fig. 4) shows that the hot ascending plume near the coordinates origin, formed by the hot fluid remaining in the inlet pipe after the flow injection has ceased, is not reproduced in this case. Also, it can be seen that, at $t$ = 1.5 s, the temperature distribution is in fact shifted toward $y$ < 0, because, in this simulation, the buoyancy term was not considered in the vertical momentum equation. Although not as pronounced, the instantaneous distribution of the simulation D1-H-I also shows this shift. In this case, the position of the centroid and the vertical size of the thermal cloud are within or very close to the variability of the 13 LES. Figure 6(c) shows that the other LES (A4-L-I, F2-L-I and G1-H-I) predict vertical sizes of the thermal cloud that fall within the variability of the 13 LES.

The axial positions and the horizontal sizes of the thermal cloud are plotted in Fig. 7. In general, the URANS-based simulations, shown in the bottom portion of each panel, predict smaller axial ranges of the thermal cloud (i.e., larger horizontal positions of the centroid) when performed with finer meshes (C1-U-I and F1-U-I, see Table II), as compared to coarser simulations (A1-U-I, A2-U-I, A3-U-I, and E1-U-I) at the three times considered. At $t$ = 1.5 s [Fig. 7(c)], the variability of the horizontal position of the centroid for the URANS-based simulations lays between 0.24 (A2-U-I) and 0.42 m (C1-U-I) with the predictions of the DNS, yielding $ c T z$ = 0.36 m and the 13 LES, yielding $ c T z$ = 0.32 m. The horizontal extension, defined as $2 \sigma T z$, ranges between 0.5 (F1-U-I) and 0.38 m (A2-U-I) with that of the DNS being 0.45 m and that of the 13 LES being 0.39 m. This indicates that, in general, at large times, the horizontal size of the thermal cloud is well captured by the URANS-based simulations, while a large scatter is observed in the prediction of the horizontal position of the centroid. This scatter depends on the combination of the grid resolution and the turbulence model used. The results suggest that, regardless of the turbulence model used (see Table II), finer meshes allow the thermal cloud to extend over each larger horizontal distances.

The estimations of the position and sizes of the thermal cloud produced by LES are plotted using blue symbols at the top of the panels of Fig. 7. At $t$ = 0.25 [Fig. 7(a)] and $t$ = 0.4 s [Fig. 7(b)], the predictions lay close or within the variability of the 13 LES, which is indicated by the gray areas. At larger times [ $t$ = 1.5 s, Fig. 7(c)], the different LES exhibit larger variability in terms of the centroid position and horizontal sizes of the thermal cloud. The horizontal extensions range between 0.39 (A4-L-I) and 0.51 m (G1-H-I). These values compare well with DNS, which yields $2 \sigma T z=$ 0.45 m and the 13 LES, which yield $2 \sigma T z=$ 0.39 m. The centroid position exhibits a significant larger variability, but its predictions are anyway close to those of DNS and the 13 LES, especially for simulations D1-H-I and B1-L-I, which lay within the variability of the 13 LES, indicated with gray areas in Fig. 7(c). According to the information shown in Table II, these differences in the axial position of the centroid of the thermal cloud at $t$ = 1.5 s can be attributed to the inherent variability of different realizations of the same turbulent flow combined with the use of different grid resolutions (note that for the LES, these range between meshes of 14.2 × 10^{6}–46 × 10^{6} elements). Overall, the results show that the variability of the different predictions, especially those relative to the horizontal position and to the extension of the thermal cloud, increases with time and simulations at larger times ( $t$ > 1.5 s) can show significant differences, up to some tens of centimeters.

### C. Particle clouds (stage II)

In this subsection, we compare and discuss the results concerning the centroid trajectories and the particle cloud size for seven different diameters (4, 8, 16, 32, 64, 128, and 256 *μ*m). Table II summarizes the different cases for stage II, which include nine simulations that consider non-evaporative particles (6 URANS and 3 LES/Hybrid) and three simulations with evaporative particles (1 URANS and 2 LES/Hybrid).

Figures 8 and 9 show, respectively, the trajectories of the centroids of the non-evaporative particle clouds for the URANS-based simulations and for the LES, together with the predictions of the DNS. The corresponding trajectories for evaporative particles are plotted in Figs. 10 and 11.

Figure 8 (URANS) shows that cloud centroid trajectories for the small particles (4, 8, and 16 *μ*m, top-row panels the figure) mostly remain afloat within the thermal puff and progressively increase their axial range, departing from the DNS result as their size increases. This is also observed for the evaporative particles, as shown in the top panel of Fig. 10 (URANS). A similar trend is also depicted in Figs. 9 and 11 corresponding to the predictions of LES.

This increased axial range for the evaporative and non-evaporative particles with diameters 8 and 16 *μ*m, with respect to the DNS, can be attributed to the reduced vertical span of the thermal puff predicted by the URANS-based and LES simulations in comparison with the DNS (see, for example, the vertical sizes of the thermal puff at $t$ = 0.4 s and $t$ = 1.5 s in Fig. 6). The reduced vertical size of the thermal puff is associated with less vertical mixing and higher axial fluid velocities that transport particles of increasing inertia (4, 8, and 16 *μ*m) to progressively larger axial positions as they remain afloat within the thermal puff.

The comparison of the top panels of Figs. 8 (URANS) and 9 (LES) shows that the axial spread of the particle centroid predicted by LES is reduced with respect to the DNS. This agrees with the enhanced vertical mixing within the thermal cloud that LES is able to reproduce in comparison with the URANS-based simulations, thus resulting in a wider thermal cloud, as shown in Fig. 6.

For the smaller particle sizes (4, 8, and 16 *μ*m), the URANS predictions (Fig. 8) do not exhibit a clear trend based on the specific turbulent model. For example, the results for cases A1-U-II-0, A2-U-II-0, and E1-U-II-0, computed with the $k$-ε models, notably differ from those reported for cases A3-U-II-0, C1-U-II-0, and F1-U-II-0, which used the $k$-ω models. Specifically, both groups of simulations predicted very similar cloud trajectories with small variability, potentially explained by differences in the mesh resolutions. Compare, for example, the DNS with cases A3-U-II-0 and C1-U-II-0 performed with the $k$-ω models but with 0.6 × 10^{6} and 6.2 × 10^{6} elements.

The trajectory of the cloud centroid for the particles with size 32 *μ*m strongly depends on the evaporation of the particles. As indicated by the DNS trajectories shown in Figs. 8–11, the 32 *μ*m particles remain afloat at $t$= 1.5 s if evaporation is considered but tend to fall if evaporation is not considered. This trend is generally well captured by URANS and LES, except for the cases F2-U-II-0 (see Fig. 9) and F1-U-II-1 (see Fig. 10).

According to the DNS, the larger particles with diameters 64, 128, and 256 *μ*m tend to fall with essentially the same trajectory irrespectively of the evaporation. For some cases (for example, A1-U-II-0, A2-U-II-0, or C1-U-II-0 in Fig. 8), the trajectories of the cloud centroid for the 256 *μ*m particles show backward trajectories (i.e., toward the plane $z$ = 0 plane) at large times. This is probably caused by the fact that, at large times, particles that have reached the bottom of the computational domain are excluded from the calculation of the centroid, which takes into account only the remaining particles injected at the end of the exhalation (which are freely falling relatively close to the inlet pipe). This effect can explain the departure from the DNS of the trajectories of the clouds for particles with the largest diameters (128 and 256 *μ*m). However, for these cases, the trajectories of the centroid at small times (namely, close to the flow injection) compare well with those provided by DNS.

Contrary to what is observed for smaller particles, the URANS-based simulations tend to underestimate the axial range of the particle cloud centroids for particles with quasi-ballistic trajectories (e.g., those with diameter equal to 64 *μ*m or larger). This may be due to the reduced vertical size of the thermal cloud predicted by the URANS-based simulations, which would explain why these relatively massive particles tend to leave the thermal cloud earlier than in DNS.

The vertical and horizontal positions of the cloud centroid as well as its size for the case of non-evaporative and evaporative particles with diameters 4 and 32 *μ*m at $t$ = 0.4 and $t$ = 1.5 s are shown in Figs. 12 and 13. For reference, the DNS predictions of the particle clouds (green lines) and the DNS predictions of the centroid and sizes of the thermal cloud (black lines) are also included. We selected these two specific particle diameters, out of the seven values considered, to illustrate the effect of the different numerical approaches on the predictions when particles remain afloat within the thermal puff for large times. In Figs. 12 and 13, the data corresponding to the non-evaporative particles are plotted on panels (a-0) to (d-0), while the data of the evaporative particles are shown in panels (a-1) to (d-1).

Figure 12 shows that the DNS predictions of the vertical position of the centroids of the particles, which tend to fall due to of gravity, are below the centroids of the thermal cloud, which tend to rise due to buoyancy. In general, the predictions of the vertical position of the centroid of the non-evaporative particle clouds (panels a-0, b-0, c-0, and d-0) are very similar among the different cases, including URANS-based simulations and LES, and in agreement with DNS, especially at short times (i.e., at $t$ = 0.4 s). However, the predictions of the vertical sizes of the particle clouds are systematically underpredicted in comparison with the DNS. A similar trend has been found and discussed in Subsection III B when comparing the vertical sizes of the thermal cloud shown in Fig. 6. The effect of evaporation on the vertical position of the centroid of the particle cloud appears to be very limited, as it can be seen by comparing the cases A1-U-II-0 with the case A1-U-II-1 corresponding to URANS-based simulations or the case F2-L-II-0 with the case F2-L-II-1 for LES. As predicted by DNS, the vertical size of the particle cloud is, in general, increased if evaporation is considered, especially for 4 *μ*m particles at $t$ = 1.5 s.

The horizontal extensions of the non-evaporative and evaporative clouds of particles with diameters 4 and 32 *μ*m are plotted in Fig. 13. As expected, the DNS results indicate that the position of the centroid is generally located at a smaller axial position compared to the centroid of the thermal cloud. At $t$ = 0.4 s, the URANS and LES predictions for the non-evaporative particles are close to the DNS results with a maximum difference, of about 0.05 m only. At larger times, the predictions exhibit a larger scatter and the maximum differences increase up to about 0.15 m. Both URANS-based simulations and LES perform similarly when compared with DNS. The horizontal size of the particle clouds is, in general, closer to the DNS results for the LES. URANS-based simulations obtained with coarser meshes (A1-U-II-0, A2-U-II-0, and A3-U-II-0) tend to overpredict the size of the cloud in comparison with simulations performed using finer meshes (C1-U-II-0, E1-U-II-0, and F1-U-II-0).

## IV. CONCLUSIONS

In this paper, we presented and discussed the results of a collaborative study, carried out by seven research teams within the framework of an international Computational Fluid Dynamics (CFD) challenge. Each participating team was asked to numerically simulate the short-term short-range flow and resulting particle dispersion generated by a prototypical violent expiratory event consisting in a short injection of hot fluid into a colder, initially quiescent, ambient. The results reported by each team have been compared against two reference cases, including a single realization provided by a Direct Numerical Simulation (DNS) and an averaged ensemble of 13 independent Large-Eddy Simulations (LES). The information obtained from these predictions is relevant for two reasons. First, to establish measures oriented to the minimization of risk infection of diseases transmitted by pathogen-laden aerosols expelled during expiratory events. Second, to generate realistic initial conditions for simulations of the long-term long-range flow and particle dispersion in indoor and outdoor scenarios. The seven research teams have performed a total of 11 simulations of the same flow configuration, but with different meshes and turbulence modeling approaches. Six simulations are based on the numerical solution of the Unsteady Reynolds Averaged Navier–Stokes (URANS) equations and five are based on the LES approach. The dispersion of both evaporative and non-evaporative particles with diameters ranging from 4 to 256 *μ*m has also been considered resulting in 12 simulations in total (7 based on the solution of the URANS equations and 5 based on the LES approach).

Overall, the URANS-based and the LES simulations correctly predict the vertical and horizontal mixing of the new jet formed during the flow injection (0 $\u2264t\u2264$ 0.4 s). At later times, of the order of 1 second after the flow injection has ceased ( $t$ ≈ 1.5 s), the vertical mixing of the thermal puff is, in general, underpredicted by both approaches with vertical sizes of the thermal cloud about 50% smaller than in DNS for URANS and 30% for LES. A larger variability in the predictions of the horizontal sizes of the thermal cloud is observed, but in general, both URANS and LES compare similarly well with the reference data. For the flow conditions considered, the specific turbulence model used for the URANS simulations has a very limited impact on the prediction of the position and extension of the frontal thermal puff. However, the use of finer meshes for the URANS-based simulations allows the thermal cloud to extend over larger horizontal distances, resulting in a better agreement with the DNS.

The trajectory and size of the particle clouds predicted by the DNS are, in general, well reproduced both by URANS-based and LES simulations. For the flow conditions simulated here, the cloud centroid exhibits wider axial ranges and lower vertical displacements than the DNS in the case of particles that remain afloat (4, 8, and 16 *μ*m). These differences are more evident for the larger (16 *μ*m) particles and for the URANS-based simulations. Their occurrence can be also connected to the under-predicted vertical mixing of the thermal puff discussed above. The 32 *μ*m particles tend to remain afloat if evaporation is considered but tend to fall otherwise. For the specific flow, thermal and humidity conditions considered in this study, this trend is generally well reproduced by the different approaches. For the larger particles (64, 128, and 256 *μ*m), which tend to follow quasi-ballistic trajectories and are less correlated with the local field compared to the smaller particles, the URANS and LES predictions are closer to the DNS.

## ACKNOWLEDGMENTS

This study was funded by the Spanish Ministerio de Ciencia, Innovación y Universidades through the Grant Nos. PID2020-113303GB-C21 and RTI2018-100907-A-I00 (MCIU/AEI/FEDER) and by the Generalitat de Catalunya through the Grant No. 2017-SGR-1234. M.Z., V.R., and N.I. acknowledge the Super-Computer Center (SCC) «Polytechnic» for providing computational resources. J.W., M.Š., and J.R. would like to thank the valuable insight given by professors Paul Steinmann and Matjaž Hriberšek and the financial support of the Deutsche Forschungsgemeinschaft, Germany under Project No. STE 544/58-2 and the Slovenian Research Agency under Project No. P2-0196.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Jordi Pallares:** Conceptualization (equal); Funding acquisition (equal); Resources (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Nikolay Ivanov:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Robert Castilla:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Pedro Javier Gamez-Montero:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Gustavo Raush:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Hadrien Calmet:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Daniel Mira:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Jana Wedel:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Mitja Štrakl:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Jure Ravnik:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Douglas Fontes:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Alexandre Fabregat:** Conceptualization (equal); Data curation (equal); Investigation (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Francisco Jose de Souza:** Data curation (equal); Investigation (equal); Visualization (equal); Writing – review & editing (equal). **Cristian Marchioli:** Conceptualization (lead); Methodology (equal); Resources (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Salvatore Cito:** Conceptualization (lead); Methodology (equal); Resources (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Akim Lavrinenko:** Conceptualization (equal); Data curation (equal); Investigation (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Hadifathul Akmal bin Norshamsudin:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Gábor Janiga:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **David F. Fletcher:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Kiao Inthavong:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Marina Zasimova:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). **Vladimir Ris:** Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: DETAILS ON LES ENSEMBLE

In this Appendix, we describe the details of the LES performed to estimate the flow variability between instantaneous fields of different realizations. We performed 13 LES of the flow dispersion, without considering the particles, with the same boundary conditions as those indicated in Sec. II but adding different small perturbations to the time evolution of the inlet velocity. Figure 14 shows the time evolutions of the inlet velocity of the 13 simulations and the corresponding intensities of the fluctuations. It can be seen that the averaged time evolution of these realizations follows well the linear increase and decrease in the inlet velocity used for the simulations of the Challenge. The fluctuations have been modulated to have similar intensities as the turbulent flow at the outlet of the mouth according to the simulations of the flow in the upper respiratory tract during a violent expiratory event (Pallares , 2022). These 13 LES were carried out up to t = 2 s (Δt = 10^{−3} s) in a reduced cylindrical computational domain (H = 50d and D = 40d, see Fig. 1) with a mesh resolution of 2.4 × 10^{6} elements using a second-order finite volume code and the WALE SGS model, as in Pallares (2022).

Figure 15 shows the time evolutions of the vertical [Fig. 15(a)] and axial [Fig. 15(b)] positions of the temperature centroid of these 13 simulations. The coordinates of the centroid of the thermal cloud ( $ c T y,\u2009 c T z$) are computed using Eq. (3). The curves corresponding to the position of the centroid plus/minus the vertical ( $ \sigma T y$) and horizontal ( $ \sigma T z$) sizes defined in Eq. (4) have been also included to indicate the evolution of the dimensions of the thermal cloud [i.e., $ c T y+ \sigma T y$ and $ c T y\u2212 \sigma T y$ in Fig. 15(a) and $ c T z+ \sigma T z$ and $ c T z\u2212 \sigma T z$ in Fig. 15(b)]. The shaded regions in Fig. 15 mark the variability of the 13 simulations within two times the standard deviation of the position of the centroid of the temperature distribution at each time. The variability of the dimensions of the cloud is also indicated with shaded regions. Consequently, if the different instantaneous positions and sizes of the temperature cloud follow a normal distribution, the shaded areas mark the region where 95% of the realizations would be. To determine the normality of the vertical and axial positions of the centroid and the sizes of the thermal cloud, we performed the Anderson–Darling normality test (Anderson and Darling, 1954) with a 1% of significance for the 13 realizations at each time. This test indicates that the different realizations satisfy the normality test for t > 0.2 s. The data at t < 0.2 fail the test probably because of the random nature of the white-noise perturbations introduced in the time evolution of the inlet velocity (see Fig. 14).

The position of the centroid of the thermal cloud and its vertical [Fig. 15(a)] and horizontal [Fig. 15(b)] sizes obtained with the DNS reported in Fabregat (2021a) are plotted in Fig. 14 with black lines. It can be seen that the vertical position and the vertical extension of the temperature cloud of the DNS is well within the variability of the 13 LES. For t < 1.2 s, the axial position and axial extension of the temperature cloud predicted by the DNS also fall within the variability of the 13 LES. At larger times (t > 1.2 s) and for the axial quantities [Fig. 15(b)], the DNS predicts a larger range, of about 5 cm, and extension for the thermal cloud. Note that for t > 1.2 s, the averaged axial position of the thermal cloud obtained with the DNS still increases with time, while the axial position for the 13 LES is essentially constant probably because of the excessive SGS viscosity at this late stage of the flow injection with a progressively reduced turbulence activity at these axial positions with a relatively coarse mesh.

## REFERENCES

*Dynamics of Meteorology and Climate*