Simulation-based study of COVID-19 outbreak associated with air-conditioning in a restaurant

The COVID-19 pandemic (caused by the SARS-CoV-2) has led to more than 69 × 10⁶ infections and 1.5 × 10⁶ deaths until now. A growing number of evidence has suggested that airborne transmission of respiratory droplets, either produced from asymptomatic or pre-symptomatic individuals or none, is a potential pathway that contributes to the wide spread of the disease.^1–4 In particular, an incident in a restaurant in Guangzhou, China, where a single asymptomatic COVID-19 patient caused the infection of eight people sitting in the same and two adjacent tables, has been widely reported in the media^5,6 and recent scientific literature¹ has suggested such a potential transmission pathway. Specifically, for airborne transmission, viruses attached to small respiratory droplets (typically <5 µm, referred to as aerosols hereafter) produced by even normal breathing or talking can stay suspended in the air for hours, and they can accumulate, particularly in indoor spaces with poor ventilation, imposing high infection risk to individuals who inhale them.^7–10 In addition, under high wind speed, large droplet size, low temperature, and high relative humidity, the airborne transmission distance can exceed the 2-m social distancing rule.^11,12 For instance, for the wind speed up to 15 km/h, the travel distance of aerosols along the wind direction can reach 6 m.¹¹ Under such circumstances, the commonly adopted social distancing rules can no longer serve as a proper preventive measure.^13,14 As a result, there is a dire need for science-driven risk assessment tools that can provide actionable and more appropriate guidelines under different practical settings. Despite numerous analytical and numerical modeling studies that provide the first order estimate of aerosol exposure and corresponding airborne infection risk,^15–20 only recently, have the computational fluid dynamics (CFD) tools been actively employed to provide a more precise risk assessment for aerosol exposure, especially its spatial and temporal variations in indoor spaces. Specifically, using Reynolds-averaged Navier–Stokes (RANS) simulation with Lagrangian particle models, Shao et al.²¹ investigated the airborne transmission risk associated with an asymptomatic patient under elevator, classroom, and grocery store settings. They introduced a local risk index and evaluated its spatial variation under different ventilation settings, which reveals indoor infection hot spots (corresponding to high risk of aerosol exposure) due to the inappropriate ventilation design. Mathai et al.²² used RANS to explore how the ventilation configuration of open and closed windows influences the airborne transmission in a passenger car cabin. The aerosols are modeled as passive scalars by an advection–diffusion equation. Abuhegazy et al.²³ investigated the airborne transmission in a classroom under various configurations of particle sizes, source locations, barrier locations, and opening windows by RANS. They found that both barriers and opening windows mitigate the airborne transmission and highlighted the importance of the filtration system on the air conditioners to remove the indoor aerosols. Shao and Li²⁴ employed both RANS and the Wells–Riley equation to evaluate the infection risk of Biden, Wallace, and the audiences from the potential sources of infection (i.e., Trump and the first lady) in the first presidential debate in 2020. On a different subject, Wang et al.²⁵ and Li et al.²⁶ explored the virus transmission during the flushing of toilets by numerical simulation. They found that the virus-laden particles can reach up to 1 m from the ground and alerted the potential massive infection through public restrooms. A novel toilet design is proposed by Wu et al.²⁷ to reduce flushing-induced plumes. Wang et al.²⁸ used RANS to study the probability of cross-infection in a hospital ward with ten ventilation configurations and found that a bottom-in-top-out ventilation configuration and short distance between the outlet and the mouths of patients benefit the reduction in the cross-infection risk. Vuorinen et al.²⁹ performed large eddy simulations (LES) and Monte Carlo simulations to understand the exposure risk of the moving people in a public indoor environment where an infected person coughs. They showed that the exposure time for getting infected can range from seconds to hours depending on the aerosol diameter and indoor turbulence level, and the spatial extent of the high-risk zone can be as large as 4 m from the infected person. In addition, using LES, Pendar and Páscoa¹⁴ studied the influence of initial velocity, aerosol size distribution, injection angle, and mouth opening area on the indoor airborne transmission and suggested that wearing masks and bending head while sneezing can effectively reduce the droplet traveling distance. Remarkably, Li et al.³⁰ conducted a simulation of airflows in the aforementioned restaurant case and used the results to infer the strong influence of air conditioning and ventilation on the reported breakout. In spite of these vigorous efforts, there have been very few CFD studies so far that can show the values and the benefits of using CFD to assess airborne infection risks in indoor spaces and to illustrate the corresponding detailed physical mechanisms that lead to airborne disease transmission in the reported real-life infections. Such work is crucial not only to demonstrate the validity of CFD to capture airborne transmission processes but also to help us develop pinpointed interventions for mitigating the disease spread.

To fill in the gap, in this paper, we present a systematic LES based investigation of indoor airflow and the associated aerosol transport in the aforementioned restaurant setting, where likely cases of airborne infection caused by asymptomatic individuals were reported and the detailed information of infection process through contact tracing are available. Our simulation employs an advanced in-house CFD solver and uses cutting-edge numerical methods [i.e., advanced immersed boundary (IB) method, stochastic modeling of Brownian motion and effect of subgrid-scale (SGS) flow structures on aerosols, and Cunningham slip correction] to resolve complex indoor processes, including turbulence, flow–aerosol interaction, thermal effect, and the filtration effect by air conditioners. Such a combination of high-fidelity numerical methods, several of which has not been adopted in the past simulation work of indoor airflow, allows us to depict a detailed picture of airborne transmission in this hallmark incident and derive a series of important insights that can benefit our understanding of this transmission pathway and the development of more effective preventive measures. We compare our results with the ones from reported contact tracing and find that they are overall consistent. This agreement provides further support for the potential airborne transmission pathway of COVID-19 in addition to the evidence provided in the literature. The rest of this paper is organized as follows. The methodology used in our simulation is described in Sec. II. Subsequently, we provide the results of simulation under different ventilation and thermal settings in conjunction with their potential connection with reported infection in Sec. III, with the conclusions given in Sec. IV.

