key: cord-0472874-g09l6sxc
authors: Korhonen, M.; Laitinen, A.; Isitman, G. E.; Jimenez, J. L.; Vuorinen, V.
title: A GPU-accelerated computational fluid dynamics solver for assessing shear-driven indoor airflow and virus transmission by scale-resolved simulations
date: 2022-04-05
journal: nan
DOI: nan
sha: 1dc84e25bbcf366af77dcbe841475359eddc5c40
doc_id: 472874
cord_uid: g09l6sxc

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.

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, consistent and mounting evidence on the airborne transmission of the SARS-CoV-2 virus has accumulated [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] . Within the aerosol physics community, the suspension time of airborne particles in air has been well-established for a century. Therefore, during the pandemic, the fluid physics research community has further revisited various factors affecting particle transport in the air. These include the impact of 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 scale-resolving 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 API, such as 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 scale-resolved 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 equa-3 tions, 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 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.

In the present work the focus is on low-speed, isothermal gas flows which can be modeled using the incompressible Navier-Stokes equations. 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 is also studied assuming one-way 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 u is the velocity, p is the pressure divided by the fluid density ρ, ν is the kinematic viscosity, b is a body force, c is the passive concentration field and α c 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: S(x, y, z) (u set − u) /τ f , and b. 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 u and p can be provided as usual. In brief, the term S(x, y, z) (u set − u) /τ f is formulated 4 Figure 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 S(x, y, z) (middle). Furthermore, the mask along the red line in the left-hand-side figure is plotted on the right.

in order to establish the desired velocity u set within the relaxation time scale τ f . By setting S(x, y, z) = 1 at the specific location, where the velocity must reach the value u set and S(x, y, z) = 0 otherwise, the coordinate dependent mask function S(x, y, z) 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

where x c , y c and z c are the volumetric center coordinates of the source region, while W x , W y and W z are the window/inlet/ventilation duct dimensions of this region. This function obtains values in the range [0,1] and B x,y,z 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 x-axis (right). Here, the width of the transition region, denoted by δ, is highly dependent on the B x 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. 5 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 using a 6th order filter. The filter is defined as follows

where the filter coefficient is chosen so that the Nyqvist frequency is zeroed in the Fourier space i.e. γ xi = ∆x 6 i π 6 . The filter resembles a standard hyperviscosity term but avoids the cross-derivatives [65] .

The pressure equation requires particular attention near the walls. In the conventional projection method, one obtains a u * from the Navier-Stokes equation without a pressure gradient and then corrects u * 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 ∇ · β(x, y, z)∇p. The Fourier 6 transform of this equation can not be directly determined. However, by adding and subtracting 1 from β, one has ∇ · β(x, y, z)∇p = ∇ · (1 + β(x, y, z) − 1)∇p. 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 (9) n times evaluating the right hand side of the equation by central differences from the previous available value of pressure (p k ) and the non-solenoidal velocity u * acquired from the momentum equation. Generally, the equation converges very quickly and a hard-coded value n = 4 is used here. The velocity field is then corrected as u = u * − ∇p, utilizing the converged solution for the pressure to yield the solenoidal field satisfying Eq. (1).

Proceeding to Lagrangian particles, the equations of motion (EoM) read

where Re p is the particle Reynolds number, τ p is the particle settling time scale, describing the delay of the particle in adjusting to altered flow conditions, and C d refers to the drag coefficient of a particle. The subscripts p and f indicate particle and fluid properties, respectively, g is the gravitational force, d 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 τ p 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. Two key quantities from the EoM above are the particle terminal velocity v * = gτ p as well as the particle sedimentation time from height h expressed as τ s = h/v * . For h = 1.5 m and d=5/10/30 µm solid particles τ s 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.

As stated earlier, we implement a numerical code to solve the governing equations (2), (3) and (9) 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 7 is illustrated in Fig. 3 . The user initializes a case and controls the subsequent simulation principally via the "SetParameters.m", "CreateFields.m", "CreateGeometry.m", "CreateSourceMasks.m" and "InitializeDrops.m" scripts. In "SetParameters.m", the user provides the necessary information regarding the simulation geometry, fluid properties and data outputting. Importantly, 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 β(x, y, z) field with primitive shapes. Then, in "CreateSourceMasks.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 "CreateGpuArrays.m". Finally, the number, size and location of Lagrangian particles is assigned in "InitializeDrops.m", after which 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 above. "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 "ContructVelocityIncrement.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 8 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 .mat file for quick access and analysis.

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

