Three-dimensional simulation of clouds of multi-disperse evaporating saliva droplets in a train cabin

Since the beginning of the COVID-19 pandemic (World Health Organization¹), CFD (computational fluid dynamics) has enjoyed a widespread use for the investigation of the inherent droplet-based mechanisms by which the virus infection can be propagated. Such a practice has flourished due to the intrinsic ability of this branch of computational physics to complement experimental and theoretical fluid dynamics by providing an alternative (extremely effective) means of simulating real processes for conditions often unavailable experimentally or not practically realizable. The intrinsic reasons of such a modus operandi (and its success) can also be directly rooted in the decision of several research groups (having different backgrounds and perspectives) to move from their traditional heartlands of applied engineering or fundamental research to new lines of inquiry aimed at a better understanding of the mechanisms of virus diffusion.

The resulting articulated numerical studies have provided disjointed glimpses of a wide variety of qualitatively and quantitatively different results in widely different parts of parameter space, essentially reflecting current interpretations or beliefs about the possible cause-and-effect relationships driving the pandemic (Mittal et al.²). This peculiar endeavor or framework has been supported by the widespread consensus that the transmission of COVID-19 occurs essentially via virus-laden “droplets,” which originate from the respiratory tract of an infected person and are expelled from the mouth and nose in a variety of circumstances. These include anomalous situations such as coughing or sneezing (where droplets can be formed by saliva and by the mucous coating of the lungs and vocal cords, Zhu et al.;³ Simha and Rao;⁴ Dbouka and Drikakis;⁵ National Research Council;⁶ Fontes et al.;⁷ and Li et al.⁸) and/or even normal activities such as breathing and talking (Jones and Brosseau;⁹ Asadi et al.;¹⁰ Bourouiba;¹¹ CDC;¹² and Smith et al.¹³).

In particular, a first analysis on which many more recent investigations have relied is the original study by Wells,^14,15 where it was shown that the evolution in time of these droplets is governed by the triadic relationship among inertia, gravity, and evaporation. It is worth highlighting that this realization has naturally led to the introduction of a kind of dichotomy in this category of studies, that is, a net distinction between “large” and “small” droplets; namely, cases where the size of droplets is such that the time they take to settle is smaller than that required for their evaporation, or, vice versa, the sedimentation process is so slow that the liquid part of the droplets is entirely consumed before they reach the ground.

These two paradigms should obviously be regarded as opposite extremes, both being severe approximations to a more complete representation of possible populations of droplets created by human beings. The former may be considered representative of practical circumstances where droplets can actually contaminate surfaces located in a certain neighborhood of the source emitting them; the latter may be used to model situations where the droplets can remain suspended in the air for relatively long times (forming an “aerosol”) and be transported accordingly at large distances from the source (Pendar and Pascoa¹⁶). The translational relevance or applicability of these realizations to the effective propagation of the infection can be immediately seen by considering that human-to-human transmission of COVID-19 is thought to occur due to some related (concurrent) mechanisms.

These mechanisms are not mutually exclusive, nor are they truly progressive. Relatively large droplets undergoing fast sedimentation can be deposited on surfaces and be transferred according to other human beings coming intentionally or inadvertently in physical contact with such surfaces (Bhardwaj and Agrawal^17,18). Aerosolized tiny droplets originating from expiratory ejecta can be transported by ambient air currents until they are inhaled by potential recipients.

A careful analysis of the existing literature, however, also indicates that there are reasons to question whether these theories are the only viable contenders for the interpretation or prediction of the infection spreading rate. As an example, the possibility for droplets to reach the nose, mouth, or conjunctiva of another human being should not be regarded as an exclusive prerogative of small droplets. Even large droplets may reach other human beings if expelled with sufficiently high momentum; that is, the “initial conditions” should also be considered as an important aspect of the overall problem (Renzi and Clarke¹⁹).

In this (already complex per se) scenario, evaporation should be seen as an additional influential factor potentially causing the transition from one mechanism to another (as the preferred mode of infection diffusion) depending on environmental conditions. This apparently innocuous observation implies that the problem also strongly depends on the “boundary conditions.” In particular, the rate of evaporation is known to be a function of the droplet surface saturation vapor pressure and the vapor pressure of the surrounding gas, which in turn depends on the degree of humidity (Dbouka and Drikakis²⁰ and Li et al.²¹). Obviously, the evaporation and movement of the droplets (after being expelled) are predominately dictated by their sizes. The evaporation rate also displays a relationship with the mass-diffusion coefficient, which in turn changes according to the droplet-to-ambient temperature difference and the velocity of the droplet with respect to the surrounding gaseous environment (Xie et al.²²). Winter conditions of low temperature and high relative humidity (RH) can cause more droplets to survive over relatively long times, which may be a possible driver of a second pandemic wave in the autumn and winter seasons (Wang et al.²³).

Superimposed on these aspects is the possibility that droplets are turned with time into small “solid residues.” Indeed, evaporation can cause quick consumption (in a few seconds) of small droplets suspended and transported in air currents, forcing them to form nuclei consisting of virions, salts that were previously dissolved in water (Vejerano and Marr²⁴ and Chaudhuri et al.²⁵), various proteins and pathogens in varying concentrations (Xie et al.²²), and a certain (small) percentage of residual water (Mezhericher et al.²⁶). These nuclei can span the range from the micrometric to the millimeter size. They can be suspended in the air for hours and are thought to play an important role in the mechanisms of diffusion of the infection as well (van Doremalen et al.;²⁷ Nicas et al.;²⁸ and Asadi et al.¹⁰). As explained before, unlike large droplets (which always require a significant amount of initial momentum in order to step away from the source), small droplets and/or related nuclei can cover large distances if properly supported by ambient flows such as those produced in indoor environments (i.e., public buildings, hospitals, homes, offices, classrooms, airplanes, trains, subways, buses) as a result of air conditioning systems and related forced flow (Tang et al.;²⁹ Zhang and Chen;³⁰ Li et al.;³¹ Eames et al.;³² Zhao et al.;³³ Balocco and Liò;³⁴ Thatiparti et al.;³⁵ Yang et al.;³⁶ Yu et al.;³⁷ Craven and Settles;³⁸ Licina et al.;³⁹ Abuhegazy et al.;⁴⁰ and Zhou and Ji⁴¹).

