2.1. Equation of state for a weakly compressible fluid

As it is composed by various gaseous components, air is compressible. For low subsonic flows with the speed of the airflow being less than a Mach number ($Ma$) of $0.1$, or a wind speed less than $100\ {\rm m}\ {\rm s}^{-1}$, we can usually treat air as an incompressible medium. This can greatly simplify the solution procedure, and the conservation equations of the flow can be solved directly by using weakly compressible SPH methods. In the explicit solution scheme, the time marching scheme needs a stable time step $\Delta t$, which depends on the Courant–Friedrichs–Lewy (CFL) condition (Morris, Fox & Zhu Reference Morris, Fox and Zhu1997). Taking into consideration the acceleration, viscous force and Courant viscous force, namely, $\Delta t\leqslant \min (\Delta t_f, \Delta t_\nu, \Delta t_c)$, they are expressed as follows:

(2.1)\begin{equation} \left.\begin{gathered} \Delta t_f \leqslant \alpha_f \min_i\sqrt{{h_I}/{\Vert\pmb{a}_i\Vert}}, \\ \Delta t_\nu \leqslant \alpha_\nu \min_i({h^2_i}/{\nu}), \\ \Delta t_c \leqslant \alpha_c \min ({h_i}/{c_0}). \end{gathered}\right\} \end{equation}

Here $\alpha _f = 0.25$, $\alpha _\nu = 0.125$ and $\alpha _c = 0.25$; $\pmb {a}_i$ represents the acceleration of particle $i$; $\nu$ denotes the kinematic viscosity coefficient; and $h$ denotes the smoothing length. Obviously, a higher sound speed will lead to a smaller stabilization time step. In order to increase the stable time step, Monaghan (Reference Monaghan1994) proposed the weakly compressible smoothed particle hydrodynamics (WCSPH) model, which used the proper artificial sound speed to replace the real sound speed. To satisfy the weakly incompressible condition, it is assumed that

(2.2)\begin{equation} Ma=\frac{\Vert\boldsymbol{v}\Vert}{c_{0}}\leqslant 0.1. \end{equation}

Subsequently, the artificial sound speed must satisfy the following condition:

(2.3)\begin{equation} c_0\geqslant 10\max(\Vert\boldsymbol{v}\Vert). \end{equation}

The artificial sound speed $c_0$ is good for ensuring an acceptable time step while ensuring the fluid flow characteristics, and a fluid with the artificial sound speed can be considered as a weakly compressible fluid. To build the relationship between pressure and density, the Tait equation is adopted as the equation of state in the WCSPH method,

(2.4)\begin{equation} p=B\left[\left(\frac{\rho}{\rho_0}\right)^\gamma-1\right]. \end{equation}

Here $\rho _0$ is the initial fluid density; $\gamma$ is the exponential coefficient, which is 1.4 for gas and 7 for water; and $B$ is related to the sound speed and mass density by

(2.5)\begin{equation} B=\rho_0 c_0^2/\gamma. \end{equation}

2.2. Governing equations and SPH formulations

In the current work, the WCSPH is adopted to perform the numerical simulations. The computational domain is discretized into particles, which have various physical properties. A kernel function is introduced as an interpolation function to approximate the physical quantities of a given point. The kernel approximation and particle approximation are two important approximations in the SPH method (Monaghan Reference Monaghan2005; Crespo et al. Reference Crespo, Domínguez, Rogers, Gómez-Gesteira, Longshaw, Canelas, Vacondio, Barreiro and García-Feal2015), through which the solutions can be obtained. In the kernel approximation, the non-local average of a field function $f(\boldsymbol {r})$ is defined as

(2.6)\begin{equation} \langle \,f(\boldsymbol{r})\rangle =\int_\varOmega f(\boldsymbol{r}')W(\boldsymbol{r}-\boldsymbol{r}',h)\,{\rm d} \boldsymbol{r}'+{O}(h^2),\end{equation}