In this study, LES is conducted, where the SGS stress is modeled by the dynamic Smagorinsky model.^31–33 The governing equations are

where u is the velocity vector, ρ is the density, p is the pressure, ν is the kinematic viscosity, S = (∇u + ∇u^T)/2 is the strain rate tensor, f_B is the buoyancy force, τ_sgs is the SGS stress tensor, f represents the forces exerted by the solid structures inside the room in the IB method, and the overbar denotes the spatial filter in LES. In the case where the turbulence in the inner boundary layer on the walls is not fully resolved, a log-law-based wall model³⁴ is used. The second-order central differencing scheme is used for the spatial discretization of both the convection and diffusion terms, and the second-order Runge–Kutta method is used for time integration. The equations are solved by the fractional-step method, and the discretized Poisson equation for pressure is solved by the Portable, Extensible Toolkit for Scientific Computation (PETSc).³⁵ 

The thermal effect is accounted for in the simulation due to the consideration that the human bodies around tables and the hot foods on the tables can lead to temperature variations, which create thermal plumes in the room. Heat transfer is simulated by an advection–diffusion equation for temperature as

where T denotes the fluid temperature and α = 22.39 × 10⁻⁶ m²/s is the thermal diffusivity. To couple the thermal effect with the airflow dynamics, the Boussinesq approximation³⁶ is applied to Eq. (2),

where g is the gravitational acceleration, β = 3400 × 10⁻⁶/C is the coefficient of thermal expansion, and T_o is the reference temperature.

The influence of the solid structures on the flow dynamics is simulated by the IB method.^37–41 This method can be used on a static Cartesian grid to comply with the structure geometry and is thus convenient to apply to structures with complex geometry. In the IB method, the velocity boundary condition is satisfied by interpolating the velocity to the grid points adjacent to the structure surfaces according to the velocity profile of the boundary layer. Linear and nonlinear interpolations are applied for cases where the boundary layer is resolved and modeled, respectively.

The transport of aerosols is simulated with a Lagrangian particle tracking algorithm. The aerosol location x_p is solved by

where u_p is the aerosol velocity obtained by

In the above equation, the first and second terms on the right-hand side represent the Stokes drag force and gravity/buoyancy force, where ρ_p and ρ_f are the density of the aerosols and the fluid, respectively, u_p is the aerosol velocity, u_f is the velocity of the fluid at the aerosol location, d is the aerosol diameter, $V= 1 6 π d 3$ is the aerosol volume, μ is the fluid dynamic viscosity, $R e p = ρ f d u f − u p /μ$ is the aerosol Reynolds number, $ϕ ( R e p ) =1+0.15R e p 0.687$⁠, and C_C is the Cunningham correction factor.⁴² Note that the Cunningham correction is applied to the Stokes drag force because the aerosols are so small that the surrounding gas cannot be modeled as a continuum medium and the non-slip boundary condition on the aerosol surface cannot be strictly applied.^43–45 The third term of the right-hand side represents the Brownian motion effect, where each component of W(t) is a Wiener process with zero mean and variance 2γk_BT, with $γ= 3 π μ d C C ϕ ( R e p )$ and k_B being the Boltzmann constant.^46,47 Both Eqs. (5) and (6) are numerically solved by an Adam–Bashforth method. To capture the aerosol deposition effect, an aerosol whose center has a distance from the wall less than its radius is immediately removed from the computation domain.

The effect of the SGS flows on the aerosol dynamics is accounted for by a stochastic model.^48,49 In this model, the fluid velocity at the aerosol location u_f is decomposed into a grid-resolvable part u_f,r and a SGS part u_f,s,

where u_f,r is obtained by interpolating the flow velocity values on the grid nodes to the aerosol location by a fourth-order interpolation scheme and u_f,s is modeled by the following Langevin equation:

where $T L s =2 σ s 2 / C 0 ε ̄$⁠, C₀ = 3, $σ s 2 =2 e s /3$⁠, $e s = ( Δ C s S ̃ ) 2 /0.3$⁠, $ε ̄ = c ε e s 3 / 2 /Δ$⁠, c_ɛ = 0.93, Δ = (ΔxΔyΔz)^1/3, C_s is the Smagorinsky coefficient, $S ̃$ is the filtered rate-of-strain tensor, and each components of dξ is a Gaussian random variable with zero mean and variance dt. Equation (8) is integrated in time with an implicit algorithm for numerical stability.

The aerosol volume fraction in the present study is ϕ = 4 × 10⁻⁷, which is much smaller than 10⁻³, the threshold of four-way coupling.⁴² Therefore, the aerosol–aerosol interaction, including aerosol collision and merging, is not considered.

Our numerical tools are validated by a forced convection case and a mixed convection case. See the  Appendix for the validation details.

The present numerical study is based on a real infection event that occurred in a restaurant in Guangzhou, China, which is illustrated in Fig. 1. On January 24, 2020, there were 89 customers dining in the restaurant during the lunch time. After this lunch, nine people sitting on table numbers one, two, and three [labeled in Fig. 1(a) as T1, T2, and T3] were infected. The other people sitting away from these three tables were not infected.³⁰ It has been believed that three people eating at T1 and two people at T3 were infected directly during this lunch.¹ Virus was transmitted from patient zero sitting at T2. In order to simplify the discussion, customers are labeled as TAPBs, as shown in Fig. 2, where A stands for the index of the tables and B stands for the index of the customers.