At this stage, to put the present work in perspective, we wish to remark that, although an additional theoretical inquiry is needed to elaborate a full accounting of all the basic forces, factors, and physical conditions governing the outbreak diffusion (and efforts are currently in progress along these lines), we now are in a situation where the CFD analysis of all the dynamics illustrated above has reached a sort of maturity, in the sense that the existing numerical techniques are in a position to convey valuable and relevant information concerning most of the questions we may ask on the problem. We wish to highlight as well that, surprisingly, the overwhelming majority of research has focused on relatively simplified configurations where the complexity, intricacy, and tortuosity of the currents transporting droplets have been filtered out in favor of a more fundamental approach (aimed at disentangling the functional relationships governing the underlying dynamics).

A relevant exception is represented by the recent work by Abuhegazy et al.⁴⁰ where the location of the “source” has been found to influence strongly the trajectory and deposition distribution of the exhaled aerosol particles and affect the effectiveness of mitigation measures such as glass barriers. According to such a study, particles larger than 20 μm entirely deposit on the ground, desks, and nearby surfaces in the room, while glass barriers can reduce the aerosol transmission of 1 μm particles by ∼92%. Moreover, by opening windows, the particle exit fraction can be increased by ∼38% compared to the case with closed windows.

In the present work, assuming the state of the art in terms of infection transmission models (that is considering a statistical population of droplets that includes both small and large variants and can evaporate), we tackle the problem considering the intricacies, which are typical of an effective (realistic) environment such as that of a train car. Following a common practice in this kind of studies (all based on hybrid Eulerian–Lagrangian techniques, which, therefore, can be said to unify the study of these subjects at least from a numerical point of view), we address the initial boundary value problem (IBVP) numerically solving the classical three-dimensional Navier–Stokes equations for the environment and additional specific equations to track the motion of particles (and their evolution in terms of temperature and mass).

As stated in the introduction, we consider evaporating droplets in the air. From a physical point of view, therefore, the problem consists essentially of a multiphase gas (including the air and the water vapor resulting from the saliva evaporation process) and the dispersed liquid droplets. We treat the gaseous phase in the framework of a variable density approach, as discussed in Sec. II A.

The governing equations for mass, momentum, and enthalpy in the classical unsteady and compressible (complete) form read

• continuity

  $∂ρ∂t+∇¯·(ρV¯)=0,$

  (1)
• species transport for vapor phase

  $∂(ρYv)∂t+∇¯·(ρYvV¯)=∇¯·(ρDv∇¯Yv)+Sv,$

  (2)
• momentum

  $∂∂tρV¯+∇¯·(ρV¯V¯)+∇¯p=∇¯·[2μ(∇¯V¯)os]+ρa¯+S¯m,$

  (3)
• enthalpy

  $ρDhDt−DpDt=∇¯·(λ∇¯T)+∇¯·(ρ∑ihiDi∇¯Yi)+2μ(∇¯V¯)os:(∇¯V¯)os+ρa¯·V¯+Se,$

  (4)

where

and $ρ, p, V¯, T$⁠, and h are the mixture (fluid) density, pressure, velocity, temperature, and enthalpy, respectively. The vector quantity $a¯$ accounts for the acceleration involved in the present problem: the component along z being the classical gravity $g=−9.81 m/s2$⁠, the other two components along x and y, representing a variation of velocity in the direction of the train motion and a centrifugal contribution in a direction perpendicular to it (due to the curvature of the train trajectory), respectively. The symbol Y_i indicates the non-dimensional mass fraction of the different components (the subscript v stands for vapor). The terms $Sv, Sm$⁠, and S_e account for the exchange of mass, momentum, and energy between the carrier fluid and dispersed droplets, respectively, which implicitly indicates that in the present work the coupling between the gas and the dispersed liquid phase is of a two-way coupling nature (as we will illustrate in detail in Sec. II B where precise expressions for these terms are provided).

Moreover, the density of the different species (only air and pure vapor here) is expressed as $ρi=ρYi$⁠; the symbol D_i denotes the species diffusion coefficient.

Accordingly, the following relationships also apply to the present problem:

where μ, λ, and C_p are the mixture dynamic viscosity, thermal conductivity, and specific heat at constant pressure, respectively. Furthermore, the gas state equation can be cast in compact form as

or in an equivalent way as

generally known as Dalton's law of partial pressures, where R_i and p_i represent the gas constant and the partial pressure of each component.

For the considered problem (relatively small value of the Mach number, i.e., M = 0.06), these equations are too complex and broad in scope. In particular, a well-known issue relating to their numerical integration is the need to keep the time integration step sufficiently small in order to capture properly the so-called acoustic wave propagation scale (a failure in doing so typically resulting in algorithm instability).

A convenient approach to circumvent this bottleneck (on which much commercial software rely, including the STAR-CCM+ v. 2021.1 platform at the root of the present study) is based on the so-called pressure splitting approach and the resulting simplifications that can be implemented in Eqs. (1)–(9) after filtering out the acoustic waves and replacing all the problem variables with equivalent series expansions in terms of the Mach number.

A rigorous justification for this modus operandi can be rooted in the study by Roller and Munz,⁴² where by means of multiple space scale asymptotic analysis it was shown that, in general, “the pressure” of a fluid can be decomposed into three parts with different physical meanings, these accounting separately for thermodynamic effects, acoustic wave propagation, and the balance of forces (dynamic pressure).

In such a context, Müller⁴³ could show that acoustics removal from the equations directly leads (with proper mathematical developments) to the so-called low-Mach-number equations, which allow for large temperature and density changes as opposed to the standard Boussinesq equations. Such mathematical developments consist of expanding in power series law of the small parameter $M2≪1$ (where, obviously, M is the reference Mach number) all the primitive variables as follows:

The reader specifically interested in this procedure may consider Beccantini et al.⁴⁴ and Benteboula and Lauriat.⁴⁵ Here, we limit ourselves to mentioning that after substituting the above expansions in the fully compressible Navier–Stokes equations, the lowest order terms in M² are collected. At the order −1, the momentum equation reduces to

Moreover, at the order zero, the following low-Mach-number governing equations are obtained for mass and momentum:

with the state equation

being required to determine the density $ρ0*(r,t)$⁠.

