DNS–URANS comparison study on cough-induced flow and particles in a large-scale circulation

$(⋅)¯$

= time-averaged value of a variable;

A

= area of the cough jet (m²);

C

= vapour concentration (–);

C_D

= drag coefficient (–);

C_μ1, C_μ2

= model constants (–);

Co

= courant number (–);

d_p

= particle diameter (m);

D

= mass diffusivity coefficient (m² s⁻¹);

e_dir

= random unit vector (–);

F_D

= aerodynamic drag force (N);

F_G

= gravitational force (N);

g

= gravitational acceleration vector (m s²);

h

= grid width (m);

I

= turbulence intensity (%);

IQR

= interquartile range (m);

k

= turbulent kinetic energy (m² s⁻²), wavenumber (m⁻¹);

m

= mass (kg);

n

= variable from a standard normal distribution N (0, 1) (–);

p

= pressure (Pa);

Pr

= prandtl number (–);

r

= average radius of the cough jet/puff (m);

r_pipe

= pipe radius (m);

Re

= reynolds number (–);

S

= maximum horizontal extend of the cough jet (m);

Sc

= schmidt number (–);

St

= stokes number (–);

t

= time (s);

t_e

= eddy lifetime (s);

t_flow

= flow-through time (s);

t_int

= particle-eddy interaction time (s);

t_r

= transit time of the particle across the eddy (s);

T

= temperature (K);

u

= velocity vector (m s⁻¹);

u

= bulk expansion velocity of the cough (m s⁻¹);

u_sa

= spatially-averaged horizontal velocity (m s^–1);

U_lid

= velocity magnitude of the lid (m s⁻¹); and

x, y, z

= Cartesian coordinates (m).

α

= entrainment coefficient (–);

β

= thermal expansion coefficient (K⁻¹);

γ

= vapour molar fraction expansion coefficient (–);

κ

= thermal diffusivity coefficient (m² s⁻¹);

λ

= wavelength (m);

λ_e

= eddy length scale (m);

ν

= kinematic viscosity (m² s⁻¹);

ρ

= density (kg m⁻³);

ε

= turbulent dissipation rate (m² s⁻³); and

η

= kolmogorov length scale (m).

0

= end of the jet phase;

a

= reference (ambient);

c

= cough;

cr

= carrier;

f

= fluid (air);

h

= horizontal;

p

= particle;

t

= turbulent; and

th

= threshold.

*

= dimensionless.

In the wake of the COVID-19 pandemic, society and the scientific community have been eager to better understand the mechanics of aerosol and droplet transmission. Violent expiratory events, such as coughing and sneezing, are critical pathways for the spread of infectious diseases, and it is crucial to include these events in simulations and models to assess the risk of airborne disease transmission (Wei and Li, 2015). A seminal study by Bourouiba et al. (2014) investigated the fluid dynamics and droplet transport of violent expiration events, demonstrating that aerosol suspension in expiration clouds prolongs the spread of respiratory pathogens. This risk is particularly high indoors, where poor ventilation and aerosol accumulation enhance the spread of pathogens (Peng et al., 2022).

In this context, computational fluid dynamics (CFD) has emerged as an essential tool to unravel the complex interactions between expiratory flows, background flows and exhaled particles. Previous numerical studies have used various techniques, such as direct numerical simulation (DNS), steady and unsteady Reynolds-averaged Navier–Stokes simulation (RANS and URANS, respectively) and large eddy simulation (LES). DNS and LES are rarely performed due to enormous computational resources required. Instead, RANS and URANS are commonly used owing to their lower computational cost. However, the accuracy of these methods in capturing airborne disease transmission remains a subject of debate. There is limited consensus on how these methods capture particle transport in ventilated environments, making robust infection risk analysis difficult.

The problem of expiratory particle transport is typically governed by the Navier–Stokes and Lagrangian particle transport equations. To solve these equations directly, DNS is the most accurate method, as it fully resolves all relevant turbulent scales in space and time. Numerous DNS studies have analysed expiratory droplet transport in quiescent environments. These studies have shown that the turbulence of the cough puff and droplet evaporation can increase the residence time of particles and significantly influence their trajectories (Chong et al., 2021; Ng et al., 2021; Rosti et al., 2021), challenging the classical Wells model (Wells, 1934). Furthermore, the DNS study by Rosti (2020) showed that the initial size of an exhaled droplet can dramatically affect its trajectory. These findings improve our understanding of short-term and short-range particle transport in quiescent environments.

Recent CFD studies have focused on the effect of expiratory flow on particle transport, with the aim of evaluating time-averaged and turbulence model-based CFD approaches in comparison with DNS in quiescent environments. Lavrinenko et al. (2022) conducted a comparative study between DNS and URANS to predict cough-induced flow and particle cloud trajectories and found that URANS struggled to accurately predict particle transport due to the unsteady nature of the cough flow, particularly for small particles. Pallares et al. (2023) performed a comprehensive investigation of the flow dynamics and particle transport of a violent expiratory event using various modelling techniques, including URANS, LES and hybrid variants (URANS–LES), and compared their results with the DNS performed by Fabregat et al. (2021a, 2021b). Irrespective of the specific turbulence model used, URANS generally underpredicted vertical turbulent mixing. In addition, both URANS and LES overpredicted the horizontal spread of small aerosol particles compared to DNS. Bi et al. (2022) compared cough-induced particle transport using URANS and LES to further validate the numerical results with experimental data. They found that while URANS and LES could reasonably predict the evolution of a human cough jet, LES showed better agreement with experiments.