Our numerical simulations are set up with the configuration shown in Fig. 1. A three-dimensional Cartesian coordinate system is used in the mathematical formulation and numerical discretization of the computational domain, in which the x-direction is set parallel to the longer wall of the room, the y-direction is set as the opposite of the direction of gravity, and the z-direction is set parallel to the shorter wall of the room. The origin is located at a corner shown in Fig. 1(a). The length of the longer wall is 17.5 m, which is located at z = 8.3 m. The lengths of the two shorter walls are 8.3 m and 6.68 m, which are located at x = 0 m and x = 17.5 m, respectively. The other wall has three segments, as shown in Fig. 1(a). The tables in this restaurant are captured in the simulation using the IB method. There are two types of table, one is round-shaped and the other is rectangular-shaped. The round-shaped tables have two sizes with diameters d = 1.2 m and d = 1.8 m, respectively, while the rectangular-shaped tables have two sizes as 0.9 × 0.9 m² and 1.2 × 0.9 m². The height of all the tables in the simulations is set as h_t = 0.75 m, based on our understanding of the typical height of the tables used in China.

There are five working HVAC (heating, ventilation, and air conditioning) systems in this room, as shown in Fig. 1(a), which are labeled in short as AC1, AC2, AC3, AC4, and AC5 and are also simulated by IB method. Four of them are installed on the wall at z = 8.3 m, and the other one is installed on the wall at x = 0 m. AC1, AC4, and AC5 have a size of 1.6 × 0.5 × 0.25 m³ (length × width × height), while AC2 and AC3 have a size of 1 × 0.5 × 0.25 m³. Each AC has both an inlet and an outlet, as illustrated in Fig. 1(a). At the inlet, air flow is sucked into the AC in the y−direction (i.e., upward vertically) with a ventilation rate q_AC. At the outlet, airflow is blown out of the AC horizontally with the same ventilation rate q_AC. On the wall at x = 0 m, two rectangular regions [colored by dark blue in Fig. 1(a)] are used to simulate the fire-door and the exhaust vent in the restroom of this restaurant, respectively. Fresh air is sucked into the restaurant room through the fire-door, while the polluted air exits through the restroom. The setup of the restroom is simplified since virus-laden droplets are mostly concentrated near AC1. The ventilation rates of the fire-door and the restroom are both set as q = 0.044 m³/s following the literature.³⁰ 

To simulate the thermal effect by the human bodies around the tables and the hot food on the tables, a simplified model is implemented, as illustrated in Fig. 1(b). Above each table, a hemispherical hot region is defined, which has a diameter the same as that of the round table below it or the diagonal length of the rectangular table below it. In each hot region, the temperature is set by a constant T = T_h. Moreover, the temperature of the air blown out of the AC outlets is set to be a constant T = T_r, which is lower than T = T_h. The temperature difference between them is denoted by ΔT = T_h − T_r. The dimensionless temperature

is used to quantify the temperature variation in the flow field. The adiabatic condition is applied at the walls, floor, and ceiling. The initial temperature in the room is set as T = T_r in the simulation.

In the simulation, the height of the mouths of all the customers is set as h_m = 1.3 m, based on the statistics of the averaged height of the sitting people.⁵⁰ Virus-laden aerosols are released from patient zero, which is marked by the red dot in Fig. 1(a) aerosols are initialized following the measurement by Shao et al.²¹ It is assumed that each breath cycle of patient zero includes a 2.37 s inhale period and a 1.58 s exhale period. During the exhale period, the aerosols are released uniformly in a circular disk with a diameter of 40 mm, which is a size based on the averaged diameter of mouths.⁵¹ The velocity of the aerosols is set as the same as the breathing flow, which is pointing to the center of T2. In the simulations, there are 44 aerosols released during each exhale period, based on the measurement result.²¹ The aerosols are set to be spherical, and the diameter of the aerosols is set as 1.5 μm following the peak value of the measured size distribution.²¹ We ignored the evaporation, splitting, and merging of droplets such that the aerosol diameter is constant. The density of the aerosols is set to be 10³ kg/m³, the same with that of water-like saliva droplets.^11,52 Based on the research by Alsved et al.,⁵³ the number of the aerosols released per second by normal talking and that by loud talking are, respectively, about two and four times as that by normal breathing. In order to consider the effect by the people normal talking or loud talking in the restaurant, cases with higher numbers of aerosols released per breath cycle have also been simulated, which produced similar qualitative results compared with the cases using the normal breathing setting. In the following sections (Secs. III and IV), only the results from the normal breathing setting are reported.

The filtration effect of the ACs and aerosol deposition is considered in our simulations. As shown in Fig. 1(a), the filtration effect is quantified by the filtration efficiency f. Once N aerosols have been sucked into an AC, the N(1 − f) of them will be re-ejected into the room randomly from the outlet of the AC. In the simulation, the filtration efficiency is set as f = 80%, based on the typical performance of HVAC systems.⁵⁴ To simulate the deposition effect, an aerosol will be labeled as deposited and deleted from the computational domain once the distance between the aerosol and the solid surface close to it (walls, tables, and ACs) is smaller than the radius of the aerosol.

The computational domain in the present study has a grid resolution of N_x × N_y × N_z = 200 × 200 × 200. The Reynolds number is defined as