These simplified equations are particularly useful as they make immediately clear that in the framework of the low-Mach-number approximation, pressure can be articulated into two components only: a thermodynamic pressure p₀ homogenous in space and allowed to vary in time, and a dynamic pressure p₁ decoupled from density and temperature fluctuations. The most remarkable implication of such a decomposition is that the dynamic pressure can be determined numerically using an approach similar to that traditionally implemented for incompressible flows in the framework of pressure-based methods; that is, the computation of the velocity can be split into three main steps. In the following, in order to illustrate this approach, we omit the asterisk and the subscripts related to the asymptotic expansions, with the zero- and first-order pressure contributions p₀ and p₁ being simply indicated for clarity as P and p′, respectively.

In an initial stage, an intermediate momentum field is computed solving an incomplete version of the momentum equation (deprived of the gradient of dynamic pressure)

In a second stage, the intermediate momentum field is formally corrected as

Where the dynamic pressure comes from the solution of an elliptic equation for pressure obtained by forcing Eq. (16) in Eq. (12),

The time derivative of the density appearing in this equation, in turn, is computed separately determining the density at different instants through the gas state equation (where the thermodynamic pressure P appears), that is,

Problem closure finally requires that the thermodynamic (constant in space) pressure is determined through a global balance of mass across the entire computational domain, that is,

for a finite-size domain (where the dotted m denotes the incoming or outgoing mass flow rate) or

for an unbounded domain.

As correctly reported by Munz et al.,⁴⁶ in general, these methods are relatively robust.

In addition to the strategies described in Sec. II A 1, a specific model is also used in this work to account for turbulence effects typically associated with the examined value of the Reynolds number (in the range between 4600 and 16 500). Given the nature of the considered carrier flow (it being steady “in mean”), in particular, a method pertaining to the general class of RANS-based technique is used (where the acronym RANS stands for Reynolds Averaged Navier–Stokes Equations).

Such a category of methods generally relies on the idea that the velocity of the fluid at any point can be decomposed into two contributions: one obtained as a time-averaged value and another formally representing an instantaneous fluctuation with respect to the mean value, that is, $V¯=V¯timeaveraged+V¯′$⁠, where obviously $V¯′=V¯−V¯timeaveraged$ and $(V¯′)timeaveraged=V¯timeaveraged−V¯timeaveraged=0$⁠.

Substituting these expressions into the original governing equations and taking a time (or ensemble) average yields the aforementioned (RANS) equations. A new term appears at the right-hand side of the momentum equation inside the divergence term, which formally plays the role of an additional stress in the flow. The extra term $−V¯′V¯′$ is known as “Reynolds stress tensor.”

Closure of the problem is then obtained by introducing the concept of turbulent viscosity by which the extra stress tensor is mathematically put in connection with the gradient of time-averaged velocity. The turbulent viscosity, in turn, is expressed as a function of the turbulent kinetic energy and the dissipation, which require two additional balance equations (to be solved together with the balance equations for the considered species, mass, momentum, and energy). In addition to the Reynolds stress tensor, other modifications are also implemented accordingly in the species and energy equations through the concepts of turbulent Prandtl and Schmidt numbers. Indeed, also the temperature and the vapor concentration can be split into a time-averaged contribution and a fluctuating part. Substitution of them into the original energy and species equations yields an additional term for each equation, which can formally be expressed as the product of a turbulent transport property and the gradient of the time-averaged version of the considered variable. In turn, these turbulent transport properties can be determined as a function of the turbulent kinematic viscosity through the concepts of turbulent Prandtl and Schmidt numbers (both assumed to be equal to 0.9 in the present work). In particular, here the Realizable $k−ε$ turbulence model (Shih et al.⁴⁷) has been adopted.

As anticipated in the introduction, we treat saliva droplets as Lagrangian particles. As the considered number of droplets is relatively limited and their size is small, interactions between the particles (collisions/agglomerations) are neglected; the presence of other droplets, or solid particles, that may be dispersed in the air is also neglected. Accordingly, each droplet is tracked individually throughout the computational domain. Although this approach may look less convenient than other methods where all coexisting phases are dealt with in the framework of a single Eulerian treatment (typically based on the introduction of a volume of fraction variable or similar concepts, see, e.g., Capobianchi et al.⁴⁸ and Lappa⁴⁹), the hybrid Eulerian–Lagrangian has distinct advantages, which make it particularly suitable (see, e.g., Capobianchi and Lappa;⁵⁰ Lappa and Burel;⁵¹ and Lappa⁵²) for the analysis of the problems like that being addressed in the present work. For each droplet, in particular, a set of 3 differential equations are solved, which describe the evolution of its position (and velocity), mass, and temperature, respectively. The additional details on the different terms appearing in these equations are provided in an ordered fashion in Subsections II B 1–II B 6.

When the conservation equation of momentum for a particle is written in the Lagrangian framework, the change in momentum is balanced by surface and body forces that act on the particle as expressed by the right-hand side of Eq. (20) (where n is the number of forces and k is an index used to indicate the generic force)

The drag force can be expressed as follows:

where $V¯s=V¯−V¯p$ is the particle slip velocity and the drag coefficient C_d is a function of the small-scale flow features around the individual particles, which can be derived with the Schiller–Naumann correlation (Schiller and Naumann⁵³).

The pressure gradient force reads

where v_p is the volume of the particle and $∇¯pstatic$ is the gradient of the static pressure in the continuous phase.

The expression for the virtual mass force can be cast in compact form as follows:

where C_vm is the virtual mass coefficient and the operator D/Dt denotes the substantial derivative.

At this stage, it should also be pointed out that, in addition to the standard forces described before, a particle in turbulent flow can also experience what is generally known as turbulent dispersion, or simply “particle dispersion.” This phenomenon is a natural consequence of turbulence itself and the extra forces that it creates at the microscopic level as a result of turbulent momentum exchange. These tend to spread discrete solid particles or droplets exhibiting inertia and a mean relative fluid–particle velocity due to gravity. Lagrangian stochastic (LS) models traditionally used to account for this process can be split into two main categories, namely, the eddy interaction model (EIM) and the random flight or walk model (RWM). They can be distinguished essentially according to the strategy implemented to statistically generate the turbulent fluid velocity in the particle surrounding [this being necessary to solve Eq. (21), i.e., the Lagrangian equation of particle motion, Huilier⁵⁴].

