Fast monotonically integrated large eddy simulation solver: validation of a new scalable tool to study and optimize indoor ventilation

Abstract

Indoor ventilation is underutilized for the control of exposure to infectious pathogens. Occupancy restrictions during the pandemic showed the acute need to control detailed airflow patterns, particularly in heavily occupied spaces, such as lecture halls or offices, and not just to focus on air changes. Displacement ventilation is increasingly considered a viable energy efficient approach. However, control of airflow patterns from displacement ventilation requires us to understand them first. The challenge in doing so is that, on the one hand, detailed numerical simulations – such as direct numerical simulations (DNSs) – enable the most accurate assessment of the flow, but they are computationally prohibitively costly, thus impractical. On the other hand, large eddy simulations (LES) use parametrizations instead of explicitly capturing small-scale flow processes critical to capturing the inhomogeneous mixing and fluid–boundary interactions. Moreover, their use for generalizable insights requires extensive validation against experiments or already validated gold-standard DNSs. In this study, we start to address this challenge by employing efficient monotonically integrated LES (MILES) to simulate airflows in large-scale geometries and benchmark against relevant gold-standard DNSs. We discuss the validity and limitations of MILES. Via its application to a lecture hall, we showcase its emerging potential as an assessment tool for indoor air mixing heterogeneity.

Impact Statement

Indoor displacement ventilation (IDV) is an energy-efficient, viable tool for reduction of exposure to infectious pathogens. In IDV, mechanically driven flow is supplemented with buoyancy-driven flow naturally generated by the occupants’ heat output. However, the detailed flow patterns in IDV are sensitive to heating, occupation patterns and room geometry. It is thus critical to have efficient tools to evaluate IDV performance with respect to airborne pathogen turbulent mixing. Yet, precise direct numerical simulations of such flows are computationally very expensive, thus impractical for complex geometries. Here, we demonstrate how a monotonically integrated large eddy simulation can be used as a cheaper alternative. With the example of a lecture hall, we showcase the potential to bridge the gap between prohibitively expensive precise flow simulations and overly simplistic models rooted in the well-mixed assumption known to inaccurately capture exposure risk in most indoor spaces that are, in fact, typically, heterogeneously mixed.

1. Introduction

Optimized indoor ventilation has to take into account multiple objectives such as air quality (Seppänen, Fisk & Mendell 1999; Fisk 2000; Erdmann, Steiner & Apte 2002; Awbi 2003; Sharpe et al. 2015; Izadyar & Miller 2022), thermal comfort (Fong et al. 2011; Frontczak & Wargocki 2011; Mora & Bean 2018; Ganesh et al. 2021) and energy consumption (Liddament & Orme 1998; Linden 1999; Pérez-Lombard, Ortiz & Pout 2008; Chenvidyakarn 2013). Since these factors influence each other, one should look at them in a joint fashion (Rackes & Waring 2014; Wei et al. 2015). Ventilation systems in educational buildings constitute a case that deserves special attention (Jia et al. 2021) due to the high population density in classrooms (Chithra & Nagendra 2012). There exist multiple reviews on the topic for further information on indoor ventilation (Liddament 2000; Chenari, Carrilho & da Silva 2016; Jia et al. 2021) including a very recent one focusing on the effect of ventilation on the spread of airborne diseases (Bhagat et al. 2024). It has become evident during the COVID-19 pandemic that inhomogeneous mixing and indoor airflow patterns also contribute to (Li et al. 2021), or mitigate, the spread of airborne pathogens (Qian et al. 2006; Li et al. 2007; Morawska et al. 2020; Morawska & Milton 2020; Bourouiba 2021a). During breathing (Arumuru et al. 2021; Yang et al. 2022), speaking (Abkarian et al. 2020), coughing and sneezing, infected individuals emit multiphase exhalation puff clouds of increasing momentum. These exhalation clouds contain pathogen-bearing droplets of mucosalivary fluid (Bourouiba, Dehandschoewercker & Bush 2014; Bourouiba 2016, 2020, 2021b; Scharfman et al. 2016). Following the comparatively high-momentum emission phase of the exhalation, these droplets are advected in the ambient airflow and may recirculate or be transported over long distances for extended duration, thus creating a risk of airborne transmission (Riley, Murphy & Riley 1978; Remington et al. 1985; Balachandar et al. 2020; World Health Organization [WHO] 2020; Bourouiba 2021b; National Center for Immunization and Respiratory Diseases, Division of Viral Diseases [NCIRD] 2021; Wang et al. 2021; Bahl et al. 2022).

To help mitigate the risk of indoor respiratory disease transmission, it is important to analyse, understand and optimize airflow and occupancy patterns, given a particular room's configuration and type of ventilation. Numerical simulations nowadays have a crucial role in this analysis. However, the temporal and spatial scales involved pose significant challenges. Respiratory flows are multiphase, containing a continuum of exhalation droplet sizes spanning submicron to millimetre diameter. These flows are emitted into rooms of size of $O(10\ {\rm m})$, in which exhalations might take $O(100\ {\rm ms})$ while full-room air changes might take $O(1\ {\rm h})$. These disparities in temporal and spatial scales make simulations of the underlying detailed fluid physics prohibitive, particularly for parametric studies requiring a large number of runs to conclude on optimal occupancy and ventilation settings for risk mitigation. To simulate the flow patterns in a room, we restrict ourselves here to one specific and typical ventilation type, namely displacement ventilation (Linden, Lane-Serff & Smeed 1990; Hunt & Linden 1999; Bhagat & Linden 2020). Such ventilation is increasingly used in modern buildings (Bjørn & Nielsen 2002; Villafruela et al. 2019) due to its energy efficiency. It is characterized by the injection of clean, cold air near or through the floor of the room. The occupants in the room produce buoyant thermal plumes with a typical heat flux of around 80–100 W per person from the body into the room (Craven & Settles 2006; Bhagat & Linden 2020). A large portion of this heat is transferred to the walls via radiative transfer (figure 17 in the Appendix). The remaining part produces an upward-directed convective flow around the occupants: the heat transfer drives the injected cold ventilation air and also the respiratory exhalations upwards toward the ceiling, where the air is removed from the room. Such a configuration leads to temperature and density stratifications, where a clean cool air zone in the lower part of the room is separated from a warmer, possibly contaminated, upper layer (Zhou et al. 2017; Bhagat et al. 2020). The thermal stratification layer height is an important indicator of the efficiency of the ventilation and has been theoretically derived, in a simple geometry, using the inlet flow rate and the total body heat output of the occupants (Linden 1999; Bhagat et al. 2020).