where $\varOmega$ denotes the compact domain, $\boldsymbol {r}'$ and ${\rm d}\boldsymbol {r}'$ represent the position in the compact domain of $\boldsymbol {r}$ and the volume, respectively, and $W(\boldsymbol {r}-\boldsymbol {r}', h)$ and $h$ are the kernel function and the smoothing length related to the compact domain range, respectively. There are several more frequently used kernel functions, such as the cubic spline function (Schoenberg Reference Schoenberg1946), the Wendland function (Dehnen & Aly Reference Dehnen and Aly2012) and the Gaussian function (Gingold & Monaghan Reference Gingold and Monaghan1977). The spatial gradient $\boldsymbol {\nabla } f(\boldsymbol {r})$ can be found as (Li & Liu Reference Li and Liu2007)

(2.7)\begin{equation} \boldsymbol{\nabla} f(\boldsymbol{r})=\int_\varOmega f(\boldsymbol{r}') \boldsymbol{\nabla} W(\boldsymbol{r}-\boldsymbol{r}',h)\,{\rm d}\boldsymbol{r}' + {O}(h^2) \approx \int_\varOmega f(\boldsymbol{r}') \boldsymbol{\nabla} W(\boldsymbol{r}-\boldsymbol{r}',h)\,{\rm d}\boldsymbol{r}'.\end{equation}

Similarly, the divergence of a vector field $\boldsymbol {\nabla }\boldsymbol {\cdot } {\boldsymbol {f}}(\boldsymbol {r})$ can be expressed as

(2.8)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{f}(\boldsymbol{r})\approx \int_\varOmega\boldsymbol{f}(\boldsymbol{r}') \boldsymbol{\cdot}\boldsymbol{\nabla} W(\boldsymbol{r}-\boldsymbol{r}',h)\,{\rm d}\boldsymbol{r}'.\end{equation}

In the particle discretization, (2.6), (2.7) and (2.8) can be expressed in the discrete form as follows:

(2.9)$$\begin{gather} f(\boldsymbol{r}_i) \approx \sum_{j=1}^{N}f(\boldsymbol{r}_j)W(\boldsymbol{r}_i-\boldsymbol{r}_j,h) \frac{m_j}{\rho_j} , \end{gather}$$

(2.10)$$\begin{gather}\boldsymbol{\nabla} f(\boldsymbol{r}_i) \approx \sum_{j=1}^{N}f(\boldsymbol{r}_j) \boldsymbol{\nabla} W(\boldsymbol{r}_i-\boldsymbol{r}_j,h)\frac{m_j}{\rho_j} , \end{gather}$$

(2.11)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{f}(\boldsymbol{r}_i) \approx \sum_{j=1}^{N}\boldsymbol{f}(\boldsymbol{r}_j)\boldsymbol{\cdot}\boldsymbol{\nabla} W (\boldsymbol{r}_i-\boldsymbol{r}_j,h)\frac{m_j}{\rho_j} . \end{gather}$$

Here $N$ denotes the particle number inside the compact support; $m_j$ denotes the mass of particle $j$; $\rho _j$ is the density of particle $j$; and $m_j/\rho _j$ represents the volume of particle $j$. In passing, we note that more accurate corrective SPH differential operators (see Li & Liu Reference Li and Liu2007) or consistent non-local particle differential operators can be found in Bergel & Li (Reference Bergel and Li2016) and Yan et al. (Reference Yan, Li, Kan, Zhang and Lai2022).

We consider the Lagrangian form of the Navier–Stokes equations in terms of mass and momentum conservation:

(2.12)$$\begin{gather} \frac{{\rm d}\rho}{{\rm d}t} ={-}\rho\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v} + \varPi , \end{gather}$$

(2.13)$$\begin{gather}\frac{{\rm d}\boldsymbol{v}}{{\rm d}t} ={-}\frac{1}{\rho}\boldsymbol{\nabla} p + \boldsymbol{g} + \varGamma . \end{gather}$$

Here $\varPi$ is the density diffusion term; $\boldsymbol {v}$ is the velocity of the fluid; $p$ denotes the hydrostatic pressure; $\boldsymbol {g}$ is the gravitational acceleration; and $\varGamma$ is the dissipative term, which can be regarded as the viscous force term (Lo & Shao Reference Lo and Shao2002), i.e.

(2.14)\begin{equation} \varGamma=(v_0\nabla^2\boldsymbol{v})_i = \sum_j m_j\left(\frac{4v_0(\boldsymbol{r}_j - \boldsymbol{r}_i)\boldsymbol{\cdot}(\boldsymbol{v}_j-\boldsymbol{v}_i)} {(\rho_j+\rho_i)(r_{ij}^2 + \eta^2)}\right)\boldsymbol{\nabla}_i W_{ij}, \end{equation}

where $v_0$ is the kinematic viscosity, and $\eta =0.1h$, which can prevent singularity when $\Vert \boldsymbol {r}_j-\boldsymbol {r}_i\Vert =0$ (Monaghan Reference Monaghan2005; Sun, Ming & Zhang Reference Sun, Ming and Zhang2015).

In computations, the conservation of mass, i.e. (2.12), is discretized as

(2.15)\begin{equation} \frac{{\rm d}\rho_i}{{\rm d}t} = \rho_i\sum\frac{m_j}{\rho_j}(\boldsymbol{v}_j - \boldsymbol{v}_i)\boldsymbol{\cdot}\boldsymbol{\nabla}_i W_{ij} + \varPi , \end{equation}

where $\varPi$ is the density diffusion term, which is expressed as

(2.16)\begin{equation} \varPi = \xi c_0 h_i \sum \frac{(\rho_j - \rho_i)(\boldsymbol{r}_j - \boldsymbol{r}_i)}{\Vert \boldsymbol{r}_j - \boldsymbol{r}_i\Vert^2}\boldsymbol{\nabla}_i W_{ij} V_j, \end{equation}

where $\xi =0.2$ is chosen for the stability of the pressure field (Antuono et al. Reference Antuono, Colagrossi, Marrone and Molteni2010), $c_0$ is the sound speed in the fluid and $V_j$ is the volume of particle $j$.

According to Lo & Shao (Reference Lo and Shao2002), the conservation of linear momentum in (2.13) can be expressed as

(2.17)\begin{equation} \frac{{\rm d}\boldsymbol{v}_i}{{\rm d}t}={-}\sum_j\left(\frac{p_j + p_i}{\rho_i}\right)V_j\boldsymbol{\nabla}_i W_{ij} + \boldsymbol{g} + \varGamma, \end{equation}

where

(2.18)\begin{equation} \varGamma = 4\nu_0\sum_j\frac{m_j(c_j + c_i)(\boldsymbol{v}_j - \boldsymbol{v}_i)\boldsymbol{\cdot}(\boldsymbol{r}_j - \boldsymbol{r}_i)}{(\rho_i + \rho_j)(\Vert\boldsymbol{r}_j - \boldsymbol{r}_i\Vert^2+0.01h^2)}\boldsymbol{\nabla}_i W_{ij}, \end{equation}