In particular, here we follow Gosman and Ioannides,⁵⁵ where this is achieved assuming that any particle passes through a sequence of turbulent eddies. The underlying concept is that the interactions between a particle and a succession of fluid eddies can be characterized by three parameters, namely, an eddy instantaneous velocity, an eddy lifetime, and an eddy size (these being functions of the considered flow Reynolds number). A Monte Carlo (MC) process is typically associated with this process. The turbulent velocity is sampled randomly from a Gaussian probability distribution function with a standard deviation, and it is kept constant for a given time of the order of the eddy lifetime. The reader interested in a more complete description of this approach, and the related assumptions and the related logical sequence of steps may consider the very recent review by Huilier.⁵⁴ 

Following a common practice in the literature, in order to mimic properly saliva droplets emitted by human beings, dispersed droplets are assumed to be composed of a water-NaCl solution. We wish to remark that, obviously, in the framework of this approximation, mucus and other viscoelastic substances potentially present in the saliva of a sick person are not taken into account. Although these substances may have an impact on the evaporation process, unfortunately reliable data on their percentages and related influence on evaporation are not available yet. The water-NaCl solution-based evaporation model traditionally assumed for saliva droplets can be found in Chaudhuri et al.,²⁵ and it will be briefly outlined in this subsection. In particular, here we limit ourselves to presenting the related governing equations where the subscripts 1, 2, and 3 denote water, air, and salt, respectively. The droplet water mass change due to evaporation can be expressed as

ρ_v is the density of the water vapor, D_v is the binary diffusivity of water vapor in air, R_s is the droplet radius, and α_g is the thermal diffusivity of the air. $BM=(Y1,s−Y1,∞)/(1−Y1,s)$ and $BT=Cp,l(Ts−T∞)/hfg$ are the Spalding mass transfer and heat transfer numbers. Y₁ is the mass fraction of water vapor, while subscripts s and $∞$ indicate the location (droplet surface and far field). $Cp,l$ and h_fg are the specific heat and specific latent heat of vaporization of the droplet liquid. Using the Raoult's law, the vapor pressure at the droplet surface for the binary solution can be determined as

where $χ1,s=1−χ3,s$ is the mole fraction at the droplet surface in the liquid phase. Considering the effect of Raoult's law and relative humidity, the vapor concentration at the droplet surface and at the far field can therefore be obtained as

where M₁ and M₂ are the molecular weights of water and air.

As for the droplet energy balance, it has been verified that the thermal gradient in the liquid phase is rather small for the conditions considered in the present work, and thus, it can be formulated as

where T_s is the droplet temperature, $ml=(4/3)πρlRs3$ and $As=4πRs2$ are the droplet mass and surface area, ρ_l and e_l are the density and specific internal energy of the mixture of salt and water, and k_g is the conductivity of the air, respectively.

Yet, in line with the existing literature and very recent efforts on this subject, we also take into account that if the concentration of solute inside the droplet exceeds a given threshold (due to water evaporation), crystallization phenomena are enabled (Ranz and Marshall⁵⁶). Considering that $Pvap(Ts,χ1,s)$ is also a function of the salt concentration in the droplet, these phenomena can be modeled through a dedicated (droplet-related) species balance equation, that is,

where Y₃ is the dissolved salt mass fraction and $ṁ3,out$ is the solute mass that leaves the solution due to crystallization. When the supersaturation ratio $S=Y3/Y3,c$⁠, with $Y3,c=0.393$ (Gregson et al.⁵⁷), is greater than 1, crystallization begins. In the present work, the growth rate of the crystal is modeled using the simplified rate equation (Naillon et al.⁵⁸ and Derluyn et al.⁵⁹)

where l is the crystal radius, $Ccr=1.14×104 m/s, Ea=58 180 J/mol$⁠, and g_cr = 1. The rate of change of the crystal mass can therefore be finally determined as⁵⁸ 

At this stage, we are in a condition to define precisely the coupling terms appearing in Eqs. (2)–(4) (by which the dispersed droplets can exert a back influence on the carrier flow). These read

where $δΩ$ denotes the volume of the generic control volume, n is the number of droplets contained in that volume, and k is an index used to indicate the generic particle. The sign minus in front of each summation obviously follows from the nature of the interphase exchange itself (an increase in the vapor concentration, carrier fluid momentum, or enthalpy must obviously correspond to a decrease in the droplet liquid mass, a decrease in the droplet momentum, and a shrinkage of the droplet temperature, respectively, and vice versa).

The validation process has been articulated in three stages. As a first step of this hierarchy, the case of a static evaporating droplet has been considered. Then, dynamic conditions (moving droplets) have been examined for both laminar and turbulent flows.

The evaporation and crystallization models illustrated in Sec. II have been implemented in the solver Star-CCM+ (see Polizio⁶⁰) and have been correlated with the experimental data by Chaudhuri et al.²⁵ The experimental setup involved a single particle suspended on an ultrasonic levitator (tec5), for which the diameter change was recorded. To mimic this case, in the numerical simulation, a cylindrical domain of 75 mm base diameter and 120 mm height has been generated, see Fig. 1. The boundary conditions are adiabatic wall on the side and pressure outlet at ambient pressure on the top and bottom faces. A polyhedral mesh of 3300 cells has been used to discretize this domain. The following environmental conditions have been assumed:

• Air temperature, 30 °C

• Relative humidity, 50%

• Air pressure, 1 atm

• No gravity effect

In order to reproduce the experiments, a single droplet has been inserted in the middle of the numerical domain, with no initial velocity. The droplet initial diameter is $338 μm$ and is composed of a water/salt solution (1% w/w). The evolution of the ratio between droplet diameter and droplet initial diameter is shown in Fig. 2.

As the reader will easily realize by inspecting this figure, the numerical results display good agreement with the experimental data. In particular, at about 130 s the equilibrium condition between droplet and air is reached, and evaporation becomes relatively weak or almost negligible. Then, crystallization occurs and the droplet diameter decreases again.

For the dynamic conditions (moving droplets), a slightly different validation strategy has been implemented. The results obtained with two different types of commercial CFD software Fluent and Star-ccm+ have been compared for a fixed (relevant) reference case.

In particular, the geometry shown in Fig. 3(a) has been implemented, discretized with a cell base size of 40 mm, locally refined around the mouth and torso areas, and simulated with the aforementioned commercial platforms (both relying on a time implicit formulation; moreover, for both cases, the convective terms have been discretized using a second-order upwind scheme). Further details are reported in the following.

As evident in the figure, the computational domain of length 6 m, width 2 m, and height 3 m includes a mannequin on the symmetry plane y = 0 at the abscissa $x=2 m$⁠. Following Dbouk and Drikakis,⁵ the mannequin average mouth print has been assumed to have a rectangular shape with an aspect ratio of $Lm/Hm=8.26$ and length L_m = 4 cm [see Fig. 3(b)].