Yang et al. (2022) have recently evaluated mechanical displacement ventilation using direct numerical simulations (DNSs) in which they studied the influence of the mechanical flow rate on the stratification layer height in a small room. Their simulations show agreement with established theory (Linden 1999; Bhagat et al. 2020) in that higher flow rates increase the height of the lower stratification layer. However, if the ventilation rate becomes too high, the layer's height no longer increases but breaks and a short-circuit flow between the in- and outlets emerges (Partridge & Linden 2017).

Such DNSs resolve fine turbulence structures and dissipation (Shishkina & Wagner 2012; Yang et al. 2022) and accurately capture the solutions to the fluid flow equations without subgrid-scale modelling of the smallest scales. However, due to their extensive computational cost, they are not scalable to parametric studies of ventilation or occupancy scenarios of realistic lecture hall dimensions, involving large numbers of simulations. Large eddy simulations (LES) used for indoor airflow (Su, Chen & Chiang 2001; Kempe & Hantsch 2017) are more amenable to parametric studies and complex geometries due to their lower computation cost gained from the subgrid-scale parameterization of the small-scale turbulence (Jiang & Chen 2001; Jiang et al. 2003; Durrani, Cook & McGuirk 2015; van Hooff, Blocken & Tominaga 2017). However, wall-bounded LES require specific validation of the modelling of the layer close to the wall (Piomelli & Balaras 2002). Moreover, the accuracy of LES for indoor air flows heavily depends on the specifics of the flow boundary conditions and scenarios (Zhang et al. 2007). Therefore, when using LES, it is critical to compare the outcome with experiments, and/or compare LES against higher resolution DNS fully resolved flow, as a first validation step (Samuel, Samtaney & Verma 2022).

Motivated by the importance of analysing and optimizing airflow patterns in lecture halls to mitigate respiratory disease transmission, we propose an efficient numerical approach to simulating airflow and tracer concentrations by employing a monotonically integrated large eddy simulation (MILES). There are existing fast LES simulations (Khan et al. 2015; Kristóf & Papp 2018; Bauweraerts & Meyers 2019; Adekanye et al. 2022; Korhonen et al. 2024). The novel aspects of the approach we present here are the techniques used for the modelling of the source terms, introduced in § A.3 of the Appendix, which can also be applied to such conventional LES simulations.

Since we introduce a newly developed numerical scheme, we devote this paper to the first step in our validation of such scheme: comparing the MILES results against those of the DNS-based study by Yang et al. (2022). Yang et al. (2022) used a fixed surface temperature of the occupant, set at 27 $^\circ$C. This was based on prior simulations, estimations and experiments by Craven & Settles (2006) and Höppe (1993), measuring the temperature of the clothing. In fact, physiologically, it is the core temperature (deep organs) that the body regulates and maintains constant; the skin organ being a peripheral sensor and effector (Morrison & Nakamura 2019). Thus, skin temperature can vary greatly across the body's surface and in response to external conditions of heat and humidity (Hardy & DuBois 1937; Houdas & Ring 2013). Hence, here, in contrast, we do not fix the surface temperature, but instead we fix the heat flux from the occupant.

Instead of modelling the true multiphase nature of the exhalations and ambient droplet-laden flows, a simplifying – although strong – assumption is often made: that CO$_2$ flow and concentration correlate with droplet/pathogen flow and concentration (Seppänen et al. 1999; Rudnick & Milton 2003; Mahyuddin & Awbi 2010). It is important to note, however, that no direct experimental verification of this assumption is available to our knowledge. The dynamics of mixing and particle–turbulence interaction indeed remain open areas of research (Yarin & Hetsroni 1994; Shaw 2003; Toschi & Bodenschatz 2009; Balachandar & Eaton 2010; Brandt & Coletti 2022; Banko & Eaton 2023). Nevertheless, for the purpose of the goal of comparing numerical schemes of MILES and DNS, we do start with the computationally affordable modelling of a gas tracer, such as CO$_2$, tracking its flow and concentration.

In § 2, we introduce the flow set-up, governing equations and numerical schemes. In § 3, we compare the results of the new MILES code against the DNS results by Yang et al. (2022), both qualitatively and quantitatively, and discuss the commonalities and differences. In § 5, we show an example/illustration of application of MILES to a realistic classroom setting. We summarize our findings in § 7 and discuss the remaining steps of validation needed to make MILES a scalable and validated decision making tool.

2. Equations, numerics and flow set-up

We solve the advection diffusion equations and the Navier–Stokes equations within the Oberbeck–Boussinesq approximation (Boussinesq 1903), together with the incompressibility constraint ${\partial u_i}/{\partial x_i} = 0$ for incompressible flows. In dimensionless form this set of equations reads

(2.1a)$$\begin{gather} \frac{\partial u_i}{\partial x_i} = 0, \end{gather}$$

(2.1b)$$\begin{gather}\frac{\partial u_i}{\partial t} + u_j\frac{\partial u_i}{\partial x_j} = {-}\frac{\partial p}{\partial x_i} + \frac{1}{Re}\frac{\partial^2 u_i}{\partial x_j \partial x_j} - \frac{1}{Fr^2}\frac{\partial z}{\partial x_i}T,\quad z = - x_3, \end{gather}$$

(2.1c)$$\begin{gather}\frac{\partial T}{\partial t} + u_j\frac{\partial T}{\partial x_j} = \frac{1}{Re\, Pr}\frac{\partial^2 T}{\partial x_j \partial x_j} + P, \end{gather}$$