The aforementioned studies were conducted over short time periods and within small domains due to computational resource limitations. In addition, these studies have been carried out in quiescent environments, neglecting ventilation effects that significantly affect long-term and long-range aerosol transport. The most common approaches used to study aerosol transport in large ventilated environments over long periods of time are RANS and URANS. Aljabair et al. (2023) studied particle transport in a ventilated room using URANS under different air conditioning settings and cough velocity profiles and found that particle transport is highly sensitive to both air conditioning and cough flow intensities. Sanada and Mat (2023) used URANS to investigate the effect of cross-ventilation on airborne particles emitted by a cough and concluded that increasing cross-ventilation intensity also increases the aerosol transfer rate due to enhanced forced convection of the air.

Wang et al. (2022) and Konstantinov et al. (2023) both used URANS to analyse aerosol transport in a train compartment, providing insight into the critical role of circulation patterns in influencing aerosol transmission. In addition, Talaat et al. (2021), Yang et al. (2018) and Yan et al. (2020) conducted URANS to investigate aerosol dynamics in aircraft cabins, with all studies agreeing that cough jets cause strong perturbations in the ambient flow, thereby increasing the movement and mixing of exhaled particles in the cabin. Schmeling et al. (2023) and Zhang et al. (2009) conducted both experimental and numerical studies of aerosol transmission in an aircraft cabin. For the numerical aspect of their work, they conducted RANS to calculate the steady-state velocity field in the cabin with transient discrete particles. Both studies highlight the challenges of matching experimental data with numerical simulations and suggest that RANS may have limitations in capturing all the details of a cough and background flows in an aircraft cabin. Building on these findings, Tamaddon Jahromi et al. (2022) provided an alternative approach by demonstrating the potential of combining the RANS and transient discrete particle approach with machine learning-based predictions to accurately predict the temporal evolution of particle dispersion in various indoor ventilation scenarios.

Although these simulations are useful for resolving flow fields and particle dynamics over longer periods and in larger spaces, their accuracy in predicting background flows, expiratory flows and their interactions with expiratory particles remains uncertain. DNS, although more accurate, is largely limited to short-term and short-range studies. This gap in evaluating the accuracy of turbulence and particle dispersion models in time-averaged simulations to predict expiratory particle transmission in ventilated public spaces or transport systems motivated us to perform a fully resolved Navier–Stokes simulation (DNS) in a large environment, typical of a ventilated domain, over an extended period of time. We then compared the results with those obtained using the more widely used URANS, under identical conditions and domain size, but with different spatial and temporal resolution and approach to turbulence.

The present study compares the results of cough-induced particles interacting with a background flow in a relatively large domain of 12.5 m³, obtained using DNS and URANS over a 30 s period. The background flow, driven by a lid, represents a large-scale circulation (LSC), typical for aircraft cabins and train compartments. The initial velocity fields used to initiate the DNS and the URANS of the cough were generated in a well-resolved Navier–Stokes simulation and a URANS of the LSC, respectively, without the presence of a cough. The time-dependent cough flow and the expiratory particle size distribution were derived from the experimental studies of Gupta et al. (2009) and Duguid (1946), respectively. Wei and Li (2017) described coughing as a two-stage process with an initial jet phase, characterised by a concentrated and directed cough flow, followed by a puff phase that occurs after the source supply ceases. Similarly, in our study the predicted flow fields, particle displacements and dispersions were analysed separately in these two phases. For statistical evaluation, the particles were additionally categorised into six clouds of different sizes.

After the introduction, the paper continues with an overview of the governing equations and the discretisation schemes in Chapter 2. Chapter 3 discusses the simulation parameters and the grid details. Chapter 4 presents the simulation details and results of the LSC simulations, together with the initial flow fields prescribed in the cough simulations. Finally, Chapter 5 analyses and compares the results of the DNS and the URANS of the cough. The paper concludes with a discussion in Chapter 6.

Although the flow is assumed to be incompressible, the differences in temperature (T) and vapour concentration (C) between the ambient air and the exhaled cough result in density variations and a non-negligible buoyancy effect. To account for these effects on the velocity field, the temperature and vapour concentration fields were coupled to the velocity field using the Boussinesq approximation (Gray and Giorgini, 1976). As the governing equations solved in the DNS and the URANS are different, they are presented separately in Sections 2.1.1 and 2.1.2.

2.1.1 Direct numerical simulation (DNS). The flow is governed by the incompressible Navier–Stokes equations and the transport equations for T and C:

Equations (1)–(3) were discretised using second-order accurate central differences to predict the flow variables at the cell centres of a Cartesian grid of hexagonal cells. The discretised equations were solved using a finite volume method provided by OpenFOAM®. To perform a DNS, we implemented Chorin’s (1968) projection method for velocity-pressure coupling and an explicit second-order accurate Euler-leapfrog time-stepping scheme. It was shown in Kath and Wagner (2016) that for the here used Cartesian grids this combination provides the accuracy needed for a DNS. The time-step for the DNS was set according to the von Neumann stability criterion derived for second-order accurate central differences by Shishkina and Wagner (2004), with an explicit time-step of Δt = 10⁻⁵ s (or Δt* = t/t₀ = 1.66 × 10⁻⁵, if made dimensionless with the jet phase duration t₀ = 0.6 s).