where $p$ represents the pressure, $V$ represents the volume of the fluid particle and $\boldsymbol {g}$ represents the gravitational acceleration. Here, $c$ represents the sound speed in the fluid, $h$ represents the smoothing length, $\boldsymbol {v}$ and $\rho$ represent the velocity and mass density, and $\boldsymbol {r}$ represents the position.

2.3. Turbulence model

A major feature of the human cough is its turbulence characteristics. For modelling and simulation of turbulent flows, the most accurate approach is DNS. However, that requires an extremely fine mesh to resolve the entire spatial and temporal scales of the turbulent flow, which usually requires expensive computational cost even at low Reynolds numbers. In the present work, we adopt the LES approach to model the turbulence via low-pass filtering of the Navier–Stokes equations rather than explicit calculation of the turbulence at every length scale. This is because the aerosol dynamics in turbulent flow is at the length scale of $100\ {\rm nm}$ to $1\ \mathrm {\mu }{\rm m}$ (Chun et al. Reference Chun, Koch, Rani, Ahluwalia and Collins2005; Rigopoulos Reference Rigopoulos2019), and it would take tremendous resources and effort to conduct DNS at such a scale.

The LES approach removes small-scale information from numerical solutions by using low-pass filtering. Then the system solver will only need to resolve the scale from the domain size $L$ to the filter size $\varDelta$, i.e. the grid size, while the unresolved scales are modelled through extra formulae. In the present work, the sub-particle-scale (SPS) model (Gotoh Reference Gotoh2001) is adopted to capture the effects of the unresolved turbulence scales,

(2.19)\begin{equation} \frac{{\rm d}\boldsymbol{v}_i}{{\rm d}t} ={-}\sum_j\left(\frac{p_j + p_i}{\rho_i} \right)V_j\boldsymbol{\nabla}_i W_{ij}+\boldsymbol{g} + \varGamma + \frac{1}{\rho _i}\boldsymbol{\nabla}\boldsymbol{\cdot}\bar{\pmb\tau}, \end{equation}

where $\bar {\pmb \tau }$ is the SPS stress tensor,

(2.20)\begin{equation} \bar{\tau}_{ij} = \rho_i (\bar{u}_i\bar{u}_j - \overline{u_i u_j}). \end{equation}

The large-eddy viscosity assumption, i.e. the Boussinesq hypothesis, is used to model the SPS stress tensor as

(2.21)\begin{equation} \bar{\tau}_{ij} = \rho_i (2\nu_t \bar{S}_{ij} - \tfrac{2}{3}\bar S_{kk}\delta_{ij}) - \tfrac{2}{3}\rho _i C_I\delta^2\delta_{ij}, \end{equation}

where $\nu _t$ is the turbulence eddy viscosity, and $\bar {S}_{ij}$ is the sub-particle scale strain tensor in indicial notation. Following Blin, Hadjadj & Vervisch (Reference Blin, Hadjadj and Vervisch2003) and Domínguez et al. (Reference Domínguez2021), we choose $C_I=0.0066$. The eddy viscosity model proposed by Smagorinsky is adopted here to obtain the dynamic viscosity of turbulent flows. In this model, the eddy viscosity is considered to proportional to the mixing length, and the characteristic turbulent velocity is based on the second invariant of the filter field deformation strain $\bar {\boldsymbol {S}}$. In this work, the mixing length is considered to be a feature of sub-particle length scale represented by the particle spacing $\Delta l$. Thus, we have

(2.22)\begin{equation} \nu_t=(C_s\Delta l)^2|\bar{\boldsymbol{S}}|, \end{equation}

where $C_s$ denotes the Smagorinsky constant and the SPS strain is given by

(2.23)\begin{equation} |\bar{\boldsymbol{S}}|=\sqrt{2\bar{S}_{ij}\bar{S}_{ij}}. \end{equation}

Rogallo & Moin (Reference Rogallo and Moin1984) suggested a range for $C_s$, which is between 0.1 and 0.24. In this work, we choose $C_s=0.12$, and the calculated results are in good agreement with the experimental data in the literature. The contribution of SPS stress to linear momentum can be calculated by using the discrete symmetric formulation (see Lo & Shao Reference Lo and Shao2002),

(2.24)\begin{equation} \frac{1}{\rho}\boldsymbol{\nabla}\boldsymbol{\cdot}\bar{\pmb\tau} = \sum_j m_j\left( \frac{\bar{\pmb\tau}_i}{\rho_i^2} + \frac{\bar{\pmb\tau}_j}{\rho_j^2}\right)\boldsymbol{\cdot}\boldsymbol{\nabla}_i W_{ij}, \end{equation}

where the sub-particle stress tensor $\bar {\pmb \tau }_i$ of particle $i$ is calculated by (2.21).

2.4. Multiphase droplet-laden cough jet model

This work is aimed at studying the interaction between saliva droplets and ambient air during coughing. The human cough jet is a mixture of multiphase fluids that is essentially a saliva droplet-laden moist gas flow (Zhu, Kato & Yang Reference Zhu, Kato and Yang2006). When the cough jet is ejected from a human's mouth, the cough jet is composed by both the gas flow and droplets, which can be regarded as a mixture of two phases, i.e. the gas phase and the solute droplet phase. Therefore, the human cough jet flow is a multiphase particle-laden fluid flow, in which the saliva droplets and respiratory air are blending, mixing, penetrating and interacting with each other discontinuously, rather than a multiphase flow with two continuous fluid phases (e.g. Yan et al. Reference Yan, Li, Zhang, Kan and Sun2019, Reference Yan, Li, Kan, Zhang and Lai2022). There has been a vast literature on modelling multiphase flow with droplets or particles, i.e. multiphase particulate flow (Crowe et al. Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2011; Podgorska Reference Podgorska2019).