The entire sidewall is treated as a symmetry plane, while the top of the room is a pressure outlet boundary. The initial temperature of the room is $T=20 °C$⁠, there is no air motion, and relative humidity is assumed to be initially zero. Droplet material is pure water.

At the initial time, a total number of 4000 particles are assumed to enter the domain through the mouth boundary shown in Fig. 3(b) with a total mass of 1.07 μg. The particle initial condition can be gathered from Table I.

Xie et al.⁶¹ conducted experimental measurements and quantified exhaled droplet's mass and size due to talking and coughing. As reported by Dbouk and Drikakis,⁵ the Rosin–Rammler distribution law, also known as Weibull distribution, with mean diameter 80 μm can be used to reproduce satisfactorily experimental measurements. In particular, here we consider an initial droplet diameter of 80 μm.

In the present subsection, the results concerning the numerical analysis of the laminar case are given in terms of the following properties of the Lagrangian phase:

• mass evaporation;

• evolution of center of gravity (CoG) coordinates.

Tables II and III summarize the numerical data and the relative percentage errors related to mass and kinematic parameters obtained performing the same droplet simulation with the two aforementioned commercial software. As also quantitatively substantiated by Figs. 4 and 5, the agreement is excellent.

The agreement also holds in terms of the shape of the particle cloud, see Fig. 6.

The adopted turbulence model is the standard $k−ε$ for both computational platforms; moreover, turbulence dispersion acts on the discrete phase. Following the same approach undertaken for the laminar case, the numerical data and relative percentage errors have been collected in Tables IV and V. The trend of the monitored parameters is shown in Figs. 7 and 8. It can be seen that the agreement between the two platforms is excellent also in this case and holds even in terms of the shape of the particle cloud, see Fig. 9.

As anticipated in the introduction, to fill a gap in the literature, in the present work we consider an interregional train passenger compartment. In order to make the outcomes of the numerical study “realistic” as much as possible, the various components and sub-systems of such a compartment are modeled in detail (as further illustrated in Subsections IV A–IV C).

Air entering the passenger cabin must be conditioned to ensure a certain degree of thermal comfort for the occupants. In general, thermal comfort can be assessed in terms of existing EN standard requirements on air velocity and temperature distribution. For the case being, the applicable standard is the EN13129.⁶² The general working principle of an HVAC (heating, ventilation, and air conditioning) system is illustrated in Fig. 10. The key elements of an HVAC system are as follows:

• INLET: air mass flow rate of fresh air coming from external ambient toward a mixing area;

• RECIRCULATION: air mass flow rate coming from the conditioned cabin toward a mixing area;

• OUTLET: air mass flow coming from conditioned cabin driven out to the ambient and replaced by fresh air;

• MIXER: it combines recirculation and fresh air;

• HVAC system: provides the mixed air with the correct characteristics (temperature and relative humidity) before entering the passenger coach;

• FAN: is the component driving the treated air into the cabin;

• CONDITIONED CABIN: is the monitored thermal comfort area.

In order to define the main characteristics of an HVAC system, an evaluation of the thermal load in the conditioned cabin has to be performed. The thermal power is composed mainly of sensible and latent heat as follows:

• SENSIBLE HEAT (Q_sens), thermal load due to temperature changes:

  • –

    Conductive heat;

  • –

    external power exchange;

  • –

    human sensible power production;

  • –

    lighting/system power;

  • –

    thermal leakage;

  • –

    refresh air.

• LATENT HEAT (Q_lat), thermal load due to humidity changes:

  • –

    people production;

  • –

    relative humidity changes due to electrical components;

  • –

    evaporation/condensing process;

  • –

    fresh air.

In particular, the power required by the HVAC system can be evaluated by applying the heat balance equation to the cabin control volume sketched in Fig. 11; this reads

where $Qsens,e$ is the sensible heat exchanged with the external ambient, $Qsens/lat,p$ is the sensible/latent heat produced by human beings, $Qsens/lat,el$ is the sensible/latent heat produced by electrical components, and $Qsens/lat,inf$ is the sensible/latent heat caused by air infiltration.

For the particular case sketched in Fig. 12, the conditioned air enters the cabin, in which all comfort parameters have to be satisfied. Then, the air is forced to pass through the grid opening and enters the vestibules area; here, the volume flow rate is split in recirculated and exhaust streams. The recirculated air enters the recirculation ducts where it is mixed with fresh air and, after being treated again by the air conditioning system to meet comfort parameters, it enters the cabin comfort area, while exhaust air is released in the external environment.

The standard EN13129 imposes quality limits for air temperature and relative humidity. As shown in Fig. 13, air relative humidity must not exceed 65% and 90% for first class quality and second class quality limits, respectively. In general, both quality limits result in a decreasing acceptable level of air relative humidity with increasing air temperature.

Other important parameters relating to the HVAC system design imposed by the EN13129 standard are typically given in terms of external design conditions and internal temperature. Typical values for a temperate zone country are reported in Table VI.

At this stage, we wish also to highlight that a “summer scenario” is not as critical as the winter one in terms of potential droplet propagation. Indeed, the droplet evaporation time is known to decrease drastically when the air temperature becomes higher. For this reason, only winter conditions are assessed in the present work as critical circumstances for droplet evolution (the reader being referred to the point B1 in Fig. 13).

As shown in Fig. 14, the computational domain for the considered interregional train passenger compartment consists of four different regions. The air conditioner placed on the top of the passenger compartment treats air that enters the computational domain through red-colored ducts. The central portion represents the passenger compartment area in which up to 80 passengers can seat, while the two blue-colored extremities correspond to the front and back vestibules areas.

The considered volumetric flow rate of treated air (⁠$V̇$⁠) and its corresponding breakdown (with regard to the regions described before) is reported in Table VII. As already explained to a certain extent in Sec. IV B, only a critical situation is investigated; that is, typical winter conditions at low temperature and high relative humidity within the comfort zone indicated by EN13129 standard and reported in Fig. 13 (point B1, see also Table VIII). Moreover, the volumetric rate of flow recirculation and exhaust airflow in each vestibule is reported in Table IX.