2.1.2 Unsteady Reynolds averaged Navier–Stokes simulation (URANS). The URANS solves the incompressible unsteady Reynolds-averaged Navier–Stokes equations, together with the corresponding transport equations for the time-averaged temperature (⁠$T¯$⁠) and vapour concentration (⁠$C¯$⁠) coupled using the Boussinesq approximation:

The turbulent viscosity (ν_t) is required to model the Reynolds stress tensor following the turbulent viscosity hypothesis (Schmitt, 2007). A two-equation renormalisation group k-ε turbulence model was used to calculate ν_t based on the solutions of additional transport equations for the turbulent kinetic energy (k) and the turbulence dissipation rate (ε) (Yakhot et al., 1992). In this model, the turbulent viscosity is expressed as $νt=Cμ1k2ε$⁠, where C_μ1 is a model constant with a value of 0.0845. The turbulent thermal diffusivity (κ_t) and mass diffusivity (D_t) are defined as κ_t = ν_t/Pr_t and D_t = ν_t/Sc_t, respectively, where the turbulent Prandtl number (Pr_t) is 0.85 and the turbulent Schmidt number (Sc_t) is 0.82.

Equations (4)–(6) were discretised using second-order accurate central differences, while the transport equations for k and ε were discretised using a first-order accurate upwind discretisation scheme. The discretised equations were solved using the OpenFOAM^® library. Unlike for the DNS, the pressure-implicit with splitting of operators algorithm was applied to obtain the pressure field. The time advancement of equations (4)–(6) was carried out using an implicit second-order accurate backward time scheme and an adaptive time-step in the range of 10⁻³ s to 10⁻² s (or 1.66 × 10⁻³ to 1.66 × 10⁻²) if made dimensionless with t₀), ensuring that the Courant number (Co.) is well below 1.

To study the particle dynamics, the expiratory particles are assumed to be discrete spherical points of constant mass. Newton’s second law of motion is applied and solved simultaneously with the flow field at each time-step, considering both the aerodynamic drag forces (F_D) and the gravitational forces (F_G), as follows:

As described in Amsden et al. (1989), the drag coefficient (C_D) is a function of the particle Reynolds number (Rep) and is given by:

where:

In the DNS, the velocity vector of the air surrounding the particle u_cr= u, i.e. the effective air velocity is equal to the instantaneous velocity at the location of the particle. In the URANS, particles are convected with the mean velocity field. To account for the effect of turbulent fluctuations on particle motion, the stochastic dispersion model is used, which adds the stochastic turbulent dispersion velocity vector:

to the mean field u_cr= ū + u_t. In equation (10), |n| is obtained to satisfy the standard normal distribution N (0, 1), and e_dir is a random unit vector. The model is active only if the particle–eddy interaction time (t_int) is greater than the time-step size (Δt), where t_int = min (t_e, t_r). Here, t_e = k/ε is the eddy lifetime, and t_r = λ_e/ε is the transit time of the particle across the eddy, with $λe=Cµ2kε$ represents the eddy length scale, where C_µ2 = 0.16432. During the interaction between the eddy and the particle, u_t remains unchanged. Afterwards, u_t is recalculated for the next eddy. A more detailed explanation of the model can be found in Gosman and Ioannides (1981).

The computational domain, shown in Figure 1, consists of a 4 m × 2.5 m × 1.25 m room and a circular pipe connected to the wall at a height of 1.6 m, representing the mouth opening. Several tests have confirmed that the domain is sufficiently large to capture the cough and the dispersed particles within the considered time period. The pipe representing the mouth opening is 6 cm long with a radius of r_pipe = 2.4 cm.

Considering that any introduced energy is dissipated at the tail of the energy spectra, the maximum wavenumber is k_max = 1/η, where η is the Kolmogorov scale. On the other hand, under a discrete grid with spacing h, the smallest resolvable wavelength is λ_min = 2 h, following the Nyquist condition. Therefore, the largest wavenumber that can be captured with this grid in a DNS is:

Thus, to resolve all relevant scales in a DNS, the condition k_DNS > k_max or:

must be satisfied, where ϵ is the unknown turbulent dissipation rate. To estimate ϵ, we follow Mazzino and Rosti (2021) and determine the scales at the end of the jet phase (t₀ = 0.6 s) to obtain $ϵ=u03/r0$⁠, where r₀ is the average radius of the cough jet at t₀, and =u₀ = r₀/t₀ represents the bulk expansion velocity of the cough (see Section 5.1 for the calculation of r₀). Thus, according to equation (12), the mean width of the grid must satisfy the criterion h < 1.8 mm (or h* = h/l₀ < 0.02). Finally, considering that the maximum Pr and Sc are both less than 1, it follows that the grid is sufficient to resolve the smallest scales in the temperature and concentration fields.