In this work, we adopt an early multiphase flow model specifically designed for two penetrating material phases (i.e. Harlow & Amsden Reference Harlow and Amsden1975a,Reference Harlow and Amsdenb; Valentine & Wohletz Reference Valentine and Wohletz1989). To capture the expiratory droplet dispersion process, we employed an SPH version of the above multiphase fluid model developed in Monaghan & Kocharyan (Reference Monaghan and Kocharyan1995), which is a Lagrangian-type mesh-free particle method. The Monaghan–Kocharyan multiphase flow model approximates the gas–droplet mixture as two interpenetrating fluids that interact by pressure and drag forces. It is assumed in this work that: (1) the droplets do not evaporate; (2) the gas particles do not condense; and (3) the stress tensor of the droplet gas may be negligible. This is a reasonable approximation because, for mild to moderate cases of COVID-19, the small-sized droplets ($d < 20\ \mathrm {\mu }{\rm m}$), which may evaporate completely in the first 2 s of their emission, are unlikely to be of any consequence in carrying infection (see Anand & Mayya Reference Anand and Mayya2020; Stadnytskyi et al. Reference Stadnytskyi, Bax, Bax and Anfinrud2020), while for large droplets ($d> 100\ \mathrm {\mu }{\rm m}$), their settling time is much shorter than their lifetime or evaporation time (see Chen Reference Chen2020; Zeng et al. Reference Zeng, Chen, Yuan, Yamamoto and Maruyama2021). For medium-sized droplets ($20\ \mathrm {\mu }{\rm m} < d < 100\ \mathrm {\mu }{\rm m}$), the evaporation effect is indeed a concern; however, the evaporation times for most such droplets are beyond the duration of the turbulent state of the coughed air $({\sim }2\ {\rm s})$, and thus the above approximation is practically acceptable.

The local interaction between the two phases can be described as the interaction among gas particles and droplet particles. As shown in figure 1, in each fluid particle compact support, the concentration or the volume fraction of the two phase fluids varies from time to time. Thus, their densities differ from their densities in the homogeneous state of the single phase under near-incompressible conditions. The current density of each phase is denoted as

(2.25a,b)\begin{equation} \hat{\rho}_p =\theta_p \rho_p,\quad p=d,g, \end{equation}

where $\rho _p$ are the equilibrium densities of the homogeneous state of the single phase, and $\theta _p$ is the volume fraction of a given phase in the mixture, which satisfies the condition (see figure 1)

(2.26)\begin{equation} \theta _g + \theta _d=1. \end{equation}

Figure 1.

Schematic illustration of the interaction scheme between gas particles and droplet particles, and the volume fractions of the droplets and ambient air.

We can then derive the continuity equations for the current gas density $\hat {\rho }_g$ and the current droplet density $\hat {\rho }_d$ as follows:

(2.27)\begin{equation} \frac{{\rm d} \hat{\rho}_g }{{\rm d}t} ={-}\hat{\rho}_g \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{v}}_g \end{equation}

and

(2.28)\begin{equation} \frac{ {\rm d} \hat{\rho}_d }{{\rm d}t} ={-}\hat{\rho}_d \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{v}}_d. \end{equation}

Since the equilibrium density $\rho _d$ of the homogeneous droplet state is known, we can determine $\theta _d=\hat {\rho }_d/\rho _d$, after we solve $\hat {\rho }_d$, and subsequently we can then find ${\theta _g=1-\theta _d}$. For a clear presentation, it is convenient to reserve the subscripts $i$ and $j$ for the gas particles, and $a$ and $b$ for the droplet particles. Thus the gas-phase continuity equation becomes

(2.29)\begin{equation} \frac{ {\rm d} \hat{\rho}_i }{{\rm d}t}= \sum _j m_j {\boldsymbol{v}}_{ij} \boldsymbol{\cdot} \boldsymbol{\nabla}_i W_{ij} \end{equation}

and the droplet continuity equation becomes

(2.30)\begin{equation} \frac{ {\rm d} \hat{\rho}_a }{{\rm d}t}= \sum _b m_b {\boldsymbol{v}}_{ab} \boldsymbol{\cdot} \boldsymbol{\nabla}_a W_{ab}. \end{equation}

Subsequently, the SPH momentum equations for gas and droplets are expressed as

(2.31)\begin{equation} \frac{ {\rm d} \hat{\boldsymbol{v}}_g }{{\rm d}t}={-} \frac{\boldsymbol{\nabla} {p}}{\rho _g}+ \frac{K}{\hat{\rho}_g}({\boldsymbol{v}}_d - {\boldsymbol{v}}_g) + \varPi+ \frac{\hat{\rho}_g - \hat{\rho}_{mix} }{\hat{\rho}_g}{\boldsymbol{g}}, \end{equation}

where $\rho _{mix} = \theta _g \hat {\rho }_g + \theta _d \hat {\rho }_d$, and

(2.32)\begin{equation} \frac{ {\rm d} \hat{\boldsymbol{v}}_d }{{\rm d}t}={-} \frac{\boldsymbol{\nabla} {p}}{\rho _d}- \frac{K}{\hat{\rho} _d} ({\boldsymbol{v}}_d - {\boldsymbol{v}}_g) + \frac{ \hat{\rho}_d - \hat{\rho}_{mix} }{\hat{\rho}_d}{\boldsymbol{g}}, \end{equation}