(2.1d)$$\begin{gather}\frac{\partial Y}{\partial t} + u_j\frac{\partial Y}{\partial x_j} = \frac{1}{Re\, Sc}\frac{\partial^2 Y}{\partial x_j \partial x_j} + S, \end{gather}$$

where $u_i$ denotes the three dimensionless velocity components, $T$ the dimensionless temperature and $Y$ the normalized passive scalar concentration. Also, $P$ and $S$ represent dimensionless source terms for heat and tracer, respectively (they are discussed in more detail in § A.3). The dimensionless control parameters are the Reynolds number $Re=U_{\infty } L_\infty /\nu$, the Prandtl number $Pr=\nu /\kappa$, the Schmidt number $Sc={\nu }/{D}$ and the Froude number $Fr={U_{\infty }}/\sqrt {gL_\infty \Delta T/T_0}$, where $L_\infty$ is a reference length, $\nu$ the kinematic viscosity, $\kappa$ the thermal diffusivity, $g$ the gravitational acceleration and $\Delta T = T_{clothing surface}-T_0 = 5$ K the temperature difference between the body and the reference temperature. As discussed in the prior section, the clothing surface temperature $T_{{clothing\ surface}} = 27\,^\circ {\rm C}$ is from Yang et al. (2022) and $T_0$ is chosen to be $295\ {\rm K} = 22\,^\circ {\rm C}$. To be consistent with the non-dimensionalization in Yang et al. (2022), we fix $Fr = 1$ and $L_\infty = 3$ m for the person in the box, resulting in a characteristic $U_{\infty } = \sqrt {gL_\infty \Delta T/T_0}$. We also fix $Pr=0.71$, which is standard for air, and $Sc = 1$ for diffusion of the passive tracer. Given the kinematic viscosity of air at 25 $^\circ$C of $\nu \approx 1.56\times 10^{-5}$ m$^2$ s$^{-1}$ and the diffusion coefficient of CO$_2$ in air $D \approx 1.6\times 10^{-5}$ m$^2$ s$^{-1}$, this is a reasonable approximation for the $Sc$ value of CO$_2$ as a passive tracer.

The approach used in this paper can be categorized as a MILES method (Boris et al. 1992; Grinstein, Margolin & Rider 2007) belonging to the larger class of implicit LES models (Margolin, Rider & Grinstein 2006). In contrast to a conventional LES method, there is no explicit subgrid-scale model. Instead, an implicit model originating from the numerical diffusion of the discretization scheme is used (Boris 1989; Boris et al. 1992; Margolin et al. 2006). To solve, we use a second-order upwind discretization for the convective terms, second-order central discretization for the viscous and diffusion terms and a first-order Euler forward time integration. The computational domain for the validation against DNS is the same as used by Yang et al. (2022) (shown in figure 15). For further, more in-depth, technical details on the method, numerical discretization and algorithm as well as details of the second-order upwinding, the geometry, boundary conditions and source terms, we encourage the reader to view Appendix, §§ A.1 to A.3.

3. Comparison of MILES with direct numerical simulation results

To validate the results of our MILES code, we compared them with DNS data by Yang et al. (2022), including temperature, velocity fields, CO$_2$ profiles and averaged values at steady state.

3.1 Qualitative comparisons

In figure 1, to compare the results qualitatively, we show averaged steady-state fields for temperature and the horizontal velocity component for three different flow rates $Q$. Note that, throughout the paper, $Q$ denotes a per-person flow rate while $\tilde {Q}$ is the total flow rate in and out of the simulation domain. Panels (a,b) show the data from DNS of Yang et al. (2022), whereas panels (c,d) show the MILES results. All plots display averaged data produced from 200 individual snapshots taken at time intervals of one characteristic time, which is equivalent to $t_{\infty } = {L_{\infty }}/{\sqrt {g L_\infty \Delta T / T_0}} \approx 4.25$ s. The overall averaging duration therefore is $200 t_{\infty }$, which is a bit more than 14 min in physical time.

Figure 1. Qualitative comparison of the steady-state fields for $Q = 0.04$, 0.09 m$^3$ s$^{-1}$; time-averaged DNS (a,b) and time-averaged MILES (c,d). (a–d left sub-panels) For each flow rate we show the temperature in a slice through the middle of the domain (i.e. around $x_2 = 1.5$ m). (a–d right sub-panels) For each flow rate we show the velocity in the horizontal direction in a slice through the middle of the domain (i.e. around $x_2 = 1.5$ m).

There is qualitative agreement between the DNS and the MILES results for different flow rates $Q$ per person, including the body-induced thermal plume and the thermal stratification layer. The fact that the DNS results are much better resolved in space is manifested by the very small features that are visible. However, the MILES code does capture most of the large-scale turbulent motion, as can be seen in the two instantaneous comparative snapshots shown in figure 2.

Figure 2. Instantaneous snapshots of DNS (a,b) and MILES (c,d) in statistically stationary steady state for a flow rate of $Q = 0.09$ m$^3$ s$^{-1}$ person$^{-1}$. (a,c) Show the temperature in the slice at $x_2 =1.5$ m and (b,d) show the horizontal velocity in the slice at $x_2 =1.5$ m.

3.2 Quantitative comparisons

For a quantitative comparison of MILES with DNS, we compare both average values taken in different regions of the median plane (see figure 3) as well as profiles for both temperature and CO$_2$ concentration (see figure 4). Essentially, we are comparing with figures 2 and 4 of Yang et al. (2022). All comparisons are performed for statistically averaged steady-state data as described above and all averages and profiles are taken from the middle slice around $x_2 = 1.5$ m.

Figure 3. Quantitative comparison between MILES and DNS. (a,c,e,g) Agreement for different temperature averages (blue, olive, green) and temperature at a single location (red) for different flow rates between the MILES code (hashed bar) and the DNS results (filled bar). (b,d,f,h) Agreement for different CO$_2$ averages (blue, olive, green) and CO$_2$ value at a single location (red) for different flow rates between the simple code (hashed bar) and the DNS results (filled bar). Note that datapoints for MILES and DNS (i.e. two bars) are displayed apart for a given identical flow rate for ease of display. Recall, however, that they are both created using the same flow rates $0.04$ and 0.09 m$^3$ s$^{-1}$.