• Avoiding the use of a supercomputer.

• Ease of use in configuring and executing a case for simple geometries with solid objects.

• 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.

• Wall models are currently not implemented in the code and one needs to rely on constant grid resolution even at the walls.

• The obstacles are presented as simple blocks.

• Thermal sources/sinks are not presently accounted for although they are known to be of high importance in indoor airflow configurations.

In the present paper, the main focus 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 9

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 Re = 5000 (window airflow 0.1m/s) and Re = 50000 (window airflow 1m/s) to better understand how Re 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 y+ < 5 − 10 are mostly reached. For Re = 50000, 88% of the computational cells have y+ < 10. For Re = 5000 all near-wall cells have y+ < 10 while ≈98% of the near-wall cells have y+ < 5.

Next, the performance of DNSLABIB is explored by introducing several simulation cases of increasing complexity. First, the code is validated in two canonical reference cases: the laminar channel flow and vortex shedding due to a rectangular body mounted into the channel. Then, we further validate the code against a scale-resolved simulation of an indoor ventilation setup performed in OpenFOAM. A number of indoor ventilation setups are then examined to study the removal of stale air from the room and assess the infection risk in these configurations in the context of COVID-19. First, a laminar channel flow driven by a pressure gradient is examined to test the IB method. The simulation parameters for this case are detailed in Tab. 1. The flow is initialized by setting a body force (acceleration) in the x-direction b =b/ρ = (b x , 0, 0). In this simple channel flow, the Navier-Stokes equations imply an analytical velocity profile across the channel U (y)

fixing the ratio b x /ν fixes the Reynolds number, which is set to Re = 500 here. As noted in Fig. 4 b) , the present simple example indicates that the analytical parabolic velocity profile is recovered in the simulation with the relative error remaining below 10 −5 .

Next, as depicted in Fig. 5 a) , a cube is placed at the centre line of the channel. The purpose of the test case is to demonstrate the performance of DNSLABIB when the wake interacts dynamically with the channel walls, resulting in vortex shedding. The cube is located at (1.2D, 0.5D) and the side length is h = D/6 while the rest of the geometry remains unaltered from the previous case discussed. Setting ∆x = ∆y we explore the case utilizing two uniform mesh densities: ∆x = h/30 and h/60. A parabolic velocity profile at the left-hand-side of the channel is forced, with no-slip conditions at the channel walls and a maximum value U max in the center of the inflow region (y = D/2). This translates to the obstacle Reynolds number of Re ob = U b · h/ν = 2 3 U max · h/ν ≈ 84 being in the 11 unsteady vortex shedding regime. The subsequent Kármán vortex street is examined in more detail in Fig. 5 , where the x-component of the velocity field is examined in steady-state in both DNSLABIB (a) and OpenFOAM (b). The reference OpenFOAM simulations are carried out with the standard incompressible PIMPLE-solver utilizing a second order accurate backward method in time, the linear (corrected) interpolation for the convection terms, and linear central differencing for the diffusion terms (see e.g. [66] for a similar approach). In contrast to DNSLABIB, OpenFOAM boundary conditions are provided in a standard manner using Dirichlet/Neumann conditions. As expected, both approaches indicate periodic vortex shedding at a distinct oscillation frequency f . 

