Minimizing the COVID-19 spread in hospitals through optimization of ventilation systems

The new COVID-19 virus pandemic has inflicted tremendous damage upon the world. In only one year since its emergence in December 2019, more than 100 million people have been diagnosed with the virus infection globally, which has led to more than 2 × 10⁶ deaths.¹ SARS-CoV-2 is a severe acute respiratory syndrome that causes COVID-19 and has a higher fatality rate in elderly patients and those with preexisting pulmonary conditions.² Although initially it was anticipated that the SARS-CoV-2 virus is transmitted through droplets from sneeze or cough (like SARS-CoV-1) and not by airborne particles and that social distancing may help break the infection chain,³ new studies suggest that the virus can indeed be transported through the air.⁴ The rapid spread of virus has overwhelmed hospitals with patients in need of intensive care, which is often limited in capacity and is generally reserved for patients with critical conditions. This has led to higher chances of infection spread to non-COVID-19 patients and healthcare workers and an overall increased probability of cross contamination. Studying the contamination spread in intensive care units is, hence, a critical step in avoiding and controlling infective diseases.^5,6

Sahu et al.⁷ studied the effects of airflow blowing angle on the stability of virus-laden particles in an ICU room in eight different cases (from 10° to 80°) numerically and found that when the airflow entered the room at angles between 30° and 50°, the time for air particles to exit through air conditioning was decreased and the particle residence time was increased for angles from 70° to 80°. Ching et al.⁸ investigated the effectiveness of hospital curtains on airborne transmitted infections numerically and experimentally and found that curtains between beds in fact reduce their spread and are most effective when they are fully open. Verma et al.⁵ used CFD (Computational Fluid Dynamics) techniques to study the particle distribution evolution in a hospital room with only one bed. The authors showed how the infection spreads and highlighted the contamination paths and concluded that higher air change rates (ACRs) cause faster removal of airborne infections and emphasized that the distance to outlet is a key parameter in controlling infections. Yuen et al.⁹ used a suction fan in an isolation room and found that airborne particles' lifetime can significantly be reduced due to the lower local mean age of air. Cho¹⁰ experimentally investigated contamination dissipation in an isolation room using various ventilation strategies and concluded that the location of supply air diffusers and outlet vents greatly affects infection spread caused by breathing in the room. The author suggests that installing the HVAC (Heating, Ventilation, and Air Conditioning) vents at lower heights results in less airborne particle dissipation.

There have been numerous studies about airborne transmitted diseases in operation and isolation rooms and processes in which the infecting agent travels from the contagious patient to others in hospitals.^11–18 Sun and Zhai presented a predicting model for airborne virus infection for indoor spaces by examining the effects of two prominent factors, viz., social distancing probability and ventilation effectiveness, and implementing them into the Wells–Riley model.¹⁹ Ge et al. investigated the exposure risk of COVID-19 virus in different rooms and sections of hospitals by collecting and testing six different air and surface samples from three hospitals for SARS-CoV-2. Their results showed that 82.6% of the samples taken from ICUs were positive and contained the DNA of SARS-CoV-2.²⁰ In another study, Jin et al. demonstrated that air samples taken from an ICU room with only one patient who was in the recovery period from COVID tested positive even days after the patient's COVID test was negative.²¹ These reports indicate the high risk of exposure to the virus in ICUs and the need for designing an optimal air conditioning system not only for high-risk areas of hospitals, but also future indoor environments.²² Verma and Sinha²³ investigated flow and temperature distributions in an ICU and showed that the stagnant regions between doctors and patients are potential contamination zones with increased chance of virus transmission. They also concluded that the arrangement and layout of patient beds are among the most critical parameters in infection control and should be considered carefully. Yu et al.²⁴ inspected several air conditioning schemes in an ICU with focus on the internal air quality (IAQ) and found that placing the air supply unit on top and exit vent at the bottom of rooms provides better IAQ performance. They also did not find any enhancement in the IAQ by adding more air supply units, and the number and size of inlets were found to depend on the dimensions of the room. The authors also conclude that increasing the air change per hour (ACH) and decreasing the inlet velocity have considerable effects on the performance of ventilation systems. Yang²⁵ studied the quality of breathing air in a four-bed hospital room with mixing and displacement ventilation methods and reported that the latter technique, which requires more air supply, provides better infection control and increases the quality of breathed air since it reduces the age of air in the room and, thus, should be the preferred method of ventilation in hospital rooms.

