A GPU-accelerated computational fluid dynamics solver for assessing shear-driven indoor airflow and virus transmission by scale-resolved simulations

We explore the applicability of MATLAB for 3D computational fluid dynamics (CFD) of shear-driven indoor airflows. A new scale-resolving, large-eddy simulation (LES) solver titled DNSLABIB is proposed for MATLAB utilizing graphics processing units (GPUs). The solver is first validated against another CFD software (OpenFOAM). Next, we demonstrate the solver performance in three isothermal indoor ventilation configurations and the results are discussed in the context of airborne transmission of COVID-19. Ventilation in these cases is studied at both low (0.1 m/s) and high (1 m/s) airflow rates corresponding to $Re=5000$ and $Re=50000$. An analysis of the indoor CO$_2$ concentration is carried out as the room is emptied from stale, high CO$_2$ content air. We estimate the air changes per hour (ACH) values for three different room geometries and show that the numerical estimates from 3D CFD simulations may differ by 80-150 % ($Re=50000$) and 75-140 % ($Re=5000$) from the theoretical ACH value based on the perfect mixing assumption. Additionally, the analysis of the CO$_2$ probability distributions (PDFs) indicates a relatively non-uniform distribution of fresh air indoors. Finally, utilizing a time-dependent Wells-Riley analysis, an example is provided on the growth of the cumulative infection risk being reduced rapidly after the ventilation is started. The average infection risk is shown to reduce by a factor of 2 for lower ventilation rates (ACH=3.4-6.3) and 10 for the higher ventilation rates (ACH=37-64). The results indicate a high potential for DNSLABIB in various future developments on airflow prediction.


Introduction
The COVID-19 pandemic has set an unprecedented demand for multidisciplinary research to comprehend the transmission mechanisms of the SARS-CoV-2 virus [1].At the onset of the pandemic, the virus was initially assumed to transmit predominantly via larger droplets and fomites present on surfaces [1].However, since the early 2020, particle size on their ability to remain airborne, the effect of relative humidity on particle shrinkage as well as particle transport over large distances in turbulent indoor airflow [15].Indeed, early scientific contributions on the airborne transmission of the virus were provided by the physics-based computational fluid dynamics (CFD) simulations, which addressed the airborne transport of small particles in different indoor settings.For instance, Ascione et al. conducted a comprehensive study on the effects of various HVAC retrofitting alternatives in a university faculty which included CFD simulations on the ventilation configurations [16].Zhang et al. performed CFD analysis of humidity and temperature distributions and ventilation performance in an indoor space using Ansys FLUENT [17].Abuhegazy et al. provided a CFD simulation of a classroom in FLUENT and a detailed discussion on how windows, glass barriers as well as aerosol size and source might effect the particle trajectories [18].Other CFD studies have considered the impact of ventilation on the distribution of aerosols from coughing using a commercial software [19] and particle trajectories in OpenFOAM [20] as well as utilizing far-UVC lightning as a virus inactivator [21].Furthermore, various elements impacting the spread of airborne particles, such as ventilation, air filters and masks, have been considered in assorted CFD publications [22][23][24][25][26][27][28].At present, the aerosol inhalation route is broadly acknowledged to be one of the key mechanisms, possibly the main mechanism, of SARS-CoV-2 transmission [1,2,29].
From the modeling perspective, Reynolds-averaged Navier-Stokes (RANS) modeling has been favored over scale-resolved simulations in indoor airflow simulations in the past despite the reduced accuracy and the evidence to its inability to capture essential turbulent phenomena in such flows [30,31].While large-eddy simulation (LES) and direct numerical simulation (DNS) can certainly mitigate these problems, these methods naturally imply more computational effort due to increased mesh sizes and level of complexity.In order to promote such scaleresolving approaches in realistic indoor airflow modeling, efficient computational approaches are therefore required.
In this context, the advances in GPU based computing may prove increasingly beneficial for these large simulations, since their architecture is well suited for performing parallel computations on large numerical systems [32] and their value specifically for CFD has also been established [33,34].In addition to incompressible Navier-Stokes solvers [35][36][37][38], successful GPU implementations have been produced for multiphase flows [39][40][41], direct numerical simulations [42][43][44] and reactive flows [45,46], for instance.Of the most recent work, GPU enabled CFD simulations based on the concept of artificial compressibility method has been demonstrated in PyFr [47][48][49], which has been since augmented with optimal Runge-Kutta schemes [50] and locally adaptive pseudo-time stepping [51] for increased performance.Additionally, in [52], a modified Chorin-Temam projection method was implemented with spectral methods and utilized in solving the Navier-Stokes equations in a periodic flow geometry on a CPU in 2D.During the pandemic, the code (DNSLAB) by Vuorinen et al. [52] was extended by the authors to 3D and rendered compatible with GPUs for periodic flows without walls.The multidisciplinary research consortium work by Vuorinen et al. in the spring of 2020 was among the first systematic CFD assessments of COVID-19 airborne transmission [22].These investigations, using several open-source CFD codes, implied that the DNSLAB runs on a GPU clearly outperformed OpenFOAM simulations on a supercomputer in terms of computational time for simple problem types, i.e. fully periodic flows.
While the computational capabilities of GPUs for CFD have certainly been recognized by many, the inherent power of these devices can be offset by the steeply increasing requirement for technical expertise as the efficient implementation of the system of equations generally involves a suitable application programming interface (API), such as Compute Unified Device Architecture (CUDA) [53].Therefore, a more simplified software environment for these GPU implementations, such as MATLAB, would be preferable for the majority of users.Indeed, MATLAB has been endorsed in many other fields of scientific computing, including neuroscience [54,55], modeling of electrical circuits and systems [56][57][58] and control and communication systems [59][60][61][62].However, similar adoption of the software in the CFD community has not been materialized and as a result, its increasing potential as an accessible computational tool may therefore be neglected.In our view, a streamlined LES simulation software with the capacity to solve very large systems rapidly in MATLAB is therefore warranted.
Hence, in an effort to bridge the research gap, we present a GPU compatible CFD solver for shear-driven airflow problems in simplified geometries.In our software, ease of use and performance are emphasized to allow scale-resolved turbulent flow simulations (similar to LES and DNS) in typical indoor environments involving ventilation airflow, for instance.The main objectives of the paper are as follows.First, to explore the possibility of performing incompressible 3D scaleresolved flow simulations on a GPU in MATLAB.Second, to implement and validate a simplified immersed boundary (IB) method in MATLAB in order to explore indoor settings with solid obstacles, walls, tables and furniture.Third, to employ the new solver to characterize indoor airflow for three different ventilation configurations and discuss the findings in the context of airborne transmission of COVID-19.
The paper is organized as follows: first, we introduce the underlying system of equations, the IB method, the concept of mask functions and sources/sinks as well as the definition of hard walls in this context.Next, we present 2 canonical reference cases and an experimental case for validating the code.Then, the performance of the newly developed DNSLABIB code is demonstrated in (1) the ventilation-induced emptying of a room of stale air, utilizing various ventilation setups, and (2) the emission of exhaled aerosols from respiratory activities such as speaking.Finally, we reiterate the main results and insights obtained in these simulations in light of the airborne infection risk.