The grid for the computational domain has been chosen after the following considerations. Three different types of grids have been considered under the same boundary conditions defined in Sec. IV C 1. The characteristic size of such coarse, middle, and fine grids is indicated in Table X (see also Fig. 15).

The outcomes of such a comparative study are reported in terms of summation of the instantaneous velocity field on local cell faces, measured in the probe positions, indicated in Fig. 16 by red points P1, P2, P3, P4, and P5, and velocity field contour on the plane section A (see again Fig. 16).

We wish to remark that, as considering the unsteady fluctuation of the velocity field is more relevant from the point of view of particle dynamics, the mesh sensitivity analysis has been based on this unsteady quantity. In particular, the aforementioned velocity summation has been recorded over a simulation time of 40 s for the just mentioned probes P1, P2, P3, P4, and P5 (the results have been reported in Figs. 17–21, respectively).

As qualitatively and quantitatively substantiated by these figures, the fluctuation of velocity in time is similar for the medium and fine grid in all the considered positions, in terms of both amplitude and frequency (with the exception of the coarse grid, an analogous distribution of peaks and valleys can be seen when the other two grids are used).

As a further demonstration of this outcome, Fig. 22 shows the mean velocity field in the plane A of Fig. 16. The significance of these figures resides in their ability to make evident that all the examined grids are able to capture the counter-rotating vortices in the upper part of the passenger cabin (in correspondence with the luggage compartment). It can be seen that the patterning behavior is essentially identical for the medium and fine grids, whereas a change in the symmetry of the main plume can be noticed when the coarse mesh is used.

On the basis of these results and related arguments, the middle mesh is therefore used for the simulations presented in Sec. IV D and IV E.

In particular, the final mesh is reported in Fig. 23. As a concluding remark for this section, we wish to highlight that, in order to capture properly boundary-layer effects, a series of 10 prism layers has been implemented, with a first layer thickness such that the wall $y+$ function is in the acceptable range of [30:80].

Droplets are released according to the initial conditions reported in Tables XI and XII for coughing and talking cases, respectively. All the passengers are assumed to hold fixed positions, that is, to be static. Moreover, four different cases are analyzed, namely:

• Case A—single passenger coughing;

• Case B—single passenger talking for 0.1 and 5 s;

• Case C—6 passengers talking;

• Case D—10 passengers talking.

For each circumstance, the released number of particles is set according to the indications provided by Chao et al.⁶³ and Stadnytskyi et al.⁶⁴ Similarly, the droplet initial diameter has been defined according to the data reported in the experimental study conducted by Xie et al.⁶¹ where a tabular CDF (cumulative distribution function) was given for the initial diameter distribution of droplets in such situations. Other relevant experimental studies (Kwon et al.⁶⁵) have shown that the cone angle of droplet ejection depends upon the gender of the person (the values used for this study are typical of the male gender).

For case B, two fundamental situations are considered, namely, B1 and B2. In the first, a talking passenger releases 1900 particles (corresponding to the amount of particles usually emitted in 1 s) “impulsively,” that is, in the same time that would be required to release them in a coughing event (0.1 s). In the second case, 2000 particles per second are released over 5 s. These situations obviously correspond to two limiting conditions, namely, a single short sentence such as “how are you?” and a more involved speech (including 5 or more complete sentences).

Droplets are assumed to enter the domain through a cone injector whose position is illustrated in Fig. 24. Specifically, the cone injector axis of case A, Fig. 24(a), has a been rotated by $20°$ and $10°$ with respect to the xz and xy planes, respectively, and the cone injector of case B, Fig. 24(b), is parallel to the x axis, while for cases C and D, Figs. 24(c) and 24(d), we have considered 6 and 10 cone injectors, respectively, parallel to x axis and directed in the speaking direction of each passenger.

The outcomes of the numerical simulations in terms of Eulerian velocity field for the conditions defined in Sec. IV C 1 are reported in Fig. 25 in terms of streamlines.

As explained before, the air enters the passenger cabin from the ceiling ducts (blue streamlines in Fig. 25), each one consisting of two rows of ventilation inlets: an external one in correspondence with the luggage compartment and an internal one located on the aisle side.

It can be seen that the airflow spans the entire passenger compartment and reaches the grille opening. At this stage, mass flow is symmetrically partitioned into two streams: 50% going in the front and 50% in the rear vestibule, respectively.

As already explained to a certain extent in Sec. IV A, airflow is then split further into a recirculation flow (green streamlines in Fig. 25) and an exhaust (red streamlines in Fig. 25) current, the first is forced in the top recirculation duct, and the second is expelled through the side ducts near the floor. Some additional insights into the flow field can be gathered from Fig. 26. EN standard requirements prescribe a maximum value for the mean velocity as a function of mean cabin temperature, to be satisfied at specific passenger positions. Figure 26 shows the behavior at three distinct heights (horizontal planes): z = 0.1, z = 0.6, and z = 1.1 m. For a cabin mean temperature of 23 °C, the maximum airspeed allowed is 0.4 m/s. As the reader will realize by inspecting this figure, on the horizontal plane at height z = 0.1 m, the mean velocity field is under the target value with the two areas near the vestibules reaching the target, in correspondence of the aisle. On the other planes, at heights z = 0.6 and z = 1.1 m, respectively, there are some areas exceeding the maximum allowed value, indicated by black lines in Figs. 26(b) and 26(c). These simple observations lead to the conclusion that the considered HVAC system satisfies the EN standard requirements except for some localized areas (which is normal). As we will show in Subsections IV E 2–IV E 5, the effective topology of the flow and its magnitude can have an influence on the particle motion.

Here, we examine the coughing case. In particular, Figs. 27 and 28 account for the time evolution of the droplet total mass and number (this figure is instrumental in showing that the evaporation process starts as soon as particles are injected into the domain).

As a fleeting glimpse into this figure would confirm, in the first 2 s, the droplet evaporation rate is very high; this stage is followed by a new phase where the rate undergoes a shrinkage due to salt concentration effects. From a quantitative standpoint, in the first 6 s about 500 droplets evaporate, while other 200 droplets are lost in the last 8 s.

These data are complemented by Fig. 29, where we have reported the droplet diameter contour at different times. Interestingly, up to a simulation time of t = 0.5 s, the ejected particle maintains the cone shape, the bigger particles being the ones that travel farther. After some time, the cloud of particles starts to feel the surrounding Eulerian (cabin) velocity field and, accordingly, undergoes some deformation. In particular, the transverse size of the cloud is maximized for t = 5 s [maximum spatial diffusion time, see the top and side views in Figs. 29(j) and 29(k), respectively].