Although the body of work regarding the study of airborne infection spread is substantial in the literature, a majority of the investigations were carried out for isolated parameters, like particle tracking, stability and life, interactions with the flow field, performance of HVAC systems, and thermal comfort, separately and there is a lack of in-depth comprehensive studies with interactive parameters in particular, as the COVID-19 virus continues to wreak havoc on the world and the number of infected patients grows beyond hospital capacities and rapid makeshift hospitals are being built. With this motivation, the current study aims to perform a numerical modeling of infection control in a five-bed hospital room that includes the aforementioned parameters and their interactions to address challenges in inhibiting risks of airborne infection and to provide strategies that mitigate unwanted hazards associated with ventilation systems. To achieve these objectives, air flow patterns are first studied, and virus-laden particles are tracked in the presence of four different ventilation strategies. Once the best system with minimal chance of infection spread is determined, an optimization process will be carried out using the response surface method (RSM) to decide on parameters with the most influence on observables such as thermal comfort, particle residence time, and air change effectiveness (ACE). The studied parameters are the air supply channel width, inlet velocity, and temperature and the results are presented in the form of four different relations that are commonly used in HVAC analyses: predicted mean vote (PMV), percentage of people dissatisfied (PPD), residence time (RT), and ACE.

The geometry of the computational domain of the ICU room and the layout of the five patient beds in this study are shown in Fig. 1. There are two inlet and two outlet vents, which are subject to change. The dimensions of the room as well as other geometrical parameters are obtained from the work of Sahu et al.⁷ and are shown in Table I. The feasibility study is carried out using Reynolds-averaged Navier–Stokes (RANS) approach with the k-ε model for turbulence transport (k is the turbulent kinetic energy and ε is the turbulent dissipation); even though this is considered a second order closure model and is less accurate than higher fidelity models like Large Eddy Simulation, it requires much less computational resources and provides adequate accuracy for the current study following the validation process. ANSYS FLUENT, a finite volume computational fluid dynamics software, is used to solve the RANS-k-ε equations. The semi-implicit pressure linked equation (SIMPLE) scheme with second order upwind spatial discretization and first order implicit time discretization for transient simulations was used to solve the Navier–Stokes equations. Particles are treated as inert with one-way interaction, which means the only effective forces influencing the trajectory of particles are the primary fluid (air) momentum and the gravity. Particle injections are performed for each patient separately. Air is injected at 0.15 m/s from vents and particles are injected from the mouth with a velocity of 0.15 m/s. The diameter of the virus is 0.001 25 mm.

The mass and momentum conservation equations for an incompressible, turbulent, and steady-state flow can be written as follows:²⁸ 

where ρ is the fluid density, u is the fluid velocity, x_i is the position, t is the time, p is the pressure, and μ is the viscosity. Since an accurate flow modeling requires a valid turbulence model, the standard k-epsilon model is utilized in this study since it has been proven to be quite effective in ventilation or low-strained flow fields. The k-epsilon equations are given as follows:²⁹ 

where μ_t is the eddy viscosity; G_k is the turbulent kinetic energy generation term; σ_k is the kinetic energy Prandtl number; σ_ε is the dissipation Prandtl number; and C_m, C_1ε, C_2ε are the model constants. Next, the movement of the Lagrangian particles is calculated by Newton's kinematic equation as follows:

where F_drag, F_saff, F_pg, and F_grav are the drag force, Saffman's lift force, pressure gradient force, gravitational force, virtual mass force, basset force, and buoyancy force, respectively. Due to the small size of the particles, the drag force has the highest influence among all the forces applied on particles by the flow.³⁰ The drag force can be calculated as follows:

where C_D is the drag coefficient for spherical particles calculated by using the correlations developed by Schiller–Newman over several ranges of particle Re. Also, the Saffman's lift force can be obtained as follows:³¹ 

where K = 2.594 and S_ij is the deformation tensor according to Li and Ahmadi³² Furthermore, the pressure gradient force F_pg is defined as follows:³³ 

Finally, the gravitational force is calculated based on the gravitational acceleration (g), as follows:

For the particle boundary conditions, the reflection condition was applied so the particles would remain active inside the domain and not just trapped by the walls as soon as they have been released. By doing so, the worst-case scenario would be considered for these particles.

The optimization process requires more explanation. First, various test cases must be executed, and the results must be analyzed prior to the optimization process.³⁴ Second, the selected parameters used for optimization must be independent of each other,³⁵ otherwise their effects cannot be studied separately. Third, it is recommended that the design points are not directly specified, and to better scatter the design points, rigorous statistical procedures such as central composite design (CCD) or Taguchi must be employed.³⁶ This is important since the DOE framework involves expanding a finite number of experiments to considerable numbers. Even though the user can select his test experiments, the considered conditions might not be qualified to represent all the possible aspects of the flow. Thus, the overall sensitivity would not be reliable.³⁷ 

In this study, the CCD method is utilized since it is considered more complex with superior performance compared to the other methods such as Taguchi or Box–Behnken design (BBD). The parameters utilized in the DOE are defined and presented in Table II.

Prior to proceeding with a comparative study, the independence of results with regard to the grid size should be verified. For this purpose, five different mesh grids with elements between 2.5 × 10⁶ and 9.1 × 10⁶ are tested with an air inlet velocity of 0.4 m/s and temperature of 293 K and the velocity profiles along the length of the room are plotted. It can be seen from Fig. 2 that increasing the number of elements beyond the grid with 6.6 × 10⁶ elements (also known as normal) does not affect the velocity profile and, hence, this grid is used for the rest of the computation since the results are independent of the grid size.

Another important step in using any numerical scheme is checking its validity; for this purpose, an experimental study of contamination spread inside an isolation room in a hospital—which was performed by Cho¹⁰—is simulated and the results are compared to both the experimental and numerical data in that study. The author used sulfur hexafluoride (SF₆) as the tracer gas in a steady-state air flow inside a 4 × 4 × 2.6 m room. SF₆ gas was injected with an exhalation velocity of 0.955 m/s, mass fraction of 0.04, and temperature of 30 °C from a source located 0.9 m above the floor. The tracer concentrations are then measured at three different locations and compared with the results obtained from CFD. A comparison between the reported data and the results obtained in this study are shown in Table III for validation purposes.

SF₆ concentrations (ppm) are measured and reported Ref. 10 at three different points 1.4 m from the ground and are shown in Table III along with the numerical results obtained in this study. The average discrepancy between the experimental results and numerical modeling is about 6.5% for all three cases, which establishes the viability of using numerical schemes for further investigations in the current study.

Since the viability of the grid and the numerical model used in the current study has been established and a level of confidence in the predictive capability of both is provided, the dispersion of airborne particles in an ICU room with five beds is simulated for four different ventilation systems, i.e., the benchmark and three different modifications (Fig. 3) and upon determining the case with the least spread, the parameters affecting infection control are studied. Subsequently, to further enhance the performance of ventilation systems with respect to conventional comfort metrics, an optimization process is carried out.

As mentioned before, the geometry of cases in the current study with the exception of ventilation vent locations is obtained from the work of Sahu et al.⁷ and, thus, for comparison the first modeled case has the exact geometry in that study and is considered the benchmark case. Two rectangular inlets and two rectangular outlets are placed in opposite walls as shown in Fig. 3. The location of vents is changed afterward: In the second case, two outlets are moved to the wall next to the two inlets, which remain unchanged. In the third case, two inlets are moved to the ceiling and the locations of the two outlets are same as the benchmark case (case I). Finally, in the last case, air is coming from one larger aspect ratio rectangular inlet in the ceiling and five outlet vents are placed on opposite side walls behind each bed.