Theoretical and numerical framework
In the present work the focus is on low-speed, isothermal gas flows which can be modeled using the incompressible Navier-Stokes equations.Such isothermal flows are relevant to e.g.window ventilation setups, where the indoor air temperature is close to the outdoor temperature and buoyancy effects are not the primary concern.This implies that heat sources, such as crowds of people or radiators, are not included in the computations.In practice, this type of a situation could emerge in the summer time in an apartment when large windows/doors are opened for creating cross-draught.Additionally, we study the transport of a passive scalar field representing the indoor CO 2 concentration.At the end of the paper, the airborne trajectories of a small number of Lagrangian particles are also studied assuming oneway coupling.The pressure-velocity coupling is based on the projection method [63], where the time integration is carried out using a 4th order explicit Runge-Kutta scheme.Additionally, the spatial derivatives in the momentum equation are discretized using 2nd order central differences in the skew-symmetric, energy conservative form (see e.g.[64]).In the projection step, the pressure equation is solved in the Fourier space using the highly efficient fast Fourier transform (fft) method in MATLAB.The fft method is considered to be a key enabler for large scale CFD simulation in MATLAB although it restricts the simulations to fully Cartesian, equispaced grids.The projection step is executed only once per time step in order to speed up the code.Based on our experience, this approximation has a negligible influence on the actual numerical solution.
The governing equations for the fluid read where  is the velocity,  is the pressure divided by the fluid density ,  is the kinematic viscosity,  is a body force,  is the passive concentration field and   is the diffusivity of the concentration.The binary mask function  is used to mark the fluid phase ( = 1) and the solid phase ( = 0) respectively.On the no-slip walls, velocities are simply multiplied by .
In Eqs. ( 1) and ( 2), two additional terms appear: (, , ) ∕  , and .The latter term is a simple body force which is needed for pressure driven flows.The former is a forcing term adjusting the velocity to a target value at the flow inlets such as windows and ventilation ducts etc.This approach is needed since the present solver is periodic in contrast to non-periodic cases where the Dirichlet/Neumann boundary conditions for  and  can be provided as usual.In brief, the term (, , ) ∕  is formulated in order to establish the desired velocity   within the relaxation time scale   .By setting (, , ) = 1 at the specific location, where the velocity must reach the value   and (, , ) = 0 otherwise, the coordinate dependent mask function (, , ) geometrically confines the momentum source to the targeted region of the geometry.A respective source term is also utilized in the scalar equation in order to investigate mixing.
A general mask function based on the hyperbolic tangent function is illustrated in Fig. 1 and reads (, , ) = 1 2 where   ,   and   are the volumetric center coordinates of the source region, while   ,   and   are the window/inlet/ventilation duct dimensions of this region.This function obtains values in the range [0,1] and  ,, defines the smoothness of the transition between the two endpoints.In the example of Fig. 1, a window is located at a wall (left) and an inflow of air with a constant velocity is forced with the mask function.The inflow velocity on a plane spanned by the blue lines smoothly decreases to zero at the windows edges (middle).The behavior of the mask function along the red line shows this continuous transition in more detail along the -axis (right).Here, the width of the transition region, denoted by , is highly dependent on the   parameter in Eq. ( 5).Furthermore, the FFT approach in the pressure equation entails a periodic solution, which requires special consideration when creating inflow/outflow boundary conditions.In general, these can be imposed as follows in our code.Inflow through an opening (e.g.window or ceiling vent) requires either (1) modeling another window on the opposite side of the room from which the flow can exit (Fig. 2 left), or (2) extending the fluid region around the room so that airflow can enter and leave the room to conserve mass (Fig. 2

right). While option
A is used in the cross-draught cases, option B is used in the ceiling ventilation case.We note that option B requires more computational resources since an extra flow passage needs to be modeled around the room.
Additionally, to model the subgrid scale effects and stabilize the flow, we utilize explicit filtering of the velocity and scalar fields at the end of each timestep.The purpose of the filtering is to dampen the small scale fluctuations in the flow field.Mathematically, we introduce a 6th order damping term using the following differential form where the filter coefficient is chosen so that the Nyquist frequency is zeroed in the Fourier space i.e.    =  6   6 .Similar explicit filtering approaches have been noted to be successful in solving convectiondominated problems in the past [65].In practice, the filtering is implemented using a second order central difference formula for the sixth derivative.In the discretized form, the differential filter above can be Fourier transformed to yield the transfer function where , ,  and  have the values 1, −6, 15 and −10.For the low wavenumbers, T → 1 while for the high wavenumbers T → 0.

Solution of the pressure equation in Fourier space
The pressure equation requires particular attention near the walls.In the conventional projection method, one obtains a  * from the Navier-Stokes equation without a pressure gradient and then corrects  * with the pressure gradient which is solved from the Poisson equation In the expression above, the mask function  is simply used to implement the Neumann boundary condition directly to the Laplacian operator ∇ ⋅ (, , )∇.The computational grid points and the mask function  are illustrated in Fig. 3.The Fourier transform of this equation cannot be directly determined.However, by adding and subtracting 1 from , one has ∇⋅(, , )∇ = ∇⋅(1+(, , )−1)∇.Finally,  the pressure equation can be recast into the following equation where the left hand side operator () now has a well-defined Fourier transform.We iterate equation ( 10)  times evaluating the right hand side of the equation by central differences from the previous available value of pressure (  ) and the non-solenoidal velocity  * acquired from the momentum equation.Generally, the equation converges very quickly and a hard-coded value  = 4 is used here.The velocity field is then corrected as  =  * − ∇, utilizing the converged solution for the pressure to yield the solenoidal field satisfying Eq. (1).We note that the solution of linear system in Eq. ( 10) could be acquired in various ways.Here, we use fast Fourier transform (fftn) in MATLAB and the standard operator relationship  ≡ − 2  −  2  −  2  .In DNSLABIB we solve Eq. ( 10) by evaluating the right hand side derivatives using finite difference method, Fourier transform both sides of the equation thereafter, and finally inverse Fourier transform in order to get pressure on iteration level .Based on our numerical tests, this approach works well herein.Fourier transform is chosen in MATLAB as a method for solving the pressure equation due to its highly efficient MATLAB implementation on GPUs.