In principle, the maximum width of the grid cells used in the DNS is 3.5 mm (or h^* ≈ 0.04) in accordance to the DNS study by Rosti et al. (2021). However, the cells are linearly refined in different subdomains as shown in Figure 2. Within the pipe subdomain, the grid is isotropic with a width of 1.2 mm which corresponds to h^* ≈ 0.0138 in dimensionless terms. The grid width in the y-direction (h_y) remains constant for all y < |r_pipe|. For y > |r_pipe|, h_y increases linearly up to the value 3.5 mm at y = 0.04 m. The grid width in the x-direction (h_x) also increases linearly from 1.2 mm at x = 0.06 m up to the value 3.5 mm at x = 0.26 m. In addition, the grid is refined from the maximum grid width of 3.5 mm normal to any wall to 0.7 mm (or h^* ≈ 0.008) such that y⁺ < 1 defined by the lid-driven circulation. Since the smallest scales that can be resolved with the grid (h/π = 0.38 mm) are smaller than the Kolmogorov scale at the end of the jet phase (η₀ = 0.64 mm), we conclude that the resolution criteria required for a DNS is fulfilled. In total, this leads to a DNS grid consisting of 372 million hexagonal cells. For comparison the URANS grid consists of 18 million cells and was generated according to previous URANS studies of cough jets in ventilated environments (Cravero and Marsano, 2022; Oh et al., 2022; Sanada and Mat, 2023). Figure 2 shows the mesh cross-sections for the DNS and the URANS near the pipe.

Figure 3 summarises the boundary conditions and illustrates a snapshot of the coughing jet carrying discrete particles and the LSC. At the inlet, a Dirichlet boundary condition was applied to the velocity, where the velocity varies according to a time-dependent cough velocity profile measured in an experiment by Gupta et al. (2009). The exhalation process starts at t = 0 s and ends at t = t₀ = 0.6 s, with a maximum jet Re of 8,500. The cough at the inlet was defined with a constant temperature of 30°C and a vapour concentration of 0.0418 (≈ 100% relative humidity [RH]).

To generate such a turbulent cough in the DNS, 16 eddy-like perturbations were added to the laminar velocity profile at the inlet of the pipe by applying the Divergence Free Synthetic Eddy Method by Poletto et al. (2013), corresponding to a turbulence intensity of 10%. In the URANS, the turbulence kinetic energy k = 1.5(I|ū|)² was prescribed at the inlet, also with a turbulence intensity (I) of 10%. Figure 4 (left) shows the prescribed time-dependent bulk cough jet velocity at the inlet.

To simulate an LSC, a lid is moved at the boundary downstream of the jet, causing the flow to move downwards near the mouth opening, as illustrated in Figure 3. Wang et al. (2008) reported that air velocities in the passenger breathing zone of an aircraft cabin range from 0 to 0.2 m/s. To achieve similar air velocities near the pipe region, a lid velocity of 0.5 m/s was prescribed. To maintain mass continuity during the exhalation process, a square outlet with an area of 0.16 m² was placed in the centre of the lower boundary with a Neumann boundary condition for velocity and a Dirichlet boundary condition for total pressure. All other boundaries are treated as non-slip walls, set to 25°C and a vapour concentration of 0.0187 (≈ 60% RH) to maintain ambient conditions. In addition, the URANS approximates the near-wall flow with a logarithmic wall function.

Simulations were initiated with resting 4,000 particles randomly positioned on the circular particle injection surface with a diameter of 2 cm (see Figure 3). Given the considerable disagreement in the literature regarding the initial particle size distribution of a cough event (Morawska, 2006), we chose an initial particle size distribution based on the experimental study of Duguid (1946), with a probability function shown in Figure 4 (right).

Particles were classified into six different size categories or “clouds” as follows: 1 µm ≤ d_p ≤ 8 µm, 8 µm < d_p ≤ 16 µm, 16 µm < d_p ≤ 32 µm, 32 µm < d_p ≤ 64 µm, 64 µm < d_p ≤ 128 µm and 128 µm < d_p ≤ 256 µm. These are referred to throughout the paper as 1–8 µm, 8–16 µm, 16–32 µm, 32–64 µm, 64–128 µm and 128–256 µm, respectively, or in combinations (i.e. 1–16 µm includes the clouds of 1–8 µm and 8–16 µm). Each cloud contains a minimum of 400 particles to ensure statistically robust results.

During the cough jet, particles were injected with a constant count rate between t = 0.015 s and t = 0.335 s. The Stokes number (St) of the particles during the cough varies between St ≈ 1 × 10⁻⁴ and St ≈ 100, depending on their diameter and cough velocity. Table 1 lists the cough, air, particle and lid properties used in our simulations.

Prior to the DNS and the URANS, the initial flow fields of a fully developed LSC were generated in a well-resolved Navier–Stokes simulation and a URANS, respectively. Both the well-resolved Navier–Stokes simulation and the URANS of the LSC solve equations (1) and (4), respectively, neglecting the buoyancy term and other transport equations for T and C, using the mesh described in Section 3.1 and discretised as discussed in Section 2.1.1. However, to generate the initial fields, the PIMPLE algorithm provided by the OpenFOAM^® library was used instead of an explicit time-marching method with a time-step size of 10⁻³ s and 10⁻² s, resulting in a maximum Co. below 4 and 1, respectively.