where the term $((\hat {\rho }_d - \hat {\rho }_{mix}) / \hat {\rho }_d){\boldsymbol {g}}$ includes the contribution of buoyancy force, and $K$ is the drag coefficient to determine the magnitude of the drag force. The latter is calculated based on the following formula:

(2.33)\begin{equation} K=\frac{\rho_g\theta_dC_D}{r_d}\Vert\boldsymbol{v}_d-\boldsymbol{v}_g\Vert, \end{equation}

where $\rho _g$ is the mass density of the ambient gas, $r_d$ is the droplet radius and $\theta _d$ is the volume fraction of the droplet phase. The drag coefficient $C_D$ in the above equation is expressed as

(2.34)\begin{equation} C_D=\begin{cases} \left[0.15(Re^*)^{0.687} +1 \right] \dfrac{24}{Re^*}, & Re^*<1000, \\ 0.44, & Re^*\geqslant 1000, \end{cases} \end{equation}

in which

(2.35)\begin{equation} Re^*=\frac{2\rho_a\Vert\boldsymbol{v}_i-\boldsymbol{v}_a\Vert\boldsymbol{r}_i}{\mu_a} \end{equation}

represents the Reynolds number of droplets, which is evaluated for each individual droplet.

As shown in figure 1, when the droplet particles are in the supports of gas particles, the droplet particle will exert a drag force on the gas particles. Inversely, when the gas particles are in the supports of droplet particles, the gas particle will exert a drag force on the droplet particles.

Finally, the balance equation of linear momentum for gas particles interacting with droplet particles can be written as follows:

(2.36)\begin{align} \frac{ {\rm d} \hat{\boldsymbol{v}}_i }{{\rm d}t} &={-} \sum _j \frac{\theta _i {p}_i + \theta _j {p}_j }{\hat{\rho} _i} V_j \boldsymbol{\nabla}_i W_{ij} - \sum _a m_a \frac{{p}_i \theta _a }{\hat{\rho} _i \hat{\rho} _a } \boldsymbol{\nabla} _i W _{ia} + \varGamma_i \nonumber\\ &\quad + \sigma \sum _a m_a {\frac{K _{ia}}{ {\hat{\rho}}_i {\hat{\rho}}_a}} \left({\frac{{\boldsymbol{v}}_{ai} \boldsymbol{\cdot} {\boldsymbol{r}}_{ai}}{{\boldsymbol{r}}_{ai}^2 + {\boldsymbol{\eta}}^2}}\right) {\boldsymbol{r}}_{ai} W_{ai} + \varPi _i + \frac{1}{\hat{\rho}_i}\boldsymbol{\nabla}\boldsymbol{\cdot}\bar{\pmb\tau}_i + \left(\frac{\hat{\rho}_i - \hat{\rho}_{mix} }{\hat{\rho}_i} \right) {\boldsymbol{g}}_i , \end{align}

where

(2.37)$$\begin{gather} \varGamma_i = 4\nu_0\sum_j \frac{m_j(c_i+c_j)(\boldsymbol{v}_j-\boldsymbol{v}_i) (\boldsymbol{r}_j-\boldsymbol{r}_i)} {(\theta_i\rho_i+\theta_j\rho_j)(\Vert\boldsymbol{r}_j-\boldsymbol{r}_i\Vert^2 +0.01h^2)}\varDelta_i W_{ij}, \end{gather}$$

(2.38)$$\begin{gather}\varPi_i = \xi c_0 h_i\sum_j \frac{(\theta_j\rho_j-\theta_i\rho_i)(\boldsymbol{r}_j-\boldsymbol{r}_i)} {\Vert\boldsymbol{r}_j-\boldsymbol{r}_i\Vert^2}\varDelta_i W_{ij}V_j, \end{gather}$$

(2.39)$$\begin{gather}\frac{1}{\hat{\rho}_i}\boldsymbol{\nabla}\boldsymbol{\cdot}\bar{\boldsymbol{\tau}}_i = \sum_j m_j\left(\frac{\bar{\boldsymbol{\tau}}_i}{(\theta_i\rho_i)^2} +\frac{\bar{\boldsymbol{\tau}}_j}{(\theta_j\rho_j)^2}\right). \end{gather}$$

Subsequently, the SPH balance equation of linear momentum for droplet particles interacting with gas particles can be written as follows:

(2.40)\begin{align} \frac{ {\rm d}\hat{\boldsymbol{v}}_a }{{\rm d}t} &={-} \sum _b \frac{\theta _a {p}_a + \theta _b {p}_b }{\hat{\rho} _a} V_b \boldsymbol{\nabla} _a W_{ab} - \sum _i m_i \frac{{p}_i \theta _a }{\hat{\rho} _a \hat{\rho} _i } \boldsymbol{\nabla}_a W _{ai} \nonumber\\ &\quad - \sigma \sum _i m_i {\frac{K _{ai}}{ {\hat{\rho}}_i {\hat{\rho}}_a}} \left({\frac{{\boldsymbol{v}}_{ai} \boldsymbol{\cdot} {\boldsymbol{r}}_{ai}}{{\boldsymbol{r}}_{ai}^2 + {\boldsymbol{\eta}} ^2 }} \right) {\boldsymbol{r}}_{ai} W_{ai}+ \left( \frac{\hat{\rho}_a - \hat{\rho}_{mix} }{\hat{\rho}_a} \right) {\boldsymbol{g}}_a. \end{align}

2.5. Solid-wall boundary condition