Lagrangian particles
Proceeding to Lagrangian particles, the equations of motion (EoM) read where   is the particle Reynolds number,   is the particle settling time scale, describing the delay of the particle in adjusting to altered flow conditions, and   refers to the drag coefficient of a particle.The subscripts  and  indicate particle and fluid properties, respectively,  is the gravitational force,  is the particle radius and  is the (bulk) density.These equations describe the trajectory of a particle which experiences the effect of gravity as well as the drag imposed by the surrounding fluid.Since the particle relaxation time scale   is generally small for particles of interest here (micrometer scale), the particle equations of motion are solved using the implicit Euler method, enabling larger time steps.Here, a uniform particle size distribution is utilized and the particle equations of motion above are solved for each individual particle.Two key quantities from the EoM above are the particle terminal velocity  * =   as well as the particle sedimentation time from height ℎ expressed as   = ℎ∕ * .For ℎ = 1.5 m and  = 5/10/30 μm solid particles   is approximately 36/9/1 min respectively if the ambient airflow is non-existent.In practice this simple analysis (see e.g.[22]) indicates that particles with sizes up to 100 μm can traverse significant distances and thus they pose a risk of also being inhaled.We further address this aspect at the end of this paper.

Overview of the DNSLABIB code
As stated earlier, we implement a numerical code to solve the governing equations ( 2), ( 3) and (10) using the MATLAB language.Accordingly, no code compilation nor external dependencies are required and the supported platforms currently include Windows, Linux and macOS.Our open-source code is freely available and the structure of the program is illustrated in Fig. 4. The user initializes a case and controls the subsequent simulation principally via the "SetParameters.m","CreateFields.m","CreateGeometry.m","Create-SourceMasks.m" and "InitializeDrops.m"scripts.In "SetParameters.m",the user provides the necessary information regarding the simulation geometry, fluid properties and data outputting.The user also specifies the maximum Courant number, as the implementation utilizes dynamic time stepping (in "AdjustDt.m").In "CreateGeometry.m", the user specifies the obstacles in the flow domain by forming the (, , ) field with primitive shapes.Then, in "Create-SourceMasks.m", the user specifies the location of the sources/sinks.
The user also decides on whether the simulation is run on CPUs or GPUs, leading to calls to either "CreateCpuArrays.m"or "Cre-ateGpuArrays.m",where the necessary arrays are converted to either GPU or CPU compliant format.An essential observation during the present study was that both the GPU and CPU execution share the same base code utilized later.Hence, the user requires no specialized knowledge on GPU programming in MATLAB to work with DNSLABIB.Based on the format of the arrays, MATLAB is capable of automatically detecting whether the subsequent operations and function calls involving these arrays are performed either on the GPU or CPU.Finally, the number, size and location of Lagrangian particles is assigned in "InitializeDrops.m".
After the computational setup and physical constants are defined as explained above, the simulation may be launched by running "NS3dLab.m".This main simulation loop calls various other functions designed to solve the fluid and particle equations which we introduced in the previous sections."SolveNavierStokes.m"("SolveScalar.m")implements the explicit RK4 time integration, updating the convective and diffusive terms in the Navier-Stokes (passive scalar transport) equations via calls to "ContructVeloc-ityIncrement.m" ("ContructScalarIncrement.m").Finally, the projection step is performed in "Project.m".In an analogous manner, "SolveDrops.m"advances the particle trajectories and velocities in time by constructing the terms in the equations of motion in "AdvanceDrops.m".Based on the settings provided by the user in "SetParameters.m"the main loop issues calls intermittently to "writeHDF5.m" to output simulation data.Since for large systems containing tens or hundreds of millions of cells and a vast number of time steps, the fluid data may routinely reach extreme sizes, and therefore, special consideration must be given to the data format.Currently, the code can be set to output the fluid data in HDF5 compliant format along with the associated XMDF metadata file or in the MATLAB native format (.mat).Additionally, the Lagrangian particle data is outputted in raw text data.In our experience, visualizing large fluid data sets in HDF5 format via external software, such as ParaView, yields superior read and rendering performance to many other options available.For standard users with more modest data outputting requirements, the flow quantities, such as the velocity field, pressure, the passive concentration (and their respective time averages) can also be forwarded to a .matfile for quick access and analysis.

Advantages and limitations of DNSLABIB
As noted, the main objective of the paper is to develop an efficient GPU-compatible code in MATLAB for performing scale-resolved 3D CFD simulations.Before proceeding to results, the expected advantages and limitations of DNSLABIB should be mentioned.The expected advantages of DNSLABIB include: • Potentially high performance, specifically for large systems, due to the GPU implementation.The primary contributors to the observed performance include avoiding loops (vectorization) as well as the efficient gpuArray structure and fft function in MATLAB • Ease of use in configuring and executing a case for simple geometries with solid objects due to MATLAB interface.• Scale-resolved simulations present the opportunity to capture physical flow features such as shear-driven turbulent flow, mixing, and flow recirculation zones.
In contrast, the main limitations of the present DNSLABIB implementation include: • Local mesh refinement is presently not possible due to the incorporation of fft for the pressure equation.The mesh resolution needs to be uniform.
• Since mesh refinement is not possible, only vents and windows with a relatively large diameter can be resolved (e.g. over 20 computational cells per diameter) at the moment.• The obstacles are presented as simple blocks.
• Considering the scope of the present paper, thermal sources/sinks are not currently accounted for although they are known to be of high importance in indoor airflow configurations.In contrast to cross-draught situations with high ventilation rates, the absence of thermal sources can be detrimental specifically in modeling interiors with lower ventilation rates as thermally-driven buoyancy flows may be more pronounced.• At present, only a single GPU can be efficiently utilized.
• Wall models are currently not implemented in the code and one needs to rely on constant grid resolution even at the walls.
Regarding the last limitation, the main focus within this study is in understanding the GPU compatibility of the present immersed boundary approach in MATLAB.Here, wall functions are not in the focus but, instead, we aim at resolving the shear-driven flows as well as possible on uniform grids.Commonly, the absence of wall models is considered detrimental to the solution accuracy in wall-turbulence driven, high Reynolds number flows, especially on coarse computational grids.Here, we carry out a sensitivity assessment on Reynolds number effects and discuss the cases  = 5000 (window airflow 0.1 m/s) and  = 50 000 (window airflow 1 m/s) to better understand how  affects the ventilation rate when the airflow velocity changes.In the present cases, near-wall air velocities are rather low speed on the order of ∼ 0.01−0.1 m∕s so that the length scales of the viscous wall layer scales + < 5 − 10 are mostly reached.For  = 50 000, 88% of the computational cells have + < 10.For  = 5000 all near-wall cells have + < 10 while ≈98% of the near-wall cells have + < 5.

General considerations
Next, the performance of DNSLABIB is explored.As background validation, the code was first validated in two canonical reference cases: the laminar channel flow and vortex shedding due to a rectangular body mounted into the channel (see Appendix A.1). Additionally, an experimental case involving a model building in a wind tunnel is simulated to test the solver in extremely high Reynolds number case in order to highlight certain advantages and drawbacks of the modeling assumptions in DNSLABIB (see Appendix A.2).
Following the basic validations, we further compare DNSLABIB against the open-source finite volume code OpenFOAM in an indoor ventilation setup at a high Reynolds number.In DNSLABIB and Open-FOAM, the grids are identical at the high shear regions of the flow, i.e. window/ventilation air jet regions, as well as at the walls.However, for computational feasibility, the OpenFOAM simulations required the use of coarser grid resolution outside these regions.A numerical comparison of the two codes is provided in the next section.
After the validation of DNSLABIB, we investigate DNSLABIB in a number of interior simulation cases to study the removal of stale air from the room and assess the infection risk in these configurations in the context of COVID-19.