and has the value Re = 6.8 × 10⁴, where ρ_a is the density of air, U_o = 1 m/s, and L_o = 1 m/s. Four cases with various ventilation rates q_AC and temperature difference ΔT are considered, as shown in Table I. The range of q_AC covers the common range of the ventilation rate of a HVAC system.⁵⁵ The temperature difference is based on the consideration of the human body temperature, the food temperature, the room temperature, and the implementation of the simplified model shown in Fig. 1(b).

The simulations were run with sufficient duration for the flow field to fully develop. Figure 3 shows the time evolution of the velocity variances $u ′ u ′$⁠, $v ′ v ′$⁠, and $w ′ w ′$ in case AC112DT10, where u′, v′, and w′ are defined in Eq. (13) and $⋅$ is the spatial–temporal average defined as

where t₀ −Δs < t < t₀ + Δs is a time window expanding from t = t₀ that is long enough to cover the temporal fluctuations and V_room is the whole volume of the room. The aerosols are released after the flow in the room has fully developed, as illustrated in Fig. 3. Then, the simulations continued until the total number of aerosols in the whole computational domain is statistically converged as a result of the balance between the number of newly released aerosols and that of deposited, filtered, and exiled aerosols.

In this subsection, the indoor airflow structures are discussed to provide a general picture of the flow field in the simulated cases. The distributions of mean flow structures and turbulent flow structures under different ventilation rates and thermal settings are illustrated, which will be connected to the aerosol exposure index and the infection risk of COVID-19 discussed in Secs. III B and III C, respectively.

To investigate the indoor flow structure, a temporal average is defined as

where t₁ is a selected time instance after the total number of aerosols reaches dynamical balance in the simulation. Based on the temporal average, the flow field is decomposed into two parts: the mean flow and the turbulence,

where $u ̄$ is the mean flow velocity and u′ is the turbulent velocity fluctuation.

Figure 4 shows the streamlines of the mean flow traced forward from the outlets of the ACs under different ventilation rates and thermal settings. The flow structures in the cases without the consideration of thermal effect are displayed in Figs. 4(a) and 4(c). It is observed that a recirculation flow with the size of the width of the room is formed in front of the ACs in cases AC14DT00 and AC112DT00, where the air ejected from the ACs first moves along the initial jet direction, turns to the floor direction once it reaches the wall, moves downward until reaching the floor, turns back toward the ACs, moves under the tables, and, in the end, is sucked into the ACs through the inlets. In the cases with the consideration of thermal effect, the flow structures are more complex, as shown in Figs. 4(b) and 4(d). It is found that the cold air blown from the ACs moves downward at an angle due to the interaction between the cold jet of a higher density and the rising thermal plume of a lower density near T3. Before reaching the floor, the cold air is heated by the interaction. Comparing cases AC14DT10 and AC112DT10, it is illustrated in Fig. 4(b) that in case AC14DT10, the heated air near the floor that originates from the jet from AC1 moves in various directions. In case AC112DT10 [Fig. 4(d)] with a stronger jet momentum, the air from AC1 moves in the initial jet direction until reaching the wall near z = 0 m. For the same reason, the jet inclination angle in case AC14DT10 is larger than that in case AC112DT10.

The exchange between the indoor and outdoor airs is investigated. To reduce the concentration of virus-laden aerosols, besides applying the filtration system, the most commonly used method is to keep exchanging indoor and outdoor airs. Fresh outdoor air can reduce the indoor aerosol concentration through the dilution effect. Figure 5 shows the streamlines of the mean flows traced backward from the restroom. As shown, the range covered by these streamlines varies with indoor conditions. The shortest range is from the tables in front of AC4 to the restroom, which occurs in cases AC14DT00 and AC112DT00, while the longest range is from the space in front of AC1 to the restroom, which happens in case AC112DT10. Figure 5 also shows that there are barely any streamlines near T1-3 that are directly connected to the restroom, and thus, it is less likely for the aerosols in this area to exit the restaurant, which is consistent with the measurement by Li et al.³⁰ 

The turbulence is quantified by the turbulent kinetic energy (TKE), defined as

and is displayed by the contours plotted near T1-3 in Fig. 6. In Figs. 6(a) and 6(c), it is illustrated that the turbulence is mostly concentrated around the jet from AC1, in the two cases without considering the thermal effect. Moreover, the TKE in Fig. 6(c) is higher than that in Fig. 6(a) in most of the displayed area, which is caused by the stronger entrainment by the jet in case AC112DT00 than in case AC14DT00. However, the distribution of turbulence is different in the cases with the consideration of thermal effect. In Figs. 6(b) and 6(d), it is found that the turbulence is mostly concentrated near the region where cold air interacts with the thermal plume, which is located above T3 in case AC14DT10 and has a broader range in case AC112DT10. Relatively weaker turbulence is also observed from Figs. 6(b) and 6(d) above the surfaces of T1 and T2, which is produced by the instability of the rising thermal plume. Comparing Figs. 6(b) and 6(d), it is shown that the turbulence in case AC112DT10 is stronger than that in case AC14DT10 in most of the area due to the stronger interaction between the cold jet and the thermal plume in case AC112DT10. Furthermore, comparing Figs. 6(b) and 6(d) with Figs. 6(a) and 6(c), it is found that the cases with the consideration of thermal effect show more intense turbulence near T1-3 than the cases without the thermal effect. In summary, different ventilation rates and thermal settings give rise to different distributions and intensity levels of the indoor turbulence. Another conclusion that is not directly shown from these results but is reasonable to expect is that the relative location between the tables and ACs also impacts the distribution of indoor turbulence.