Figure 4. Profiles of average temperature (a) and CO$_2$ concentration (b) in the middle slice around $x_2 = 1.5$ m over the height of the room for flow rates $Q = 0.04$, 0.09 m$^3$ s$^{-1}$ person$^{-1}$. The solid lines show the DNS and the dashed lines the MILES results. The average is computed at every height from the blue and yellow regions in figure 3. More details about the averaging regions are given in the Appendix in § A.5.

3.2.1 Steady-state-averaged values

Figure 3 shows averaged values for two different flow rates, $Q$, taken in different areas of the middle slice shown in the inset. The error bars denote the minimum and maximum values of the 200 snapshots taken into account for the temporal averaging. The details of the averaging regions are shown in the Appendix in figures 18 and 19. For both $Q = 0.04$ m$^3$ s$^{-1}$ person$^{-1}$ and $Q = 0.09$ m$^3$ s$^{-1}$ person$^{-1}$ very good agreement between DNS and MILES is seen. The fact that the CO$_2$ concentrations do not agree as well for $Q = 0.04$ m$^3$ s$^{-1}$ person$^{-1}$ might be related to the bigger differences in the way the breathing and the heat transfer are implemented in DNS and MILES.

3.2.2 Steady-state profiles

An important outcome of the simulations is the so-called cleaning height that denotes the height of separation of the lower, cool and presumably cleaner air versus the upper, warmer and presumably contaminated air. At that transition, both temperature and CO$_2$ concentration profiles show steep gradients. To investigate how well the MILES code captures the cleaning height, we compare steady-state profiles for both temperature and concentration. Figure 4 shows profiles for $Q = 0.04$ and 0.09 m$^3$ s$^{-1}$ person$^{-1}$ for the DNS reference (solid) and MILES (dashed). As already seen from the average values, the agreement for the temperature profile is good for both $Q = 0.09$ m$^3$ s$^{-1}$ person$^{-1}$ and $Q = 0.04$ m$^3$ s$^{-1}$ person$^{-1}$, while the agreement for the CO$_2$ concentration is again worse than that for temperature for both flow rates. To get a good agreement for the temperature profile for $Q = 0.04$ m$^3$ s$^{-1}$ person$^{-1}$, we used a higher grid resolution of $300^3$ in contrast to $150^3$ for $Q = 0.09$ m$^3$ s$^{-1}$ person$^{-1}$. More details about mesh convergence can be found in the Appendix in § A.7. In § A.10 of the Appendix, we also included profiles for a low flow rate of $Q = 0.01$ m$^3$ s$^{-1}$ person$^{-1}$. The agreement between MILES and DNS for this low flow rate, however, is worse, even when using the higher grid resolution of $300^3$ for the MILES code. The reason that we did not include this case in this discussion is the fact that the discrepancy could be related to the ground truth DNS profile. As we show in § A.9 of the Appendix, the time required to reach the steady state for low flow rates can be up to multiple hours. Therefore, it is possible that the DNS profile in Yang et al. (2022) might not have been obtained in fully steady state.

3.3 Multiple-occupant cases: four occupants

As an extension of the one-occupant case and as an intermediate step towards full occupancy cases, such as lecture halls, we compared MILES with DNS data with multiple occupants. The set-up remains identical to the one-occupant case, except that only one total flow rate of $\tilde {Q} = 0.08$ m$^3$ s$^{-1}$ is considered, i.e. four $Q = 0.02$ m$^3$ s$^{-1}$ person$^{-1}$. Four occupants are placed within the box at the corners of a square with side lengths $d = 1$ and $0.5$ m.

3.3.1 Qualitative agreement

Figure 5 shows averaged steady-state data obtained with DNS (a,b) and with MILES (c,d). Temperature and vertical velocity are shown for both separation distances. The agreement is good in all cases. When increasing the number of occupants to four, the forcing terms in the equations are increased. In such case, we observe that phenomena such as the layer height are less dominated by local equilibria and the match between MILES and DNS is improved despite the low flow rate of $Q = 0.02$ m$^3$ s$^{-1}$ per occupant. This may indicate that the discrepancy between MILES and DNS is smaller with multiple occupants in the room. Thus, the MILES method might be fruitfully applied to large lecture halls with tens or hundreds of people.

Figure 5. Comparisons of simulations with multiple occupants. (a,b) Show DNS results for the separation distances $d = 0.5$ m and $d = 1$ m. (c,d) Show the corresponding MILES results. For both distances (a–d left sub-panels) show the average steady-state temperature and (a–d right sub-panels) the vertical velocity component in slices through the occupants at around $x_2 = 1.25$ m and at $x_2 = 1.0$ m for $d = 0.5$ m and $d = 1.0$ m, respectively.

3.3.2 Quantitative agreement

A quantitative comparison for the four-occupant case is now discussed. Figure 6 reports the vertical velocity component on a slice through two of the four people at a height of $x_3 = 1.8$ m, i.e. above the heads. Comparison between DNS and MILES results show good agreement regarding peak location and vertical velocity magnitude, for both $d=0.5$ m and $d=1.0$ m.

Figure 6. Vertical velocity profile above the occupant's heads at height of $h = 1.8$ m for the two cases with four occupants.

4. The CPU time comparison between DNS and MILES

In § 3 we compared DNS and MILES in terms of reliability of the results. We find good agreement for medium and high flow rates with one occupant and with multiple occupants. The MILES method also has a big advantage which is its time efficiency. Table 1 shows a comparison between the runtimes and costs of DNS and MILES for a case with one occupant in the room.

Table 1. Comparison of runtime and cost for DNS and MILES. According to Walker (2009), we have assumed a cost of 0.1 USD per CPU hour. The DNS code was run on a cluster using 384 cores. The MILES code was run on the Sherlock cluster at Stanford University using a single NVIDIA A100 SXM GPU. We assumed a price of 2.5 USD per hour for such a GPU (Jolt ML Ltd 2024). The runtimes show some variability and are taken from the slowest simulations available.