3D turbulent flow indoors
In the following, airflow is simulated in a more challenging validation case at  = 5000 involving a large indoor space (8 m × 8 m × 3 m along -, -and -axis, respectively) with open windows.The space is divided into two larger rooms and a corridor (3 rooms in total) which are connected by doorways.A cross-draught is then generated as air enters and exits through the windows, ventilating the room.In addition to the  = 5000 with window airflow peak velocity   = 0.1 m∕s, we also investigate a higher inflow velocity of   = 1 m/s corresponding to  = 50 000.The main validation is based on the  = 5000 case which is also rather challenging in terms of fluid dynamics while the  = 50 000 case is investigated to better understand the  sensitivity of DNSLABIB and to show that the code stays numerically stable at extreme Reynolds numbers as well.Both of these cases are turbulent and hence non-trivial.In OpenFOAM, the standard pimpleFoam solver was employed in conjunction with the Gamma scheme for the convection term in the momentum equation to yield stable simulations in these cases.A similar approach has been validated and utilized earlier in the literature in various studies by the authors using scale-resolving simulations (see e.g.[66]).agreement.The key physical phenomena of the turbulent flows are marked in (a) as follows.I: The shear layer generated turbulence exhibits a Kelvin-Helmholtz instability a few window widths from the window (see also Fig. 9).II: Negative velocities and formation of recirculation zones are observed next to the walls aligned with the flow direction.III: The flow accelerates at the more narrow doorways.These physical phenomena are well known and expected and it is therefore crucial that they are correctly captured by DNSLABIB and in qualitative agreement with the OpenFOAM code.
Fig. 6(a) details the location of three sampling lines on the midcut plane, along which the x-velocity component is interpolated and compared between the two computational tools in (b) for the  = 5000 case.The comparison indicates good agreement between DNSLABIB and the OpenFOAM simulations.Notably, the comparison in the present flow setup is computationally rather demanding because the flow is highly transitional and the window width is large compared to the room dimensions so that the walls are relatively close to the shear layers.Considering these aspects, the present results for DNSLABIB at  = 5000 can be considered to be in very good agreement with OpenFOAM.For  = 50 000, the agreement is still satisfactory.We remind that for the lower Reynolds number all near-wall cells have + < 10 while 98% of them are below + < 5.For the higher Reynolds number approximately 88% of the cells have + < 10.Hence, the present results provide numerical evidence that, in particular for the  = 5000 case, the mean velocity gradients are well resolved all the way to the walls.

DNSLABIB performance on a GPU
Next, we discuss the simulation time for the three-room configuration emphasizing the performance of DNSLABIB in comparison to OpenFOAM.The number of total computational cells for the DNSLAB case is approximately 33.3 ⋅ 10 6 while the OpenFOAM case contains ∼ 16.2⋅10 6 cells.As noted earlier, the grids for both codes are consistent at the high shear regions of the flow, i.e. window/ventilation air jet regions, as well as at the walls (identical grid size ).However, due to computational limitations, the OpenFOAM simulations required utilizing a coarser grid resolution (grid size ≈ 2) outside these regions.The DNSLABIB run for a reference room simulation takes approximately 48 h on a single NVIDIA Tesla A100 GPU.In contrast, the parallel OpenFOAM simulation, executed on 1040 CPU cores (Intel Xeon Gold 6230), which was noted to be optimal for this case, consumes approximately 155 h of computational time for the respective simulated physical time.Not only is the DNSLABIB run by a factor of 3x faster than the OpenFOAM run with almost double the computational cell count, DNSLABIB can run on a more lightweight platform containing a single GPU, avoiding using a supercomputer.
A systematic DNSLABIB scaling test for the reference room case is shown in Fig. 7(a) detailing the computational time as a function of the mesh size as the mesh is refined.The computational time is defined as the execution time of a single (NVIDIA Tesla V-100) GPU to reach a simulation time of 120 s (starting from  0 = 0 s).The relationship is linear which is an ideal result indicating no superfluous overhead is generated in simulations involving a larger number of computational cells.This result is reasonable, since the simulation data is completely contained within the VRAM of the GPU in our benchmark cases and thereby any degree of extraneous communication overhead should be avoided.Additionally, Fig. 7(b) portrays the ratio of the execution time per simulation time step while executing DNSLABIB on both a single CPU core (Intel Xeon W-2123) and a more modest GPU (here, a NVIDIA GP104GL) as the system size  is increased.At low , executing DNSLABIB on a CPU yields superior performance as  CPU ∕ GPU < 1.However, the GPU implementation yields similar performance at  ≈ 4 ⋅ 10 4 for this case and with system sizes extending beyond this, the GPU implementation rapidly outpaces the CPU variant in terms of the required execution time.At  = 2 ⋅ 10 6 , the CPU execution time per time loop is 15-fold in this particular instance.However, the actual performance benefit is naturally dependent on the specific hardware and studied case.

Room setups and overview of the airflow
In addition to the previously validated case, Fig. 8 displays two additional ventilation setups with the white cloud portraying the entering fresh airflow which is assumed to be isothermal.Here, the studied rooms are empty even though simple block shaped furniture can be easily included in DNSLABIB.Fig. 8(a) corresponds to the original 3room case while (b) shows a case where the room dividers have been removed and the resulting space becomes one single open room.In panel (c), the room dividers remain removed but the ventilation airflow is now generated by vents located at the ceiling as the windows are closed.In reality, the incoming airflow is commonly directed by grids to enhance mixing, break strong air currents and thereby enhance the indoor airflow comfort as well.Here, however, the inflow is simply modeled as plain air jets.The parameters applied in these cases are detailed in Table 1 and in the following discussion, the cases are referred to as 3RW (3 rooms with windows), 1RW (1 room with windows) and 1RV (1 room with vents), respectively.For consistency, we have matched the inflow rate of air in these cases.
Panels (d)-(e) in Fig. 8 also display the instantaneous spatial patterns in the CO 2 concentration field for the three setups at the midcut plane at  = 20 s for 3RW (d), 1RW (e) and 1RV (f) cases respectively ( = 50 000).Here, the black (white) color designates areas of stale (fresh) air.For instance, in (e) and (f), numerous pockets of stale air appear which remain stagnant and poorly ventilated.Some of these pockets are expected and intuitive (room corners) while some are less intuitive such as the boundary of the wall separating the two Fig. 7.The computational time required as a function of the system (grid) size  plotted (a) for the 3RW case at  = 50 000 (high ventilation rate).The relationship seems linear, which is ideal, indicating there is no extraneous communication overhead during the simulation.The device memory (NVIDIA Tesla V100) inhibits further benchmarking beyond  = 2 ⋅ 10 7 .The execution time per a single simulation timestep can be additionally evaluated for the present code.Executing DNSLABIB on both a single CPU core (Intel Xeon W-2123) and a GPU (NVIDIA GP104GL), the ratio of the execution times on different hardware can be plotted as a function of  (b).