Before using the simulated data to predict the COVID-19 infection risk through airborne transmission, first, in this section, we illustrate and analyze the distributions of aerosols under different ventilation rates and thermal settings. Explanation on the aerosol distribution is also given in this section, based on the flow structure illustrated in Sec. III A.

The indoor spatial–temporal distribution of the aerosols under various ventilation rates and thermal settings is quantified by the aerosol concentration C(x, t), which is defined as the number of aerosols per unit volume. The mean value of C during the statistical steady state, which is also called aerosol exposure index in this study, is used to quantify the temporally averaged distribution of aerosols or the degree of being exposed to the aerosols in the restaurant, which is numerically calculated in the present simulations as

where n_i,j,k is the total number of aerosols that have been in the grid cell (i, j, k) at each time step during the time period t₁ < t < t₂, ΔV is the volume of the grid cell, and n₁₂ is the total time steps for the simulations running from t₁ to t₂.

In the simulations, the aerosols are transported by three forces on the right-hand side of Eq. (6), i.e., the Maxey–Riley equation, including the Stokes drag force, the gravitational force, and the Brownian motion effect force. It is found from the results that in the simulated cases, the Stokes drag force f_stokes is dominant over the other two forces in most of the computational domain. It can be further decomposed into four parts by substituting Eq. (13) into Eq. (6),

where $u ̃$ is the ensemble averaged velocity of the aerosol, f_m is the force due to the mean flow velocity, f_t is the force due to the turbulent fluctuation, f_pe is the force due to the mean velocity of the aerosol, and f_pt is the force due to the fluctuating velocity of the aerosol. Based on the above decomposition, an aerosol, in general, experiences less Stokes drag if its mean velocity $u ̃ p$ is closer to the local mean flow velocity $u ̄ f$ because of the opposite effects of f_pe and f_m. Therefore, the mean flow tends to drag an aerosol to move along with its streamline by the combined effect of f_m and f_pe. The local turbulence tends to randomly disperse a group of aerosols^56–60 through the fluctuating shear forces f_t + f_pt. The aerosol exposure index in the statistically steady state results from the combination of all these effects.

Figure 7 shows the iso-surfaces of the aerosol exposure index $C ̄$ for $C ̄ =0.1 C ̄ max$ (in black) and $C ̄ =0.01 C ̄ max$ (in yellow) under different ventilation rates and thermal settings, where $C ̄ max$ is the maximum of $C ̄$⁠. Comparing the iso-surfaces of $0.1 C ̄ max$ and $0.01 C ̄ max$⁠, we can see that the former is closer to patient zero and located around T2, while the latter reaches T1 and T3 in cases AC14DT00, AC112DT00, and AC14DT10. In case AC112DT10, it is observed from Fig. 7(d) that the iso-surface of $0.01 C ̄ max$ is preferentially concentrated around T3, in contrast to T1. The values of C_max in the different cases are C_max = 2740.9 in case AC14DT00, C_max = 783.6 in case AC14DT10, C_max = 276.9 in case AC112DT00, and C_max = 279.3 in case AC112DT10. Therefore, there are more aerosols concentrated inside the iso-surface of $0.01 C ̄ max$ in case AC14DT00 than the other three cases because the stronger airflows due to the higher q_AC and thermal effect in cases AC14DT10, AC112DT00, and AC112DT10 drive more aerosols to collide and deposit on the walls near T1-3 so that the total numbers of aerosols have lower values in the latter three cases. Comparing the range reached by the iso-surfaces of $0.01 C ̄ max$ in the x−direction with that reached in the z−direction, we can observe that in all the cases, the aerosol exposure index decays faster along the x−direction than in the z−direction since the mean flow drag f_m is dominant in the z−direction over in the x−direction due to the jet flow from AC1. Moreover, the topography of the iso-surface of $0.01 C ̄ max$ varies with the ventilation rates and thermal settings. Figures 7(a) and 7(c) display L-shaped contours of $C ̄$ inside the iso-surface of $0.01 C ̄ max$ in cases AC14DT00 and AC112DT00. The lower aerosol exposure index in the region between the two sides of the L-shape is due to the turbulent dispersion around the jet flow. Figures 7(b) and 7(d) show a more complex distribution of the aerosol exposure index with a rougher topography of the iso-surface of $0.01 C ̄ max$ in cases AC14DT10 and AC112DT10 because the rising plumes due to the thermal effect shown in Figs. 4(b) and 4(d) create a more complex drag force field. Moreover, Fig. 7(d) shows a preferential concentration on the T3 side than on the T1 side due to the transport by the stronger mean flow in case AC112DT10 compared with case AC14DT10.

Since the infections at T1 and T3 are believed to be directly related with the airborne transmission,³⁰ much research attention is paid to the regions close to T1 and T3. The iso-surface of $C ̄ / C ̄ max =−1.8$ (the threshold is selected as a representative one based on the isolines concentrated near T1 and T3, see Figs. 9 and 11) and the mean flow streamlines near T1 and T3 are shown in Figs. 8 and 10, respectively. The isolines of $C ̄ / C ̄ max$ near T1 and T3 at the height of mouth h_m are also displayed in Figs. 9 and 11, respectively. The correlations between the contours of the aerosol exposure index and the mean flow streamlines are evident.