Since SPH is a non-local particle method, each particle has a finite-size compact support that contains its interacting neighbouring points. However, for SPH particles near the boundaries, they do not have a complete domain support. This could introduce a non-physical response, and it may cause lower accuracy, even significant errors, in computational results. In this work, the ghost particle technique (Adami, Hu & Adams Reference Adami, Hu and Adams2012) is adopted to approximate the interface between the boundary and the fluid phase, as shown in figure 2. Compared with other techniques to impose a solid-wall boundary condition, such as mirror particles, the main advantage of ghost particles is its simplicity to enforce the boundary conditions for complex boundary geometries, and, once the particles are initialized, the boundary is well described throughout the simulation.

Figure 2.

Fluid particles $(\bullet )$ interact with ghost particles denoting the wall $(\circ )$ to guarantee the whole support for kernel interpolation.

In figure 2, fluid particles $(\bullet )$ near the solid wall will interact with ghost particles representing the wall $(\circ )$ within the support region. When the fluid particles are governed by the equation of motion, the position and velocity of the ghost particles are fixed still. By choosing the wall velocity for viscous interaction, the free-slip or no-slip boundary condition can be imposed at the wall. The free-slip boundary condition is adopted by ignoring the viscosity between the solid boundary and the fluid phase. By imposing the no-slip condition, one can obtain the smoothed velocity field in the fluid phase through the ghost particle velocity (Adami et al. Reference Adami, Hu and Adams2012) as follows:

(2.41)\begin{equation} \tilde{\boldsymbol{v}}_i = \frac{\sum\limits_j\boldsymbol{v}_j W_{ij}}{\sum\limits_j W_{ij}}. \end{equation}

Subsequently, the wall velocity ${\boldsymbol {v}}_w$ is obtained as

(2.42)\begin{equation} \boldsymbol{v}_w = 2\boldsymbol{v}_i - \tilde{\boldsymbol{v}}_i, \end{equation}

where $\boldsymbol {v}_i$ denotes the prescribed solid-wall velocity. Similarly, the wall pressure $p_w$ can also be obtained by the summation of neighbouring fluid particles,

(2.43)\begin{equation} p_w=\frac{\sum\limits_fp_fW_{wf}+(\boldsymbol{g}-\boldsymbol{a}_w) \boldsymbol{\cdot}\sum\limits_f\rho_f\boldsymbol{r}_{wf}W_{wf}}{\sum\limits_f W_{wf}}, \end{equation}

where $\boldsymbol {a}_w$ is the prescribed acceleration of ghost particles and ${\boldsymbol {a}_w} = {\textbf{0}}$ is applied to the fixed boundary. Furthermore, the density of ghost particles can be inversely solved through the Tait equation as

(2.44)\begin{equation} \rho_w=\rho_{0}\left(\frac{p_w}{B}+1\right)^{{1}/{\gamma}}. \end{equation}

2.6. Open-boundary condition

In this work, to model the cough jet inflow as well as the outflow at the other side of the simulation domain, the boundaries of the computational domain need to enforce specific flow conditions. For cough jets, the inlet boundary conditions can be achieved by specifying the inflow velocity for buffer particles in the inlet region, whereas other velocity or pressure conditions need to be specified or imposed at the outlet boundary. To impose these boundary conditions, an open-boundary condition algorithm (Tafuni et al. Reference Tafuni, Domínguez, Vacondio and Crespo2018) is utilized here. The working principle of the open-boundary condition algorithm is briefly summarized in figure 3. The blue circles represent the fluid particles in the computational domain, and then a permeable threshold is set in front of the edge of the domain, which is depicted by a dashed line. Buffer particles are created and maintained during the computation process to be the same as the fluid particles except that their physical information is passed from the fluid region by ghost particles (green squares). The position of the ghost particles is then obtained by mirroring the normal distance of the boundary particles from the permeable surface into the fluid. After that, the physical quantities are interpolated and extracted from the neighbouring fluid particles to the ghost particles; the latter then passes those properties back to their corresponding boundary particles. If a fluid particle passes through the permeable surface, it will be transferred to a buffer particle, which will also cause the mirroring of a ghost particle back into the fluid. If a boundary particle passes through the edge of the domain from the buffer region, it will be removed from the computational domain. Conversely, if a buffer particle travels through the permeable surface, it will become a fluid particle such that its ghost particle will be discarded, and a new boundary particle will be inserted near the other edge of the buffer zone along the flow direction of the buffer particles.

Figure 3.

Schematic illustration of the open-boundary condition.

The physical information of the ghost particles is extracted by the standard particle interpolation using their neighbouring particles that reside within their smooth length, which is expressed as

(2.45)\begin{equation} f^g(r_i)=\frac{\displaystyle\sum f^g(r_j)W^g_{ij}m_j/\rho_j} {\displaystyle\sum W^g_{ij}m_j/\rho_j}, \end{equation}

where $f^g(r_i)$ is the physical quantity at the ghost particle, the summation is performed on the neighbouring fluid particles adjacent to the ghost particles, and $W^g_{ij}$ is the weight function between the ghost particle $i$ and its fluid neighbour $j$.

2.7. Shifting algorithm

The particle spacing in the SPH method is an important factor for the problem of numerical stability, especially in simulations of turbulent flow, because in this case the particles cannot keep a uniform distribution. As a result, noise is introduced in the velocity and pressure fields, and in some cases unrealistic voids could be generated inside the fluid field. To resolve this issue, we adopt a so-called shifting algorithm that can adaptively redistribute SPH particles and make them quasi-uniform. We first assume that the number of particles passing through unit surface per unit time is proportional to the particle velocity, and then the particle moving velocity and subsequently the particle moving distance can be found. Based on the particle concentration, a particle shifting distance $\delta \boldsymbol {r}$ (Lind et al. Reference Lind, Xu, Stansby and Rogers2012) can be obtained as