In the following, airflow is simulated in a more challenging validation case at Re = 5000 involving a large indoor space (8 m x 8 m x 3 m along x-, y-and z-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 Re = 5000 with window airflow peak velocity U in = 0.1m/s, we also investigate a higher inflow velocity of U in =1m/s corresponding to Re = 50000. The main validation is based on the Re = 5000 case which is also rather challenging in terms of fluid dynamics while the Re = 50000 case is investigated to better understand the Re 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. Fig. 6 illustrates the setup in more detail for Re = 50000. The Re = 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 z = 1. 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. 7 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 Re = 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 Re = 5000 can be considered to be in very good agreement with OpenFOAM. For Re = 50000, the agreement is still satisfactory. We remind that for the lower Reynolds number all near-wall cells have y+ < 10 while 98% of them are below y+ < 5. For the higher Reynolds number approximately 88% of the cells have y+ < 10. Hence, the present results provide numerical evidence that, in particular for the Re = 5000 case, the mean velocity gradients are well resolved all the way to the walls.

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. The DNSLABIB run for a reference room simulation takes approximately 48 hours on a single NVIDIA Tesla A100 GPU. In contrast, the parallel OpenFOAM simulation, executed on 1040 CPU cores (Intel Xeon Gold 6230), consumes approximately 155 hours of computational time for the respective simulated physical time. Not only is the DNSLABIB run by a factor of 3x faster than the Open-FOAM 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. 8 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 seconds (starting from t 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. 

In addition to the previously validated case, Fig. 9 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. 9 a) corresponds to the original 3-room 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 Tab. 3 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. 9 also display the instantaneous spatial patterns in the CO 2 concentration field for the three setups at the midcut plane at t = 20 s for 3RW (d), 1RW (e) and 1RV (f) cases respectively (Re = 50000). 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 windows in e). The formation of stagnation zones is Figure 9 : 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. 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 Re = 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 Property 1 room (windows) 3 rooms (windows) 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.

where V ([V ]=m 3 ) is the room volume andV 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. 10 demonstrates the CO 2 concentration in the 3RW simulation (Re = 50000), 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.

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. 11 , the mean CO 2 content as a function of simulation time t is plotted for the studied cases. Panel a) corresponds to simulations at Re = 50000 and b) to Re = 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 · t) + 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 Re = 50000 and Re = 5000 as ACH = 2 · 1.2m 2 · 1.0m/s / [8m · 8m · 3m] = 45 1/h and ACH = 2 · 1.2m 2 · 0.1m/s / [8m · 8m · 3m] = 4.5 1/h, respectively. The computational ACH values for the various cases can be acquired by imposing a fit of the form f (t) = C 0 exp(−ACH f it · t) + C 1 to the data presented in Fig. 11 . These values, displayed in Tab. 4, vary between 82% -150% (Re = 50000) and 75% -139% (Re = 5000) of the theoretical mixing ventilation value. For Re = 50000, 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 Re = 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 Figure 11 : The time-evolution of the mean carbon dioxide content in each of the ventilated room scenarios presented earlier for the high ventilation rate, Re = 50000 (a) and the low ventilation rate, Re = 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. The ACH values estimated for Re = 5000/50000 simulations via curve fitting f (t) = C 0 exp(−ACH · t) + C 1 to each CO 2 displayed in Fig. 11 . The last two columns detail the range of ACH values based on the error estimates (see Fig. 11 ) and their values normalized by the theoretical ACH value.

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. 10 and 6. 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, the present numerical findings at Re = 5000 and Re = 50000 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. 10 and the CO 2 profiles they entail, the distribution of CO 2 content in these cases is examined in Fig. 12 , where the normalized distribution function f (CO 2 ) is plotted for different instances of time, where t = 0 s denotes the start of the simulation. The respective DNSLABIB and OpenFOAM results for 3RW are displayed in a) and b) (Re = 5000) and c) and d) (Re = 50000) while the profiles from the 1RW and 1RV simulations obtained with DNSLABIB are plotted in c) and d) (Re = 50000). 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. 11 , 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. Additionally, the distribution is displayed at t = 20 s, t = 60 s, t = 160 for these cases at Re = 50000 (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).

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 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 (P inf ) can be calculated as follows