Near T1, the mean flow streamlines in Figs. 8(a) and 8(c) indicate that there exists a recirculation zone above T1 in cases AC14DT00 and AC112DT00. Inside the recirculation zone, the streamlines originate from the region around T2, first move toward the negative z−direction, and then gradually turn to T1 with a spiral-like trajectory. The recirculation zone above T1 is generated by the interaction among the walls, the surface of T1, and the local airflows. It is illustrated in Figs. 9(a) and 9(c) that from $C ̄ / C ̄ max =−0.5$ to −1.8, the isolines of $C ̄ / C ̄ max$ near T1 tend to be elongated while turning to T1. This pattern of spreading is due to the convection by the mean flow in the recirculation zone observed in Figs. 8(a) and 8(c). In Fig. 8(b), two peaks on the iso-surface of $C ̄ / C ̄ max =−1.8$ are observed around T1 in case AC14DT10. The left peak is associated with a recirculation zone attached to the wall, which is formed by the thermal plume rising from T1 indicated by the red streamlines. Another peak is located between T1 and T2, where hot air passes through, indicated by the streamlines together with the green arrows in Fig. 8(b). Corresponding to the two peaks observed in Fig. 8(b), it is illustrated in Fig. 9(b) that from $C ̄ / C ̄ max =−0.5$ to −2, the isolines of $C ̄ / C ̄ max$ are preferentially distributed in the region between T1 and the wall, where there is the recirculation zone, and in the region between T1 and T2, where the hot airflow caused by the thermal passes by. The spreading of aerosols to the first region is due to the convection by the thermal plume that generates the recirculation, and the spreading to the latter region is due to the convection by the hot airflow passing by. As for case AC112DT10, it is found in Fig. 9(d) that the magnitude of the gradients of the aerosol exposure index $|∇ C ̄ |$ around the isoline of $C ̄ / C ̄ max =−1.7$ is significantly larger than those in the other cases, which indicates the less spreading of the aerosols in the region colored by dark blue in Fig. 9(d). The less spreading is consistent with the observation in Fig. 8(d). Near T1, the streamlines are directly connected to AC1, indicating that the aerosols in the dark blue region in Fig. 9(d) are mostly transported by the mean flow along these streamlines. However, there are less aerosols concentrated in front of AC1, as shown in Fig. 7(d). Therefore, less aerosols are transported to the region near T1 in case AC112DT10. The irregular distribution of the isolines of $C ̄ / C ̄ max$ in Fig. 9(d) also indicates that the turbulence due to the thermal effect around T1, which is also illustrated in Fig. 6(d), plays a role in dispersing the aerosols by exerting the turbulent shear force f_t + f_pt on the aerosols.

To illustrate the aerosol exposure index near T3, Figs. 10(a) and 10(c) show the streamlines of the mean flow inside the iso-surface of $C ̄ / C ̄ max =−1.8$ in cases AC14DT00 and AC112DT00. Along these streamlines, the air moves from the region close to patient zero to the wall and the inlet of AC1. The isolines of $C ̄ / C ̄ max$ displayed in Figs. 11(a) and 11(c) are elongated in the same direction around T3 from $C ̄ / C ̄ max =−1.2$ to $C ̄ / C ̄ max =−1.9$⁠. The elongation indicates the aerosol spreading, which is due to the transport by the mean flow shown in Figs. 11(a) and 11(c). Comparing Fig. 11(a) with (c), the isolines in the latter figure are sharper near the spot below AC1. This is because the stronger suction in case AC112DT00 exerts more forcing on the aerosols, which attracts more aerosols to AC1 compared with case AC14DT00. As shown in Figs. 11(b) and 11(d), in both the cases AC14DT10 and AC112DT10, there exists a peak of the iso-surface of $C ̄ / C ̄ max =−1.8$⁠, which is located above the surface of T3. The isoline of $C ̄ / C ̄ max =−2.1$ in Fig. 11(b) indicates a local high aerosol exposure index region, which corresponds to the peak observed in Fig. 10(b) in case AC14DT10. In addition, in Fig. 11(d), from $C ̄ / C ̄ max =−1.4$ to $C ̄ / C ̄ max =−2.2$⁠, the elongation of the isolines above T2 corresponds to the peak observed in Fig. 10(d). The streamlines in parallel with the green arrows in Figs. 10(b) and 10(d) indicate that the aerosol spreading in this region is due to the mean flow convection associated with the heated air to interact with the cold air in the jet above T3. Comparing the peak observed in case AC14DT10 with the one observed in case AC112DT10, the latter is located lower because in the latter case, the cold jet from AC1 has stronger momentum, which makes the converging point closer to the surface of T3.

In this section, first, a discussion about using the simulated results to predict the infection risk of COVID-19 is presented. Then, the simulation-based distribution of the infection risk is illustrated and explained by linking it to the analysis in Sec. III B.

There are various factors that can determine the total infection risk of COVID-19 associated with airborne transmission. According to Mittal et al.,⁶² the total infection risk R_in is the product of the risk functions of a series of complex factors and processes, which can be expressed as

where R_fm is the infection risk due to the fluid mechanical factors (e.g., ventilation rate, temperature difference, and initial velocity of droplets) and R_ph is the infection risk due to the physiological factors (e.g., the health condition of different people and the survival rate of the virus during transmission). It is important to clarify that in the present study, only R_fm, which is the risk due to the fluid mechanical factors, is analyzed based on the CFD results. The infection risk R_fm is connected to the simulated aerosol exposure index distributions by an infection model,⁶² which is illustrated in Fig. 12. An exposure region in front of a person with volume 2000 ml is defined, and R_fm is defined by averaging the mean virus-laden aerosol exposure index inside the volume as

where V_e is the volume of the exposure region.