In all the cases, a steady-state airflow inside the room is established prior to particle injections. The pathlines of airflow for each case are shown in Fig. 4.

It can be seen from the flow pathlines that the ventilation system introduces a bulk flow with noticable eddy size insdie the room. To quantify the size of the largest eddies produced by the HVAC system, the turbulent length scale is calculated according to Ref. 26 as follows:

in which u' is the turbluence intensity and ε is the turbluence dissipation rate. Contours of the turbulence length scale along with streamlines are plotted in the planes of breathing for both sides of the rooms in Fig. 5.

The significance of the size of the largest eddies in the room manifests itself in the “entrapment” of airborne particles, which in turn leads to increased local particle residence time and consequently higher chance of transmission. As shown in Fig. 5, the integral length scale, which is on the order of the largest eddy, has a maximum of roughly 3 m for all cases; however, the largest eddy is furthest away from patient beds in mod III whereas in the other cases, they are located in the proximity of patients. The location and size of eddies smaller than the integral scale are also important; as can be seen from streamlines, in all cases except mod III, there are significant vortices in the bulk-flow in the domain, which essentially shows how particles will travel from one patient inside the room.

To better understand the evolution of particle pathways in time, injections are done for each patient separately for all cases. Figure 6 shows these injections for all the cases. It is evident that the released particles experience a completely different pathway for each patient for a given airflow despite having the same injection mass flow rate. For instance, the maximum and minimum pathway lengths for particles from injection to the outlet vent in the benchmark case are 6.35 m for patient b and 2.13 for patient c, respectively. For case IV, these numbers are 0.9 m and 0.95 m respectively. Longer pathways correspond to larger eddy sizes which in turn means further traveling distance of particles trapped in the eddy prior to the exit from the room.

Although it is evident from Figs. 5 and 6 that the ventilation system design affects infection spread significantly, the results should be quantified for accurate comparison. For this purpose, various parameters are defined for particles as well as the ventilation system. The first parameter is the dimensionless timescale, T, which is defined as

where $τparticle$ is the particle residence time (s) and $τroom$ is the bulk-flow residence time (s) inside the room (which is the ratio of the room volume to the total volume flow rate, i.e., $τroom=VV̇$ = 548.1 s). Three values of the minimum, maximum, and mean dimensionless timescale of particles averaged for all patients in each case are calculated and plotted in Fig. 7.

The larger the value of T, the longer it takes for the particles to exit the room and consequently the higher the chance of airborne exposure for the staff and non-infected patients. It is evident from Fig. 7 that a poorly designed ventilation system, the benchmark case for instance, can prolong the particle residence time by an order of magnitude. A smaller value of T, in contrast, minimizes the risk of contamination spread. Mod III, for example, has a mean T of 0.01, meaning that the average time a particle spends in the room is less than 6 s [compared to ∼6000 (s) in the benchmark case]. Another point worth noting is the large residence time distribution, which shows the degree to which contamination spread is dependent on the location. Ideally, there should not be any correlation between the location of infected patients and the level of spread of the infectious agent and this occurs only in mod III since the difference between the maximum and minimum T is negligible. For other cases, however, there is at least one order of magnitude difference in T, which implies a high level of dependency of spread to location.

Another parameter is the extracted mass of injected particles (m_e) normalized by the total injection mass (M), i.e.,

This normalized extracted mass of particles (m_n) is 0 at the time of injection and will be 1 when the mass is completely extracted from room. m_n is measured over time for each patient and then averaged for each case and the results are plotted in Fig. 8.