where I is the number of infectious people in the modeled setting, q is the rate of generation of the infectious units termed "quanta", p is the respiratory rate of a person, t is the exposure time andV 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 ([c q ]=1/m 3 ) of quanta in the room air which (I=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 Q(t) is the effective dose a person has accumulated during time t and c q (t) is the time-dependent average concentration of quanta in the room air. We assume that c q is directly proportional to the previously discussed mean CO 2 content, i.e. c q (t) = C 0 exp(−ACH · t), where C 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. 13 a) presents the infection risk based on Eq. (18) with the initial values of (a) 100 quanta and (b) 500 quanta homogeneously spread in the room volume, C 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 hour 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 Re = 50000 and Re = 5000 plotted with solid and transparent curves respectively. Furthermore, as the respiration rate of a person, the value of p = 1.2m 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 minutes compared to the probability of 0.9-1.9 % calculated for the well-ventilated 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 (Re = 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, 22 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.

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. 14 presents two time snapshots of an exhaling person 0.2 seconds (a) and 0.8 seconds (b) after start of exhalation. The air pulse is modelled as a particle laden jet with a diameter of D = 3 cm. The Reynolds number of the exhaled jet is Re = 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 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 the particle sedimentation time in still air with τ s being 36 min, 9 min and 1 min for particle sizes of 5 µm, 10 µm and 30 µm, respectively.

Finally, Fig. 15 illustrates the aerosol emitting person placed in a room with a crossdraught and room dividers (3RW). The person emitting the infectious aerosols is located in the left upper corner of the room at a) t = 0 s, b) t = 20 s, c) t = 40 s and d) t = 100 Figure 14 : 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] . s respectively. The simulations suggest that the smallest aerosols (≤ 30 µm) can rapidly travel significant distances following the airflow, bypassing the solid walls.

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.

In the present work, we endeavored to create a CFD tool for scale-resolved simulations with reduced computational effort. Therefore, the DNSLABIB software was programmed in the MATLAB language and implemented on a 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.

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. 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 25 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 scale-resolved 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 Q 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 seconds 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 near-by 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 Re = 5000, this aspect may not be completely critical for low-speed indoor airflows if most of the near-wall y+ 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.

The DNSLABIB package is freely available at https://github.com/Aalto-CFD/DNSLABIB.

Airborne transmission of respiratory viruses

COVID-19: the case for aerosol transmission

High-resolution large-eddy simulation of indoor turbulence and its effect on airborne transmission of respiratory pathogens-model validation and infection probability analysis

Consideration of the aerosol transmission for COVID-19 and public health

Aerosol transmission of SARS-CoV-2? evidence, prevention and control

Transmission of COVID-19 virus by droplets and aerosols: A critical review on the unresolved dichotomy

The flow physics of COVID-19

Persistence of severe acute respiratory syndrome coronavirus 2 in aerosol suspensions

Aerosol and surface stability of SARS-CoV-2 as compared with SARS-CoV-1

Identifying airborne transmission as the dominant route for the spread of COVID-19

Airborne transmission of severe acute respiratory syndrome coronavirus-2 to healthcare workers: a narrative review

COVID-19 vulnerability: the potential impact of genetic susceptibility and airborne transmission

Probable airborne transmission of SARS-CoV-2 in a poorly ventilated restaurant

Modelling airborne transmission of SARS-CoV-2 using CARA: risk assessment for enclosed spaces

Spread of infectious agents through the air in complex spaces

The design of safe classrooms of educational buildings for facing contagions and transmission of diseases: A novel approach combining audits, calibrated energy models, building performance (BPS) and computational fluid dynamic (CFD) simulations

Simulation study on indoor air distribution and indoor humidity distribution of three ventilation patterns using computational fluid dynamics

Numerical investigation of aerosol transport in a classroom with relevance to COVID-19

The role of air conditioning in the diffusion of Sars-CoV-2 in indoor environments: A first computational fluid dynamic model, based on investigations performed at the vatican state children's hospital

Modeling transient particle transport in transient indoor airflow by fast fluid dynamics with the markov chain method

Predicting airborne coronavirus inactivation by far-UVC in populated rooms using a high-fidelity coupled radiation-CFD model

Modelling aerosol transport and virus exposure with numerical simulations in relation to SARS-CoV-2 transmission by inhalation indoors

Numerical study of three ventilation strategies in a prefabricated COVID-19 inpatient ward

On airborne virus transmission in elevators and confined spaces

Modeling airborne pathogen transport and transmission risks of SARS-CoV-2

Modelling airborne transmission and ventilation impacts of a COVID-19 outbreak in a restaurant in guangzhou, china

Investigating the influences of ventilation on the fate of particles generated by patient and medical staff in operating room

Fluid dynamics simulations show that facial masks can suppress the spread of COVID-19 in indoor environments

Coronavirus disease (COVID-19): How is it

Fifty years of CFD for room air distribution

LES over RANS in building simulation for outdoor and indoor applications: a foregone conclusion?, in: Building Simulation

GPU computing in medical physics: A review

Recent progress and challenges in exploiting graphics processors in computational fluid dynamics