The flow-through time (t_flow) is calculated from the lid length (L = 2.5 m) and the lid velocity (U_lid = 0.5 m/s), giving t_flow = L/U_lid = 5 s. Figure 5 demonstrates for the well-resolved Navier–Stokes simulation of the LSC that the velocity magnitude at the three points along the pipe centreline, i.e. 0.1 m, 1 m and 3 m from the pipe inlet, fluctuates around a stable mean value after $tflow* ≈ 160$⁠.

The well-resolved Navier–Stokes simulation and the URANS of the LSC were both run up to $tflow*=360$⁠, and the resulting flow fields were used as initial conditions for the DNS and the URANS of the cough, as shown in Figure 6. Figure 7 further illustrates the time-averaged velocity magnitude distributions along the lines y = 0 m and x = 2 m in the mid-xy plane, averaged over $tflow*=160$ to $tflow*=360$⁠. The time-averaged velocity profiles predicted by the well-resolved Navier–Stokes simulation and the URANS are in close agreement. Any differences can be attributed to the different lateral and spatial resolutions of the simulations and turbulence treatments.

The jet phase refers to the first 0.6 s of the cough simulation during the cough event (see Figure 4). The horizontal (axial direction of the pipe) velocity of the cough is typically the most dominant velocity component and represents the primary direction of flow during the jet phase. Figure 8 shows the horizontal velocity fields u_x and ū_x of the cough jet at four different times along the mid-xy plane predicted in the DNS and the URANS, respectively.

Overall, the mixing layer of the jet in the URANS appears smoother compared to the fluctuating mixing layer of the DNS. At t = 0.1 s, just after the jet reaches its maximum bulk flow velocity at t ≈ 0.07 s, the horizontal extent and overall shape of the jet are similar in both simulations. At t = 0.2 s, this agreement persists, although the shape predicted by the URANS is more symmetrical than that of the DNS. By t = 0.35 s, the URANS jet does not extend as far as the DNS jet. Finally, at the end of the jet phase, the non-smooth nose of the DNS jet is not visible in the URANS jet, which is about 18% shorter and has a more symmetrical shape than the DNS jet on the left.

To practically interpret the differences in the horizontal velocity fields, we calculated the spatially averaged horizontal velocity:

in a subdomain defined by 0 m ≤ x ≤ 1.1 m, –0.2 m ≤ y ≤ 0.2 m and –0.2 m ≤ z ≤ 0.2 m, where the total volume V_t = 0.044 m³, and N is the total number of cells in the subdomains. In equation (13), u_h,i represents the horizontal velocity of the i-th cell, where u_h = u_x for the DNS and u_h = ū_x for the URANS.

The time series of u_sa during the jet phase were calculated for both the DNS and the URANS and are shown in Figure 9 (left). At the very beginning, i.e. before the start of the cough, both simulations predicted a negative u_sa, indicating that the background flow in both simulations is directed in the opposite direction to the jet flow.

The URANS overpredicts u_sa by 0.011 m/s at t = 0 s compared to the DNS. To properly compare u_sa of the jets during the jet phase, we shifted ū_x by adding 0.011 m/s to account for the initial discrepancy caused by the background flow on u_sa (assuming this initial error remains constant throughout the jet phase). In Figure 9 (left), the results of u_sa from the DNS (u_x) and the offset curve from the URANS (ū_x + 0.011) closely follow each other, indicating that the URANS accurately predicted u_sa.

To quantify the rate at which ambient air is entrained into the jet, we further analysed and compared the entrainment coefficient, α = s/r, during the jet phase, as shown in Figure 9 (right). Here, s and r represent the maximum horizontal extent and the average radius of the cough jet, respectively. The outer edge of the jet is defined based on the threshold vapour concentration $Cth=C0+Cc 1000$⁠, which is slightly larger than the ambient vapour concentration, about 0.1% of the maximum vapour concentration of the cough. The average radius r is calculated as r = (r_y + r_z)/2, where r_y = 0.5 A_y/s represents the average radius of the jet in the mid-xy plane and r_z = 0.5 A_z/s represents the average radius of the jet in the mid-xz plane. Here A_y and A_z denote the areas of the jet in the mid-xy and mid-xz planes, respectively, as defined by the outer edge of the jet.

Up to t = 0.1 s, α is similar for both simulations, showing a rapid decrease during the initial phase of the jet. From 0.1 to 0.6 s, the entrainment observed in the URANS exceeds that in the DNS, with time-averaged α of 0.13 and 0.11, respectively. This overall 18% overprediction of jet entrainment in the URANS corresponds closely to the 18% shorter horizontal extent of the jet also predicted by the URANS, as shown in Figure 8 at t = 0.6 s. The larger entrainment volume of ambient air results in a larger jet volume (i.e. greater spatial dispersion) and consequently a reduced horizontal extent of the jet.

Figure 10 shows the spatial distribution of the exhaled particle clouds of the DNS (top) and the URANS (middle) at the end of the jet phase from the lateral- and rear-views.