It can be seen from Fig. 8 that m_n for all cases goes to 1 ultimately, which means that all the injected mass exits the room eventually. However, the time over which maximum extraction occurs determines the effectiveness of the ventilation system. For the benchmark case, it takes about 6.7 times the bulk-flow residence time (3600 s) for the ventilation system to empty the room from particles. This would be done in less than 10 s (less than 0.02 of $τroom$⁠) in mod III when the air enters from the ceiling and exits from both sides of the room. Figure 8 shows that mod I, in which the air comes from the ceiling again, performs better than the benchmark case and mod II, but is still not as good as mod III. Results from Fig. 8 could be used to define another timescale associated with the extraction of injected material for each ventilation system. By applying a simple mass balance, In + Generation − Out = Accumulation, on the particles after injection, we have

in which c is the concentration, $v̇$ is the volumetric flow rate, and V is the total volume of the room. Equation (7) is an ordinary differential equation, which can be solved through separation of variables and integration as follows:

Therefore, if the logarithm of the effluent particle concentration normalized by the initial concentration (ln(c/c_out)) is plotted with respect to time, the extraction timescale (⁠$τ̃$⁠) is easily obtained by taking the slope of the plot,

To find this time for each case, graphs of Ln(m_e/M) averaged for all patients are plotted and the slopes are calculated using linear fitting as shown in Fig. 9. It should be noted that when $τ̃$ = $τroom$⁠, Eq. (16) shows the performance of a perfectly stirred reactor (PSR); so for comparison, the PSR predictions are also plotted.

It is evident from Fig. 9 that the extraction time for each case is −1/(dy/dx) of the fitted line. Therefore, $τ̃$ is 1250, 769.2, 555.5, and 7.2 s for the benchmark (case 1), mod II, mod I, and mod III respectively. This further proves that the ventilation system in mod III outperforms all the other cases as it can completely extract the particles up to four orders of magnitude faster than the other cases. Hence, this is proposed as the most effective ventilation system to control infection spread inside the patient rooms. By comparing the results with the PSR predictions, it can be concluded that except mod III, all modeled cases act similar to a PSR, which as suggested by its name and by definition guarantees perfect homogeneity of the mixture inside a control volume, which can be catastrophic when it comes to mixing of a contagious agent such as SARS-CoV-2 virus in a room. In fact, in places like ICUs or isolation rooms, the mixing characteristic of air conditioning systems should be far from perfect and the only case in which this is achieved is mod III.

After it is stablished that in the considered scenarios, mod III—in which air is coming from the ceiling and the outlet vents are placed on both walls behind the patient beds—performs the best, the effects of three different parameters, viz., inlet channel width (W), inlet air velocity (V), and inlet air temperature (T), on the system performance are studied in a systematic way and an optimization scheme is carried out to enhance the performance even further. The optimization scheme of the RSM is used and 20 design points are considered using the central composite design method, and three different response surfaces are created: 1—predicted mean vote (PMV), 2—percentage of people dissatisfied (PPD), and 3—air change effectiveness (ACE). These are standard parameters introduced by ASHRAE and are used to determine the thermal comfort of ventilation systems. PMV is calculated as follows:

where M is the rate of metabolic generation and L is the thermal load, which is calculated based on numerous parameters such as metabolic heat loss, surface temperature of the body, and effective thermal resistance of clothing material.²⁷ 

The PPD is then calculated using PMV from Ref. 27 as follows:

For the optimization process, the width is varied from 0.048 to 0.15 m, the velocity is varied from 0.22 to 0.58 m/s, and temperature is varied from 294.1 to 297.1 K. Contour maps of PMV and PPD are shown in Fig. 10.