Table 1
The simulation parameters for the three room configurations for the two Reynolds numbers.windows in (e).The formation of stagnation zones is expected based on known fluid dynamics and flow recirculation near the room corners.Furthermore, turbulent airflow affects the mixing of the fresh and stale air.We note that qualitatively very similar results are observed for the cases with  = 5000 but the flow simply evolves 10 times slower due to the lower window airflow velocity (not shown herein for brevity).From the viewpoint of infection risk, strong variance in the CO 2 content of a room also indicates a potential spatial variance in the infection risk of an airborne disease.While indoor CO 2 measurements have recently been employed as a proxy to monitoring the virus concentration during the COVID-19 pandemic [67][68][69][70], one could argue that the information provided by CO 2 sensors can yield misleading output if such local variances are not quantified while designing the measurement setup.The observations highlight the importance of accounting for the uncertainty resulting from the geometric features of indoor spaces.Gaining a complete three dimensional description of the airflow characteristics is typically only accessible via scale-resolved CFD simulations.In practice, personal CO 2 meters offer real time air quality monitoring at the location of any individual.

Ventilation characteristics
A commonly used engineering metric for the ventilation rate is the air changes per hour ([ACH] = 1/h) defined as where  ([ ] = m 3 ) is the room volume and V is the volumetric airflow into the room ([ V ] = m 3 /h).In practice, ACH depends on the room airflow details, heat sources and geometrical features.ACH value can be measured by CO 2 measurements [71].Here, we estimate ACH as follows [72].A room is initially saturated with a relatively high CO 2 content stale air (here: 1000 ppm).Then, fresh outdoor air with a low CO 2 concentration (here: 400 ppm) is released into the room via windows or ventilation ducts.The mean CO 2 concentration as a function of time can then be monitored to yield the ACH value.Fig. 9 demonstrates the CO 2 concentration in the 3RW simulation ( = 50 000), where the stale air is gradually displaced by clean air over time.The window jets are noted to either impinge on the opposite wall in the back room or alternatively exit almost directly through the doorway to the corridor.Regions of lower ventilation performance can emerge near wall corners within re-circulation areas (e.g. the region in close proximity of the window at lower left corner of the frame).For the 1RW case, the absence of the additional walls imply that the window jets exit through the opposite windows with relatively less mixing of stale and clean air.2, vary between 82%-150% ( = 50 000) and 75%-139% ( = 5000) of the theoretical mixing ventilation value.For  = 50 000, the displayed ACH values are high  yet the values are consistent with some of the reported values obtained in experimental settings involving natural ventilation [73,74].For  = 5000, the ACH values (3.4-6.2) correspond to typical values observed in mechanical ventilation setups.For both Reynolds numbers, the ventilation generated by the ceiling vents (1RV) outperforms the various cross-draught setups (3RW, 1RW) and among the two cross-draught simulations, the space with room dividers exhibits enhanced ventilation performance.This is due to the improved mixing of the air masses via a combination of jet impingement, turbulence and flow re-circulation within the back room, which were already discussed in conjunction with Figs. 9 and 5.This suggests that the net impact of the solid obstacles on the ACH value is more ambiguous than one might anticipate as it may depend on the exact details of the airflow and airflow-obstacle interactions.However,

Table 2
The ACH values estimated for  = 5000/50 000 simulations via curve fitting  () =  0 exp(−ACH ⋅ ) +  1 to each ⟨CO 2 ⟩ displayed in Fig. 10.The last two columns detail the range of ACH values based on the error estimates (see Fig. 10) and their values normalized by the theoretical ACH value.the present numerical findings at  = 5000 and  = 50 000 based on full 3D numerical data imply that the theoretical ACH value may be highly inaccurate and off-set by a factor of 0.75-1.51.
In order to explore the fine structures in Fig. 9 and the CO 2 profiles they entail, the distribution of CO 2 content in these cases is examined in Fig. 11, where the normalized distribution function  (CO 2 ) is plotted for different instances of time, where  = 0 s denotes the start of the simulation.The respective DNSLABIB and OpenFOAM results for 3RW are displayed in (a) and (b) ( = 5000) and (c) and (d) ( = 50 000) while the profiles from the 1RW and 1RV simulations obtained with DNSLABIB are plotted in (c) and (d) ( = 50 000).Initially, the probability distribution peaks at CO 2 = 1000 ppm.As the simulation and the state of ventilation progresses, the distribution veers towards the lower end of the CO 2 spectrum as anticipated.Furthermore, the distribution functions clearly imply non-homogeneous mixing of the air, perceived as flat and uniform distributions, supporting the observations made earlier.
Therefore, based on the present numerical findings, we find that the airflow significantly affects the mixing patterns and dilution of the CO 2 concentration in the configurations considered here.The routinely employed definition of ACH is limited to conditions of perfect and extremely rapid mixing, rarely encountered in realistic indoor ventilation setups.In reality, as exemplified by the results in Fig. 10, the standard deviation in CO 2 concentration levels can be in the order of 10-20% of the mean value or higher, also pointing to non-homogeneous mixing of the air masses.However, the numerical results indicate that the main difference between the present empty rooms stems from the flow geometry and air supply configuration while the local variation of ACH can be relatively important as well.In the presence of large pieces of furniture (e.g.book shelves) or room dividers, the local variation of ACH may become more prominent which could be considered more in the future.