Such figures are particularly useful as they can be used to get quantitative information on the particle spreading process (covering 6 m and affecting 6 rows of seats). For this case (with the hypotheses defined in Sec. IV D), the droplets would hit directly only 7 passengers.

The companion (talking) case is depicted in the present section. As explained before, for this case two different scenarios have been simulated (see also Table XII):

• Case B1—1900 droplets ejected impulsively (short sentence such as “hello guys” or “how are you?”);

• Case B2—2000 droplets per second ejected over 5 s (continuous speech including several complete sentences).

The time evolution of droplets number is reported in Figs. 30 and 31 for cases B1 and B2, respectively. It is worth noticing that for the latter case the total number of particles never achieves the expected value of 10 000 due to the evaporation process, which (as time passes) causes some droplets (the smaller ones) to disappear quickly. Notably, in both situations, a rapid decrease in the droplet number can be seen just after the completion of the injection (droplet release) phase. This is yet due to the evaporation of the droplets with smaller size. After this stage (which lasts approximately 5 or 6 s regardless of whether a short sentence or continuous speech is considered), a less steep branch is obtained, which physically corresponds to the timeframe in which slower evaporation of bigger droplets occurs.

Interestingly, as evident in Fig. 32, if the non-dimensional mass is reported as a function of time, the decrease rate (i.e., the angular coefficient of the straight lines mimicking the quasi linear behavior of the considered trend) is independent from the injected mass (see Table XIII).

An explanation/justification for this finding can be elaborated in its simplest form on the basis of the argument that, regardless of the number of droplets being released in the cabin (1900 and 10 000 for cases B1 and B2, respectively), the amount of water evaporating is so small (in terms of volume fraction) that it is not able to change the air humidity content appreciably [on which the evaporation rate depends, as mathematically expressed by Eq. (25)]. Another remarkable implication of this observation is that the maximum distance the droplets can travel away from the source emitting them does not depend on their number (as substantiated by the 3D views in Fig. 33); rather, it is a function of the droplet diameter (distribution) only and of the ambient temperature.

This is the reason for which cases C and D being discussed in Secs. IV E 4 and IV E 5, respectively, are simulated considering the impulsive injection condition (with 1900 droplets) only (this condition leading to notable computational savings and being sufficient to evaluate the maximum possible extension attained by the cloud).

Before moving to the cases with multiple sources, however, meaningful insights also follow from a comparison of Fig. 33 with the analogous ones for the coughing testbed Fig. 29.

Although the behaviors of particle total mass and number are relatively similar in terms of trends and amplitudes, for the talking case, the droplet evolution in space is seemingly much more limited and, accordingly, the passenger exposition risk is lower.

Notable differences can also be identified in the morphology of the cloud shape and the related spatiotemporal evolution.

As a concluding remark for this section, we wish to report that, for the sake of completeness, we have also simulated the case in which the train undergoes a variation of velocity along its main direction of motion (i.e., $ax≠0$⁠, a_y = 0, all the cases presented previously being obtained for $ax=ay=0, az=−g$⁠). More precisely, we have examined the interesting scenario corresponding to an emergency braking of the considered interregional train (V_max = 160 km/h). For such a situation, a typical deceleration value is of the order of 1.3 $m/s2$⁠. In order to mimic these circumstances, $ax=−1.3 m/s2$ has been considered in the governing equations in addition to a_y = 0 and $az=−g$ and the case B1 has been simulated one more time. The related results, along with the companion case B1 without deceleration, are reported in Fig. 34. As witnessed by this figure, the effect of the considered inertial force can be considered almost negligible, which explains why the remaining cases (see Secs. IV E 4 and IV E 5) have been yet simulated assuming that the train moves at a constant speed.

A more involved scenario is considered in this subsection where the number of talking heads is increased to six [uniformly distributed in the cabin as shown in Fig. 24(c)]. The corresponding time evolution of the droplets is reported in Fig. 35 at different simulation times; in this figure, we have also indicated the possible “area of influence” related to each talking passenger, clearly showing that 6 talking passengers may infect up to 9 passengers (in a statistical sense, 7.5% of the passengers may infect up to 11% of the remaining passengers).

In view of results reported in the above two subsections, here we discuss the even more realistic scenario in which a series of 10 talking passengers are disposed (evenly spaced) as shown in Fig. 24(d). The total mass and number of droplets are 10 times the corresponding values for the single talking passenger case; moreover, the related temporal trends (not shown) are rather similar to those already discussed for case B1.

Figure 36 shows the time evolution of the droplets, specifically each different color is representative of a specific position of the considered talking passenger. In this case (for the sake of brevity), only a global view of the complete scenario inside the train cabin is reported (diagonal and top views only).

In line with the findings already illustrated in Subsections IV E 2 and IV E 3, these results confirm that for a talking passenger near the centerline of the train cabin, the maximum droplet diffusion area (or “generated risk” area) is attained after 5 s after the emission of the droplets. Interestingly, however, a completely different situation can be noticed with regard to the passengers located near the vestibules and at intermediate positions.

This can be appreciated in Fig. 37 where dark lines are used to highlight the aforementioned “areas of influence.” In particular, for the passengers near the vestibules, the cloud shape (or risk region) tends to be stretched toward the vestibules following the flow going in that direction. As a result of the acceleration undergone by the air (as it moves toward the vestibules), the maximum extension of the cloud of particles is no longer attained after t = 5 s (rather it changes with time continuously).

Other meaningful information, following natural from an inspection of Fig. 38, concerns the evident relatively high risk, which passenger temporarily located inside the vestibules would be exposed to. In those regions, the percentage of recirculating flow is relatively high (42% of total volumetric flow rate, as reported in Table IX). As a matter of fact, particles released from nearby passengers reach these areas in just a few seconds.

As a concluding remark, it is worth highlighting that this simulation also provides potentially useful statistical data (to be considered together with those already obtained for the 6-talking-heads case). Ten passengers (evenly spaced inside the cabin) talking without a face mask can produce droplets able to hit about 20 passengers (12% of passengers could theoretically infect up to 25% of the remaining cabin occupants).

The diffusion of the COVID-19 pandemic depends on a multitude of influential factors, many of which are of a purely thermo-dynamic or fluid-dynamic nature. These, in turn, require different levels of analysis, which range from the study of the intrinsic (fundamental) physical mechanisms by which nature operates to more practical aspects and the specific intricacies connected with the complexity of the environment where human beings leave and operate.