The reason for the fast runtime of the MILES code is the fact that all operations on a high level are implemented as element-wise tensor additions, subtractions, multiplications, divisions or $\min \max$ operations using tensors in PyTorch (Paszke et al. 2019) or use the high-performance stencil computations in ParallelStencil.jl (Omlin & Räss 2022). For approximations where neighbouring elements are relevant, such as gradients, the tensors are shifted against each other using tensor slicing, which results in the same element-wise structure. The computations are performed and the tensors are stored on GPUs which yields a drastic improvement in runtime compared with using arrays on CPUs. The aim of the current study, however, is not at all for it to be a competitive comparison between the DNS and MILES codes. The implementation can be done in many ways, here, we used PyTorch (Paszke et al. 2019) for an early version of the code and Julia (Bezanson et al. 2017) together with ParallelStencil.jl (Omlin & Räss 2022) for the current version.

5. Application to lecture hall scenarios

The designated use case for the newly developed MILES code, especially due to its computational efficiency, is the investigation of flow, temperature and CO$_2$ and other tracer distributions throughout rooms on the scale of an entire lecture hall. Here, we consider one specific lecture hall at ETH Zürich (Lecture Hall ML E 12). Typical of larger lecture halls, the room has a slanted (amphitheater) layout with displacement ventilation. We expect the simulation results for this lecture hall to carry over qualitatively to other lecture halls which share its main characteristics. For different room types (e.g. seminar rooms, meeting rooms, offices or laboratory spaces), one would need to perform similar studies with the applicable room geometries, seating layout and ventilation types.

5.0.2.1 Geometry and boundary conditions

Figure 7 shows a three-dimensional visualization of the geometry used for the simulations. The uniform resolution for the lecture hall case was chosen to be $\Delta x_1 = \Delta x_2 = \Delta x_3 = 0.05$ m. All solid objects including people are modelled using staircase geometries. The inlets are located below the seats, whereas the outflow is taken to be uniform through the ceiling except for where the lamps are. The combined inlets and outlets honour a prescribed volumetric flow rate of $\tilde {Q} = 4500$ m$^3$ h$^{-1}$ which is equivalent to approximately $Q = 0.02$ m$^3$ s$^{-1}$ person$^{-1}$. Uniform velocity is enforced at both inlet and outlet surfaces. We fix $Fr = 1$, $T_0 = 295$ K, $\Delta T = 5$ K and $L_\infty = (L_x+L_y+L_z)/3$, resulting in a characteristic $U_{\infty } = \sqrt {gL_\infty \Delta T/T_0}$. We also fix $Pr=0.71$ and $Sc = 1$. The walls are assumed to be adiabatic and impenetrable with a no-slip condition for the velocity. Occupants are modelled similarly as in the previous validation cases, featuring an impenetrable adiabatic core surrounded by a layer of cells, where in total a heat input of $P_{body}$ per occupant is applied. Two additional heat sources are modelled for this lecture hall: the projector in the back with a heat output of $P_{projector} = 1200$ W and a media rack with an estimated heat output of $P_{mediastation} = 500$ W. The simulated occupants’ mouths are modelled with a small number of cells representing a constant (in time) scalar CO$_2$ source term equivalent to a physical source term of 7.5 l min$^{-1}$ with a CO$_2$ concentration of 4 %. In addition to the CO$_2$ exhaled by all simulated occupants in the lecture hall, six additional tracer fields were included in our analysis to distinguish exhaled air coming from different sources. The mouths are, again, shown in red in figure 7.

Figure 7. Set-up for the ‘ML E 12’ lecture hall featuring a typical displacement ventilation set-up. The inlets are located below the seats (green) and the outlet is assumed to be uniformly distributed over the ceiling (transparent red). The golden parts represent the regions around the occupants, where the heat source term is introduced. The mouths (red) denote the region where the passive tracer (CO$_2$) is introduced. All walls, tables and other solid objects are adiabatic and impenetrable (grey). The projector (violet) in the back and the media rack (yellow) are additional heat sources. The mouths of the occupant spreading an additional tracer in addition to CO$_2$ are highlighted in blue. Note that the $x_3 = 0$ plane does not coincide with the floor in the front part of the lecture hall. This is because we included an additional staircase in the far right corner.

5.0.2.2 Base case scenario

The investigated scenario is a case with 50 % occupancy where people are sitting uniformly distributed in the lecture hall with empty seats to the left and right of each of the simulated occupants. For the base case configuration the heat power per occupant was set to $P_{body} = 42$ W distributed in layers of one cell thickness around each occupant.

5.1 Results

Snapshots and averaged steady-state temperature and CO$_2$ concentration maps are shown in figure 8. The snapshots show only the turbulent fluctuations and demonstrate that not only an averaged steady-state field is computed but also the transient evolution of the quantities. The slices on the bottom of figure 8 show the time-averaged temperature and CO$_2$ concentration. In the front part of the lecture hall one can see the characteristic separation layer seen in the one- and four-occupant simulations. Also, we observe that the CO$_2$ concentration is highest in the top rear part of the lecture hall.

Figure 8. The MILES results for the temperature (a) and for the CO$_2$ concentration (b) for the base case scenario of the lecture hall. The top row shows instantaneous snapshots and the bottom row time-averaged steady-state fields.

5.2 Spreading pattern

Figure 9 shows time-averaged steady-state exhaled air concentration maps for this scenario in the approximate breathing zone. The red dots denote the locations of selected emitters. The exhaled air from emitters sitting in the upper parts of the lecture hall (top end of the pictures) has a smaller domain of influence which is mainly limited to the upper region. The two emitters located on the side mainly spread towards the corners, whereas the exhaled air by the spreader in the middle is sucked up by the heat-induced airflow of the projector. The exhaled air from emitters sitting in the front part of the lecture hall, however, can spread much wider. Again, for the two spreaders at the sides, the exhaled air is mainly transported to the back and towards the corners of the lecture hall. The air exhaled in the middle of the front row is transported to both sides of the lecture hall with a preference for the right side. Streamlines for two of the locations are shown in figure 13.