(2.46)\begin{equation} \delta\boldsymbol{r} ={-}D\boldsymbol{\nabla} C, \end{equation}

where $C$ represents the particle concentration and $D$ denotes the diffusion coefficient, which controls the shifting magnitude. Then, the gradient of particle concentration can be obtained by using the SPH gradient operator:

(2.47)\begin{equation} \boldsymbol{\nabla} C = \sum_j\frac{m_j}{\rho_j}\boldsymbol{\nabla} W_{ij}. \end{equation}

The proportionality coefficient $D$ should be set large enough to provide an effective particle shifting without introducing obvious errors or instabilities. This is achieved by enforcing the von Neumann stability condition for a convection–diffusion equation,

(2.48)\begin{equation} D\leqslant\frac{h^2}{2\Delta t_{max}}, \end{equation}

where $\Delta t_{max}$ represents the maximum time step determined by the CFL condition. Combining the equation above with the CFL condition, one can find the shifting coefficient $D$ by

(2.49)\begin{equation} D=Ah\Vert\boldsymbol{v}\Vert \,{\rm d}t, \end{equation}

where $A$ is a dimensionless constant, which depends on the problem set-up and the specific discretization, and ${\rm d}t$ is the time step. The value of $A$ is in the range of 1–6, and we used $A= 2.0$ in this work.

2.8. Similarity protocol in experimental data

In this work, we shall compare the numerical results with the experimental data obtained by Wei & Li (Reference Wei and Li2017). Since their experiments were carried out by using a water tank and glass beads, a similarity protocol is utilized in their work to quantitatively map the water tank experimental data to the real scenario in air. By neglecting the nonlinearity of the drag force and ignoring the force due to fluid acceleration, the added-mass force and the Basset history force, under the Stokes’ region ($Re_p={| u_f-u_p| d_p}/{\nu }<1$), the equation of particle motion is written as

(2.50)\begin{equation} \frac{{\rm d}\boldsymbol{v}_p}{{\rm d}t} =\frac{1}{St}(\boldsymbol{v}_f-\boldsymbol{v}_p - \boldsymbol{v}_\tau), \end{equation}

where ${\boldsymbol {v}}_f$ and ${\boldsymbol {v}}_p$ are the velocities of the fluid phase and the particle phase, $\boldsymbol {v}_\tau$ represents the particles’ terminal settling velocity, and $\tau$ represents the particles’ relaxation time. They are respectively defined as follows:

(2.51a,b)\begin{equation} \boldsymbol{v}_\tau = \frac{(\rho_p-\rho_f)d^2_p\boldsymbol{g}}{18\mu},\quad \tau=\frac{\rho_pd^2_p}{18\mu}. \end{equation}

Normalizing (2.50) by the characteristic length scale $D$ and characteristic velocity scale $U_c$ will give the following expression:

(2.52)\begin{equation} \frac{{\rm d}\bar{\boldsymbol{v}}_p}{{\rm d}\bar t}=\frac{1}{St}(\bar{\boldsymbol{v}}_f -\bar{\boldsymbol{v}}_p - \bar{\boldsymbol{v}}_\tau), \end{equation}

where $St={U_c\tau }/{D}$ is the Stokes number. The average injection velocity is chosen as the characteristic velocity $U_c$, which is given by

(2.53)\begin{equation} U_c=\frac{1}{t_{inj}}\int_{0}^{t_{inj}}U(t)\,{\rm d}t, \end{equation}

where $t_{inj}$ is the duration of the cough, or the duration of the initiating air-breathing phase. In the water tank modelling, the same Reynolds number, Stokes number and ratio of terminal settling velocity to characteristic velocity scale are used in order to make the mapping consistent, the corresponding formulae being given as

(2.54)$$\begin{gather} Re_c = \frac{U_{c,g}D_g}{v_g}=\frac{U_{c,w}D_w}{v_w}, \end{gather}$$

(2.55)$$\begin{gather}St_{c,g} = St_{c,w} , \end{gather}$$

(2.56)$$\begin{gather}\frac{u_{\tau,g}}{U_{c,g}} = \frac{u_{\tau,w}}{U_{c,w}}. \end{gather}$$

Combining the equations above leads to

(2.57)$$\begin{gather} \frac{U_{c,w}}{U_{c,g}} = \sqrt[3]{\frac{\nu_w}{\nu_g}\frac{\rho_{p,g}}{\rho_{p,g}-\rho_g}\frac{\rho_{p,w} -\rho_w}{\rho_{p,w}}} , \end{gather}$$

(2.58)$$\begin{gather}\frac{D_w}{D_g} = \sqrt[3]{\left(\frac{\nu_w}{\nu_g}\right)^2 \left(\frac{\rho_{p,g}}{\rho_{p,g}-\rho_g}\frac{\rho_{p,w}-\rho_w}{\rho_{p,w}}\right)^{{-}1}}. \end{gather}$$

The above two equations define the geometric and boundary conditions for particle experiments in water tanks, e.g. the nozzle diameter and exit velocity. More details may be found in Wei & Li (Reference Wei and Li2017).

2.9. Numerical modelling and parameters