In the lateral-views of Figure 10, particles smaller than 64 µm are predominantly trapped within the jet in both simulations. In contrast, larger particles tend to separate from the cough jet, with significant differences between the predictions of the simulations. For example, 24% of particles in the 64–128 µm range have already left the cough jet in the DNS, whereas only 14% have done so in the URANS. For particles in the 128–256 µm range, the DNS predicts that 92% have already left the cough jet, compared to 82% in the URANS. It is worth noting that the particles in the range of 64–128 µm and 128–256 µm resolved in the URANS show a 10% greater entrapment than those in the DNS. This discrepancy can be attributed to the greater entrainment observed in the URANS, which results in a larger jet volume encompassing more particles.

Further analysis of the rear-views of the particle distributions in Figure 10 reveals an uneven lateral dispersion of particles within the cough jet in the DNS, whereas in the URANS the dispersion pattern is more uniform in the radial direction, resulting in a more circular shape.

The statistical comparison of the lateral-views in Figure 10 (bottom) provides information about the horizontal transport of the particles (along the x-axis), while the comparison of the rear-views provides information on their spatial transport (along the z-axis). In this analysis we used the median as an indicator of particle cloud displacement and the interquartile range (Q3 – Q1) as a measure of particle cloud dispersion. In addition, the horizontal displacement and the lateral dispersion of the particle clouds are summarised in Table 2.

The horizontal displacements of particles smaller than 32 µm resolved in the URANS are on average about 23–25% shorter than those in the DNS due to the shorter extent of the jet itself (with these particles being predominantly trapped within the jet). This difference decreases with increasing particle size and becomes almost negligible for particles larger than 128 µm.

The most significant overprediction of lateral dispersion is observed for the particle cloud with the smallest particle size (1–8 µm) by 47%, while the most accurate lateral dispersion predictions occur for particles in the 32–64 µm range by 8%. In contrast, the lateral dispersion of the particle cloud with the largest particle size is underpredicted by 27%.

The jet phase, which ends at t = 0.6 s, is followed by the puff phase which lasts until the end of the simulation (t = 30 s). To analyse and compare the flow fields, the horizontal (u_x and ū_x) and vertical (u_y and ū_y) velocity fields are shown in Figure 11 at four different times: t = 2 s, 5 s, 15 and 30 s, along the mid-xy plane. The contour line of the outer edge of the cough puff corresponds to the vapour concentration threshold C_th (see Section 5.1 for the detailed definition of C_th).

Examining the horizontal velocity component shown in Figure 11, it is clear that the cough puff and the LSC of the URANS are smoother than those of the DNS, which resolves all the details of the fluctuations. At t = 2 s, the URANS underpredicts the horizontal extent of the cough puff by 20%, similar to what was observed at the end of the jet phase (t = 0.6 s). This underprediction becomes more pronounced over time. By t = 30 s, the puff extends horizontally up to 3 m in the DNS, whereas it is 33% shorter in the URANS. One reason for this discrepancy is that the puff in the URANS almost stops propagating horizontally after 15 s, whereas in the DNS the frontal portion continues to propagate within the room. The DNS shows the puff splitting into many large and small segments, whereas in the URANS it remains smooth and connected, with a broader span. In addition, the upper part of the cough puff resolved in the URANS disturbs the upper boundary layer of the LSC more than in the DNS (clearest at t = 15 s and t = 30 s), indicating a more intense interaction with the strong backward flow of the LSC, which could also contribute to the weakened puff propagation.

The vertical velocity fields in Figure 11 reveal that the URANS smooths out the vertical fluctuations of both the puff and the background flow. Notably, the vertical extent of the puff predicted by the URANS is much larger than that predicted by the DNS. This could be due to a stronger interaction between the downward background flow near the mouth and the puff, causing the cough puff to almost reach the ground by the end of the simulation. This behaviour is not observed in the DNS, possibly because the faster propagation of the puff (due to smaller entrainment) allows it to escape the strong downward background flow near the mouth.

To evaluate the displacements of the exhaled particles, we analysed the time series of the horizontal and vertical displacements of the particle clouds, as shown in the top two graphs of Figure 12. Horizontal and vertical displacements are defined by the x and y distances between the median particle cloud position and the centre of the mouth inlet, respectively. The lower plots in Figure 12 show the differences in horizontal and vertical displacements resolved in the DNS and URANS for each particle cloud.

In the DNS, particles smaller than 32 µm show similar horizontal displacements up to 15 s, reaching their maximum horizontal displacement during the simulation. After this point, they begin to move slightly backward, with smaller particles showing greater backward displacement. This behaviour can be attributed to the backward background flow, which was particularly strong near the top wall. Notably, particles in the 32–64 µm range already diverge from the smaller particles after 2 s, resulting in shorter displacements until t = 22 s. However, these particles eventually surpass the smaller ones and reach the largest horizontal displacement of all particle clouds by the end of the simulation, with a final displacement of 1.55 m. Particles larger than 64 µm stop their horizontal propagation within 10 s, with final displacements ranging from 0.7 to 0.95 m.

In the URANS, particles smaller than 64 µm show similar horizontal displacements during the early puff phase (up to the 10-second mark). Subsequently, particles in the 16–64 µm range surpass the smaller ones, reaching the largest final horizontal displacement of 1.4–1.6 m by the end of the simulation, while the final displacement for the particle clouds with smaller particles (1–16 µm) is 1.15–1.2 m. It is worth noting that while the puff in the URANS may not propagate further after 15 s, the horizontal displacement of particles in the 16–64 µm range continues to increase. Particles larger than 64 µm stop propagating horizontally within 10 s, with final displacements ranging from 0.7 to 0.85 m for particles in the 64–256 µm range, with the shortest horizontal propagation observed for the largest particles (128–256 µm).