Figure 9. Lecture hall base case: spread of exhaled air for different emitter locations in the lecture hall. The red dot in each panel depicts the location of the emitter. Note: the depicted slice is not perpendicular to the $x_3$ axis but slanted because of the slanted room geometry.

To gain more confidence in these findings, the robustness of these results was examined by varying different simulation input parameters and by comparing the resulting spreading maps to the baseline case shown in figure 9.

5.3 Robustness of spreading patterns

The main modelling assumptions required for these simulations is the introduction of the source terms. We therefore study the robustness of the tracer spreading with respect to the heat input per occupant as well as the heat source-term distribution.

5.3.1 Heat input – qualitative effect

The heat transfer by the occupants is a major driver of the flow in such lecture halls. The heat output of one occupant in the box of the DNS was approximately $12$ to $15$ W depending on the flow rate considered. However, typical values used for realistic studies are $\approx$80–130 W (Spitzer, Hettinger & Kaminsky 1982) per occupant of which approximately 50 W are transferred via convection (Höppe 1993). Here, we assumed that loss due to radiation gets absorbed by the walls where the added heat contributes very little to natural convection, thus it is ignored. The default spreading pattern in figure 9 is based on a heat flux of 42 W per occupant. Figure 10 shows the same spreading pattern for a heat flux of 15 W per person.

Figure 10. Lecture hall with different heat input: spread of exhaled air for different emitter locations in the lecture hall. Compared with the base case in figure 8, and figure 9 the heat input per occupant here is 15 W. The red dot in each panel depicts the location of the emitter. Note: the depicted slice is not perpendicular to the $x_3$ axis but slanted because of the slanted room geometry.

For the rear part of the lecture hall, the spreading patterns in the top row of figure 10 look comparable to the ones for 42 W. The spreading patterns for the three locations in the front (bottom row of figure), however, are different compared with the baseline case. Especially further away from the source the magnitude is much reduced for the front left and front centre spreader. The spreader in the front right shows a unique feature of very high concentration towards the upper right corner similar to the pattern observed for the same spreader location in figure 9.

5.3.2 Heat input – quantitative effect

As the heat input per occupant is increased from 15 to 70 W, the room temperature increases (figure 11a,c,e) as well as the temperature around each simulated occupant (figure 11a,e). Additionally, the temperature gradient from the front to the rear of the lecture hall also increases (figure 11e). The large temperature values towards $x_1=0$ are a result of the high temperatures established just below the ceiling in the back of the lecture hall where the projector provides a strong heat source (figure 11e). At the back of the room (figure 11c), individual simulated occupants can no longer be discerned from the spatially averaged steady-state temperature profiles, suggesting that the rear row analysed is above the mixing height of the room.

Figure 11. Average steady-state temperature (a,c,e) and C$O_2$ concentration (b,d,f) for different heat input values (15,42,70) W. All profiles are temporally averaged steady-state data spatially averaged over different regions: (a,b) averaged over slices $x_1\times x_3$, [7.5 m, 8 m] $\times$ [2.8 m, 3.2 m] for each value along $x_2$; (c,d) averaged over slices $x_1\times x_3$, [2.7 m, 3.2 m] $\times$ [4.45 m, 4.95 m] for each value along $x_2$; (e,f) averaged over slices along $x_2 \times x_3$ [2 m, 2.5 m] $\times$ (0.5 m region above head) for each value along $x_1$. More details about the region above the head can be found in § A.5 in the Appendix.

The concentration profiles (figure 11b,d,f) show less clear a trend than the corresponding temperature profiles. However, the simulation results suggest that the peak tracer concentrations are highest for the lowest heat input per person (15 W) explored. The lower heat input per person results in a lower plume velocity and hence higher residual concentrations, although the relative effect seems small.

The layering height in the lecture hall is estimated by revisiting the (1.1) in Yang et al. (2022) derived from Schmidt (1941), Morton, Taylor & Turner (1956) and Linden et al. (1990)

(5.1)\begin{equation} \tilde{Q} = cn^{2/3}B^{1/3}(h-h_v)^{5/3}, \end{equation}

where $n$ is the number of occupants, $B$ is the buoyancy flux per occupant, $h_v$ is the virtual origin and $c$ is an empirical constant ($c = 0.105$ according to Morton et al. 1956; Linden 1999). Rearranging (5.1), we obtain the dependency of the cleaning height $h$ as a function of the number of occupants $n$ and the buoyancy flux per occupant $B$, as

(5.2)\begin{equation} h = \left(\frac{\tilde{Q}}{cn^{2/3}B^{1/3}}\right)^{3/5} + h_v. \end{equation}

We compare our simulated layering heights with the theoretically predicted values by (5.2). Figure 12 shows average steady-state profiles of temperature (a) and concentration (b) in the front part of the lecture hall.

Figure 12. Simulated average steady-state profile and predicted cleaning height of temperature (a) and CO$_2$ concentration (b). The simulated profiles (solid lines) are obtained from a slice at around $x_2 = 7.4$ m and between approximately $10.75\ m\leq x_1 \leq 13.27\ m$. The dot in the left panel denotes the cleaning height defined by the steepest temperature gradient. The dashed lines represent the predicted cleaning heights from (5.2). Spatially averaged profiles of temperature in the front part of the lecture hall for two different resolutions (c). The profiles are obtained from a slice around $x_2 = 7.4$ m and for $x_1$ between approximately 10.75 and 13.27 m.

In our simulations, we define the layering height as the point of steepest gradient in temperature, as depicted as a dot in the temperature profile in figure 12. The theoretically predicted layering height is shown by the dashed horizontal lines. The latter were obtained from (5.2) with $n= 63, \tilde {Q} ={\frac {4500}{3600}}$ m$^3$ s$^{-1}$, $B = \{15+\frac {1200+500}{63},42+\frac {1200+500}{63},70+\frac {1200+500}{63}\}$ W. The constant contribution of the projector and mediastation, 1200 and 500 W, respectively, were added proportionally to $B$, thus increasing the heat flux per occupant by $(1200+500)/63$ W. The virtual origin height $h_v$ was chosen to be $h_v = -1.214$ m to match the cleaning height for a heat input of 42 W per occupant. The same value was then used for the cases of $15$ and 70 W.