In the present study, the dispersion of evaporating saliva droplets into a train cabin has been analyzed from both the traditional coarse-grained Eulerian (i.e., continuum) perspective and from a fine-grained micromechanical level in which all the saliva droplets have been tracked individually together with the related content of water, momentum, and energy. Special care has been devoted to the description of the required mathematical models and numerical methods, with the explicit intent to create a theoretical framework on which other future studies may rely. In doing so, the state of the art has been considered in terms of existing paradigms that have already proven to successfully deal with most (if not all) of the above-mentioned influential factors. Following experimental studies on the subject, dilute dispersions of multi-dispersed droplets have been considered. On the one hand, the multi-dispersed nature of these distributions has allowed the implementation of realistic settings, where the proper differences (in terms of droplet number and distribution) affecting talking and coughing events have been taken into account. On the other hand, the dilute nature of the considered multi-phase flow has represented the necessary pre-requisite for the application of a two-way coupling strategy (thereby alleviating us from the burden to account for particle–particle interactions in the frame of four-way coupled or similar numerical approaches). The evaporating nature of the droplets has been also modeled, taking into account both concentration and temperature gradients. Moreover, the formation of solid nuclei due to the crystallization of salt dissolved in the droplets has also been considered, as this can yet have a non-negligible influence on the evaporation process. The high-fidelity representation of boundaries has not been limited to the environment, but has been applied also to the involved human beings (modeled as static mannequins). A part of our meticulous description of such methodological aspects has also resulted from the realization that existing studies are rare and sparse, which has hindered to a certain extent the development of general criteria and consensus about the most relevant approach for the analysis of this class of problems. In the majority of cases, these can be considered weakly compressible and can be therefore treated in the frame of pressure-based solvers and related extensions to low-Mach-number flows.

To make the outcomes of such a study as much realistic as possible, all the intricacies (in terms of the forced topology of the flow and morphology of the boundaries limiting it) have been considered and implemented. This has required a priori detailed analysis of the HVAC (heating, ventilation, and air conditioning) system for interregional trains, the effective standards used by engineers to design such systems and of the typical droplet emission processes associated with human beings. Moreover, fully three-dimensional (extremely demanding) numerical simulations have been used to adequately account for all these aspects.

This framework has led us to identify the average behavior of particles by revealing their spatiotemporal evolution. We have connected such statistics to the evaporation process, giving deeper insights into the particle transport mechanisms in conjunction with the considered complex flow topology and typical physiological or natural events such as talking and coughing.

It has been shown that with all the relevant effects taken into account, the paths of the incompressible liquid particles are intertwined and connected with the intricacies of the flow established inside the train cabin as a result of the air conditioning system. Nevertheless, the initial conditions (nature of the considered physiological event) have also a remarkable impact on the droplet cloud evolution in terms of symmetry and droplet transport rate, which finally result in a different risk level for the occupants of the cabin.

In a nutshell, the main outcomes of the present numerical study can be summarized as follows.

A single coughing passenger seating in proximity to the compartment centerline can produce droplets potentially impinging on up to 7 distinct human beings. In this specific case, particles are faster than those released in an equivalent talking event. Therefore, as expected, a coughing event leads to a more dangerous scenario.

More realistic circumstances, however, are represented by the situation in which a relatively high number of passengers are talking. In order to obtain statistically meaningful data for these cases, we have considered talking mannequins periodically positioned along with the entire extension of the compartment. The main outcomes of the related simulations are summarized in Table XIV.

Some additional insights stem naturally from Table XV where we have collected specific data extracted from the case with 10 talking passengers (Fig. 37). Under a slightly different perspective, the companion Table XVI indicates the maximum, minimum, and mean number of non-talking passengers exposed to the risk of infection as a function of a “discrete subset” R of talking heads (the variation in the numbers being due to the specific talking passengers considered as source of droplets).

Finally, Fig. 39 shows the relationship between the number of passengers “at risk” (hit by droplets) and the number of droplet emitting heads. The significance of this final figure lies essentially in the evidence it provides about the quasi-linear relationship between the considered quantities and the related ratio (which we loosely define as the “train cabin infection spreading ratio” η). As the reader will realize by inspecting this figure, for the considered conditions η takes a value $≈2$ (obviously, this value being valid only in the considered range of the talking heads number; a decrease in η should be expected due to saturation effects as the number of talking passengers increases).

Another interesting way to think about the present results is to consider that these simulations have revealed that if evaporation is properly taken into account, the size of the droplet clouds attains an asymptotic (maximum) value that does not depend on the time during which the talking events occur. This is because droplets evaporate completely as a certain threshold distance from the source is exceeded. Despite this “limitation” affecting the droplet lifetime, the droplet clouds can interfere with other passengers and lead to an intricate matrix of cause-and-effect relationships, as shown in Tables XV and XVI.

As a concluding remark, we wish to point out that, although it has been applied to a train cabin, this model could be used to track the particle spreading process in all those circumstances where some flow control is implemented (through air conditioning) and the Reynolds numbers take values comparable to that considered here, for example, the cabin of airplanes, offices, and restaurants.

For the sake of completeness (and in order to get the additional insights into these processes), we have finally assessed the relative importance of molecular and turbulent diffusion through comparison of relevant contour maps showing the spatial distribution of the kinematic (molecular) viscosity (satisfying the well-known Sutherland's law for gases) and the turbulent viscosity. This has revealed that turbulent diffusion is dominant (the turbulent viscosity being three orders of magnitude larger than the molecular counterpart, see Fig. 40). A similar concept can also be applied to the other transport properties (thermal diffusivity and species diffusion coefficient) as both the molecular and turbulent Prandtl and Schmidt numbers take (approximately) unit values in the present work.

An exciting prospect for the future is to extend such analysis to non-nominal situations such as those resulting from an anomaly in the conditioning system. Moreover, more sophisticated evaporation models shall be elaborated that can account for the presence of mucus in the saliva droplets (especially with regard to coughing events).

The authors would like to thank M. Eid and M. Gorlani (respectively, General Manager and Technical Director at Blue Engineering) for supporting this research.

Moreover, they would like to thank A. Gallina (Aerothermal Engineer at Blue Engineering) for some useful discussions about the performances of HVAC systems and thermal comfort requirements.

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