The lower left plot in Figure 12 shows that the URANS underpredicts the horizontal displacement of particles smaller than 64 µm during the early puff phase compared to the DNS, with this underprediction being more pronounced for smaller particles. This behaviour is consistent with the shorter puff propagation observed in the URANS during the early puff phase. After this phase, the magnitude of the underprediction decreases, especially for particles in the 16–64 µm range, which finish the simulation with a very small error of 0.05 m. We conclude that the horizontal displacements of the particles in the URANS are slower in the early puff phase but faster thereafter compared to those in the DNS. Although the frontal portion of the puff continues to propagate in the DNS after t = 15 s, the displacement of smaller particles in the URANS still exceeds that in the DNS overall. Particles in the 64–128 µm range are underpredicted by about 0.08 m during the early puff phase and maintain this difference throughout the simulation. Despite these discrepancies, the URANS accurately predicts the horizontal displacement of particles larger than 128 µm with almost zero error, where gravitational forces are dominant.

Examining the vertical displacements of the particle clouds, we observed that in both simulations, the larger the particle, the faster it loses height due to the increasing influence of gravitational forces. In the DNS, particles smaller than 16 µm showed a slight increase in height over the simulation time, ending up about 10 cm higher than the height of the mouth. Particles in the 16–32 µm range experienced a slight descent, losing about 20 cm in height by the end. Particles in the 32–64 µm range almost reached the ground by the end of the simulation, while those in the 64–128 µm and 128–256 µm ranges settled completely within 9.5 and 3 s, respectively.

In contrast, the URANS predicted no ascent above the mouth for any particle cloud throughout the simulation. Particles smaller than 16 µm descended by about 0.20 m, while those in the 16–32 µm range descended more, about twice as much. Particles in the 32–64 µm range almost settled by the end of the simulation, reaching a height of just 0.1 m above the ground. In addition, particles in the 64–128 µm and 128–256 µm ranges settled completely within 9.5 and 3 s, respectively.

The lower right plot in Figure 12 shows that the URANS predicts the vertical displacements of larger particles more accurately than smaller particles. In the URANS, particles smaller than 32 µm consistently overpredicted the height loss, ending up 0.2–0.25 m lower than those in the DNS, which are mostly transported within the cough puff. As mentioned above, the cough puff in the URANS almost reaches the ground at the end of the simulation, and the smaller particles mainly transported by this puff end up at lower heights.

In addition, the accuracy of the URANS in predicting the vertical displacement of the 32–64 µm particle cloud shows a distinct pattern throughout the simulation. Initially, it has the same accuracy as the smaller particles in the early puff phase, but then the vertical difference between the simulations decreases in the mid puff phase (10–20 s) and becomes negligible in the late puff phase (20–30 s). However, at the end of the simulation, these particles are slightly higher in the URANS compared to those in the DNS, with a difference of almost 0.1 m. It is worth noting that these particles settle almost completely in both simulations, suggesting that this difference is likely to disappear shortly after t = 30 s. During the early puff phase, particles in the 64–256 µm size range settle completely in both simulations, resulting in no vertical difference at all.

To evaluate the dispersion of the exhaled particles, we analysed the time series of the interquartile range (IQR = Q3−Q1) of the particle clouds, as shown in the top graph of Figure 13. The IQRs were calculated by averaging over all three dimensions, where IQR = (IQR_x + IQR_y + IQR_z)/3. The IQR serves as an important indicator of particle dispersion around the median position of the particle cloud. In addition, we calculated the dispersion differences for each particle cloud, shown in the bottom plot of Figure 13.

The dispersions of particles smaller than 32 µm resolved in the DNS follow logarithmic growth trends, reaching magnitudes of 1.1–1.4 m by the end of the simulation. Overall, we observed that smaller particles show larger dispersions. Particles in the 32–64 µm, 64–128 µm and 128–256 µm ranges show their highest dispersion after 2.5, 7 and 22 s with magnitudes of 1.05 m, 0.81 m and 0.75 m, respectively. These values indicate the start of particle deposition on the ground. Particles in the 64–128 µm and 128–256 µm ranges settle completely within 9.5 and 3 s, respectively, with final dispersions of 0.4 m and 0.5 m (where IQR_y = 0 m due to complete deposition on the ground). The dispersion trends of the particle clouds resolved in the URANS show similar characteristics to those in the DNS, but with different magnitudes and timing.

The lower plot in Figure 13 illustrates the differences in the particle cloud dispersions resolved in the DNS and URANS. It can be seen that in the early puff phase, the URANS overpredicts the dispersion of particles smaller than 32 µm by 0.05–0.1 m compared to the DNS. This overprediction can be attributed to the larger entrainment of the cough flow resolved in the URANS. These overpredictions decrease with time, forming a concave downward curve. Between 5 and 7 s, the overprediction decreases and is replaced by underprediction. The underprediction of the dispersion of these particles increases with time, resulting in a dispersion approximately 0.5 m smaller than that observed in the DNS.