5.4 Slightly perturbed seating scenarios

To check the robustness of the spreading patterns with respect to slightly perturbed seating configurations, we performed two additional simulations with 42 W per person. We compared the base case seating configuration shown in figure 32(a) with the two scenarios in figure 32(b,c). From the results in figure 33, it becomes clear that these slight perturbations in seating configuration and number of occupants do not have a strong effect on the resulting magnitude of the average profiles of temperature and CO$_2$ concentration. However, when moving all occupants positions by one seat as in scenario 1 in figure 32(b), we can clearly see the local changes of areas of high temperature and CO$_2$ concentration.

5.5 Projector

We also compared results of scenarios with and without the projector. The results are presented in figure 34 in § A.13. Surprisingly, despite the fact that the projector is by far the strongest heat source in the lecture hall, the profiles of temperature and CO$_2$ are even less affected than due to the slight perturbations discussed above. The largest difference is visible at the very back of the lecture hall (figure 34e), where the projector provides a strong heat source. The profiles in figure 34, however, do not cover much of the uppermost part of the lecture hall, where we would expect a stronger effect.

5.6 Resolution

We also show that, with respect to the mean temperature profile, the resolution of $h = \Delta x = \Delta y = \Delta z = 0.05$ m is sufficient for this lecture hall simulation case. Figure 12(c) shows the spatially averaged temperature profiles in the front part of the lecture hall for resolutions of both $h = 0.0$ and 0.025 m. Both lines are averaged from 114 snapshots between 30 and 45 min of physical time.

5.7 Streamlines and flow structures

For a better understanding of how the spreading patterns above emerge, we show streamlines of the time-averaged steady-state velocity field, initialized at two of the 6 spreader locations in figure 9. We would like to point out that these are streamlines and not trajectories obtained by particle tracking. However, using frequent snapshots of the MILES solution for particle tracking would be an ideal extension and usage of the method. The streamlines in figure 13(a) show the spread of the exhaled air from the middle position of the front row. They are first going up but do not go all the way to the ceiling, indicating that the buoyancy force is not strong enough to push all the way through the warm upper layer. They are then mostly going towards the front of the lecture hall. These streamlines going to the front are then deflected to both sides of the lecture hall. The ones deflected towards the mediastation are sucked upwards by the flow due to the heat released by the mediastation in the front left corner. The other streamlines deflected to the other side return from the front wall and propagate towards the seating area in the front right corner, leading to an elevated concentration, as shown in bottom centre of figure 9. Figure 13(b) shows the spread of exhaled air from a second front row release position on the far right (bottom right corner in figure 9). Different from the release position in figure 13(a), the streamlines mostly occupy only the right side of the lecture hall, which explains to first order the higher concentration in exhaled air in the right part of the lecture hall in figure 9.

Figure 13. (a) Streamlines based on the averaged steady-state flow field. The streamlines are initialized at the location of the spreader in the middle of the front row shown in figure 7. (b) Streamlines based on the averaged steady-state flow field. The streamlines are initialized at the location of the spreader on the right end of the front row shown in figure 7.

These streamlines are based on the time-averaged velocity field, it is therefore not clear how stable the streamline path is over time.

6. Limitations

We provide a summary of the limitations of the present study regarding the physical assumptions and numerical discretizations used.

First, the physical limitations are listed. These include the fact that no active exhalation is taken into account; the multiphase nature of the exhalation is also not modelled with exhalations treated as passive tracers; the over-temperature and the momentum of the exhaled air are neglected. Compared with heat/temperature, the CO$_2$ exhalation is simplified more drastically in the MILES model, which could explain the larger differences in the agreement between DNS and MILES for the CO$_2$ concentration values and profiles observed in figures 3 and 4. Also, the radiation is not modelled explicitly but simply accounted for via an empirical reduction in the convective heat loss term (Hardy & DuBois 1937; Höppe 1993). The fact that in reality the exhaled air is multiphase in nature, at core body temperature, which is significantly higher than the clothing temperature or temperature of the surrounding, and carries forward momentum, all significantly influence the tracer and pathogen distribution (Bourouiba et al. 2014). All this can be modelled with an Euler–Lagrangian approach as done in Chong et al. (2021) and Ng et al. (2021) but this is computationally much more costly, so that at most seconds of a single respiratory event can be modelled.

Second, we list the numerical scheme limitations. The simulations we show have an implicit turbulence model tightly coupled to the spatial discretization and therefore grid resolution. We had to use a higher grid resolution of $300^3$ for the flow rate of $\tilde {Q}=0.04$ m$^3$ s$^{-1}$ compared with $150^3$ for $\tilde {Q}=0.09$ m$^3$ s$^{-1}$ to get better agreement with the DNS data. Moreover, for some simulations, we observe spurious spatial oscillations towards the ceiling of the lecture hall. From figure 23 in the Appendix, we have some indication that these oscillations become smaller with increasing grid resolution.

7. Conclusion and outlook