Infection risk
Having extracted the ACH values from each simulation case by numerical fitting in the previous discussion, we next proceed to the more practical implications of these results in terms of the infection risk.Therefore, a framework for relating the ACH and infection risk is required.During the COVID-19 pandemic, the classical Wells-Riley Fig. 10.The time-evolution of the mean carbon dioxide content in each of the ventilated room scenarios presented earlier for the high ventilation rate,  = 50 000 (a) and the low ventilation rate,  = 5000 (b).The concentration in the room, which is initially saturated with stale air, decreases in an almost exponential manner.Clear deviations from the theoretically derived behavior is observed.equation has been utilized extensively for infection risk assessment [3,[75][76][77][78][79][80].According to the Wells-Riley model, the infection probability (  ) can be calculated as follows where  is the number of infectious people in the modeled setting,  is the rate of generation of the infectious units termed ''quanta'',  is the respiratory rate of a person,  is the exposure time and V is the air exchange rate (in units of [m 3 /h]).This form of the model assumes (1) a steady state situation reached over a longer period of time during which the infectious emit virus to the air, (2) immediate and uniform mixing in the room so that distance from the source is not taken into account, and (3) constant removal of airborne particles by the ventilation.As a remark, under steady state conditions, the argument in the exponential function is simply the inhaled dose which is proportional to the average concentration ([  ] = 1/m 3 ) of quanta in the room air which ( = 1) can be calculated simply as follows [22] Here, we address the infection risk in the three indoor settings, assuming that an infectious person has occupied the space for a period of time, saturated the room with exhaled air (high CO 2 content) and dispersed infectious quanta to the space which remain infectious.Then, the infected person leaves the room, ventilation is commenced and the infection risk for a person entering the room starts accumulating.We therefore rewrite the Wells-Riley model as follows where () is the effective dose a person has accumulated during time  and   () is the time-dependent average concentration of quanta in the room air.We assume that   is directly proportional to the previously discussed mean CO 2 content, i.e.   () =  0 exp(−ACH ⋅ ), where  0 represents the initial, homogeneous concentration of quanta in the room saturated with stale air.This assumption considers only the small aerosols which remain airborne for very long periods of time (see next section).Fig. 12(a) presents the infection risk based on Eq. ( 19) with the initial values of (a) 100 quanta and (b) 500 quanta homogeneously spread in the room volume,  0 = {100, 500}∕(8 ⋅ 8 ⋅ 3) ≈ {0.52, 2.60} quanta/m 3 .In the COVID-19 context [81,82] such values could be representative to a person performing activities in an indoor setting over a 1 h period, releasing pathogens either at a moderate (medium vocal activities, such as talking) or very high rate (singing), respectively.The present demonstrative cases are displayed for  = 50 000 and  = 5000 plotted with solid and transparent curves respectively.Furthermore, as the respiration rate of a person, the value of  = 1.2 m 3 ∕h is applied herein [22].The key observation from panel (a) is that a high enough ventilation rate reduces the average infection risk extremely efficiently.
The average infection risk approaches ≈ 24% probability level if ventilation is switched off for the considered time window of 30 min compared to the probability of 0.9-1.9%calculated for the wellventilated cases.Even in the worst-performing ventilation setup (1RW), the infection risk reduces by a factor of 10.Similar reductions are also apparent in the results presented in panel (b), where the absolute reduction in the infection risk is even greater (≈ 80%).For the lower ventilation rate ( = 5000), the risk is typically reduced by a factor of 2 for each setup.For instance, in (a), the ventilation setups at low inflow rate reduce the infection risk from 24% to 9%-13% and in (b), from 80% to around 37%-53%.Yet, this can be considered to be a significant gain.As demonstrated herein, window ventilation and enhancement of mechanical ventilation could be considered to be powerful complementary tools in keeping the society open and reducing infection risks in public places such as schools, choir practices, shops and bars.As an additional note, HEPA filters deliver clean air at a certain volumetric flow rate and in analogy with ACH definition, an effective air change eACH (volume of filtered air/room volume) can be defined.HEPA filters offer an energy efficient way in reducing indoor virus concentration to increase the effective ACH value ACH * = ACH + eACH and correspondingly lower the virus concentration indoors.In the future, DNSLABIB could be used to model air filtration devices positioned at different indoor locations via the volumetric source terms.

Dispersion of airborne particles
In the previous section, we discussed the reduced infection risk associated with various ventilation setups.The estimates considered only the fraction of the airborne particles which are small enough so that they can be directly correlated with the indoor CO 2 concentration.But what size particles are small enough to follow the airflow?During the COVID-19 pandemic there has been a major scientific debate on the distinction between an aerosol and a droplet.In the aerosol community it is well understood that particles up to sizes of 50-100 μm are able to easily remain airborne over extended times and distances because they evaporate quickly [22].However, until COVID-19, the medical literature, including WHO in their early guidance in 3/2020, adhered to an erroneous 5 μm cutoff.Such an unfortunate misconception biased the early attention towards surface transmission which was a major error corrected later on in the pandemic when COVID-19 was noted to be airborne by WHO [29].Presently, there is a broad scientific consensus on the airborne route as a major driver for the ongoing pandemic [5,83].Next, we discuss how solid particles between 1-100 μm travel in the air using DNSLABIB.Fig. 13 presents two time snapshots of an exhaling person 0.2 s (a) and 0.8 s (b) after start of exhalation.The air pulse is modeled as a particle laden jet with a diameter of  = 3 cm.The Reynolds number of the exhaled jet is  = 7000.The gray cloud represents the passive scalar field generated in this expulsion of air posing high CO 2 concentration while the red/green/blue droplets present large (90-100 μm), medium (10-30 μm) and small (1-10 μm) particles emitted Fig. 12.The average infection risk accumulates as a function of time.Here, a scenario is illustrated where an infectious person has released virus to the air and a susceptible person arrives to the room at time  = 0 s.Initial conditions (a)  0 = 100 quanta∕  and (b)  0 = 500 quanta∕  .All configurations are highly effective in reducing the risk of infection at higher ventilation rates.At lower ventilation rates (transparent lines), the risk reduction is also significant.from the airways, respectively.These size classes are chosen in order to comply with the clinical evidence from human expired aerosol measurements [84].The distribution of the particle sizes in our simulation is 20% (90-100 μm), 53.3% (10-30 μm) and 26.7% (1-10 μm) of the total particle count (200), respectively.Here, particle evaporation is neglected so that the observed picture corresponds to a conservative situation where particles sink significantly faster than they would sink in reality e.g. at a RH = 30%-40% indoor relative humidity [22].
Notably, the smaller particles up to 30 μm remain airborne and they are transported easily through the air by the spreading air jet for considerable distances, the trajectories having a clear correlation to the exhaled CO 2 plume.This is consistent with our previous findings [22] where particles of size 20 μm were shown to be transported over shelves of a supermarket between two aisles.However, the particle behavior also simply follows from the particle sedimentation time in still air with   being 36 min, 9 min and 1 min for particle sizes of 5 μm, 10 μm and 30 μm, respectively.The results displayed here also suggest revisiting social distancing guidelines and protection measures in combating infectious respiratory disease.Within one minute, the two individuals separated from the initial aerosol source by two large walls have been exposed to the smallest aerosols.Since the aerosols travel from room to room and further to the corridor, this example clearly supports the notion that in general, room dividers, plexi-glasses or visors do not provide sufficient protection from aerosols.

Conclusions
In the present work, we endeavored to create a CFD tool for scaleresolved simulations with reduced computational effort.Therefore, the DNSLABIB software was programmed in the MATLAB language and implemented on a single GPU.Currently, DNSLABIB performance is mainly limited by the GPU memory.The largest system sizes that Here, an exhalation is modeled as a particle laden jet.In this particular example, solid and non-evaporating particles of size 30 μm are noted to remain airborne and easily travel horizontally to reach the airways of another person.For practical mucus droplets, evaporation shifts the critical size to a much larger value, up to 100 μm [5,22].we have explored are currently in the order of 100M computational cells which is possible with a single NVIDIA A100 (40Gb RAM) GPU.Extending the capabilities of MATLAB to the realm of CFD, DNSLABIB provides the end-user with prior CFD experience a starting point for further exploration.Considering the scope of the paper, the focus herein was on isothermal flows which are relevant to e.g.window ventilation setups, where the indoor air temperature is close to the outdoor temperature and buoyancy effects are not the primary concern.We note that, in the 2020's, a major portion of the world population still lives in circumstances with relevance to these conditions.DNSLABIB implements solid obstacles in a simplified manner via the Immersed Boundary (IB) method and the simulations are executed on a GPU.DNSLABIB was validated in two canonical reference cases, the pressure-driven channel flow and a channel flow past a bluff body.Additionally, validation was conducted against experimental data of a model building in a wind tunnel.The code was utilized to study three separate indoor ventilation configurations using scale-resolving simulations.In one of these cases, a comparison against a respective OpenFOAM simulation was performed as well.A superior performance of DNSLABIB over the corresponding OpenFOAM implementation was discovered, the speed-up factor being in the order of 3 while avoiding usage of a supercomputer which is considered as a major step advocating the usage of GPUs for indoor airflow assessment.
The ventilation characteristics in the three cases were studied by monitoring the CO 2 concentration as it is routinely adopted as a measured proxy for airborne transmission risk.The CO 2 distributions and mean value over time revealed highly inhomogeneous mixing of air, contrasting the common assumption of ideally mixed air employed to determine the ACH value.The actual ACH for each case was determined and significant deviations from the theoretical value were noted.Collectively, the results indicated the presence of strong local variations in the CO 2 concentration indoors.This further emphasizes the need to better understand the full 3D airflow features indoors using scaleresolved simulations to ensure high air quality both locally and on average.
The extracted ACH parameters were further applied in the Wells-Riley model to assess the infection risk associated with each ventilation setup in the context of COVID-19.The provided examples indicate significantly reduced infection risk if windows are opened immediately after entering a room with stale, virus-contaminated air.However, it should be highlighted that the actual benefits of ventilation emerge only if the inhalation dose  remains small enough which can be achieved by (1) minimizing exposure time and (2) minimizing airborne virus concentration via enhanced ventilation and/or air filtration or masks.Even in the most conservative scenarios with the lowest ventilation rates (ACH = 3.4-6.3),we provided numerical examples on cases where the infection risk was reduced by a factor of 2, which is very significant per se.For the well ventilated cases (ACH = 37-67), over a 10-fold decrease in the infection risk was observed.
Finally, the transportation of respiratory particles in an exhaled jet was investigated.The studied condition resembles a case with a relative humidity of RH = 100% where particle size reduction due to evaporation does not occur [22].In such a conservative setup, solid particles up to 30 μm in size were witnessed to remain airborne for prolonged periods of time and they were therefore shown to affect observers in the considered premises over extended distances as well.It is clear that the traditional 5 μm cutoff is simply erroneous.However, the numerical example in question indicates that the large particles with sedimentation time on the order of ∼1-10 s will exit the exhaled air jet region quickly after which they will settle on the floor and surfaces.For large particles (here: 90-100 μm), it is clear that the direction and strength of the exhaled and the ambient airflow will play a major role in how well those particles may reach other people's airways nearby before settling down.The medium particles (here: 10-30 μm) are clearly able to stay airborne for a long duration of time and travel over extended distances.The smallest particles (here: <10 μm) are able to remain airborne for very long periods of time.The inhalation of small and medium aerosols is the most likely route of virus transmission as they can be unconditionally inhaled over short and long distances and their viral content [1] as well as their number-concentration is abundant [84].
The complexity underlying the ventilation cases studied here also suggests numerous venues for future research.A logical continuation could be the development of a wall model to more accurately consider boundary layers at high shear surfaces.However, as noted here in particular for  = 5000, this aspect may not be completely critical for low-speed indoor airflows if most of the near-wall + values remain small enough.Additionally, since indoor airflow can be considerably affected by buoyancy effects, such as heat generated by the occupants, adding an appropriate coupling (e.g.Boussinesq approximation) between the heat sources (sinks) and the momentum equation would be reasonable as well.

A.2. Model building in a wind tunnel
In the work by Jiang et al. [86], the authors designed an experimental setup where a hollow model room structure with an open doorway which was placed in a wind tunnel and then subjected to turbulent flow with constant (average) inflow velocity.The authors additionally performed LES simulations to compare with the experiments.Due to the meticulous and extensive gathering of experimental data on multiple line sources by the authors, this ventilation study of a building serves as a benchmark for DNSLABIB.This high  benchmark is challenging because DNSLABIB does not allow local grid refinement near the walls.As described in [86], the door of the model building is placed towards the flow.The dimensions of the model setup are detailed in the upper left figure in Fig. A.17 while the computational details are displayed in Table A.5.The computational domain size is 8 by 4 by 4 in the streamwise () and spanwise directions (, ), respectively.At the inlet, the freestream velocity is chosen to produce the Reynolds number   = 1.4 ⋅ 10 5 , where  refers to the model building edge length.The computational domain is bounded by solid walls to emulate A model geometry of a building with a windward opening was set in a wind tunnel in [86].The dimensions of this miniature is detailed in the upper left figure of the panel.The upper right figure shows the instantaneous streamwise velocity and the locations of the line sources, at which the velocity profiles are compared with DNSLABIB and the experiments in [86].Here,  refers to the streamwise component and  to wall-normal components of the velocity, respectively.The mean quantities are plotted over the normalized height coordinate Ẑ.Additionally, the forcing regions for the velocity and the schematic representation of the associated mask functions are displayed.

Table A.5
The comparison between the simulation parameters in this work and the LES simulations as reported by Jiang et al. a wind tunnel environment and no-slip boundary conditions are applied for the velocity at solid walls.Additionally, periodic boundary conditions are imposed on the inlet and outlet, where the velocity forcing region is modeled with a smoothly varying hyperbolic function (tanh) approximating a flat-top profile in the  and  directions, respectively.In the forcing region, the streamwise velocity is forced to a profile  (ℎ) =  0 exp(−12.0⋅ (ℎ)) in the -and -directions, where ℎ denotes the height from each wall.This yields a profile which is always within 5% of the logarithmic velocity profile witnessed in the experiments in Ref. [86].The forcing regions at the inlet and the outlet have a width of 0.6 in the streamwise () direction.
The windward opening allows air to enter the model building and the velocity profiles across line sources (upper right figure) are compared between DNSLABIB and the experiments at the centerline of the model.Qualitatively, the results show the correct flow patterns: (1) a coherent vortex structure on top of the building and (2) strong recirculation zones both on the front and the back sides of the model building.The streamwise (lower row, Figs. 1 and 3 from the left) and the spanwise components of the velocity (lower row, Figs. 2 and 4 from the left) display a good quantitative agreement inside the model building while the highest discrepancies are observed outside of the model at the vortex forming above the roof of the building as expected.
Here, the size of the eddy is predicted correctly by DNSLABIB but the steep relative shifts between the minima and maxima in the streamwise velocity are not captured.This is an expected result, however, since the effects of the high Reynolds number employed in the experiment are not as pronounced inside the building as they are outside of the building.DNSLABIB does not allow wall refinement and thereby the wall shear stress is not correctly captured for such high  flow conditions.This could be mostly likely remedied by introducing a suitable wall function to the DNSLABIB simulation which is beyond the scope of the present work, however.The present high  test case complements the previous 2D validation study showing the advantages and limitations of DNSLABIB under such challenging flow conditions.

Fig. 1 .
Fig. 1.In the demonstration above, a window is located at a wall (left) and air is entering through the window.A plane, spanned by the blue lines, displays the inflow velocity field, constrained in place with the mask (, , ) (middle).Furthermore, the mask along the red line in the left-hand-side figure is plotted on the right.  + ∇ ⋅ () =   ∇ ⋅ [(, , )∇] + (, , ) ⋅   −    (3)

Fig. 2 .
Fig.2.The spectral approach employed in solving the pressure equation in this work necessitates maintaining periodicity in the solution.This may be achieved naturally by a suitable inflow/outflow setup (left) or extending the fluid region beyond the system of interest (right).

Fig. 3 .
Fig. 3. Illustration of the computational grid points in the solid walls (black symbols) and the fluid region (gray symbols) for the present immersed boundary method.Here,  = 0∕1 inside the wall/fluid parts.

Fig. 4 .
Fig. 4. A flow chart illustrating the code execution pipeline.Importantly, the GPU/CPU executions utilize the same base code.The arrays are initially converted into a GPU or CPU compliant format and are then passed to the functions implemented in MATLAB.Depending on the format of the arrays, the functions are automatically executed on either a GPU or CPU.

Fig. 5 .
Fig. 5. Instantaneous -component of velocity at the midcut plane using DNSLABIB (a) and OpenFOAM (b).The time-averaged velocity fields obtained in DNSLABIB and OpenFOAM are visualized in (c) and (d) respectively.I Kelvin-Helmholtz instability, II recirculation zones and III flow acceleration.

Fig. 5
illustrates the setup in more detail for  = 50 000.The  = 5000 case would remain qualitatively very similar with 10× slower timescales due to the 10-fold lower velocity scales.The midcut plane displaying the cross-section of the room is portrayed here at  = 1.5 m.The -component of the instantaneous velocity field is presented comparing DNSLABIB (a) and OpenFOAM (b) at the midcut.The panels (c) and (d) present the temporally averaged -component of velocity for DNSLABIB and OpenFOAM respectively.Here, red color implies positive value for the velocity component, while blue suggests negative values.First, the mean velocity data in (c) and (d) indicate a good

Fig. 6 .
Fig. 6.The time averaged -component of velocity obtained in DNSLABIB and OpenFOAM along three sampling lines (right) for  = 5000 (low ventilation rate).The sampling lines are displayed on the left at the midcut of the simulation geometry and the arrows indicate the plotting direction (from 0 to 8 m).The agreement appears very good.

Fig. 8 .
Fig. 8. Three ventilation setups with an identical airflow rate.(a) Cross-draught via windows and the space is divided to 3 smaller rooms (3RW).(b) Cross-draught via windows without room-dividers (1RW).(c) Vents located at the ceiling without room-dividers (1RV).(d)-(e) display slices of the instantaneous CO 2 fields for the respective cases.
The mean CO 2 concentration can be monitored to yield the actual ACH value of each ventilation setup which differs from the theoretical ACH value.In Fig.10, the mean CO 2 content as a function of simulation time  is plotted for the studied cases.Panel (a) corresponds to simulations at  = 50 000 and (b) to  = 5000.The standard deviation is plotted as well with the error bars to indicate the uncertainty of the CO 2 distributions for each case.Simultaneously, the analytical expression 400 ⋅ exp(−ACH ⋅ ) + 600 based on the ventilation theory of perfectly mixed air is displayed.The theoretical ACH coefficient is determined as the ratio of the flow mass fluxes involved and the room volume for both  = 50 000 and  = 5000 as ACH = ( 2 ⋅ 1.2 m 2 ⋅ 1.0 m∕s ) ∕ [8 m ⋅ 8 m ⋅ 3 m] = 45 1/h and ACH = ( 2 ⋅ 1.2 m 2 ⋅ 0.1 m∕s ) ∕ [8 m ⋅ 8 m ⋅ 3 m] = 4.5 1/h, respectively.The computational ACH values for the various cases can be acquired by imposing a fit of the form  () =  0 exp(−ACH   ⋅ ) +  1 to the data presented in Fig. 10.These values, displayed in Table

Fig. 11 .
Fig. 11.The distribution function  (CO 2 ) of the carbon dioxide content at  = 200 s,  = 600 s,  = 1600 s for the 3RW case at  = 5000 (low ventilation rate) in DNSLABIB (a) and OpenFOAM (b).Additionally, the distribution is displayed at  = 20 s,  = 60 s,  = 160 for these cases at  = 50 000 (high ventilation rate) in (c) and (d).Finally, the profiles are also plotted for the DNSLAB simulations of the 1RW case (e) and 1RV case (f).

Finally
, Fig. 14 illustrates the aerosol emitting person placed in a room with a cross-draught and room dividers (3RW).The person emitting the infectious aerosols is located in the left upper corner of the room at (a)  = 0 s, (b)  = 20 s, (c)  = 40 s and (d)  = 100 s respectively.The simulations suggest that the smallest aerosols (≤30 μm) can rapidly travel significant distances following the airflow, circumventing the solid walls.

Fig. 13 .
Fig. 13.Here, an exhalation is modeled as a particle laden jet.In this particular example, solid and non-evaporating particles of size 30 μm are noted to remain airborne and easily travel horizontally to reach the airways of another person.For practical mucus droplets, evaporation shifts the critical size to a much larger value, up to 100 μm[5,22].

Fig. 14 .
Fig. 14.A simulation of a room with the person emitting infectious aerosols located in the top left corner (3RW case at  = 50 000, high ventilation rate).The (a), (b), (c) and (d) panels correspond to the simulation times  = 0 s,  = 20 s,  = 40 s and  = 100 s, respectively after the onset of aerosol emission by the person.

Fig. A. 15 .
Fig. A.15. Velocity in the laminar channel flow (a) and matching the analytic, parabolic velocity profile (b).

Fig. A. 16 .
Fig. A.16.A cube in a channel flow.The -component of velocity obtained in DNSLABIB (a) and OpenFOAM (b).

Fig
Fig. A.17.A model geometry of a building with a windward opening was set in a wind tunnel in[86].The dimensions of this miniature is detailed in the upper left figure of the panel.The upper right figure shows the instantaneous streamwise velocity and the locations of the line sources, at which the velocity profiles are compared with DNSLABIB and the experiments in[86].Here,  refers to the streamwise component and  to wall-normal components of the velocity, respectively.The mean quantities are plotted over the normalized height coordinate Ẑ.Additionally, the forcing regions for the velocity and the schematic representation of the associated mask functions are displayed.

Table A .
4 A cube in a channel flow: comparison of the Strouhal number and the drag coefficient obtained by DNSLABIB and OpenFOAM on two different grid resolutions.