The dispersion of particles in the 32–64 µm range shows a distinct behaviour compared to other particle clouds, being underpredicted in the URANS throughout the simulation. 20 s after the initiation of the cough, the dispersion error for these particles becomes smaller relative to the smaller particles, corresponding to only 30% of the dispersion observed for particles smaller than 32 µm. The best accuracy in dispersion prediction is observed for particles larger than 64 µm, where the URANS underpredicts the IQRs by only about 0.08 m at the end of the simulation, corresponding to about 16% of the final dispersion observed for particles smaller than 32 µm.

Overall, the dispersions of particles smaller than 32 µm resolved in the URANS are overpredicted during the early puff phase and then underpredicted for the remainder of the puff phase, with the worst accuracy. In contrast, the dispersion of medium sized particles (32–64 µm) is underpredicted throughout the simulation but with better accuracy than that of smaller particles. The dispersion of larger particles (64–256 µm) is predicted with the highest accuracy.

This study aimed to evaluate the performance of URANS compared to DNS in predicting cough-induced flow and particle behaviour within a background flow defined by a lid-driven LSC. The initial flow fields for the cough simulations were generated separately for the DNS and the URANS by running well-resolved Navier–Stokes and URANS of the LSC, respectively. In the well-resolved Navier–Stokes simulation of the LSC, we captured detailed turbulent structures of the background flow, providing a much more detailed representation of the flow patterns, while the URANS smoothed out these details. Although the time-averaged velocity magnitudes were similar, small differences were observed, likely due to the lower lateral and spatial resolution of the URANS and its treatment of turbulence. Although these differences may seem small, they could have a significant impact on particle transport.

Qualitative differences in the cough jets were observed during the jet phase. For example, the horizontal propagation was underpredicted by the URANS shortly after the cough was initiated. The spatially averaged horizontal velocities in the jet region indicated that the jets were dynamically similar, with comparable bulk velocities. However, the URANS showed a more intense backward flow in the background, suggesting that the bulk velocity measurements of the background flow along the horizontal axis were poorly predicted in the jet region. Further evaluation of the entrainment indicated that the shorter horizontal propagation of the jet in the URANS was due to an overall overprediction of entrainment, resulting in a shorter and wider jet compared to the DNS. Analysis of particle distributions at the end of the jet phase showed that particles smaller than 64 µm were mostly trapped within the cough jet and had similar lateral and horizontal expansions. Larger particles tended to detach from the cloud, with slightly fewer doing so in the URANS compared to the DNS, possibly due to the greater entrainment of ambient air.

At the end of the puff phase (t = 30 s), the cough puff in the DNS propagated almost 3 m, whereas in the URANS it was limited to 2 m, with the maximum horizontal propagation almost stopping after 15 s. The URANS overpredicted the vertical extent of the puff, with the puff almost reaching the ground due to the strong downward background flow. In contrast, the puff in the DNS interacted less with this downward flow due to lower entrainment and faster horizontal propagation. These differences in the horizontal and vertical velocity fields significantly affected particle cloud displacements and dispersions: horizontal particle displacements were underpredicted in the URANS during the early puff phase, especially for smaller particles, but this underprediction decreased with time. The vertical position of particles smaller than 32 µm was consistently underpredicted in the URANS, probably due to the greater descent of the cough puff. Vertical displacements of particles larger than 64 µm were well captured by the URANS, while medium sized particles (32–64 µm) were underpredicted during the early puff phase but overpredicted during the late puff phase, showing a different trend compared to other particle clouds.

Overall, the dispersion accuracy in the URANS deteriorated with decreasing particle size. The URANS overpredicted the dispersion of particles smaller than 32 µm during the early puff phase (up to 10 s) due to greater entrainment (especially in the lateral direction) but underpredicted it afterwards with a larger error, mainly due to the smaller horizontal dispersion observed in the URANS. In contrast, medium sized particles were consistently underpredicted throughout the URANS, even during the early puff phase.

These results highlight the need for a better representation of the time-averaged flow field and cough-induced particle dynamics in such scenarios. We conclude that fully resolving the turbulence of both the cough and background flows dramatically affects the resulting fields compared to a widely used time-averaged turbulence model-based simulation such as URANS. Our study highlights the weaknesses of this popular approach, particularly in predicting cough-induced flow. Another challenge is to accurately represent the particle dispersion effects of turbulent fluctuations, which are often calculated using particle dispersion models. In this study, we used a widely used stochastic particle dispersion model, but it was difficult to evaluate the model explicitly because the time-averaged flow field was not accurately represented by the URANS. Consequently, the dispersion error could be attributed to several factors. A possible approach to evaluate the particle dispersion model is to extract the spatially filtered and time-averaged flow fields of a cough puff directly from the DNS, which is a future goal of this study.

The work was supported by the Initiative and Networking Fund of the Helmholtz Association of German Research Centres (HGF) under the CORAERO project (KA1-Co-06). The authors gratefully acknowledge the scientific support and HPC resources provided by the German Aerospace Center (DLR). The HPC system CARO is partially funded by “Ministry of Science and Culture of Lower Saxony” and “Federal Ministry for Economic Affairs and Climate Action”.

Funding: Helmholtz-Gemeinschaft, No. KA1-Co-06.