Motivated by the COVID-19 pandemic and the goal to help mitigate the risk of respiratory infectious disease transmission, a number of numerical studies on indoor ventilation have emerged recently (Vuorinen et al. 2020; Liu et al. 2021; Talaat et al. 2021; Oksanen et al. 2022; Wang & Hong 2023; Korhonen et al. 2024). Detailed LES simulation results (Liu et al. 2021) for the Guangzhou restaurant scenario (Bostock 2020; Chang 2020) have shown a connection between regions of high aerosol exposure due to airflow inhomogeneity and reported infections, as well as the usefulness of computational fluid dynamics as an instrument for the estimation of airborne infection risk. Reynolds-averaged Navier–Stokes (RANS) results have been presented for a very similar case as the validation case used here with a single mannequin in a box (Wang & Hong 2023). The saturation of the interface height as a function of the ventilation rate agree very well with the results from Yang et al. (2022). However, both Liu et al. (2021) and Wang & Hong (2023) use a fixed temperature boundary condition for the people in the simulations, which does not allow for a control of the heat flux compared with the fixed heat flux condition used in this work. Different Dirichlet temperature values were used in Auvinen et al. (2022) for the head and the clothing estimated from infrared images. The highly resolved LES simulations of a big room ($V=170$ m$^3$, 286 million cells) are combined with aerosol concentration measurements to successfully validate the LES model and study two different strategies to reduce the transmission risk (Auvinen et al. 2022). However, to simulate 1 h of physical time, the simulations take 6 days on more than 800 CPUs. Our method, using a coarser resolution as compared with Auvinen et al. (2022), can simulate a 1 hour lecture hall scenario of the same order of magnitude in size ($V=580$ m$^3$, 7.1 million cells) on a single GPU in less than 24 h using the latest GPU hardware. There are other existing fast LES simulations (Khan et al. 2015; Kristóf & Papp 2018; Bauweraerts & Meyers 2019; Adekanye et al. 2022). Most of them rely either on a conventional LES approach using subgrid-scale models, make use of a different method such as the lattice Boltzmann method or require more involved computations. Also, none of these methods uses the source-term approach employed here to model the heat loss of the occupants. RANS has also been used for efficient simulation of the flow inside an aircraft cabin (Talaat et al. 2021). However, the authors do not mention temperature being taken into account, which would be important to capture buoyancy effects driven by the heat loss of the people. In Vuorinen et al. (2020), coughing is modelled and the momentum transfer by the cough is accounted for. Long-range transport $O$(10 m) of aerosols is shown. For the OpenFOAM simulations, an implicit LES approach was used as a subgrid-scale model (Vuorinen et al. 2020), similar to the MILES approach used here. Recently, a highly efficient GPU-based LES simulation tool that uses binary masks very similar to the ones used here has been reported (Korhonen et al. 2024). Despite the high efficiency, no cases with physical simulation duration greater than 30 min are presented. Here, we have presented a newly developed MILES code for ventilation simulation. This simple yet reasonably accurate numerical scheme, and especially its powerful implementation with the use of PyTorch (Paszke et al. 2019) (earlier version) or Julia (Bezanson et al. 2017; Omlin & Räss 2022) (current version) for GPUs, make the new code a potentially valuable tool for flow simulation indoors. The use of source terms for heat and scalar concentration in the form of binary masks was shown to be a possible way to avoid the high computational cost related to boundary layers with Dirichlet conditions for temperature and scalar concentration. The technique of the modelling via source terms is not restricted to a MILES method. It could be similarly used together with conventional LES approaches. We compared the MILES results with DNS data by Yang et al. (2022) for displacement ventilation in a room with single and multiple occupants. Good qualitative agreement can be observed for temperature, CO$_2$ and velocity profiles. Moreover, we illustrate the application of MILES to a lecture hall.

The sensitivity of tracer spreading to the choice of values for various modelling parameters, such as heat input per occupant and heat layer thickness around the occupants, was investigated. The effects on temperature and CO$_2$ concentration values resulting from an increased heat layer thickness from $h=0.05$ to 0.1 m for the lecture hall case are small, where the larger thickness tends to smooth out peaks. This can be explained by the fact that the larger heat release layer has a more spread out heat source term which is smaller in numerical magnitude per cell because the source term is applied in more cells. For the two different cases we found that the effect of the heat layer thickness on the spreading pattern in the lecture hall is negligible. A change in the heat input of between 15 and 70 W showed a direct effect of increased temperature values in the lecture hall.

Generally, a potentially infected host, i.e. the source, sitting in the back typically influences a small region, mostly bounded by the rear of the lecture hall. The exhaled air of a source sitting in the front is transported up and to the back, resulting in predominantly increased concentrations toward the back. If one were to extrapolate the flow and concentration of the tracer to the flow and concentration of pathogen-bearing exhalation droplets or droplet residues, the simulations would suggest that people sitting in the back of the lecture hall are exposed to higher concentrations of exhaled air and therefore experience a higher risk of infection, compared with people in the front. Thus, occupants in the back would have a higher risk of exposure compared with the well-mixed limit, and people seated at the front a lower risk of exposure, again compared with the well-mixed limit. These results, from MILES, are supported qualitatively by recent experiments by Rusch & Rösgen (2022), who used direct measurements of CO$_2$ with human subjects present in lecture halls of similar geometry and ventilation type as examined here. Their measurements also show increased scalar concentrations in the back of slanted auditoriums, further demonstrating concentration inhomogeneity, providing further evidence that the well-mixed limit is not a suitable approximation for risk assessment and mitigation of airborne contaminants and pathogens in key gathering locations such as auditoriums, movie theatres, etc.

Further, on the fluid dynamics front, a very interesting aspect and one of general importance is the amount of heat an occupant introduces to the flow. The values in the DNS resulting from isothermal surfaces at 27 $^\circ$C range from 12 to 15 W per occupant, which is much lower than expected based on experience (Craven & Settles 2006; Bhagat & Linden 2020). Our sensitivity studies show that the exact heat flux value might not have a huge effect on concentration maps obtained for the entire lecture hall. They also show that a change in heat input has a direct effect on temperature values and via buoyancy also on CO$_2$ concentrations. It would be of interest to obtain the DNS results using a prescribed heat flux condition. One could then study both the effect of the boundary condition itself but also further validate the MILES approach for these applications.

Finally, studies using our MILES method should be performed to explore a wider parameter space, e.g. the geometry of the lecture hall, the ventilation type, the movement of bodies (Mingotti et al. 2020) and different respiratory activities. One direct extension of the work could involve producing frequent MILES solution snapshots and using the resulting flow fields for a particle tracking study. Furthermore, additional studies are necessary to investigate to which degree CO$_2$ concentration and the concentration of exhaled droplet residues actually correlate, which remains a major unknown. In addition, validation of simulation results against measured CO$_2$ or droplet/particle concentrations is also necessary to ensure that the simulations correctly incorporate the essential physics to capture the flow and transport of exhalations. These factors are crucial for better understanding building ventilation systems in general and how they can be leveraged for improved infection risk mitigation.