**3. Numerical methodologies**

The CFD software ANSYS Fluent [19] is used to solve the Reynolds-Averaged Navier-Stokes (RANS) flow equations for the conservation of mass and momentum. The continuity equation is:

$$\frac{\partial \rho}{\partial t} + \vec{\nabla} \cdot \left(\rho \vec{V}\right) = 0 \tag{1}$$

where *t* is time, *ρ* is density and *V* is the velocity vector. The governing equations for conservation of momentum are:

$$\frac{\partial \left(\rho \vec{V}\right)}{\partial t} + \rho \vec{V} \left(\vec{\nabla} \cdot \vec{\nabla}\right) = -\vec{\nabla}p + \vec{\nabla} \cdot \overline{\vec{\tau}} + \rho \vec{g} \tag{2}$$

where *p* is pressure, *τ* is fluid stress tensor and *g* ! is gravitational body force. The energy equation is:

$$\frac{\partial}{\partial t}(\rho E) + \overrightarrow{\nabla} \cdot \left[\overrightarrow{V}(\rho E + p)\right] = \overrightarrow{\nabla} \cdot k\_{\text{eff}} \overrightarrow{\nabla} T + \overrightarrow{\nabla} \cdot \left(\overrightarrow{\overline{\tau}}\_{\text{eff}} \cdot \overrightarrow{V}\right) \tag{3}$$

where *E* is total energy and *keff* is effective conductivity for the total fluid thermal conductivity and turbulent thermal conductivity.

The realizable *k*-*ε* turbulence model was employed in this study, which contains a formulation for the turbulent viscosity for better predictions for the flows involving boundary layers in strong adverse pressure gradients, separation and re-circulation zones. The transport equations for turbulence kinetic energy *k* and turbulence dissipation *ε* are:

$$\frac{\partial}{\partial t}(\rho k) + \vec{\nabla} \cdot \left(\rho k \vec{V}\right) = \vec{\nabla} \cdot \left[\left(\mu + \frac{\mu\_t}{\sigma\_t}\right) \vec{\nabla} k\right] + \mathbf{G}\_k + \mathbf{G}\_b - \rho \varepsilon \tag{4}$$

and

$$\frac{\partial}{\partial t}(\rho \varepsilon) + \vec{\nabla} \cdot \left(\rho \varepsilon \vec{V}\right) = \vec{\nabla} \cdot \left[\left(\mu + \frac{\mu\_t}{\sigma\_\varepsilon}\right) \vec{\nabla} \varepsilon\right] - \rho \mathbf{C}\_2 \frac{\varepsilon^2}{k + \sqrt{\nu \varepsilon}} + \mathbf{C}\_{1\varepsilon} \frac{\varepsilon}{k} \mathbf{C}\_{3\varepsilon} \mathbf{G}\_b \tag{5}$$

where *μ* is molecular viscosity of the fluid, *μ<sup>t</sup>* is the turbulence viscosity, *Gk* is the generation of turbulent kinetic energy due to mean velocity gradients and *Gb* is the generation of turbulent kinetic energy due to buoyancy. The constants are *C*1*<sup>ε</sup>* ¼ 1*:*44, *C*<sup>2</sup> ¼ 1*:*9, *σε* ¼ 1*:*2, and *σ<sup>k</sup>* ¼ 1*:*0. The constant *C*3*<sup>ε</sup>* determines the degree to which turbulence dissipation *ε* is affected by buoyancy and is calculated using the flow velocity component parallel and perpendicular to the gravitational vector.

The release of respiration particles uses the discrete phase model in ANSYS Fluent. The model employs an Euler-Lagrange approach where the fluid flow in the continuous phase is solved using the Navier-Stokes equations, while trajectory of the particles is predicted by:

$$\frac{du\_p}{dt} = F\_D \left( u - u\_p \right) + \mathbf{g} \frac{\left( \rho\_p - \rho \right)}{\rho\_p} \tag{6}$$

where *u* is the fluid phase velocity, *up* is particle velocity, *FD u* � *up* is the drag force per unit particle mass:

$$F\_D = \frac{18\mu}{\rho\_p d\_p^{\frac{2}{2}}} \frac{C\_D Re}{24} \tag{7}$$

where *ρ<sup>p</sup>* is density of the particle and *dp* is particle diameter. The Reynolds number for the particle motion is:

$$Re = \frac{\rho d\_p |u\_p - u|}{\mu} \tag{8}$$

The pressure-based Navier-Stokes (PBNS) solver is used in this study. The pressure and velocity are coupled using the semi-implicit method for pressure-linked equation (SIMPLE) algorithm. The least squares cell based and PRESTO! schemes are used to discretize gradient and pressure, respectively. Momentum is discretized using the QUICK scheme, and second-order upwind is used to discretize the energy, turbulent kinetic energy and dissipation rates.

The time step, Δ*t* for the simulations was selected such that the Courant-Fredrichs-Levy (CFL) number is close to one:

$$\text{CFL} = V \frac{\Delta t}{\Delta \mathbf{x}} \tag{9}$$

where *V* is the inlet velocity and Δ*x* is the cell size. A grid resolution study was performed, as detailed in **Appendix A**, to estimate error and accuracy of the mesh.