In this work, large-scale numerical simulations are performed to investigate the cough jet development and droplet dispersion. The computational domain and enlarged details of the inlet are illustrated in figure 4, with its dimension size set as 3.8 m in length, 2.2 m in height and 1.4 m in width, in which the ambient air is contained in an air chamber surrounded by rigid walls. The density and viscosity of the ambient air and jet air in the computational domain are $\rho _g=1.181\ \mathrm {kg}\ \mathrm {m}^{-3}$ and $\nu _g=1.86\times 10^{-5} \mathrm {Pa}\ \mathrm {s}$, respectively. A cough is modelled as the injection of a limited amount of air into the stationary environment. The transient characteristics of the cough flow are achieved by applying the flow velocity at the inlet area as a boundary condition. The domain is initially discretized into 98 872 692 particles including the solid-wall boundaries with the particle spacing $\Delta x=0.005\ \mathrm {m}$. An open inlet is put at one side of the box with 1.6 m in height from the bottom of the box to mimic the real cough position, as suggested by Dbouk & Drikakis (Reference Dbouk and Drikakis2020). The area of the inlet is set as $0.0005\ \mathrm {m}^{2}$, which is approximately composed by a square section with width $r=0.02\ {\rm m}$ and height $d=0.025\ {\rm m}$, and the corresponding flow velocity to be applied as inlet boundary condition can be reproduced from Gupta et al. Reference Gupta, Lin and Chen2009, as illustrated in figure 5.

Figure 4.

Geometrical set-up of the computational domain.

Figure 5.

Velocity profile of the inlet flow.

The lognormal distribution has been used when we dynamically generate droplets in the jet at the injection region following Duguid (Reference Duguid1946), which obeys the following expression:

(2.59)\begin{equation} y=f(x\mid\mu,\sigma)=\frac{1}{x\sigma\sqrt{2{\rm \pi}}} \exp\left(\frac{-(\log x-\mu)^2}{2\sigma^2}\right),\quad \mathrm{for}\ x>0. \end{equation}

The observed statistical distribution for the size of the droplet particles in cough jets may be approximated as a lognormal distribution (Johnson et al. Reference Johnson2011; Van Sciver et al. Reference Van Sciver, Miller and Hertzberg2011; Han, Weng & Huang Reference Han, Weng and Huang2013; Wang, Xu & Huang Reference Wang, Xu and Huang2020). The two parameters for the lognormal distribution equation (2.59) can be calculated as $\mu =-11.1503$ and $\sigma =0.9660$. Droplet injection is achieved in the following steps: (i) The desired number of droplets to be injected in each size range is specified according to the experimental data. (ii) Each time air phase particles are generated and pushed into the inlet, a random drawing is performed to determine whether or not droplet particles should be created. (iii) If new droplet particles are to be injected into the ejection flow, a random drawing is performed again to choose the droplet radius, while the total number of droplets must obey the lognormal distribution. During a fixed ejection period with a fixed number of ejection particles, which is about 125 000 particles in our simulation, and the droplet occurrence probability will ensure that among 125 000 particle ejections there will be 10 000 times when we eject a droplet particle. Then each time that a droplet particle is ejected, we shall draw a random number to determine the size of the droplet particle based on the lognormal distribution.

In order to explore the effects of the injected droplets in the exhaled air and measure the statistics of the data, we performed 10 numerical simulations, whose details are listed in table 1. They are classified into two groups. One is performed with the two-way coupling considered between the air phase and the droplets. In the other group, there are no droplets being injected into the jet flow. All the simulation cases share the same geometry and discretization as well as the boundary conditions. They only vary little by small velocity perturbations at ejection. That is, random perturbations of 1 % amplitude are applied to the original ejection velocity profile to form five different velocity profiles in total, which are then assigned to the inlet boundary condition. Those profiles are then adapted in five pairs of numerical simulations to study the turbulence statistics of the ambient air phase. Evaporation and condensation of the droplets are not considered in this work, as our focus is concentrated on the influence of the droplets on the motion of the turbulent puff, and its possible effects inversely on droplet dynamics and dispersion. As predicted by the constant-temperature model (Wells Reference Wells1934; Langmuir Reference Langmuir1918), the evaporation time for a droplet to be fully evaporated to its dry nucleus is proportional to the initial diameter squared, which is the so-called $d^2$ law. As for the Lagrangian statistics method that is commonly used in one-way coupling simulations, the motion of a droplet is dictated by the drag force from its surrounding gas phase, which varies with the surface area affected by the convective effects.

Table 1.

Details of the numerical simulations. Here $v^*$ is the peak ejection velocity, and $t_{inj}$ is the ejection duration. Simulation cases CVD20V22a–e will have droplets generated in the injection flow, while no droplets are created in cases CVD20N22a–e. The cases in the two groups sharing the same alphabetical indexing share the same ejection velocity profile.

In this regard, one would expect non-evaporating/condensing droplets to travel a shorter distance for drops with small diameters, or a slightly longer distance for those with large diameters. Such estimation is based on the fact (which is suggested by numerical investigations that considered humidity) that the supersaturated region right after the injection flow may help to increase the volume of the droplets for a short time, while a humid puff may help to reduce the evaporation rate of the droplets inside. In addition, the rate of evaporation/condensation is dependent not only on the relative velocity of droplets, but also on the temperature and humidity of the surrounding air. Since our work is mainly concerned with the study of droplet dispersion in indoor environments, i.e. under room temperature, condensation could be neglected even with 90 % relative humidity according to the results suggested by Ng et al. (Reference Ng, Chong, Yang, Li, Verzicco and Lohse2021). By neglecting evaporation, our prediction may give optimized results with respect to the realistic situation.

In the simulations, both respiratory air and droplets will come out from the inlet into the open space based on the prescribed velocity profile. The duration of the cough simulation is 4 s, and at the end of the simulation, the total number of particles, both air and droplet particles, will exceed 125 million.