Real-time 3D fluid simulation on GPU with complex obstacles

Practical CFD simulations on programmable graphics hardware using SMAC

Implementation of a semi-implicit pressure-based multigrid fluid flow algorithm on a graphics processing unit

CUDA implementation of a Navier-Stokes solver on multi-GPU desktop platforms for incompressible flows

An accelerated 3D Navier-Stokes solver for flows in turbomachines

A multi-GPU accelerated solver for the three-dimensional two-phase incompressible Navier-Stokes equations

Solving incompressible two-phase flows on multi-GPU clusters

Numerical solution of the two-phase incompressible Navier-Stokes equations using a GPU-accelerated meshless method

Direct numerical simulation of turbulent flow in a square duct using a graphics processing unit (GPU)

GPU accelerated flow solver for direct numerical simulation of turbulent flows

Direct numerical simulation of turbulence using GPU accelerated supercomputers

Accelerating multi-dimensional combustion simulations using GPU and hybrid explicit/implicit ODE integration

Accelerating S3D: a GPGPU case study

PyFR: An open source framework for solving advection-diffusion type problems on streaming architectures using the flux reconstruction approach

On the utility of GPU accelerated high-order methods for unsteady flow simulations: A comparison with industry-standard tools

A high-order cross-platform incompressible Navier-Stokes solver via artificial compressibility with application to a turbulent jet

Optimal Runge-Kutta schemes for pseudo time-stepping with high-order unstructured methods

Locally adaptive pseudo-time stepping for high-order flux reconstruction

Dnslab: A gateway to turbulent flow simulation in matlab

Scalable parallel programming with cuda: Is cuda the parallel programming model that application developers have been waiting for?

Electromagnetoencephalography software: overview and integration with other EEG/MEG toolboxes

The brian simulator

Passive and active battery balancing comparison based on MATLAB simulation

MATLAB based modeling to study the performance of different MPPT techniques used for solar PV system under various operating conditions

Modelling of electric and parallel-hybrid electric vehicle using Matlab/Simulink environment and planning of charging stations through a geographic information system and genetic algorithms

Model predictive control system design and implementation using MATLAB®

Contemporary communication systems using MATLAB, Cengage Learning

MIMO-OFDM wireless communications with MAT-LAB

Spectral methods: evolution to complex geometries and applications to fluid dynamics

On the implementation of low-dissipative Runge-Kutta projection methods for time dependent flows using OpenFOAM®

Dust diffusion in protoplanetary disks by magnetorotational turbulence

Large-eddy simulation of local heat transfer in plate and pin fin heat exchangers confined in a pipe flow

Assessment of CO2 and aerosol (pm2. 5, pm10, ufp) concentrations during the reopening of schools in the COVID-19 pandemic: The case of a metropolitan area in central-southern spain

CO2 concentration visualization for COVID-19 infection prevention in concert halls

Indoor air quality in naturally ventilated classrooms. lessons learned from a case study in a COVID-19 scenario

Recommendations for ventilation of indoor spaces to reduce COVID-19 transmission

High energy efficiency ventilation to limit COVID-19 contagion in school environments

Analysis of airflow through experimental rural buildings: Sensitivity to turbulence models

Natural ventilation for the prevention of airborne contagion

Natural ventilation for reducing airborne infection in hospitals

Association of the infection probability of COVID-19 with ventilation rates in confined spaces

The efficacy of social distance and ventilation effectiveness in preventing COVID-19 transmission

Estimating COVID-19 exposure in a classroom setting: A comparison between mathematical and numerical models

Natural ventilation strategy and related issues to prevent coronavirus disease 2019 (COVID-19) airborne transmission in a school building

Exhaled CO2 as a COVID-19 infection risk proxy for different indoor environments and activities

A coupled computational fluid dynamics and Wells-Riley model to predict COVID-19 infection probability for passengers on long-distance trains

Estimation of airborne viral emission: Quanta emission rate of SARS-CoV-2 for infection risk assessment

Quantitative assessment of the risk of airborne transmission of SARS-CoV-2 infection: prospective and retrospective applications

Ten scientific reasons in support of airborne transmission of SARS-CoV-2

Modality of human expired aerosol size distributions

We thank the Academy of Finland for their financial support (grant No. 335516) and the Aalto Science-IT project for the high-performance computational resources.