Effects of varying the mentioned parameters on comfort conditions are visible in Fig. 10. According to ASHRAE, a PMV of zero represents thermal neutrality and changes between −0.5 and 0.5 are considered within the comfort zone. The comfort zone would normally have a PPD of less than 20%, which means the occupant satisfaction rate is more than 80% for a given thermal condition.²⁷ For a constant channel width, as the velocity is increased the PMV is decreased to larger negative values, i.e., less comfortable conditions. Likewise, as the velocity is increased for a constant width, the PPD increases. It should be noted even for small values of velocity, PMV and PPD exit the comfort zone conditions as the width is increased beyond certain values. Increasing the temperature for a constant channel width on the other hand increases the PMV from larger negative numbers, for example, −2 at T = 294.16 K toward the confront zone limit, −0.5 at 297.14 K. Similarly, the PPD gets closer to the comfort conditions as the temperature is increased at a constant width. Finally, both PMV and PPD approach the comfort limit as the velocity is decreased at constant temperature. It is worth noting that the iso-values of PMV and PPD in Fig. 10 show that for all three parameters, a certain value can be obtained through various choices of velocity, channel width, and temperature. In other words, there is a range of parameters in which thermal conditions are considered desirable.

The ACE is a measure of air displacement within an enclosed environment. It is the ratio of air residence time to the average age of air at a reference height (typically breathing height); according to ASHRAE criteria for modern buildings, it should be larger than 0.95. An ACE of unity represents perfect mixing of the air distribution; hence, to avoid local stagnation of air ACE should be close to unity across the room. Since ACE is defined based on the ventilation system's residence time, it can be rewritten using T (dimensionless timescale) in Eq. (13) as follows:

Here, $τavg$ is the average age of air and is calculated utilizing a user-defined scalar inside the FLUENT software.

Like PMV and PPD, for a constant width, as the velocity decreases the ACE approaches the comfort condition. Using the contours in Fig. 10, three different quadratic fittings are proposed for calculating PMV, PPD, and ACE as a function of width, temperature, and velocity as follows:

Using the above equations, we eliminate all the other variables that may not be readily available during the design process and make the calculations easy.

As mentioned before, the response surfaces in Fig. 10 show that thermal comfort condition is obtained within a range of values for the three studied parameters, hence a desirability function is defined based on PMV, PPD, and ACE criterion that shows the limits for variable selection. The ideal thermal conditions correspond to the desirability of 1, and by plotting the response surfaces one can easily select the parameters that provide comfortable conditions. The desirability surfaces are plotted in Fig. 11.

The maximum desirability obtained in the current study is 0.865 as shown in Fig. 11. This value can be obtained through various design selections. For example, for a velocity of 0.22 m/s and temperature of 297.137 K any value of width between 0.104 and 0.108 m results in a desirability of 0.865. To show the performance of the ventilation system and the distribution of thermal criterion, another set of computation is performed with fixed values of channel width, air inlet velocity, and temperature and the contours of PMV, PPD, ACE, and age of air are plotted in Fig. 12.

Even though the predicted values for standard thermal comfort parameters fall in the suggested range, the distribution of these parameters are also significant as the ideal condition must have a uniform distribution. Figure 12 shows that not only the predicted values for PMV, PPD, and ACE are close to ideal conditions, the distribution of parameters specifically in proximity of patient beds is in fact uniform, which further proves the ability of the proposed design to provide thermal comfort.

The effects of design parameters on the performance of ventilation systems to control the spread of airborne particles in an ICU are studied numerically. Four different cases with various locations of inlet and outlet vents are considered and the spread of particles is studied. Two new criteria for the ventilation system are defined—viz., dimensionless timescale (T) and extraction timescale $τ̃$—and the performance of all the cases are compared. The results show that when the air inlets are placed in the ceiling and outlets are placed behind the patient beds (mod III), the infection spread is least probable since the particles exit the room 2 orders of magnitude faster than the other cases. Once it is established that this system provides a significant enhancement in infection control, an optimization process is performed to understand the effects of design variables (air inlet width, inlet air velocity, and inlet air temperature) on the thermal comfort conditions (PMV, PPD, ACE) according to the suggested standard values, and the relations for calculating each of these parameters based on the design variables are proposed. Desirability functions that are comprised of all three thermal condition parameters are used to determine the range of variables that result in thermally comfortable conditions and a maximum desirability of 0.865 is obtained.

The authors declare no potential conflict of interest of any kind throughout conducting of present study.

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