The distributions of R_fm of the customers sitting around T1 are shown in Fig. 13. For the cases without the consideration of thermal effect, it is observed from Figs. 13(a) and 13(c) that in cases AC14DT00 and AC112DT00, T1P1 has the highest R_fm, T1P3 has the lowest R_fm, and T1P2 and T1P4 are in between. As illustrated by the aerosol exposure index and streamlines shown in Sec. III B, this distribution is due to the recirculation zone above T1, which creates a negative gradient of $C ̄$⁠, $−∇ C ̄$⁠, roughly pointing from T1P1 to T1P3. For the cases with the consideration of thermal effect, Fig. 13(c) shows the similar trend in case AC14DT10 compared with the observations from Figs. 13(a) and 13(b). However, the mechanism that leads to this distribution is different. T1P2 is located in a region where $−∇ C ̄$ is roughly pointing from T1P1 to T1P2 due to the local thermal plume-formed recirculation zone, while T1P4 is located in a region where $−∇ C ̄$ is roughly pointing from T1P1 to T1P4 due to the hot airflows passing by. T1P3 is located the furthest to T1P1. Moreover, comparing Fig. 13(d) with Figs. 13(a), 13(c), and 13(d), it is observed that the R_fm value in case AC112DT10 is around 10% of those in cases AC14DT10 and AC112DT00 and 1% of that in case AC14DT00 because T1P1–T1P4 are all located in the region where aerosols are less concentrated and most of the local airflows are directly from AC1. In the real infection event that occurred in this restaurant, T1P3 is the only person who was not infected during the lunch.³⁰ Without considering the factors beyond fluid mechanics, the present predictions around T1 in cases AC14DT00, AC14DT10, and AC112DT00 show a same trend as in the real infection.

The R_fm of the customers sitting around T3 is displayed in Fig. 14. It is observed from Figs. 14(a)–14(c) that the distributions of R_fm in cases AC14DT00, AC14DT10, and AC112DT00 have an approximate V-shape, where the two ends have higher R_fm than those in between. The difference is that the ratio between the left end in the V-shape distribution to the right is lower in cases AC14DT00 and AC112DT00 than in case AC14DT10. The V-shape distributions observed in cases AC14DT00 and AC112DT00 are due to the airflow sucked into AC1, which creates a distribution of $−∇ C ̄$ above T3 roughly pointing from T3P7 to T3P4 and from T3P6 to T3P3. However, the V-shape distribution observed in case AC14DT10 is due to a more complex mechanism. Around T3P1 and T3P2, $−∇ C ̄$ is created by the hot airflows interacting with the cold jet, while around T3P7 and T3P6, $∇ C ̄$ is created by the hot airflow being sucked into AC1. Moreover, a distribution of a different shape is observed in case AC112DT10 in Fig. 14(d), where the R_fm values of T3P6 and P3P7 are much larger than those of the other five people (the ratios are higher than 2). Similar to the mechanism in case AC14DT10, the distribution of $−∇ C ̄$ is due to the combined effect of the hot airflows being sucked by AC1 (around T3T7 and T3T6) and the hot airflows interacting with the cold jet (T3P2). The observed much higher risks of T3P6 and T3P7 are due to the stronger cold jet in case AC112DT10 that drives the concentrated region close to the surface of T3. In the real infection event, it is believed that T3P3 and T3P6 were infected during the lunch directly.³⁰ The prediction on the distributions of R_fm does not show a higher R_fm of T3P3 and T3P6 than other people. However, it shows the approximate V-shape distributions in cases AC14DT00, AC112DT00, and AC14DT10, which indicate that the people sitting at the two ends of the table have more risks than the ones in between, which is consistent with the real infection at T3.

In this section, further discussions are provided based on the simulated results.

First, we would like to offer some remarks on the applicability of our CFD for infection risk prediction. Specifically, our results have shown that CFD can provide detailed information on the spatial variation of aerosol concentration and aerosol exposure time as well as the aerosol transport process to illustrate the airborne transmission process and show remarkable agreement with infection pattern reported in the real scenario. Moreover, the CFD results can be used to systematically evaluate the effect of complex environmental factors such as ventilation design, thermal plumes, and turbulence on indoor airborne transmission risks. Such information may assist the development of more accurate reduced order models for risk assessment and the implementation of preventive measures under different indoor settings. Nevertheless, we acknowledge that the infection risk evaluation from our current CFD is only derived from the aerosol exposure index. To yield a more substantiated metric of infection risk, a relevant infection-dose model, currently not available for SARS-CoV-2, is needed. In addition, our prediction of the aerosol exposure index can be further improved if the more detailed information about the geometry of the space, ventilation, air conditioning, thermal condition, and human behaviors is accessible.

Second, beyond a quantification of the aerosol exposure index, we would like to highlight the ability of our CFD to determine the potential airborne transmission pathways, which can directly lead to actionable preventive measures. Here, we use a reverse-time tracing method to determine the origin of aerosols that cause the possible infection exposure of an individual. Specifically, for example, Fig. 15 depicts two intriguing pathways that lead to the infection of T3P3. One, marked by the blue and green aerosols in the figure, shows the potential infection exposure associated with the aerosol transport under the tables. This result suggests the need to shield the space under a table, in addition to the space above, when the interplay between ventilation flow and thermal plumes induces the complex indoor airflows. The other pathway, highlighted with red aerosols, points to the infection potentially caused by the returning aerosols from the AC due to the limited filtration efficiency. To further quantify the relative significance of such a transmission pathway, we introduce an AC exposure fraction factor, defined as

