UNIVERSIDADE FEDERAL DE SÃO CARLOS CENTRO DE CIÊNCIAS EXATAS E DE TECNOLOGIA PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA QUÍMICA NUHU AYUBA Application of CFD to Study the Dynamics of Droplets Expelled During Coughing Aplicação de CFD no Estudo da Dinâmica de Gotículas Expelidas Durante a Tosse SÃO CARLOS - SP AUGUST, 2024 ii UNIVERSIDADE FEDERAL DE SÃO CARLOS CENTRO DE CIÊNCIAS EXATAS E DE TECNOLOGIA PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA QUÍMICA NUHU AYUBA Doctoral Thesis presented to the Postgraduate Program in Chemical Engineering at the Federal University of São Carlos as a requirement to obtain the degree of Doctor in Chemical Engineering. Supervisor: Gabriela Cantarelli Lopes Co supervisor: Gabriel Henrique Justi SÃO CARLOS – SP AUGUST, 2024 iii iv v ACKNOWLEDGMENT First and foremost, I would like to express my deepest gratitude to my supervisor, Dr. Gabriela Cantarelli Lopes, for her invaluable guidance, continuous support, and encouragement throughout the development of this study. Her expertise and mentorship were fundamental to the success of this research. I also extend my sincere thanks to my co-supervisor, Dr. Gabriel Henrique Justi, for his insightful suggestions and constructive feedback, which significantly enhanced the quality of this work. I am grateful to the members of the examination committee—Prof. Dr. Liliana de Luca Xavier Augusto, Prof. Dr. Marcela Kotsuka da Silva Câmara Bastos, and Prof. Dr. Vadila Giovana Guerra Bettega—for their valuable time, contributions, and critical evaluation of this dissertation. My heartfelt appreciation goes to my wife, Jady, and my daughter, Victória, for their unwavering support, patience, and understanding throughout this academic journey. I also extend special thanks to my brothers Minighom, Suleiman, and Seth in Nigeria, whose enduring encouragement and belief in my potential have always motivated me to strive for excellence in my academic pursuits. A special acknowledgment is due to Dr. Toni Jefferson Lopes for his immeasurable support before, during, and after the completion of this research. I am also particularly thankful to my friend Habila Thomas Yusuf, who played a pivotal role in assisting with my doctoral application to São Carlos. Your support has always gone beyond friendship. I would also like to thank Dr. Antonio Carlos Daltro de Freitas and Fernando Pedro Dias for their support at the onset of this research. My sincere thanks to the African community in the cities of São Luís and São Carlos for enriching my social experience during this period—indeed, all work and no play makes Jack a dull boy. vi Special thanks to my fellow doctoral colleagues Ricardo Arbach Fernandes de Oliveira, Victor Ferreira, and Karla Raphaela Braga de Melo for their technical assistance and support, especially in handling hardware and computer systems essential to this research. I am also grateful to my friends Anacleto Domingos, Osmilde Miranda, Joseph Osei, Ojubanire Babatunde, Osires Bideu, Francisco Malungo, Kiabanga Zamba, and Elvis Kasanga, whose kind words and encouragement helped me through difficult times. I would like to thank the coordination and staff of the Department of Chemical Engineering at UFSCar for providing the necessary infrastructure and a conducive environment for conducting this study. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. vii ABSTRACT This study investigates the behavior of respiratory droplets in indoor environments, applying the Eulerian-Lagrangian method to improve understanding of droplet transmission in the context of Infectious Respiratory Diseases (IRDs). The research is divided into four key sections. The first part evaluates the effects of different turbulence models—SST k-ω, standard k-ε, and Reynolds Stress—on the evaporation of a single droplet, aiming to identify a model that is both precise and computationally efficient. The second part compares the use of graded tetrahedral and hexahedral meshes in simulating droplet behavior, assessing their relative effectiveness and computational cost. In the third part, the influence of phase coupling, droplet size, and relative humidity on the evaporation process is examined, while the fourth part explores the impact of wind boundary conditions on droplet dispersion. To further validate the findings, the models developed were applied under the climatic conditions of São Carlos, São Paulo, Brazil, to observe how the droplets behave in a real-world setting. The results demonstrate that hexahedral meshes provide more accurate droplet behavior representation compared to tetrahedral meshes, which showed discrepancies with previous studies. The effect of phase coupling was found to be minimal, with smaller droplets evaporating rapidly irrespective of ambient humidity, while larger droplets fell to the ground. The study also reveals that higher humidity levels slow down evaporation, with medium-sized droplets (around 50 µm) persisting longer and traveling farther—up to approximately 1.40 meters. This suggests a minimum social distancing distance in indoor, quiescent environments to minimize the transmission of IRDs. Additionally, smaller droplets evaporated quickly, and their dispersion increased at lower wind speeds. The research concludes that the no-slip boundary condition is more cost-effective for indoor simulations. These findings contribute valuable insights into understanding droplet behavior and optimizing models for IRD transmission analysis in indoor spaces. Keywords: Eulerian-Lagrangian; Respiratory droplets; Tetrahedral and hexahedral mesh; Turbulence model; droplet evaporation; Relative humidity. viii RESUMO Este estudo investiga o comportamento das gotículas respiratórias em ambientes internos, aplicando o método Euleriano-Lagrangeano para melhorar a compreensão da transmissão de gotículas no contexto de Doenças Respiratórias Infecciosas (DRIs). A pesquisa está dividida em quatro seções principais. A primeira parte avalia os efeitos de diferentes modelos de turbulência—SST k-ω, k-ε padrão e Tensões de Reynolds—sobre a evaporação de uma única gotícula, com o objetivo de identificar um modelo que seja preciso e computacionalmente eficiente. A segunda parte compara o uso de malhas tetraédricas e hexaédricas graduadas na simulação do comportamento das gotículas, avaliando sua eficácia relativa e o custo computacional. Na terceira parte, é examinado o efeito do acoplamento de fase, o tamanho das gotículas e a umidade relativa sobre o processo de evaporação, enquanto a quarta parte explora o impacto das condições de contorno do vento na dispersão das gotículas. Para validar ainda mais os resultados, os modelos desenvolvidos foram aplicados nas condições climáticas de São Carlos, São Paulo, Brasil, para observar o comportamento das gotículas em um ambiente real. Os resultados demonstram que as malhas hexaédricas fornecem uma representação mais precisa do comportamento das gotículas em comparação com as malhas tetraédricas, que apresentaram discrepâncias com estudos anteriores. O efeito do acoplamento de fase foi encontrado como mínimo, com as gotículas menores evaporando rapidamente, independentemente da umidade ambiente, enquanto as gotículas maiores caíam ao chão. O estudo também revela que níveis mais altos de umidade desaceleram a evaporação, com gotículas de tamanho médio (cerca de 50 µm) persistindo por mais tempo e viajando mais longe—até aproximadamente 1,40 metros. Isso sugere uma distância mínima de distanciamento social em ambientes internos e tranquilos para minimizar a transmissão de DRIs. Além disso, gotículas menores evaporaram rapidamente e sua dispersão aumentou com velocidades de vento mais baixas. A pesquisa conclui que a condição de contorno de não-deslizamento é mais econômica para simulações em ambientes internos. Esses resultados contribuem com valiosas informações para entender o comportamento das gotículas e otimizar modelos para análise da transmissão de DRIs em espaços fechados. Palavras-chave: Euleriano-Lagrangiano; Gotículas respiratórias; Malha tetraédrica e hexaédrica; Modelo de turbulência; Evaporação de gotículas; Umidade relativa. ix LIST OF FIGURES Figure 1- Airborne transmission of IRDs. (Self-elaborated diagram with safe distance according to results of this work). .............................................................................................. 3 Figure 2- The direction of the injection of the respiratory droplets ........................................... 6 Figure 3- Disperse phase boundary conditions......................................................................... 45 Figure 4- Mesh levels ............................................................................................................... 40 Figure 5- Graded mesh with quadrilateral elements (2D) or hexahedral elements 3D .............. 9 Figure 6- Diagrams of ideal (a) tetrahedral and (b) hexahedral meshes. ................................. 12 Figure 7- Coupling between the continuous phase and discrete particle. ................................ 34 Figure 8- Cough velocity profile. ............................................................................................. 36 Figure 9- Relationship between vapor mass fraction of pure water and relative humidity at 293.15 K and 1 atm. ................................................................................................................. 37 Figure 10- Computational respiratory chamber.. ...................................................................... 39 Figure 11- Structured 2D quadrilateral grid ............................................................................. 41 Figure 12- Grid convergence test. ............................................................................................ 47 Figure 13- Droplet evaporation convergence test. .................................................................... 48 Figure 14- Respiratory droplet evaporation (a) 0% relative humidity (b) 80% relative humidity. ................................................................................................................................... 50 Figure 15- Evaporation of cough droplet particle using (SST) k-ω, k-ε and Reynolds Stress (RSM) turbulence models (a) 1 μm and (b) 100 μm. ............................................................... 51 Figure 16- Single drop evaporation under various humidity and droplet sizes. ....................... 53 Figure 17- Geometry of the respiratory chamber. .................................................................... 54 Figure 18- Projections of the isomeric, side, and cutout views of a tetrahedral mesh with (a) 40 thousand elements (b) 150 thousand elements (c) 600 thousand elements, and (d) 2.4 million elements. ...................................................................................................................... 56 Figure 19- Projections of the isomeric, side, and cutout views of a hexahedral mesh with (a) 40 thousand elements (b) 150 thousand elements (c) 600 thousand elements, and (d) 2.4 million elements. ...................................................................................................................... 57 Figure 20- Richardson Extrapolation for (a) Tetrahedral mesh (b) Hexahedral mesh (c) Tetrahedral and Hexahedral mesh at 0.1 s. ............................................................................... 62 Figure 21- The evaporation of respiratory droplets with different sizes under 0% relative humidity condition in a quiescent space; comparison with the first group of authors. ............ 65 x Figure 22- The evaporation of respiratory droplets of different sizes under 0% relative humidity condition in a quiescent space; comparison with the second group of authors. ....... 66 Figure 23- Visual representation of respiratory particles 10 µm using tetrahedral and hexahedral meshes at the specified interval of time under relative humidity conditions of (a) at 0% and (b) at 80%. ................................................................................................................... 67 Figure 24- Visual representation of respiratory particles 50 µm using tetrahedral and hexahedral meshes at the specified interval of time under relative humidity conditions of (a) at 0% and (b) at 80%. ................................................................................................................... 68 Figure 25- Visual representation of respiratory particles of 100 µm in diameter using tetrahedral and hexahedral meshes at the specified interval of time under a relative humidity condition of 0%. ....................................................................................................................... 69 Figure 26- Visual representation of respiratory particles of 100 µm in diameter using tetrahedral and hexahedral meshes at the specified interval of time under a relative humidity condition of 80%. ..................................................................................................................... 70 Figure 27- (a) Tetrahedral mesh and (b) hexahedral with the grading regions. ....................... 71 Figure 28- Respiratory droplets with a diameter of 100 µm laying on the bottom of the domain floor of a closed space with a relative humidity of 80% calculated by a hexahedral mesh. ......................................................................................................................................... 73 Figure 29- Cough chamber. ...................................................................................................... 76 Figure 30- Mesh with its boundaries. ....................................................................................... 77 Figure 31- Mean droplet diameter as a function of time for different meshes. ........................ 80 Figure 32- G3 aspect ratio distribution. .................................................................................... 82 Figure 33- Droplet evaporation of cough droplets in a quiescent room at 0% relative humidity. .................................................................................................................................................. 83 Figure 34- Models 1,2 and 3 using Cases 3 and 4 (a) horizontal direction (b) vertical direction (c) M2 Case 4; dispersed particle cloud formed at t = 0.10 s. ... Erro! Indicador não definido. Figure 35- Models 1,2 and 3 using Cases 5 and 6 (a) horizontal direction (b) vertical direction (c) M2 Case 6; dispersed particle cloud formed at t= 0.10 s. ................................................... 87 Figure 36- Models 1, 2 and 3 using Cases 7 and 8 (a) horizontal direction (b) vertical direction (c) M2 Case 8; compacted particle cloud formed at; Left t= 0.10 s and Right 10.00 s. ........... 89 Figure 37- Respiratory flow domain. ....................................................................................... 92 Figure 38- Transverse section of the numerical mesh built with hexahedral elements. ........... 93 Figure 39- Droplet cloud at different velocities. ...................................................................... 97 Figure 40- Droplet cloud at different velocities. ...................................................................... 98 xi Figure 41-YZ plane velocity profile under a wind velocity of 3 m/s (a) no-slip (b) zero shear stress and (c) pressure outlet. .................................................................................................... 99 Figure 42- Reynolds profile of wind at high velocity (4/m) (a) no-slip (b) zero-shear stress (c) pressure outlet. ........................................................................................................................ 100 Figure 43- 50 μm droplet behavior under (a) 0.5 m/s (b) 1 m/s and (c) 3 m/s. ...................... 103 xii LIST OF TABLES Table 1- Eulerian boundary conditions. .................................................................................... 43 Table 2- Lagrangian boundary conditions ................................................................................ 44 Table 3- Mesh quality measurement ......................................................................................... 11 Table 4- Some turbulence models applied in the CFD studies of respiratory droplets. ........... 14 Table 5- Mesh information. ...................................................................................................... 42 Table 6- Boundary condition of the Euler-Lagrangian model. ................................................. 45 Table 7- Grids applied in this study. ........................................................................................ 47 Table 8- Simulation conditions used by Li et al. (2020). ......................................................... 49 Table 9- Time taken by 100 μm droplet to evaporate at different injection velocities and 0% humidity. ................................................................................................................................... 52 Table 10- Details of (a) tetrahedral and (b) hexahedral meshes. .............................................. 55 Table 11- Boundary condition for mesh independence test and result calculations. ................ 58 Table 12- Respiratory particle diameter at a specific time (a) Tetrahedral mesh (b) hexahedral mesh .......................................................................................................................................... 60 Table 13- Richardson extrapolation for tetrahedral mesh. ....................................................... 61 Table 14- Richardson extrapolation for hexahedral mesh. ....................................................... 61 Table 15- GCI at using different time intervals. ....................................................................... 63 Table 16- asymptotic range of convergence check ................................................................... 64 Table 17- Droplet evaporation time estimated by the tetrahedral and hexahedral meshes respectively at different humidity conditions. .......................................................................... 72 Table 18- Total CPU time measured at t = 0.52 s for the simulation. ...................................... 74 Table 19- Models 1, 2, and 3 and their respective physical configurations. ............................ 78 Table 20- Simulation cases and their respective characteristics............................................... 78 Table 21- Simulation conditions. ............................................................................................. 79 Table 22- Mesh parameters. ...................................................................................................... 80 Table 23- Horizontal distance traveled by cough droplets. ...................................................... 84 Table 24- Details of the mesh generated using hexahedral elements. ...................................... 92 Table 25- Boundary condition of the Euler-Lagrangian model. ............................................... 94 Table 26. Wall conditions used in this research. ....................................................................... 95 Table 27- Details of computational time at different air velocity applied in this work. ......... 101 xiii LIST OF ABBREVIATIONS AND ACRONYMS 2D ..................................................................................................................... Two dimensional 3D ................................................................................................................... Three dimensional BC ................................................................................................................. Boundary condition CFL ....................................................................................................... Courant-Fiedrichs-Lewy DNAs ....................................................................................................... Deoxyribonucleic acid DRW ............................................................................................. Discrete Random Walk model FSM ......................................................................................................... Fractional Step method GCI ....................................................................................................... Grid Convergence Index IRDs ............................................................................................. Infectious respiratory diseases KHRT............................................................................ Kelvin-Helmholtz, and Rayleigh-Taylor LES .......................................................................................................... Large Eddy Simulation PISO ..................................................................... Pressure-Implicit with Splitting of Operators RANS .................................................................................... Reynolds-averaged Navier-Stokes RNG k-ε ................................................................................................ Re-Normalisation Group RSM ........................................................................................................ Reynolds-stress models SARS ................................................................................... Severe acute respiratory syndromes SIMPLE .................................................. Semi-Implicit Method for Pressure-Linked Equations SIMPLEC .................................................................................................... SIMPLE-Consistent SST k-ω .................................................................................................... Shear Stress Transport TAB ...................................................................................................... Taylor Analogy Breakup xiv TB ............................................................................................................................ Tuberculosis UDF ......................................................................................................... User Defined Function xv LIST OF SYMBOLS Latin symbols Ad ..................................................................................................... Surface area of droplet [m2] Cw,∞ .................................................................... Water concentration in the bulk gas [kmol/m3] Cw,d ........................................................... Water concentration at the droplet surface [kmol/m3] C1ε ........................................................................... Production term of the the energy equation C2 ................................................................................................................................... Constant C2ε ................................................................................. Dissipation term of the epsilon equation Cc ....................................................................... Cunningham correction factor (dimensionless) cd ....................................................................................... Heat capacity of the droplet [J/kg K] CD ........................................................................................................................ Drag coefficient Cij ..Convective transport term. Represents the change in Reynolds stress due to the mean flow [kg/(m·s³] Cε .................................................................................................................................... Constant Cμ ..................................................................................................... Turbulent viscosity constant Disk ..................................................................................................................... Dissipation of k DL,ij ....................................................................................... Molecular diffusion term [kg/m·s³] DT Turbulent diffusion term [kg/m·s³] Dv ....................................................................................... Diffusion coefficient of vapor [m2/s] Dω ...................................................................... Cross-diffusion term, Diffusive term of omega dp .................................................................................................... Damping coefficient [N·s/m] Echild ................................................................................................ Energy of the child droplets en .......................................................................................................... Coefficient of restitution EParent..................................................................................... Energy equation of parent droplet F⃗ . .................................................................................................................... Force vector [m/s²] F1 ..................................................................................................................... Blending function Fij .................................................................................... Additional source/sink terms [kg/m·s³] Fother ............................................................Other external forces acting on the mass [kg/m·s³] fβ ............................................................................................................................. Function of B G̃k ................................... Production of turbulent kinetic energy due to mean velocity gradients Gb ....................................................... Generation of turbulence kinetic energy due to buoyancy xvi Gij ....................................................................................... Buoyancy production term [kg/m·s³] Gk ........................... generation of turbulence kinetic energy due to the mean velocity gradients Gω ....................................................................... Production of omega due to velocity gradients h ........................................................................................ coefficient of heat transfer [W m2/ K] hd,g .......................................................................................... latent heat of vaporization (kJ/kg) hfg ................................................................................................................... Latent heat [kJ/kg] hg .............................................................................................................. Enthalpy of air [kJ/kg] k ....................................................................... Turbulent kinetic energy, Spring constant [N/m] K ............................................................................................................ Proportionality constant kc ................................................................................................. Mass transfer coefficient [m/s] L3 .......................................................................................................... Cell length at level 3 [m] Le ............................................................................................................... Eddy length scale [m] L2 ................................................................................................................. Cell length at level 2 m ............................................................................................................ Mass of the particle [kg] ṁ ........................................................................................................................ Mass flow [kg/s] ṁd→g. ..................................................................................................... Mass transfer rate [kg/s] ma .......................................................................................................... Component source term Mt0 ................................................................................................ Initial turbulent Mach number MWw ............................................................................. Molecular weight of species [kg/ kmol] n ..................................................................................................... Total number of cells applied nd ............................................................................................................... Density number [m-3] P ............................................................................................................................... Pressure [Pa] Pg .......................................................................................................... Air pressire Pressure [Pa] Pij ........................................................................................................ Production term [kg/m·s³] Poutlet ............................................................................................................ Outlet pressure [Pa] Psat ................................................................................................. Saturated vapor pressure [Pa] Pstatic .............................................................................................................. Static Pressure [Pa] Ptotal ................................................................................................................ Total pressure [Pa] r ......................................................................................................................... Mesh growth rate R ........................................................................................ Universal gas constant [Pa L/mol K] Rβ ......................................................................................... Reynolds number for quantity beta Rideal ................................................................................................................... Ideal radius [m] xvii Rk ............................................................................................... Reynolds number for k quantity Rβ .................................................................................................................................... Constant Rω ............................................................................................. Reynolds number for ω quantity Re ....................................................................................................... Relative Reynolds number S .......................................................... Mean rate-of-strain tensor, Shape factor (dimensionless) Sc ....................................................................................................................... Schmidt number Sh .................................................................................................................... Sherwood number t ....................................................................................................................................... Time [s] T ................................................................................................................. Integral time scale [s] T∞ .................................................................................. Temperature of the surrounding air [K] Tb ........................................................................................................... Boiling temperature [K] tcross ................................................................... Time needed for a particle to cross an eddy [s] Td ........................................................................................ Respiratory droplet temperature [K] tinteraction ...................................................................................................... Interaction time [s] TL ............................................................................................. fluid Lagrangian integral time [s] Tv .............................................................................................................. Vapor temperature [K] v ............................................................................................................................... Velocity [m] v' ................... Fluctuating component of the continuous phase velocity due to turbulence [m/s] v̅ ............................................................................ Mean velocity of the continuous phase [m/s] v⃗ d ........................................................................................ Velocity of respiratory droples [m/s] v⃗ g ....................................................................................... Velocity of the dispersed phase [m/s] v1,n ........................................................................................ of the particle befor collision [m/s] v2,n ........................................................................... Velocity of the particle after collision [m/s] vd .......................................................................................... Respiratory droplets velocity [m/s] Vi ....................................................................................................................... Cell volume [m3] Videal ................................................................................................................ Ideal volume [m3] vinlet ............................................................................................................... Inlet velocity [m/s] Vreal ........................................................................................................................... Real volume vrel ........................................................................................................... Relative velocity [m/s] We ......................................................................................................................... Weber number x ....................................................................... Displacement of the mass from equilibrium [m] Ya,vap ................................................................................. The respective mass fraction of vapor Yk ......................................................................................................... Turbulent dissipation of k xviii Yω ........................................................................................................ Turbulent dissipation of ω Greek symbols α0 ............................................................................. Constant related to turbulent kinetic energy α∞ ........... Constant related to the rate of production and dissipation of turbulent kinetic energy α∞ * ........... Constant related to the rate of production and dissipation of turbulent kinetic energy β* ........................................................................................................Turbulence model constant β∞ * ................................................................................................................................... Constant Γk,ω .............................................................................................. Effective diffusivity of k and ω ε ............................................................................................................................ dissipation rate εij ....................................................................................................... Dissipation term [kg/m·s³] ζ ................................... Normally distributed random number with zero mean and unit variance λ ............................................................................................ Cough droplets mean free path [m] μg ...................................................................... The molecular viscosity of the fluid [kg m-1 s-2] μt ................................................................................................. Turbulent viscosity [kg m-1.s-1] ξ* .......................................................................................................... Constant, Model constant σ .................................................................................................... Particle surface tension [N/m] ρa,vap ......................................................................................................... Vapor density [kg m-3] ρg ............................................................................. Density of the respiratory droplets [kg m-3] σk ........................................................................................................... Turbulente Prandtl for k σε ............................................................................................................. Turbulent Prandtl for ε τ ......................................................................................................... Particle relaxation time [s] τ̿g ........................................................................................................................... Viscous stress τe........................................................................ Characteristic lifetime Eddy [s], Time scale [s] φij ...............................................................................Pressure-strain correlation term [kg/m·s³] xix SUMMARY 1. INTRODUCTION ................................................................................................................ 1 2. LITERATURE REVIEW ..................................................................................................... 3 2.1 Airborne transmission ................................................................................................... 3 2.2 CFD studies of respiratory events ................................................................................ 4 2.3 Boundary characterization ............................................................................................ 5 2.4 Droplet size and distribution ......................................................................................... 6 2.5 Computational domain and meshing ............................................................................ 7 2.6 Discrete Particle Model (DPM) .................................................................................. 13 2.7 Turbulence model ....................................................................................................... 13 3. STUDY OBJECTIVE ......................................................................................................... 16 4. MATHEMATICAL MODEL ............................................................................................ 17 4.1 Eulerian model ............................................................................................................ 17 4.2 Lagrangian formulation .............................................................................................. 25 4.3 Coupling ..................................................................................................................... 33 5. SIMULATION .................................................................................................................... 35 5.1 Work organization ...................................................................................................... 35 5.2 Processing capacity ..................................................................................................... 35 5.3 Common to all parts of the studies ............................................................................. 35 5.3.1 Parameters and boundary conditions .......................................................................... 36 5.3.2 Simulation configuration ............................................................................................ 37 5.3.3 Run interruption .......................................................................................................... 38 6. PART 1 – TURBULENCE MODEL STUDIES ............................................................. 39 6.6.1 Mesh independence test .............................................................................................. 46 6.6.2 Droplet independence test .......................................................................................... 48 6.6.3 Validation 49 6.6.4 Droplet evaporation .................................................................................................... 50 6.6.5 Turbulence model comparison ................................................................................... 51 7. PART 2 – MESH STUDIES .............................................................................................. 54 7.7.1 Particle cloud .............................................................................................................. 66 xx 7.7.2 Particle evaporation/trapping time .............................................................................. 72 7.7.3 Total cpu time ............................................................................................................. 74 8. PART 3 - INFLUENCE PHASE COUPLING, DROPLET DIAMETER AND RELATIVE HUMIDITY ON HUMAN COUGH DROPLET EVAPORATION ........................... 75 8.1 Computational domain & mesh .................................................................................. 75 8.2 Numerical process ...................................................................................................... 77 8.2.1 Simulated cases ........................................................................................................... 77 8.2.2 Mesh independence test .............................................................................................. 79 8.2.3 Droplet evaporation .................................................................................................... 80 8.2.4 Mesh choice ................................................................................................................ 81 8.2.5 Validation of the results .............................................................................................. 82 8.2.6 Model, humidity, and diameter influence ................................................................... 83 9. PART 4 - APPLICATION ................................................................................................. 91 10. CONCLUSIONS ............................................................................................................... 107 1 1. INTRODUCTION Infectious respiratory diseases (IRDs) such as severe acute respiratory syndromes (SARS), tuberculosis (TB), influenzas (avian and swine) have been interest of study for a long time (Gupta et al., 2009 and Morawska, 2006). The rapid spread of SARS-CoV-2 virus, better known as the new Corona Virus disease (COVID-19) in the early 2020 (Kucharski et al., 2020) has intensified the need to elaborate more studies about IRD. Currently, IRD studies are not just limited to the medical fields, but have also extended to others. For example, in Engineering, pathogens may not directly be studied. Instead, parameters that influence transportation of contaminated respiratory droplets are studied. Good examples of parameters that influence the spread of IRDs are the particle size and velocity. For instance, large droplets plunge to the ground when they are released, on the contrary, smaller droplets tend to travel larger horizontal distances. Additionally, climatic conditions such as relative humidity, temperature and wind velocity and direction contribute greatly to the spread of IRDs. In this study, sometimes, simplifying the flow domains and the physical properties of the respiratory droplets are needed to reduce computational cost, yet providing results that are similar to real events. For instance, geometries used to carry out the flow studies are made in such a way to permit the application of coarser meshes. This is done by avoiding the use of curves to build the flow domains. To simplify the respiratory droplets, contents such as mass composition of virus (Branche et al., 2014; Honkinen et al., 2012; Wolfel et al., 2020) proteins, lipids, carbohydrates, DNAs and salts (Redrow et al., 2011) are neglected. Hence, respiratory droplets with the physical properties of pure water are applied. To execute the CFD studies of droplets emitted during respiratory events, the application of the appropriate numerical models is important. Considerations in the transient formulation of the flow governing equations, turbulent models, multiphase flow formulation which depends on the number of phases applied, droplet evaporation, particle break (where applicable), droplet dispersion due to turbulence, droplet evaporation and boundary conditions should be made. These ample combination makes this study open for new findings. 2 This analysis is conducted to contribute to the existent results available in the literature. To this end, the CFD of respiratory droplets is studied using the ANSYS Fluent commercial software adopting the Eulerian-Lagrangian model. It is divided in four parts. In the first part, the objective is to study turbulence effect on the evaporation of a single respiratory droplet. The (SST) k-ω, k-ε and Reynolds Stress (RSM) turbulence models are compared applying a 2D simulation. In the second part of the research, the objective is to build a 3D mesh, comparing two different types: the trihedral and hexahedral. In the third part of the research, the objective is to study the influence of air relative humidity on human cough droplet evaporation with the aim to use a 3D tetrahedral mesh to evaluate the evolution of the droplet at different relative humidity conditions. The last part has the quest to apply the tests carried out in previous parts of the research with the objective to study the behavior of the respiratory droplets under the climatic conditions of São Carlos, São Paulo, Brazil. In addition to, a complimentary study on behavior of the flow at different boundary conditions is carried out in this part of the research. In summary, this study is expected to contribute to the understanding of the spread of viruses and bacteria carried by droplets, enabling the definition of standards and recommendations to reduce the transmission process of diseases caused by this transmission mechanism. Additionally, this study is also expected to aid in choice and application of the numerical models and meshes that can be applied in the CFD studies of respiratory droplets in closed and open environments. 3 2. LITERATURE REVIEW 2.1 Airborne transmission Airborne transmission is a primary route for the spread of IRDs, where pathogens are expelled from an infected person through coughing, sneezing, talking, or breathing. They are expelled with the respiratory droplets. The larger droplets fall to the ground, medium droplets travel larger horizontal distances and small droplets evaporate instantly when they are released. Even though the liquid part of the droplets may evaporate totally, parts remain suspended in the air for extended periods in form of aerosols (Morawska and Cao, 2020) increasing the risk of inhalation by susceptible individuals. The image in Figure 1 shows how the droplets are classified when they are released into the surrounding. Figure 1- Airborne transmission of IRDs Pessoa 1 Pessoa 2 Sick person Healthy person Aerosols persist in the air Larger droplets Evaporation of smaller droplets (Self-elaborated diagram with safe distance according to results of this work) Apart from airborne transmission, contaminations can take place through fomites (Van Doremalen et al., 2020) and high-touch surfaces (Ong et al., 2020) which are surface contaminations. Contaminations can also take place through direct contact (Kwok et al., 2015; COVID, C., 2020). 4 2.2 CFD studies of respiratory events CFD has become an important tool for studying the behavior of the droplets (Sedighi et al., 2023). One of the reasons for its regular application by the scientific community is its lower cost when compared to experimental methods (Chen et al., 2022 and Cao et al. 2022). The CFD studies of respiratory events have been carried out in two forms. The first involves the study of the behavior of the respiratory parameters within the human respiratory system. Classical examples are the studies carried out within the lungs and trachea of a human. For instance, authors like Faizal et al. (2020), Augusto et al. (2016), Tsega et al. (2018), Qi et al. (2014), Tsega et al. (2022) and Azarnoosh et al. (2016) have used CFD to carry out investigations that range from particle deposition on the respiratory tract to respiratory flow patterns (studied with specific flow parameters such as velocity, pressure, and shear stress). For external respiratory events, Oh et al. (2022), Chillón et al. (2023), Zee et al. (2021) and Li et al. (2018) have applied CFD to study the particles expelled from events such as coughing, sneezing, breathing, talking, etc. To carry out these studies with more precision, the availability of experimental data (Zee et al., 2021) even though they are limited (Cao et al. 2022), the choice and application of adequate physical model (Li et al., 2018), solver (Walsh and Boyle, 2020; Li et al., 2021), boundary conditions (Dao et al., 2022 and Aliabadi et al., 2010), and adequate simplification and meshing of the flow domain (D'Alessandro et al., 2021 and Jiang, 2020) is important. Also, an adequate numerical model needs to be applied to calculate the evolution of the respiratory droplets under specific environmental conditions. Eulerian-Lagrangian model has been applied by many authors like D'Alessandro et al. (2021), Dbouk and Drikakis (2020), Yan et al. (2019), Li et al (2020) and Li et al. (2018). Eulerian model solves the corresponding transport equations related to the continuous phase (Chen et al., 2022, Oh et al., 2022, D’Alessandro et al., 2021) while the Lagrangian model tracks and calculates each particle’s position and velocity and updates them as they move through the flow field. (Xiao et al., 2023, Chéron et al., 2023, Zhang and Chen, 2007). In the CFD of respiratory events, the air 5 and water vapor make up the continuous phase (Eulerian phase) which is solved by the Eulerian model, while the droplets which are largely water form the Lagrangian phase (Chillón et al., 2023). Regarding the computational tool applied to study the CFD of respiratory droplets, ANSYS Fluent, which is the commercial software is applied in this work, provides pressure and density-based solvers. The respiratory events fit into incompressible flows; therefore, the pressure- based solver can be applied. This solver offers coupled and segregated algorithms (ANSYS, 2013). The coupled algorithm solves the momentum and continuity equations in a coupled form while the segregated algorithm solves the flow equations individually. The application of the coupled algorithm improves the convergence speed when compared to the segregated counterpart. However, the later algorithm requires lesser computational memory. Therefore, the time-cost factor puts the coupled scheme at an advantage over the segregated one. In the coupled algorithm, the SIMPLE, SIMPLEC, PISO, coupled, and FSM are the pressure-based coupling schemes available in the commercial software adopted in this study. Some of these have been reported in literature; SIMPLE by Cao et al. (2022) and Oh et al. (2022). while D'Alessandro et al used PISO, which is recommended for transient calculations which permits the application of large time-steps, and highly distorted mesh (ANSYS, 2013). 2.3 Boundary characterization In general studies of CFD of respiratory events, setting up boundary conditions is crucial to obtaining qualitative results (Chen et al.,2022 and Li et al., 2023). The boundary conditions include the characteristics attributed to the velocity inlets, pressure outlets, stationary or moving walls and the flow space. Gupta et al. (2009) have highlighted flow direction, respiratory surface area, temperature, and respiratory droplet size distribution as the specific boundary conditions. The respiratory particles can be emitted in the direction normal to the respiratory surface area, as applied in the works of Redrow et al. (2011) and Bourouiba et al. (2014), or inclined at a specific angle as discussed by Oh et al. (2022) and Gupta et al. (2009). The inclined surface inlets may approximate the reality since exhalation from mouth region is sometimes not normal to the y-axis (Figure 2). However, many researchers use injections normal to the y-axis due to its 6 economic advantage since they reduce the complexity of the geometry at the inlet, thereby reducing the number of elements applied in the meshing process. Figure 2- The direction of the injection of the respiratory droplets Normal injectionϴ Inclined injection (Self-elaborated) The temperature of the respiratory droplets injected into the flow domain is analogous to the human body temperature. For instance, Redrow et al. (2011) and Chillón et al. (2023) have applied 310.15K while Oh et al., (2021) have applied 307.15K. The first group of authors have applied the standard body temperature while the second have applied the temperature of the body at rest. As for the ambient temperature and humidity applied in the CFD studies, different values have been applied by researchers; specifically, characterized by where the research was conducted. 2.4 Droplet size and distribution Even though there are studies on drops size distribution, the inconsistencies in the values presented by authors make it difficult to elaborate concise studies about fluid dynamics of respiratory events. Undoubtedly, the inconsistencies arise from different factors such as physiological state of people, environmental factors and stages of contamination (for instance at advanced stages of IRDs, coughs and sneezes are normally more violent). To explain the physiological factor, Han et al. (2013) presented results showing that heights, weights and gender influenced the forced vital capacities (FVCs measured in ml) of sneezes. That is to say; a respiratory activity may present different particle sizes depending on the aforementioned factors. Additionally, patients at advanced stages of contamination tend to release larger and higher concentrations of 7 droplets. In terms of particle sizes emitted, no standards have been established (Scharfman et al., 2016). Therefore, Li et al. (2018) have used sizes ranging from 10 to 100 𝜇𝑚 for cough droplets. Gralton et al. (2011) and Perella et al. (2020) have studied 1-500 µm. Wei and Li (2015) discussed cough particle sizes that range from 1 to 1000 𝜇𝑚. Weber and Stilianakis (2008) have applied a large range of 1-2000 µm; this range covers most of the rage analyzed by the CFD of IRDs community. Other authors that have contributed contributed significantly include Yang et al. (2007), Zhang et al. (2019) and Lee et al. (2019. In general, larger respiratory droplets fall to the ground immediately after they have been released (Redrow et al., 2011; LI et al., 2018 and Chillón et al., 2023). On the other hand, smaller droplets persist in the air for a longer period (Bourouiba et al., 2014). In this study, droplets of 1 µm, 10 µm, 50 µm, 100 µm, 200 µm and 500 µm, which represent a range of particles were used. The range of droplet sizes are distributed so that effect of relative humidity, droplet size and wind velocity and other parameters can be studied within a wide range of particle size distribution. 2.5 Computational domain and meshing The domain strategic reduction is also a crucial step towards achieving modeling and simulation with lower computational cost. The domain reduction involves the removal of unnecessary sections such as the human parts of the body that have no direct impact on the simulation results. Such a strategy reduces computational cost (Cao et al., 2022 and Duan et al., 2015). In most of the indoor studies of respiratory events, the whole computation domain (which is a representation of closed surroundings like a room) is represented by a rectangular geometry. The inlet of such geometry (which in reality is the human mouth) is represented by simple geometries in the literature. For instance, a rectangular inlet has been applied in the work of Chillón et al., (2023) and D’Alessandro et al. (2021), a circular one has been used in the work of Oh et al. (2022), Li et al. (2018) and Redrow et al. (2011), a squared inlet has been used by Gupta et al. (2009) and Pan et al. (2022), and lastly, ellipse has been applied in the work of Oh et al. (2021). The use of simple geometries makes meshing more convenient at the domain inlet, since relatively 8 coarser cells can be applied (Cao et al., 2023). In general, it has not been reported that the inlet geometry influences in a large scale, however, it reflects in the computational cost. After a proper domain reduction, numerical meshing is the next step. It has been observed that in the literature, authors have applied different kinds of meshes to carry out the CFD of respiratory events. These meshes have been built with tetrahedral, hexahedral, and polyhedral elements. At this point it is important to highlight that mesh building is not limited to the application of the aforementioned elements, but can also be carried out with prismatic and pyramidal elements (Baldan et al., 2023). Even though the later elements are not frequently applied, they appear in negligible proportions during meshing. Faleiros et al. (2022), Yan et al. (2019) and Redrow et al. (2011) have used tetrahedrons to discretize their flow domains to study the CFD of respiratory events. On the other hand, Chillón et al. (2023), Cao et al. (2022), D'Alessandro et al. (2021), Dbouk et al. (2020) and Li et al. (2018) have applied hexahedral mesh in their studies. Lastly, Oh et al. (2022), Zee et al. (2021) and Li et al. (2020) have used polyhedral mesh. All the aforementioned researchers have meshed their flow domains with different numbers of elements ranging from hundreds of thousands to units of millions of mesh cells. This means they might be working at different ranges of precision, depending on whatever criteria each researcher must have adopted. It has been noticed that authors have applied a different element-domain volume relationship. Redrow et al. (2011) have used 500,000 elements to mesh a 17.94 m3 domain (~ 31,000 elements/m3); Chillón et al. (2023) have used 877,000 elements to mesh their 50.4 m3 flow domain (~ 17,000 elements/m3); D'Alessandro et al. (2021) applied 6 million cells to mesh their 12 m3 flow domain (~ 500,000 elements/m3), and Oh et al. (2022) have applied about 500,000 elements to mesh a 14.4 m3 flow domain (~ 35,000 elements/m3). The element-domain volume relationships mentioned above reflect the complexity of each scenario. For instance, a domain containing one human is simpler than one with two, and an empty room is less complex than an airplane cabin. The various studies by the mentioned authors were conducted under different conditions. Consequently, it is expected that more complex domains will have higher mesh densities regardless of their dimensions. 9 In the CFD studies of respiratory activities, the choice of cell geometry to build a mesh is influenced by on the complexity of the flow geometry or the cost, above all, to obtain an accurate result. For instance, it is difficult to apply hexahedrons (which are orthogonal cells) to build mesh on complex geometries (Wang, 2021), even though they have a reduced error approximation to the surface normal gradient, when compared to tetrahedrons and polyhedrons (Baudouin et al., 2014 and Benzle et al., 1995), in other words, the hexahedral meshes yield more accurate solutions than their tetrahedral counterparts for the same number of elements (Biswas, 1998). However, since most of the studies in this field are carried out with simplified geometries that can be meshed up with almost all forms of geometries available, the choice of application is more influenced by computational cost. In some literature, the inlet region (which represents the human mouth in reality) of the computational domain is meshed with the highest cell density. A gradual increment in average cell size is observed while moving away from that region (Figure 5). Therefore, the coarsest elements are found in the farthest regions away from the inlet. Figure 3- Graded mesh with quadrilateral elements (2D) or hexahedral elements 3D Inlet region (self-elaborated) The reason for applying this form of mesh grading is to achieve a more accurate discretization of the droplets and transport equations (D’alessandro et al., 2021). When the Discrete 10 Particle Model (DPM) is applied to track the motion of each respiratory droplet using a surface injection in commercial software like Fluent, each face cell at the inlet region releases a single particle per DPM flow time (Fluent, 2013). The measure of quality is eminent in every numerical mesh generated to solve any CFD problem. In the software used to elaborate this work, mesh quality measurement can be carried out through measuring metric factors such as orthogonality, skewness, aspect ratio, warping, element quality among other prominently factors. Table 3 presents different parameters to evaluate the quality of mesh. Accessing mesh details is important because it is an important criterion of precision (Baudouin et al., 2014) that reflects in the modeling and simulation results. In this research, the skewness has been mostly used to measure the quality of mesh. As a complement to Table 3, this measurement of mesh quality is explained thus. The measure of the quality of mesh was carried out through the evaluation of the skewness which compares a cell geometry applied in meshing to an ideal cell geometry. For instance, the ideal tetrahedral and hexahedral cells are respectively formed with equilateral triangular and orthogonal quadrilateral faces. Therefore, the skewness for each cell geometry is measured differently. For a tetrahedral cell, the skewness is calculated from Eq. 2.11 which is the ratio of the real cell volume to the ideal volume. 𝑆𝑘𝑒𝑤𝑛𝑒𝑠𝑠 = 1 − 𝑉𝑟𝑒𝑎𝑙 𝑉𝑖𝑑𝑒𝑎𝑙 2.11 The real cell volume always satisfies the limit 𝑉𝑟𝑒𝑎𝑙 ≤ 𝑉𝑖𝑑𝑒𝑎𝑙. The volumes of the ideal and real cells are respectively calculated by the relation 𝑉𝑖𝑑𝑒𝑎𝑙 = 0.5132𝑅𝑖𝑑𝑒𝑎𝑙 3 , where 𝑅𝑖𝑑𝑒𝑎𝑙 is the radius of the external sphere formed around the tetrahedral cell as shown in Figure 3 (a) ii. To calculate the real volume of (i), the same relation is used. 11 Table 1- Mesh quality measurement Mesh metric Description Accepted range Orthogonality Measures the angle between the cell faces and the vector connecting the cell centers. Good orthogonality means that these vectors are perpendicular (or close to it), which helps in minimizing numerical diffusion and ensuring accurate flow calculations. Values above 0.7 are generally considered good, but the acceptable range may vary depending on the specific application and solver settings Aspect ratio Expresses ratio of the longest edge length to the shortest edge length of a cell. High aspect ratio cells can lead to inaccuracies in the solution, especially in regions with high gradients. Values below 5-10 are generally acceptable, but this can vary depending on the flow conditions and solver. Cell quality Often based on the deviation of the cell shape from an ideal shape (e.g., an equilateral tetrahedron for 3D tetrahedral meshes). The quality metric can be calculated using different methods, such as the equiangular skew or the ratio of the cell volume to the volume of an ideal cell with the same circumscribing sphere. It ranges from 0 (worst) to 1 (best). Skewness Measures how much a cell deviates from an ideal shape (e.g., a perfect square in 2D or a cube in 3D). High skewness can lead to numerical errors and convergence issues Zero (0) is the ideal value. Values below 0.5 are generally considered good; values above 0.9 are often problematic. (Adapted from ANSYS, 2013). The skewness of a hexahedral cell measures the deviation of the angles between face pairs from the ideal 90° vertex angles of a perfect hexahedron. The ideal and distorted hexahedral 12 cell geometries are depicted in Figure Figure 6 (b) i and ii, respectively. For hexahedral cells, the value of skewness is calculated by the relation in Eq. 2.12. 𝑠𝑘𝑒𝑤𝑛𝑒𝑠𝑠 = 𝑚𝑎𝑥 ( 𝜃𝑚𝑎𝑥 − 𝜃𝑒 180 − 𝜃𝑒 , 𝜃𝑒 − 𝜃𝑚𝑖𝑛 𝜃𝑒 ) 2.12 where 𝜃𝑚𝑎𝑥 and 𝜃𝑚𝑖𝑛 are the maximum and minimum angles formed by the distorted cells, which are non-orthogonal in this case. 𝜃𝑚𝑖𝑛 is the angle for an equiangular cell (60 for a triangle, 90 for a square) (ANSYS, 2012). Figure 4- Diagrams of skewed and ideal (a) tetrahedral and (b) hexahedral elements that are used to form numerical meshes ϴmin Rideal 90 ° (Adapted from ANSYS, 2012) For both tetrahedral and hexahedral meshes, the skewness value ranges from 0 to 1. Generally, a mesh is considered to be of higher quality when its skewness values are closer to 0. It (b) (a) 13 is typically recommended that simulations are performed using meshes with an average skewness value below 0.5 (ANSYS, 2012). 2.6 Discrete Particle Model (DPM) Discrete Particle Model in Computational Fluid Dynamics is a numerical method used to track individual particles in a fluid flow. It is used to simulate the behavior of particles, such as their motion, interaction with each other, and interaction with the surrounding fluid (Li et al., 2020). In DPM surface injection, each mesh-cell at the surface of the inlet represents a particle. Therefore, each mesh-cell injects one particle per unit time into the flow domain. When a large number of particles is injected into the respiratory domain, the computational cost increases, because an explicit numerical scheme may require a lower time step to keep the Courant-Fiedrichs-Lewy (CFL) number low enough to avoid numerical inaccuracy. This particular dimensionless number is related to the velocity 𝑢, timestep (t) and the element size (𝑥) (Eq. 2.13). 𝐶𝐹𝐿 = 𝑣𝐷𝑡 𝐷𝑥 2.13 CFL sets restrictions on fluid particles so that they move within one grid length (Ferziger et al., 2019) to avoid inaccurate solutions in explicit integral schemes. However, the implicit scheme permits the application of higher CLF numbers (Blazek, 2015, and Miettinen and Siikonen, 2015). For instance, CLF value of 25 has been applied by Faleiros et al. (2022). However, the application of CFL below the unit value is recommended. In the work of Dbouk and Drikakis (2020), it was observed that CLFmax > 0.2 produced slight discrepancies in particle cloud properties trajectory despite applying an implicit formulation. 2.7 Turbulence model Among other important considerations to make in the CFD studies of respiratory droplets emitted during respiratory droplets, selecting a suitable turbulence model is also important 14 (Liu et al., 2013 and Cao et al., 2022). The Large Eddy Simulation (LES) and Reynolds-averaged Navier-Stokes (RANS) are prominently applied in this type of study. In the work of Cao et al. (2022), the two methods have been studied and the authors concluded that even though the LES model gives a more accurate result, it has a higher computational cost when compared to the RANS models. It is known that the computational cost is an important factor when carrying out CFD studies of respiratory events. However, such consideration becomes secondary when it risks the precision of the results. The commercial software applied in this study offers several RANS turbulence models, that extend from the one-equation Spalart Allmaras model to the most elaborate second-order closure model; the Reynolds Stress Model (RSM). Nevertheless, the k-ε and k-ω have been frequently applied. The authors in Table 4 have used different kinds of RANS turbulence model to carry out their research. Table 2- Some turbulence models applied in the CFD studies of respiratory droplets Turbulence models (SST) k-ω RNG k-ε Realizable k-ε Standard k-ε Authors -Cao et al. (2022) -Faleiros et al. (2022) -Dbouk and Drikakis (2020) -Cao et al. (2022) -Li et al. (2018) and Redrow et al. (2011) -Yan et al. (2019) -Cao et al. (2022) -Oh et al. (2022) -Zee et al. (2021) -Cao et al. (2022) (self-elaborated) Even though all the models present satisfactory results, the (SST-k-ω) has been recommended by Faleiros et al. (2022). The reason behind the recommendation is the fact that it ensures proper zone selection and modeling without any user interference. In such a model, the turbulence is solved using k-ε in the near-wall regions, while k-ω is activated in regions far from 15 the wall (Menter, 2003). Such a model could increase the precision at the inlet (mouth) region since it has a high cell density. 16 3. STUDY OBJECTIVE The main objective of this research is to apply the Eulerian-Lagrangian to carry out the CFD studies of droplets emitted during respiratory events using the commercial software FLUENT 14.5. Specifically, the objective of the research is divided into four parts: First Part (Section 6): I. Objective: Study the effect of turbulence on the evaporation of a single respiratory droplet. II. Details: Compare the SST k-ω, k-ε, and Reynolds Stress (RSM) turbulence models using a 2D simulation. Second Part (Section 7): III. Objective: Conduct a 3D mesh analysis. IV. Details: Compare two prominent meshes (tetrahedral and hexahedral) applied in CFD studies of droplets emitted during human respiratory events. Third Part (Section 8): V. Objective: Investigate the influence of phase coupling, droplet diameter and relative humidity on human cough droplet evaporation. VI. Details: Use a 3D tetrahedral mesh to calculate the evolution of droplets at different relative humidity conditions. Fourth Part (Section 9): VII. Objective: Apply the tests results from previous research (Firs Part, Second Part and Third part). VIII. Details: Study the behavior of respiratory droplets under the climatic conditions of São Carlos city, São Paulo State, Brazil, and conduct a complementary study on flow behavior at different boundary conditions. 17 4. MATHEMATICAL MODEL 4.1 Eulerian model In this research, a multiphase-phase flow is conducted, therefore the mathematical model is built to calculate the conservation of mass, momentum and energy for each component of the flow. The air and water vapor make up the continuous phase in this research. Therefore, a single conservation of mass is given by Eq. 4.1. 𝜕(𝜌𝑔) 𝜕𝑡 + 𝛻 ⋅ (𝜌𝑔𝑣𝑔) = �̇�𝑑→𝑔 4.1 where 𝜌𝑔 denotes the density (kg m-3) of the continuous phase, 𝑡 represents the time (s) and 𝑣𝑔 represents the mass-averaged velocity (m s-1). The term �̇�𝑑→𝑔 stands for the mass transfer rate through evaporation from the discrete phase (which is composed of the respiratory droplets) to the continuous phase which is as a result of coupling between the continuous and discrete phases (kg m-3 s-1). As the continuous phase is a multicomponent mixture of air and water vapor, the density, 𝜌𝑔, varies according to the mixture rule defined in Eq. 4.2. 𝜌𝑔 =∑𝑌𝑖 𝜌𝑖 4.2 where 𝑌𝑖 and 𝜌𝑖 are the mass fraction and the density of specie i in the mixture, respectively. The momentum conservation of the Eulerian phase can be represented by Eq. 4.3. 𝜕(𝜌𝑔𝑣𝑔) 𝜕𝑡 + ∇ ⋅ (𝜌𝑔𝑣𝑔𝑣𝑔) = − ∇𝑃𝑔 + ∇ ⋅ �̿�𝑔 + 𝐹𝑔 + 𝑆𝑔 4.3 18 The terms 𝑃𝑔, �̿�𝑔 and 𝐹𝑔 respectively represent pressure (Pa), viscous stress, volumetric force, source term for momentum exchange between the continuous and discrete phase. Given that the respiratory case study in this work is turbulent, the phase stress tensor (τ̿𝑔) is defined by Eq. 4.4. τ̿𝑔 = 𝜇𝑔 [(∇𝑣𝑔 + ∇𝑣𝑔 𝑇) − 2 3 ∇ ⋅ 𝑣𝑔𝐼]̿ 4.4 where (∇𝑣𝑔 + ∇𝑣𝑔 𝑇) term represents the rate of strain tensor, which is symmetric and accounts for the deformation (both shear and elongation) of the fluid. 𝐼 ̿is the identity tensor, a unit matrix that preserves the dimensional integrity of the terms in the equation. The second term on the right-hand side represents the effect of volume dilation. The energy conservation of the continuous phase is given by Eq. 4.5, accounting for the internal energy change due to pressure, viscous effects, and heat flux. 𝜕(𝜌𝑔ℎ𝑔) 𝜕𝑡 + ∇ ⋅ (ρ𝑔𝑣𝑔ℎ𝑔) = ∇ ⋅ (τ̿𝑔 ∙ 𝑣𝑔) − ∇ ⋅ 𝑞𝑔 − �̇�𝑑→𝑔ℎ𝑑,𝑔 4.5 where the first and the second terms of Equation 15 represent the local time rate of change of the enthalpy per unit volume and the convective transport of enthalpy, respectively. The first and the second terms on the right-hand side of the equation represent the rate of work done by viscous stresses and the divergence of the heat flux, which can include conductive and possibly viscous dissipation components. �̇�𝑑→𝑔 and ℎ𝑑,𝑔 are the mass transfer rate from the droplet to vapor and the latent heat of vaporization. This is also the source term which results from the exchange of energy between the discrete phase and the continuous phase. Finally, the Eulerian phase component mass balance is given by Eq. 4.6. 𝜕(𝑌𝑖𝜌𝑔) 𝜕𝑡 + ∇ ⋅ (𝜌𝑔𝑣𝑔𝑌𝑖) = − 𝛻 ∙ 𝐽𝑖 + 𝑆𝑖 4.6 where 𝑌𝑖 is the mass fraction of the specie i, 𝑆𝑖 is the component source term, and 𝐽𝑖 is the diffusion flux of specie i, which arises due to gradient of concentration. This is given by Eq. 4.7. 19 𝐽𝑖 = −( 𝜇𝑔 𝑆𝑐𝑔 + 𝜇𝑔,𝑡 𝑆𝑐𝑔,𝑡 )∇ ∙ 𝑌𝑖 4.7 where 𝜇𝑔 and 𝜇𝑔,𝑡 are respectively air viscosity and its turbulent contribution. 𝑆𝑐𝑔 and 𝑆𝑐𝑔,𝑡 are the respective molecular Schmidt number and its turbulent contribution. While there are several methods to conduct turbulent CFD simulations, each has its application advantages and limitations. These approaches are Reynolds-Averaged Navier-Stokes (RANS), Large Eddy Simulation (LES), and Direct Numerical Simulation (DNS). RANS is the simplest and the most popular method to solve turbulence in fluid flows. The idea about this method is that; since turbulent flows are unpredictable and occur over large- scale range of length and time, it becomes necessary to work with time-average of the equations governing the flow. For real fluids, the instabilities observed in the flows can be measured through interaction of non-linear inertial terms and viscous term contained in the flow equations. Through averaging, the aforementioned complexities add-up to RANS transport equations to form unknown variables or second moments. Therefore, higher-order correlations are needed to solve the unknown flow variables. However, the variables need to be supplemented before any solution can be obtained. Hence, algebraic or differential relations (also known as turbulent closures) are needed through defining a number of auxiliary equations. RANS closure usually focuses on modeling Reynolds stress tensor (Perot and Wang, 1999). The method offers the most economic means of solving complex turbulent industrial flow problems. In this approach, turbulence modeling is achieved through Reynolds’ decomposition of flow variables into average and fluctuating parts, respectively. Further substitution of the Reynolds- decomposed variables in the Navier-Stokes equation and taking an average result to an unknown Reynold stress tensor. This unknown force is modelled in order to obtain a solution (Afonsi, 2009). Typical examples of closure models that can be applied to solve RANS equations include Spalart-Allmaras One-Equation, k-ε, k-ω, Reynolds Stress Model (RSM), Laminar- Turbulent Transition Models, among others. Although all of them can be applied specifically to 20 solve turbulent flow equations, the approach of algebraic closure adopted in each method to solve Reynolds’s tensor vary. 4.1.1 Standard 𝒌 − 𝜺 turbulence model The 𝑘 − 𝜀 is semi-empirical, therefore, its derivation depends on phenomenological or empirical considerations. It is modeled based on transport equations for turbulent kinetic energy (𝑘) and its dissipation rate (𝜀). Its robustness, considerable accuracy and economy in solving a wide range of engineering flows and heat transfer are some of the reasons why it is applied frequently in simulations. In the derivation of standard 𝑘 − 𝜀, it is assumed that the flow is fully turbulent, and the effects of molecular viscosity are negligible. This method is only valid for fully turbulent flows. Thus, the turbulent kinetic energy (𝑘) and dissipation rate (𝜀) are obtained from Eqs. 4.8 and 4.9. 𝜕 𝜕𝑡 (𝜌𝑘) + ∇. (𝜌𝑘𝑣) = ∇. [(𝜇 + 𝜇𝑡 𝜎𝑘 ) ∇𝑘] + 𝐺𝑘 + 𝐺𝑏 − 𝜌𝜀 + 𝑆𝑘 4.8 𝜕 𝜕𝑡 (𝜌𝜀) + ∇. (𝜌𝜀𝑣) = ∇. [(𝜇 + 𝜇𝑡 𝜎𝜀 )∇𝜀] + 𝐶1𝜀 𝜀 𝑘 (𝐺𝑘+𝐶3𝜀𝐺𝑏) − 𝐶2𝜀𝜌 𝜀2 𝑘 + 𝑆𝜀 4.9 At high Reynolds number, the viscosity term is solved by combining 𝑘 and 𝜀 as presented in Eq. 4.10. 𝜇𝑡 = 𝜌𝐶𝜇 𝑘2 𝜀 4.10 where 𝐺𝑘 and 𝐺𝑏 represent the generation of turbulence kinetic energy due to the mean velocity gradients and generation of turbulence kinetic energy due to buoyancy, calculated as described in Modeling Turbulent Production in the k- ε Models of Fluent (2013). 𝜎𝑘 and 𝜎𝜀 are the turbulent Prandtl numbers for 𝑘 and 𝜀 , respectively. 𝑆𝑘 and 𝑆𝜀 are user-defined source terms. The 21 constants 𝐶𝜇 = 0.09, 𝐶1𝜀 = 1.44, 𝐶2𝜀 = 1.92, 𝐶𝜇 = 0.08, 𝜎𝑘 = 1.0 and 𝜎𝜀 = 1.3 have been determined from experiments for fundamental turbulent flows, and are widely accepted and used as default values for viscous models. 4.1.2 Shear stress transport (SST) k-ω The (SST) k-ω was developed by Menter et al. (2003) in order to eliminate the boundary condition dependency, which the model has with respect to free fluid flow. It is similar to the standard k-ε model. The only peculiarity lies in some adjustments the model incorporates. In (SST) k-ω, standard k-ε and k-ω are multiplied by a unit blending function which is designed to activate either of the models depending on the region where the turbulence would be solved near the walls or away from it (Menter et al., 2003). Transport equations for SST k-ω model are represented by Eqs. 4.11 and 4.12. 𝜕 𝜕𝑡 (𝜌𝑘) + ∇. (𝜌𝑘𝑣) = ∇. (𝛤𝑘∇𝑘) + �̃�𝑘 + 𝑌𝑘 4.11 𝜕 𝜕𝑡 (𝜌𝜔) + ∇. (𝜌𝜔𝑣) = ∇. (𝛤𝜔∇𝜔) + 𝐺𝜔 + 𝑌𝜔 4.12 where 𝑣, 𝛤𝑘 and 𝛤𝜔 are the velocity and effective diffusivity of 𝑘 and 𝜔. �̃�𝑘 , and 𝐺𝜔 generation of turbulence due to mean velocity gradients in 𝑘 and 𝜔. 𝑌𝑘 and 𝑌𝜔 represent turbulent dissipation of 𝑘 and 𝜔. The term 𝐷𝜔 represents a cross-diffusion term. The effective diffusivity of 𝑘 and 𝜔 as are given by Eqs. 4.13 and 4.14. 𝛤𝑘 = 𝜇𝑔 + 𝜇𝑔,𝑡 𝜎𝑘 4.13 𝛤𝜔 = 𝜇𝑔 + 𝜇𝑔,𝑡 𝜎𝜔 4.14 22 The 𝜎𝑘 and 𝜎𝜔 are the turbulent Prandtl numbers for 𝑘 and 𝜔. The term 𝜇𝑔,𝑡 stands for the turbulent viscosity and is computed by Eq. 4.15. 𝜇𝑡 = 𝛼∞ ∗ ( 0.024 + 𝜌𝑘 6𝜇𝜔 1 + 𝜌𝑘 6𝜇𝜔 ⁄ ) 4.15 For turbulent flows, 𝛼∞ ∗ = 1 (Eqs. 4.16 to 4.22 are all functions of 𝜇𝑡). 𝜎𝑘 = 1 𝐹1 𝜎𝑘,1 + 1 − 𝐹1 𝜎𝑘,2 4.16 𝜎𝜔 = 1 𝐹1 𝜎𝜔,1 + 1 − 𝐹1 𝜎𝜔,2 4.17 𝐹1 = 𝑡𝑎𝑛ℎ {𝑚𝑖𝑛 [𝑚𝑎𝑥 ( √𝑘 0.09𝜔𝑦 , 500𝜇 𝜌𝑦2𝜔 ) , 4𝜌𝑘 𝜎𝜔,2𝐷𝜔 +𝑦2 ]} 4.18 𝐷𝜔 + = 𝑚𝑎𝑥 (2𝜌 1 𝜎𝜔,2 1 𝜔 𝜕𝑘 𝜕𝑥𝑗 𝜕𝜔 𝜕𝑥𝑗 , 10−10) 4.18 �̃�𝑘 = min (𝐺𝑘 , 10𝜌𝛽 ∗𝑘𝜔) 4.19 In a manner of Boussinesq hypothesis, 𝐺𝑘 = 𝜇𝑡𝑆 2 4.20 Where S is mean rate-of-strain tensor. And 𝐺𝜔 = �̃�𝑘 𝛼 𝑣𝑡 4.21 23 For SST 𝑘 − 𝜔, 𝛼∞ = 𝐹1 ( 𝛽𝑖,1 𝛽∞ ∗ − 𝑘2 𝜎𝑤,1√𝛽∞ ∗ ) + (1 − 𝐹1) ( 𝛽𝑖,2 𝛽∞ ∗ − 𝑘2 𝜎𝑤,2√𝛽∞ ∗ ) 4.22 Equations 4.24 to 4.27 are functions of k dissipation (𝐷𝑖𝑠𝑘). 𝐷𝑖𝑠𝑘 = 𝜌𝛽 ∗𝑓𝛽𝑘𝜔 4.23 𝑓𝛽 = { 1 𝑥𝑘 ≤ 0 1+680𝑥𝑘 2 1+400𝑥𝑘 2 𝑥𝑘 ≥ 0 4.24 𝑥𝑘 = 1 3 𝜕𝑘 𝜕𝑥 𝜕𝑤 𝜕𝑥 4.25 𝛽∗ = 𝛽𝑖 ∗[1 + 𝜉∗𝐹(𝑀𝑡)] 4.26 𝛽𝑖 ∗ = [ 0.27(𝑅𝑒𝑡/𝑅𝑒𝛽) 4 1 + ( 𝑅𝑒𝑡 𝑅𝑒𝛽 )4 ] 4.27 where 𝑅𝛽 = 8, 𝛽∞ ∗ = 0.09 and 𝜉∗ = 1.5 Equations 4.29 to 4.34 are functions of dissipation of ω (𝐷𝑖𝑠𝜔). 𝐷𝑖𝑠𝜔 = 𝜌𝜔 2𝛽𝑓𝛽 4.28 𝛽 = 𝛽𝑖 [1 − 𝛽𝑖 ∗ 𝛽𝑖 𝜉𝐹(𝑀𝑡)] 4.29 24 𝐹(𝑀𝑡) = { 0 𝑀𝑡 ≤ 𝑀𝑡0 𝑀𝑡 2 −𝑀𝑡0 2 𝑀𝑡 > 𝑀𝑡0 4.30 𝑓𝛽 = 1 + 70𝑥𝜔 1 + 80𝑥𝜔 4.31 𝑥𝜔 = | 𝑄𝑖𝑗Ω𝑗𝑘𝑆𝑘𝑖 (𝛽∞∗ 𝜔)3 | 4.32 𝑀𝑡 2 = 2𝑘 𝛾𝑅𝑇 4.33 𝑄𝑖𝑗 = 1 2 ( 𝜕𝑢𝑖 𝜕𝑥𝑗 − 𝜕𝑢𝑗 𝜕𝑥𝑖 ) 4.34 where at high Reynolds numbers, 𝛽𝑖 ∗ = 𝛽∞ ∗ = 0.09, 𝑀𝑡0 = 0.25, 𝛼∞ ∗ = 1, 𝛼∞ = 0.52, 𝛼0 = 1 9 , 𝛽𝑖 = 0.072, 𝑅𝛽, 𝑅𝑘 = 6, 𝑅𝜔 = 2.95, 𝜉 ∗ = 1.5, 𝜎𝑘 = 2.0, 𝜎𝜔 = 2021. Since SST k-𝜔 is based on both standard k-𝜔 and 𝑘 − 𝜀, in order to blend the two models together, 𝑘 − 𝜀 has been transformed into k and -𝜔 equations which resulted in the inclusion of diffusive term 𝐷𝜔 as represented by Eq. 4.35. 𝐷𝜔 = 2𝜌(1 − 𝐹) 𝜔𝜎𝜔,2 𝜕𝑘 𝜕𝑥𝑗 𝜕𝜔 𝜕𝑥𝑗 4.35 4.1.3 Reynolds stress model (RSM) These models do not use Eddy Viscosity hypothesis. However, they model transport equations for all components of the Reynolds stress tensor and dissipation rate in the fluid. They are suitable for complex flows since it includes effects of streamline curvature, sudden changes in 25 the strain rate, buoyancy flows compared to (SST) k-ω and Realizable 𝑘 − 𝜀. Specifically, they model solves 5 viscous Equations. Nevertheless, researches have shown that Reynolds Stress model is not superior to two-equation models, especially in anisotropic flow. The momentum equation is given by Eq. 4.36. 𝜕𝜌(𝑣𝑖′𝑣𝑗′̅̅ ̅̅ ̅̅ ̅) 𝑡 + 𝐶𝑖𝑗 = −𝐷𝑇,𝑖𝑗 + 𝐷𝐿,𝑖𝑗 + 𝜑𝑖𝑗 − 𝜀𝑖𝑗 4.36 where 𝜕𝜌𝑣𝑖′𝑣𝑗 ′ 𝜕𝑡⁄ is a local derivative, 𝐶𝑖𝑗 is the convective term, −𝐷𝑇,𝑖𝑗 is the turbulent diffusion, 𝐷𝐿,𝑖𝑗 is the molecular diffusion, , 𝜑𝑖𝑗 is the pressure strain, −𝜀𝑖𝑗 is dissipation. In the model, the terms 𝐶𝑖𝑗, 𝐷𝐿,𝑖𝑗, −𝑃𝑖𝑗 and −𝐹𝑖𝑗 are calculated. On the other hand, −𝐷𝑇,𝑖𝑗, −𝐺𝑖𝑗, 𝜑𝑖𝑗 and −𝜀𝑖𝑗 are modeled. Details of the closure equations are available in ANSYS (2013). 4.2 Lagrangian formulation The organization of this stage is carried out in the following steps; General Equations, droplet dispersion, droplet evaporation, droplet break, and coupling. 4.2.1 Force balance equation In this formulation, the particle movement is calculated applying force balance using the Newton’s second law of motion. Therefore, the trajectory of any discrete particle is predicted by integrating the force balance on the Lagrangian reference frame (Zahari et al., 2018). The force balance is equated to the particle inertial as written in Eq. 4.37. 𝑑𝑣𝑑 𝑑𝑡 = 𝐹𝐷(𝑢𝑔 − 𝑣𝑑) + 𝑔 (𝜌𝑑 − 𝜌𝑔) 𝜌𝑑 4.37 where 𝜌𝑑 represents the density of the droplet (kg m-3) and 𝐹𝐷 stands for the fluid drag force per unit mass given by Eqs. 4.38. 26 𝐹𝐷 = 0.75 𝜇𝐶𝐷𝑅𝑒 𝜌𝑑𝑑𝑑 3 4.38 where 𝐶𝐷 is the drag coefficient explained in Liu et al. (1993), and 𝑑𝑑 is the particle diameter (m). The 𝑅𝑒 in Eq. 37 is the relative Reynolds number, which is given by Eq. 4.39. 𝑅𝑒 = 𝜌𝑑𝑑|𝑣 𝑑 − 𝑣 𝑔| 𝜇𝑔 4.39 4.2.2 Droplet dispersion The dispersion of droplets caused by turbulence can be calculated using a stochastic model. Specifically, the Discrete Random Walk (DRW) model was employed in this study. This stochastic model is used to compute particle dispersion within a turbulent flow field (Gosman and Loannides, 1983). The droplet velocity is, then, recalculated as a sum of its mean velocity (�̅�) and a fluctuating contribution (𝑣′), as presented in Eq. 4.40. 𝑣𝑔 = �̅� + 𝑣 ′ 4.40 The fluctuating (𝑣′) is treated as a random variable and simulated using Gaussian random value with zero mean and a unit variance. Each of the particle eddy is described by a Gaussian distributed random velocity fluctuation 𝑣′, as shown in Eq. 4.41. 𝑣′ = √2𝑘/3 𝜁 4.41 Where k and 𝜁 are the turbulent kinetic energy (m2 s−2) and the normally distributed random number with zero mean and unit variance. The fluctuating component can be computed based on local turbulence characteristics, typically derived from the turbulence kinetic energy and the turbulence dissipation rate. 27 The dispersion of the respiratory droplets is calculated based on the concept of integral time scale (T), as shown in Eq. 4.42. For smaller tracer particles that dislocate with zero drift velocity with reference of the transporting fluid, the integral time is given as Eq. 4.43, which is the fluid Lagrangian integral time (𝑇𝐿). 𝑇 = ∫ 𝑣′𝑑(𝑡)𝑣 ′ 𝑑(𝑡 − 𝜏) 𝑣′2̅̅ ̅̅ 𝑑𝜏 ∞ 0 4.42 𝑇𝐿 = 𝐶𝐿 𝑘 𝜀 4.43 where 𝜏 represents particle relaxation time (s), which is proportional to the dispersion rate, and means; the larger integral time, the larger the turbulence in the flow. The 𝐶𝐿 is an empirical constant that ranges from 0.1 to 0.3 depending on the turbulent model, for example 0.15 for the k-ε turbulent models, 0.30 for the RSM. For the k-ω models, 𝜔 = 𝜀/𝑘 in Eq. 4.43. In DRW model, the interaction of a droplet with discrete stylized fluid phase turbulent eddies given by a Gaussian distribution of random velocity fluctuation (Eq. 4.44) and the characteristics lifetime Eddy (𝜏𝑒) (Eq. 4.45) are simulated. 𝜏𝑒 = 2𝑇𝐿 4.44 It is considered that the particle interacts with the fluid phase over a minimum Eddy lifetime (4.45) 𝑡𝑖𝑛𝑡𝑒𝑟𝑎𝑐𝑡𝑖𝑜𝑛 = 𝑓 (𝜏, 𝑡𝑐𝑟𝑜𝑠𝑠) 4.45 where tcross is the time needed for a particle to cross an eddy, and it is expressed by Eq. 4.46. 𝑡𝑐𝑟𝑜𝑠𝑠 = − 𝜏 𝑙𝑛 [1 − ( 𝐿𝑒 𝜏|𝑣𝑔 − 𝑣𝑑| )] 4.46 28 where 𝐿𝑒 and |𝑣𝑔 − 𝑣𝑑| are the eddy length scale, and the magnitude of the relative velocity, respectively. The time scale of the particle 𝜏𝑑 and the eddy length scale 𝐿𝑒 are given by Eqs. 4.47 and 4.48. 𝜏𝑑 = 1 18 𝜌𝑑𝑆𝑑𝑑 2𝐶𝑐 µ𝑑 4.47 𝐿𝑒 = 𝜏𝑑√ 2𝑘 3 4.48 where S and Cc are shape factor and Cunningham slip correction factor (given by Eq. 4.49). 𝐶𝑐 = 1 + 2𝜆 𝑑𝑑 (0.4𝑒 −( 1.1𝑑𝑑 2𝜆 ) + 1.257) 4.49 where λ is the cough droplets mean free path. When the interaction time (Eq. 4.45) is reached, a new value of instantaneous velocity is calculated by using a new value of 𝜁 in Eq. 4.44. 4.2.3 Droplet evaporation The evaporation of cough droplets is governed by the diffusive flux of droplet vapor in air (Li et al., 2020). In the simulations presented in this study, it was considered that the evaporation initiates when the droplet temperature, 𝑇𝑑, reaches the vaporization temperature, 𝑇𝑣. The evaporation continues until the droplet completely disappears or when it reaches the boiling temperature, 𝑇𝑏, as represented by Eq. 4.50. 𝑇𝑣 ≤ 𝑇𝑑 < 𝑇𝑏 4.50 When the evaporation rate is slow, it can be assumed to be primarily governed by gradient diffusion. In this case, the vapor transitioning from the droplet into the gas phase is 29 assumed to be driven by the difference in vapor concentration between the droplet surface and the bulk gas. The droplet mass variation is, then, given by Eq. 4.51). 𝐷𝑚𝑑 𝐷𝑡 = −𝑘𝑐(𝐶𝑣,𝑑 − 𝐶𝑑,∞) 𝐴𝑑𝑀𝑊𝑤 4.51 where 𝑘𝑐, 𝐶𝑤,𝑑 and 𝐶𝑤,∞ are, respectively, the mass transfer constant (m s-1), the water concentration at the droplet surface (kmol m-3) and the water concentration in the bulk gas (kmol m-3), respectively. 𝑀𝑊𝑤 and 𝐴𝑑 are the molecular weight of water (kg kmol-1), and surface area of droplet (m2), respectively. The mass transfer coefficient 𝑘𝑐 is calculated from Froessling Equation (Eq. 4.52) through correlating Sherwood (𝑆ℎ) and Schmidt (𝑆𝑐). 𝑆ℎ = 𝑘𝑐𝑑𝑑 𝐷𝑣 = 2 + 3𝑅𝑒𝑑 1/2 𝑆𝑐1/3 5 4.52 where 𝐷𝑣 and 𝑑𝑝 are, respectively, the diffusion coefficient of water vapor in air (m2 s-1), and the Schmidt number. Eq. 4.52 correlates the data for Reynolds numbers ranging from 2 to 800 and Schmidt numbers from 0.6 to 2.7. The data from Evnochides and Thodos correlate for Reynolds numbers between 1,500 and 12,000 and Schmidt numbers between 0.6 and 1.85. This work was conducted within this range. At the droplet surface, where the temperature is 𝑇𝑑, the vapor concentration is calculated by assuming that the partial pressure of vapor at the surface is equivalent to the saturated vapor pressure (Psat). 𝐶𝑣,𝑑 = 𝑃𝑠𝑎𝑡 𝑅𝑇𝑑 4.54 where R is the universal gas constant (atm L mol-1 K-1). Additionally, the concentration of the vapor present in the bulk gas is calculated from Eq. 4.55. 30 𝐶𝑣,∞ = 𝑝 𝑅𝑇∞ 𝑌𝑣 4.55 where 𝑝 (Pa) is the absolute pressure onto the domain, 𝑇∞ (K) is the temperature of the surrounding air and 𝑌𝑣 is the local bulk mole fraction of the water vapor. For high evaporation rates, the effect of convection during droplet evaporation is relevant, thus modeled based on that principle. Li et al. (2018) and Lu et al. (2020) presented the evaporation rate from the droplet surface into the bulk gas phase (Eq. 4.56). 𝑑𝑚𝑑 𝑑𝑡 = 𝑘𝑐𝐴𝑝𝜌∞ln (1 + 𝑌𝑣,𝑑−𝑌𝑣,∞ 1 − 𝑌𝑣,𝑑 ) 4.56 where 𝜌∞ is the density of bulk gas (air, in this case). 𝑌𝑣,𝑑 and 𝑌𝑣,∞ are the respective vapor mass fraction at the droplet surface and vapor mass fraction in the bulk gas (air). During the process of evaporation, the temperature of the droplet is updated according to the heat balance equation between the droplet and the surrounding (continuous phase) (Eq. 4.57). 𝑚𝑑𝑐𝑑 𝑑𝑇𝑑 𝑑𝑡 + 𝑑𝑚𝑑 𝑑𝑡 ℎ𝑓𝑔 = ℎ𝐴𝑑(𝑇∞ − 𝑇𝑑) 4.57 where 𝑐𝑑 is the heat capacity of the droplet (J kg-1 K-1), ℎ𝑓𝑔 is the latent heat (J kg-1), ℎ is the coefficient of heat transfer (W m2 K-1). 4.2.4 Droplet Break- Taylor Analogy Breakup (TAB) Collisions and breakup of particles occur rampantly in nature. For example, they are seen during the formation of raindrops, activities of nuclear reactors, spray combustion, firefighting, and material spraying processes (Pan et al., 2009). Specifically in this work, when the respiratory particles expelled travel through the air, it is believed that they are subjected to forces that may lead to their breakup. Therefore, to conduct such analysis, the choice for the model will 31 depend on the flow characteristics, most especially indicated by the Weber number (We). The Weber number is dimensionless number compares the relevance of inertial forces when compared to surface tension as given by (Eq. 4.58) (Campbell, 2015). 𝑊𝑒 = 𝐷𝑟𝑎𝑔 𝐹𝑜𝑟𝑐𝑒 𝐶𝑜ℎ𝑒𝑠𝑖𝑜𝑛 𝐹𝑜𝑟𝑐𝑒 = 𝜌𝑑𝑟𝑑𝑣𝑑 2 𝜎 4.58 where 𝜎 is the particle surface tension (N m-1). When We = 1, there are virtually no unbalanced forces acting on the particles (Campbell, 2015). Therefore, it may be assumed that flows are free from surface turbulence. However, We >100 is a value that is considered high (Pan et al., 2009). At such Weber magnitude, the surface turbulence becomes visible and particles can be seen splashing (Campbell, 2015). Among the particle breakup models, the TAB is recommended for low Weber numbers (We ≤ 100). At the transitional stage (We ≥ 100), KHRT is recommended, while Wave and SSD are recommended for high Weber. In this work, the TAB model was applied since Weber ranged between minimum and maximum values of 0.378 and 78.39, which is within the recommended range. The TAB relates the oscillating and distortion of a droplet to the spring-mass system. For instance, it relates the restoring force of a spring to fluid surface tension force, external force to droplet drag force, and damping force to droplet viscosity force. It is assumed that the oscillations attain a critical value which leads to the breakup of parent droplets into several child droplets. According to O’Rourke and Amsden (1987), the damped forced droplet oscillation is governed by Eq. 4.59. 𝐹𝑜𝑡ℎ𝑒𝑟 𝑚 − 𝑘𝑥 𝑚 − 𝑑𝑝 𝑚 𝑑𝑥 𝑑𝑡 = 𝑑2𝑥 𝑑𝑡2 4.59 where 𝐹𝑜𝑡ℎ𝑒𝑟: Other external forces acting on the mass (N). m: Mass of the particle (kg). 32 k: Spring constant (N/m). x: Displacement of the mass from equilibrium (m). dp: Damping coefficient (N·s/m). 𝑑𝑥 𝑑𝑡 : Velocity of the mass (m/s). 𝑑2𝑥 𝑑𝑡2 : Acceleration of the mass (m/s²). From Taylor’s analogy, the coefficients of Eq. 4.59 are given by Eqs. 4.60, 4.61, and 4.62. 𝐹𝑜𝑡ℎ𝑒𝑟 𝑚 = 𝐶𝐹 𝜌𝑔𝑣 2 𝑑 𝜌𝑑𝑟 4.60 𝑘 𝑚 = 𝐶𝑘 𝜎 𝜌𝑑𝑟3 4.61 𝑑𝑝 𝑚 = 𝐶𝑑 µ𝑑 𝜌𝑑𝑟2 4.62 In nondimensional form, 𝑦 = 𝑥/𝐶𝑏𝑟 is substituted in Eq. 4.59 to form Eq. 4.63. 𝐶𝐹𝜌𝑔𝑣 2 𝑟𝑒𝑙 𝐶𝑏𝜌𝑑𝑟 − 𝐶𝑘𝜎 𝜌𝑑𝑟 3 𝑦 − 𝐶𝑑µ𝑑 𝜌𝑑𝑟 2 𝑑𝑦 𝑑𝑡 = 𝑑2𝑦 𝑑𝑡2 4.63 where 𝑟𝑑 is the droplet radius. The dimensionless constants 𝐶𝑏, 𝐶𝐹, 𝐶𝑘 and 𝐶𝑑 are given by 0.5, 0.33, 8 and 5, respectively. Droplet breakup occurs when 𝑦 > 1, To determine the child droplet size, the energy equations of the parent and child are respectively equated. The energy equation (𝐸𝑃𝑎𝑟𝑒𝑛𝑡) of the parent droplet is given by Eq. 4.64. 𝐸𝑃𝑎𝑟𝑒𝑛𝑡 = 4𝜋𝑟 2𝜎 + 𝑘 𝜋 5 𝜌𝑑𝑟 5 [( 𝑑𝑦 𝑑𝑡 ) 2 + 8𝜎 𝜌𝑑𝑟3 𝑦2] 4.64 33 where 𝜋 5 is geometric factor, ( 𝑑𝑦 𝑑𝑡 ) 2 is the square of the rate of change of y with respect to time, representing kinetic energy, 8𝜎 𝜌𝑑𝑟 3 𝑦 2 is the potential energy term related to the su