where n_e and n_eac are the total number of aerosols and the number of aerosols from AC1 that enter the exposure region during the averaging time. From Fig. 16(a), it is shown that f_ac can reach as high as 30% for the individuals sitting around T1 in case AC112DT10 due to the directional transport of aerosols through the airflows generated straight from AC1, as shown in Fig. 8(d). For the same reason, in cases AC14DT00 and AC112DT00, T1P3 yields the peak value of f_ac at about 15%. Among all the cases considered, the lowest f_ac takes place in AC14DT10 because the airflow from AC1 is strongly disturbed by the rising thermal plumes causing more omni-directional dispersion of aerosols in this case. In contrast, the T3P1 in case AC112DT10 has the maximum f_ac among all the cases due to the directional transport of aerosols through the jet flow from AC1. Finally, at T3, the f_ac peaks at T3P3 due to its location adjacent to AC1 with values of about 20% and 11% for the cases AC14DT00 and AC112DT00, respectively. These results indicate that the aerosols returning from ACs is a potentially important factor when evaluating infection risks, and an increase in the AC filtration efficiency is highly desirable for mitigating indoor infection risks.

In this paper, we presented a systematic CFD-based investigation of indoor airflow and the associated aerosol transport in a restaurant setting, where likely cases of airborne infection of COVID-19 caused by asymptomatic individuals were reported and detailed information of infection process through contact tracing is available. We employed an advanced in-house LES solver and other cutting-edge numerical methods (i.e., IB method, stochastic modeling of Brownian motion and effect of SGS flow structures on aerosols, and Cunningham slip correction) to resolve complex indoor processes simultaneously, including turbulence, flow–aerosol interplay, thermal effect, and the filtration effect by air conditioners. Using the aerosol exposure index derived from our simulation, we are able to provide a spatial map of the airborne infection risk in the restaurant under different air conditioning and thermal settings. Our results have shown a remarkable direct linkage between regions of the high aerosol exposure index and the reported infection patterns in the restaurant, providing strong support to the airborne transmission occurring in this widely-reported incident. Moreover, using flow structure analysis based on the mean flow streamlines and the distribution of turbulent kinetic energy and reverse-time tracing of aerosol trajectories in the simulation, we are able to further pinpoint the influence of different environmental parameters (i.e., thermal, ventilation rate, and filtration efficiency) on the infection risks and various potential airborne transmission pathways that can lead to the infection. Specifically, the thermal plumes caused by the temperature difference between ambient air and human/table are shown to cause local recirculation flows in low-ceiling/confined spaces that can substantially increase infection risks. In addition, the interplay between thermal plumes and ventilation flow from air conditioners can result in complex flow patterns that can transport the aerosol circumventing regular shielding and lead to a potential exposure of an individual behind a shield. Moreover, the returning aerosols from air conditioning/ventilation systems due to the limited filtration efficiency can also cause aerosol exposure of individuals adjacent to or facing the ventilation outlets. Our results show that these aerosols can contribute as high as 30% to the total infection risk in some ventilation and thermal settings. These results highlight the need for a proper shield design or placement of shielding according to the local flow patterns as well as the need for improving the filtration efficiency of our air conditioning/ventilation system. Overall, our research has demonstrated the capability and value of high-fidelity CFD tools for airborne infection risk assessment and the development of effective preventive measures, which can be further strengthened if proper infection-dose models and more detailed information of specific settings are available.

H.L. and S.H. contributed equally to this work.

J.H. would like to acknowledge the support of the University of Minnesota Rapid Response Grant from the Office for Vice President of Research (OVPR). The authors would also like to acknowledge the Minnesota Supercomputing Institute (MSI, http://www.msi.umn.edu) at the University of Minnesota for providing resources that contributed to the research results reported within this paper.

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

Two benchmark cases from the literature were simulated to validate our simulation algorithm. Case 1 is about forced convection, which was measured by Posner et al.⁶³ The setup of this case is shown in Fig. 17. The dimension of the room model is 0.914 × 0.305 × 0.457 m³. A square inlet and a square outlet vent are located on the ceiling of the room and the length of their sides is 0.1 m. The inlet velocity is 0.235 m/s. A partition of height 0.15 m is located at the center of the room. The grid number used in LES is 120 × 40 × 60. Besides the inlet and outlet boundary conditions at the vents, no-slip boundary condition is applied on the walls and the partition. We compared the measured and simulated vertical velocity along two lines, as shown in Fig. 17, where line 1 is parallel to the central axis of the inlet vent and line 2 goes through the center of the partition. The comparison is plotted in Fig. 18. As shown, the present LES result is very close to the experiment result.

We have also conducted convergence tests on our numerical algorithm. We performed LES for case 1 with the grid numbers of 80 × 30 × 45, 120 × 40 × 60, and 180 × 60 × 90. The comparison of the vertical velocity along line 1 is plotted in Fig. 21. As shown, our simulation algorithm obtains good convergence performance with different grid resolutions.

Case 2 is about mixed convection, which was measured by Blay et al.⁶⁴ The setup of this case is shown in Fig. 19. The dimension of the room model is 1.04 × 1.04 × 0.7 m³. An inlet vent and an outlet vent span along the z–direction, and their widths are 0.018 m and 0.024 m, respectively. The inlet velocity is 0.57 m/s. The grid number is 120 × 120 × 60. The boundary conditions for the flow velocity are the same as those in case 1. For the setup of heat transfer, the floor temperature is 35 °C, the temperature on all the walls and the ceiling is 15 °C, and the temperature of the inlet flow is also 15 °C. We compared our velocity and temperature results along line 1 (x = 0.52 m and z = 0.35 m) and line 2 (y = 0.52 m and z = 0.35 m) with the previous experiment.⁶⁴ As shown in Fig. 20, our simulation agrees with the measurement. Grid convergence test for validation case 2 has also been performed in a way similar to validation case 1 and is not repeated here.