key: cord-0987655-y849jdkv authors: Gimenez, Juan M.; Idelsohn, Sergio R.; Oñate, Eugenio; Löhner, Rainald title: A Multiscale Approach for the Numerical Simulation of Turbulent Flows with Droplets date: 2021-06-26 journal: Arch Comput Methods Eng DOI: 10.1007/s11831-021-09614-6 sha: 55f9fe14d606a962afe7af4f50788e4f853feda6 doc_id: 987655 cord_uid: y849jdkv A multiscale approach for the detailed simulation of water droplets dispersed in a turbulent airflow is presented. The multiscale procedure combines a novel representative volume element (RVE) with the Pseudo Direct Numerical Simulation (P-DNS) method. The solution at the coarse-scale relies on a synthetic model, constructed using precomputed offline RVE simulations and an alternating digital tree, to characterize the non-linear dynamic response at the fine-scale. A set of numerical experiments for a wide range of volume fractions, particle distribution sizes, and external shear forces in the RVE are carried out. Quantitative results of the statistically stationary turbulent state are obtained, and the turbulence modulation phenomenon due to the presence of droplets is discussed. The developed synthetic model is then employed to solve global scale simulations of flows with airborne droplets via the P-DNS method. Improved predictions are obtained for flow conditions where turbulence modulation is noticeable. The current crisis due to the SARS-CoV-2 virus, which by May 2021 has infected over 165 million people and killed millions worldwide, has led to a significant increase of research in the area of modeling and simulation of infectious diseases in air. In particular, the simulation of aerosols propagation with available numerical approaches has been widely employed to understand and prevent the propagation of infectious diseases. Several recent works can be found that consider a wide range of situations: the propagation of aerosols in the built environment to model the reopening of closed spaces after the crisis [1, 2] , verifying and planning social safety distances [3, 4] , aerosol transport and deposition in a human respiratory tract [5] , and, in general, to understand and prevent the propagation of pathogens in air considering realistic situations [6, 7] . However, and maybe due to the need to provide immediate results, there have not been substantial developments to validate and improve the standard numerical approaches. In particular, there are some known weaknesses of well known techniques for the turbulence modeling of such kind of multiphase flows. In the present work, a novel approach for the detailed simulation of the turbulent dispersed flow of water droplets in air is presented and discussed. The proposed methodology may be applied to a broad range of multiphase flows, and as such can be used to enhance aerosol modeling for the prediction of pathogen transport and transmission in air. Multiphase flows are defined as the simultaneous flow of materials with two or more phases. Here, the continuous (or carrier) phase, that occupies a continually connected region of space, and may be either gaseous or liquid. The dispersed phase, that occupies disconnected (discrete) regions of space, can consist of either solid (particles), liquid (droplets) or gas (bubbles). In the following, the term 'particles' will be used regardless of whether their state is solid, fluid or gas. The multiscale and nonlinear interactions between the carrier and the dispersed phases lead to complex flow physics and pose unique modeling challenges. Also, many of these flows involve turbulence. The simultaneous presence of two of the most challenging topics in fluid mechanics, namely multiphase flows and turbulence, is still an unsolved problem. The stochastic behavior of the turbulent nature of the carrier flow is altered due to the presence of the dispersed phase, which is known as turbulence modulation. Phenomena of enhancement (production) and attenuation of the turbulent kinetic energy of the carrier phase occur simultaneously because different mechanisms can be active at different length (and time) scales. Attenuation in a particle-laden flow can be explained by a higher inertia, the increased dissipation arising from particle drag, and a larger effective viscosity. On the other hand, production is a consequence of a larger velocity fluctuation due to wake dynamics and self-induced vortex shedding and the buoyancy-induced instabilities due to density variations [8] . The overall modulation depends on the relative strength of the mechanisms, which is a result of the characteristics of the particles, such as volume fraction p , mean diameter d p , mass loading p , also on the carrier phase flow properties, the inertial effects measured by the flow (Re) and the particle ( Re p ) Reynolds numbers, and the response effects quantified by the Stokes S t number. Some models have been proposed to consider the relative influence of these phenomena and predict the turbulence modulation by effect of the dispersed phase [9] [10] [11] . Nevertheless, each of these models has shown only a restricted validity range [12, 13] . Detailed numerical simulations have shown some potential in analyzing the turbulence modulation in specific multiphase configurations [14, 15] . However, the possibilities granted by the currently available computational power are not supported by the experimental data needed to validate the numerical results [16] . Therefore, there still is a lack of a universal indicator for turbulence modulation. This is mainly due to the wide range of length and time scales in the particulate flow systems and the complexity of studying them-either numerically or experimentally. Multiphase flows are intrinsically multiscale in nature. Considering a simplified separation of scales, particle-fluid and particle-particle interactions that occur at the microscale (order of some particle diameters d p ) lead to organized mesoscale structures (tens of d p ) that affect macroscale flow behavior [17] . In the last years several proposals using multiscale numerical strategies have been put forward. A variety of modeling and simulation methods for gas-liquid systems considering all the scales was proposed in [18] . Other approaches focus on the behavior of the thin boundary layer at the discrete phase boundaries, dealing with the mass and heat transfer phenomena that occur there [19] , or the collision of fluids droplets [20] . In the context of gas-solid flows, the particle-resolved direct numerical simulation of the microscale governing equations has been employed for understanding physics and obtaining quantitative information for meso/macroscale developments [21] . Although these developments have shown significant advances in their particular application scope, none has contributed towards the improvement of the turbulence modulation modeling in dispersed flows. In fact, just a few proposals on this topic are found in the literature. These are applied to energy-based turbulence models that account for the effect of turbulence modulation through additional source terms which depend on the Stokes number only [22] . Because the addition of these extra source terms usually leads to a prediction of nonrealistic rates of attenuation and production [23] , they do not generally improve the simulation results. This suggests that there is a necessity to develop a turbulence model suitable for particle-laden flows [24] , and that a systematic study of the results of high fidelity simulations of particle-laden flows can provide further insights [25] . The methodology here proposed to manage multiscale and turbulent multiphase flows is based on the novel Pseudo-Direct Numerical Simulation (P-DNS) method [26] . This multiscale numerical strategy relies on accepting as a premise that the solution on a fine enough mesh is a reliable result. The computation is carried out at two levels or scales, termed coarse scale and fine scale, respectively. The most expensive part of the computations, the solution of the fine-scale, is performed offline on representative volume elements (RVE). The outcome of the RVE simulations provides a database of homogenized stress tensors as a function of several dimensionless numbers which characterize the flow. The P-DNS approach was initially developed for homogeneous flows with massive instabilities [26] . It was later extended to solve general transport equations [27] , and improved to deal with turbulent flows that are out of equilibrium [28] . The results reported show an increased accuracy in comparison to LESbased approaches when employing similar coarse discretizations to solve standard turbulent flow benchmarks. The main aim of this work is to employ the P-DNS approach to define a set of multiphase RVEs whose simulation results allow us to develop a synthetic model that summarizes the behavior of the carrier phase turbulence modulation phenomena due to the presence of the discrete phase. The methodology is applied to construct a database for the specific case of a dispersed flow of water droplets into air, mimicking the presence of saliva droplets on air during aerosol propagation. Making use of high performance computing facilities, a set of numerical experiments for a wide range of volume fractions, particle distribution sizes, and inertial forces introduced as shear loads are carried out. Quantitative results of the statistically stationary turbulent state are obtained. This brings insights into the turbulence modulation phenomenon. The database of RVE results is stored in an alternating digital tree scheme to facilitate its use in general flow solvers via a fast interpolation technique. Global case studies of airborne droplet flows are analyzed, highlighting the differences and similarities found using the fine-scale model developed vis-a-vis a standard turbulence modeling. The P-DNS method is a numerical approach to solve multiscale fluid dynamic phenomena overcoming the computational burden associated with direct numerical resolution of the governing equations.This framework is based on four main concepts (see Fig. 1 for a schematic summary): (a) The unknowns fields are solved separately in two meshes, the coarse mesh and the fine mesh. (b) The fine mesh solution is obtained by direct numerical simulation of many separated subproblems in small domains (the RVEs) under consistent boundary conditions emanating from the coarse mesh fields. (c) The information to be transmitted between the coarse and the fine mesh is written in a non-dimensional way. In this manner, the RVE database is built via offline simulations for all the possible configurations. (d) The database is used as input to construct a synthetic model to emulate the fine scale response. The set of synthetic models for each different multiphase flow is called CFD vademecum. The P-DNS method can be seen as an adaptation of the variational multiscale (VMS) method [29, 30] , where the fine solution is solved numerically instead of analytically. Also, from the homogenization methods' point of view, it can be seen as an evolution of Finite Element square (FE 2 ) methods [31] , where the most expensive part of the computations is performed offline. The P-DNS procedure, which is described in Algorithm 1, has been successfully employed for solving general transport equations including convective, diffusive, and reactive phenomena with time-varying solutions [27] , and also for homogeneous turbulent flows [26] . The P-DNS method relies on the assumption that the RVE represents the physics of the problem solved in a sub-domain without conforming exactly to the same geometry nor the physical parameters of the related problem in the coarse mesh. This can be done because the governing equations of most problems in mechanics can be written in a dimensionless form. For example, in the advection-diffusion case, the dimensionless number that characterizes the problem is the Peclet number. As long as the same Peclet number is preserved in the coarse mesh and in the RVE, and the boundary conditions are properly chosen, the dimensionless results in both scales will be the same [27] . For turbulent flows, it was found that the viscous incompressible flow equations with an added stress tensor have to be solved in the coarse scale, while the entry parameter to the fine scale is a dimensionless tensor, also called instability index (Id). Here, the RVE defines a statistically anisotropic but homogeneous turbulent flow where an equilibrium stress tensor is computed [28] . Due to the complex nature of multiphase flows, the application of the P-DNS methodology to numerically solve the equations governing these flows is a challenging task. A general RVE configuration is established to properly represent any sub-domain of the coarse scale considering the local features of the carrier flow and the presence of the dispersed phase. In this context, the dimensionless parameters (RVE inputs) and the dimensionless results (RVE outputs) must be defined. The particular physics of water droplets in air is selected as a case study. Thereafter, a large number of RVE simulations are performed, covering a wide range of values for the inputs. In a multiphase flow, two critical parameters define the level of interaction between the phases: the volume fraction occupied by the dispersed phase ( p ) and the mass loading ( p ), defined as the ratio of mass of the dispersed to the carrier phase. For small p and p the turbulent carrier flow dominates the dynamics of the discrete phase, a regime which is known as one-way coupling. The (back) influence of the dispersed phase on the carrier phase dynamics can not be neglected when p ≈ 1 , leading to the so-called two-way coupling regime. When p > 0.001 the interactions between particles become important and collisions, agglomeration, and break-up phenomena must be considered, leading to the so-called four-way coupling regime [32] . The discrete phase is described using a Lagrangian reference frame for each particle, while the carrier-phase is considered as an incompressible flow with an Eulerian technique. The study of the dispersed phase behavior makes use of a standard Discrete Element Method (DEM) [33] . The dynamics of a particle p with mass m p and inertia moment I p are devoted to the Newton's second law leading to: where p and p are the linear and angular velocities of particle p, respectively, p is the particle-fluid interaction force acting on particle p, while g p is the gravitational force. Also, c pj and pj are, respectively, the contact force and contact torque acting on particle p by its n c p contacts (which could be either particles or walls), while nc pk is the non-contact force acting on particle p from particle k or from other sources (for instance electromagnetic or biological forces in some active particulate flows), with a total of n nc p non-contact sources. For the sake of completeness, contact and noncontact forces between particles were included in Eqs. (2.1a) and (2.1b), but they are not considered in the range of experiments carried out in this work. The multiphase systems studied here are assumed to be diluted and the collision between particles is not taken into account. Under these assumptions, the model is known as the Lagrangian pointparticle approach which perhaps has the longest history in turbulent dispersed multiphase flow computations [34] . The gravitational and buoyancy forces are computed as one total force acting on a spherical particle p as: where p and d p are, respectively, the density and diameter of particle p while is the gravity acceleration vector. The contact between the particles and the surrounding fluid is modeled by the interaction force p which, in its general form, can be written as: where d is the drag force, ∇p the pressure gradient force, ∇⋅ the viscous force, vm the added mass force, and are other forces such as the Basset force, the Saffman force for the lift, the Magnus force, etc. In gas-liquid systems with large density differences ( f ≪ p ) drag force dominates and the remaining forces can be neglected [35] . The drag force depends on the relative velocity − p , where is the fluid velocity, and is calculated as: where Re p is t he par ticle Reynolds number Re p = ( f | − p |d p )∕ f , being f is the fluid dynamic viscosity. Finally, the drag coefficient C d is an empirically measured function of Re p , given by [36] In the case of gas-liquid systems, the assumption of spherical particles used by Eq. (2.5) holds for Weber numbers less than one. The governing equations for an incompressible fluid phase in the presence of a secondary particulate phase are the volume-averaged continuity and momentum equations, i.e. where = f f + (1 − f ) p with f the volume fraction occupied by the fluid phase, p the pressure, and the tensor is the viscous stress of the fluid phase. The momentum transfer term , which enforces the two-way coupling between the particles and the fluid-phase, is defined in each computational cell as where N p is the number of particles located in a cell of volume V c . The fluid-phase volume fraction is computed for a given cell as The governing equations presented above cover a wider range of applications of concentration and mass loading. Note that for diluted concentrations one can approximate f ≃ 1 , simplifying several expressions. In the context of the P-DNS framework, the solution at the fine scale is obtained by subdividing the global domain into small representative volume elements (RVEs) that are simulated separately by imposing consistent boundary conditions in accordance with adequated dimensionless parameters. By taking previous developments as a basis [26] , the RVE is defined as a cubic domain of volume R and boundary R centered at the origin of coordinates. Without any loss of generality, the velocity and pressure fields in the RVE can be expressed as where f and p f are the field fluctuations. The upper-index f denotes the fine scale field following the P-DNS nomenclature. The mean velocity and pressure gradients, and , are defined as The mean velocity and pressure fields, ̃ and p , are required to satisfy where ̃ = 0 and p = 0 are imposed as the RVE reference frame is set to move with the coarse scale velocity [26] . Equations (2.11a), (2.11b), (2.12a), and (2.12b) imply that null spatial average is required for the fields fluctuations and their gradients, it is: In the RVE cubic domain, conditions Eqs. (2.11a) and (2.12a) are satisfied shifting the solution by its instantaneous spatial average. Additionally, by imposing periodic boundary conditions for f and p f , Eqs. (2.11b) and (2.12b) are satisfied. This results in the following periodic jump conditions for and p, Similarly to the homogeneous flow case, the target of RVE simulations is to obtain statistically steady fluid inertial stresses, which are stored in a database. The spatial mean of the inertial stresses is defined as: For a statistically steady quantity Q(t), its equilibrium value is defined as assuming that statistical equilibrium is reached at time t i . In equilibrium state the (dimensionless) inertial stresses of a homogeneous fluid depend solely on the instability indices Id , defined as where S is the symmetric part of , as the instabilities are not influenced by the solid body rotation part. Moreover, with the appropriate tensor rotations [26] , tensor Id can be expressed in a coordinate system where the axial strain velocities (diagonal terms) vanish. Then, this tensor can be written as (2.14) which depends on only two parameters, Id 1 and Id 2 , with Id 1 ≥Id 2 . The strategy to construct the database is as follows: first the values of H, f and f are set equal to 1. Then, the range of the desired Id 1 and Id 2 values is swept by varying the velocity gradient. The presence of particles must also be managed through dimensionless parameters. As the RVE is thought to represent the physics of any sub-domain of the coarse scale, a wide range of droplet concentrations, sizes, and fluid/particles density ratios needs to be analyzed in advance. In order to take into account the polydispersion of particle sizes, the Rosin-Rammler distribution is employed to seed an initial cloud of particles randomly located in the RVE. The mass fraction Y of droplets of diameter lower than d p is given by: where d and n are the representative dimensionless diameter and the spread parameter respectively. The non-dimensionalization of the particle diameters is performed considering the coarse scale subdomain reference size H c , i.e. d =d c ∕H c with the upper-index c denoting for coarse scale data. For large values of n, the distribution tends to a uniform particle size. On the other hand, the value of n is between 1.5 and 8 for the initial droplet size distribution of most aerosols [4, 37] . Also, the Rosin-Rammler distribution is easy to be reconstructed from a set of particles by processing their diameters. This feature enables us to supply the coarse mesh simulations with the database results in the P-DNS multiscale approach. As the same droplet distribution can be satisfied with a different number of particles, it is also required to impose the initial dispersed phase volume fraction p . Note that the concentration is already a dimensionless parameter. As the objective is to compute equilibrium statistics for a particular configuration, the initial value of p has to be preserved Finally, Fig. 2 presents a conceptual scheme of the configuration of a typical RVE. The solution of the discrete phase [Eqs. (2.1a), (2.1b)] is performed using a semi-implicit integration approach [38] . This guarantees numerical stability using moderate timesteps when managing the wide range of timescales introduced by the polydispersion of particle sizes. As the computing requirements are proportional to the number of particles, the strategy of computational parcels is employed. A parcel p agglomerates a set of N̂p particles. The main hypothesis is that nearby particles experience the same changes in their properties. The extensive quantities p are computed as: Thus, the governing equations are computed for an individual particle p and the effect of that particle is transferred back to the fluid N̂p times. Note that there is a restriction on the largest particle diameter when the particle-point approach is employed. The largest particle diameter should be smaller than the mesh size h. In LES simulations, d p < 0.2 , where is the smallest resolved eddy size, while for DNS simulations d p < 0.2 , where is the Kolmogorov scale. The equations for the continuum phase [Eqs (2.6a), (2.6b)] are solved with the standard finite volume method. The timestep is restricted such as the Courant number does not exceed the unity, while the pressure-velocity coupling is solved with the PIMPLE algorithm [39] imposing a convergence of three orders for the velocity and pressure residuals. Second-order operators are used for both the spatial and time discretizations. As the particles motion is integrated independently from the flow solver, it is not difficult to envision situations where for the extreme cases of very light or very heavy particles physically meaningless or unstable results may be obtained. In order to avoid these undesired numerical artifacts, the limiting procedures taken from [40] are employed during particle updates. Multiphase flows are extremely complex and the spectrum of possible configurations is so large that a full analysis of the turbulence modulation phenomena over the entire range of cases is unaffordable at present. In this work, and aligned to the current urgencies for fast and accurate numerical predictions to help alleviate the pandemic crisis, a synthetic model from pre-computed (offline) simulations of representative volume elements is built to characterize the nonlinear dynamic response of the multiphase system composed by water droplets in air. This synthetic model will represent a specific entry in a general CFD vademecum for multiphase flows. In the sequel, air is considered for the carrier phase, and water droplets for the discrete phase (in order to mimic human saliva) with a density ratio of p ∕ f = 1000 . The twoway coupling approach is selected according to the expected mass loading during the flow of droplets on air. The concentration range has been restricted to 10 −6 ⪅ p ⪅ 10 −3 . As the saliva could originate from coughing, sneezing, respiration or talking, and also each droplet could evaporate, a wide range of particle sizes has to be considered for every potential sub-domain of the coarse scale. Therefore, the coarse scale particle diameters d c p could vary between approximately 1 μ m and 500 μ m [1] , while typical element sizes H c on meshes employed for the coarse scale simulations are between 1 and 50 mm, leading to dimensionless particle diameters d c p ∕H c ranging from 10 −5 to 10 −2 . As the total number of parcels for a given RVE simulation is restricted due to computational requirements, the initial seeding procedure clusters particles with similar diameters into parcels with an average diameter that satisfies a prescribed minimum volume per parcel. This volume is usually set to the ratio between the total particles volume and the maximum number of parcels. In the set of simulations carried out up to 50,000 parcels were considered. On the other hand, the energy inflow to the system via the shear forces applied on boundaries (which are modeled by prescribing the Id 1 and Id 2 parameters), is taken from [28] . Table 1 presents the range selected for the input parameters of the multiphase RVE. The combination of these, restricted to Id 1 ≥Id 2 , leads to a set of 810 RVE configurations, which were simulated using an RVE discretization of 41 3 equal-size cells. The initial state of an RVE with a particular configuration is shown in Fig. 3 . As stated above, the timestep is restricted to a maximum Courant number of one for the carrier phase. The total simulation time is set according to an a priori estimation of the time required to achieve converged statistics. As an example, Fig. 4 presents the temporal evolution of the instantaneous and the time-averaged values of some components of the inertial stress tensor for a particular configuration. As shown in Fig. 4a , the turbulent behavior of the flow in the RVE leads to a chaotic evolution of the velocity fields which makes it impossible to predict the instantaneous spatial averages of the stresses ( ). However, due to the equilibrium conditions imposed to the RVE, i.e. production and dissipation rates of the kinetic energy are balanced, the time-averaged statistics eventually converge to an equilibrium value ( ̄ ) (Fig. 4b) . Using Eq. (2.17), this tensor is non-dimensionalized, which leads to ̃ . This latter tensor is stored in the database. Each simulation requires between 25 and 75 core-hours. To evaluate the entire set of simulations in a reasonable time slot, the total computational load of over 40K core hours was distributed into several computing nodes of an available high-performance computing facility [41] . The next subsection deals with a summary of the results obtained in this study. For completeness, the full database, i.e. the entire set of stress tensors computed for each RVE configuration, can be found in [42] . Some definitions are introduced to facilitate the post-processing and discussion of the results. For a given time t the kinetic energy of the carrier phase is computed as which is a statistically steady quantity whose equilibrium value k can be computed using Eq. (2.16). Henceforth, this scalar quantity is used to characterize the turbulent flow in the RVE as it summarizes in a unique and positive value some features of the tensor ̄ . Although the computed dataset is large and it can be examined from several points of view, here the presentation of the RVE results is focused on analyzing the effect of the particles on the kinetic energy of the carrier flow. Figure 5 presents the time-averaged kinetic energy obtained by 126 RVE simulations. Results are clustered in three subfigures according to the particle mean diameter: large particles ( d = 0.005 ), medium particles ( d = 0.0005 ) and small particles ( d = 0.00005 ). The cases with Id 2 = 0 and varying Id 1 , p and n are analyzed, while reference results previously obtained solving unladen RVEs under same shear stresses, i.e. homogeneous flows with the same Ids, are also presented for comparison purposes [28] . Turbulence modulation can be seen as the change of intensity of the time-averaged kinetic energy of the carrier flow due to the presence of the particles regarding the particle-free flow. The unladen cases show, as expected, that an increase in the shear stresses leads to an increase of the kinetic energy. When the presence of the discrete phase is included in the analysis, Fig. 5a shows that the presence of large particles increases the carrier flow turbulence if the particle concentration is above some threshold. The value of this threshold increases with the shear stress (i.e. the value of Id). For example, large particles in a concentration of p ≈ 10 −4 noticeably enhance the carrier fluid turbulence when Id 1 = 100 . However, a value of p ⪆ 5 × 10 −4 is required to modify the unladen situation when Id 1 = 500 . From the analysis of individual simulations, the turbulence augmentation is obtained as a consequence of higher magnitudes of the averaged velocity field of the carrier phase. They are induced by the higher inertia of the particles and their constructive interaction with chaotic carrier flow patterns. On the other hand, small particles attenuate the turbulence for almost all the configurations analyzed, with the exception of low concentrations with strong imposed shear stress (Fig. 5c) . For medium size particles (Fig. 5b) there is again a critical minimum concentration over which turbulence modulation occurs. In this case, and similarly to small particles, high concentrations tend to attenuate the turbulent intensity of the carrier flow. No apparent generalized dependence of the variation of kinetic energy on the spread parameter n (filled and dashed lines compare same RVE configuration but different n) is discernible. However, although the particle mean diameter d p apparently is the most sensitive parameter to determine the turbulence modulation behavior, in some particular cases the dispersion could modify the results. For example, in the cases of kinetic energy enhancement considering large particles, the simultaneous presence of smaller particles, i.e. when n = 1.5 is chosen, reduces the rate of turbulence increase. Also, for configurations with medium size particles near the critical concentration threshold, using small values of n (large dispersion) produces a kinetic energy increase, which may be associated with the inclusion of large particles. This phenomenon is negligible when the concentration rises and the global results lead to turbulence attenuation due to the massive presence of small particles. As expected, as the concentration parameter p diminishes the results obtained in laden RVEs tend to the unladen cases for the same Id. However, even a very diluted presence of small particles is enough to completely modify the flow behavior: for example, with Id 1 = 100 and p = 1.5 × 10 −6 , a strong attenuation is observed, and a lower volume fraction is needed to neglect the effects of those small particles. For completeness, it is necessary to analyze solutions where three-dimensional shear stresses are imposed. Figure 6 presents the results when Id 1 is fixed to 200 and Id 2 varies. As previously presented in an RVE-based analysis of homogeneous flows [28] , the imposition of shear stress over the third dimension does not ensure an increase of the turbulent intensity. However, the presence of particles modifies this behavior. For large particles (Fig. 6a) , turbulence augmentation is generally found above a certain critical concentration ( p ⪆ 5 × 10 −4 ), with exception of cases with large values of Id 2 where some turbulence attenuation is obtained. In the case of small particles (Fig. 6c) the turbulence attenuation is greater when Id 2 is lower. Moreover, when Id 1 = Id 2 a slight turbulence increase is observed even for medium and high particle concentrations. The cases of medium size particles present strong attenuation for small values of Id 2 and concentrations above ( p ⪆ 5 × 10 −4 ). This effect reduces when Id 2 increases. Also, an undefined behavior for p ≈ 10 −4 is observed. There a dispersed and a non-dispersed diameters distribution induces turbulence enhancement and attenuation respectively. For the three particle sizes, there is a critical concentration value that separates negligible from perceptible modulation, but this critical value of p depends on the particle size. As shown here, for a given RVE configuration the flow behavior can be unpredictable. Several mechanisms coexist in the RVE, whose dependence on the dimensionless input parameters can not be determined in advance. Spectra analysis Selecting Id 1 =100 and Id 2 =0 as a reference case, in Fig. 7 the dimensionless energy spectra is analyzed for the case without particles (Fig. 7a) , and compared with the spectra of other configurations: with turbulence augmentation (Fig. 7b) , with turbulence attenuation (Fig. 7d) , and with negligible modulation (Fig. 7c) . There, the energy spectra corresponding to a large number of timesteps (gray lines) are averaged to obtain a representative energy spectrum for the specific RVE configuration (bold and colored lines). The unladen case is employed here to compute physical and numerical scales to evaluate the adequacy of the mesh discretization chosen to compute a direct numerical solution. The parameters governing the small-scale motion, which is known as Kolmogorov scale, include the dissipation rate per unit mass and the kinematic viscosity . With these parameters, one can define length, time, and velocity scales as follows: under the consideration for the RVE simulation, = 1 and the dissipation rate can be approximated as = u 3 ∕l , with u and l the reference velocity and length respectively. In view of the results, a reasonable choice for a basic analysis is selecting u = u rms , i.e. the root mean square of the velocity field, and l = L , i.e. the integral length scale. For example, for the analyzed case in Fig. 7a , it is easy to see that the integral length scale is L ≈ 1∕(2 ) , and computing u rms = √ 2k ≈200, = u 3 ∕L ≈ 2 200 3 can be estimated, leading to ≈ 0.012 . This latter is slightly smaller than the grid size chosen to discretize the RVE, x = 1∕41 . This fact leads to a proper representation of every turbulent scale. One can observed in the spectra that the largest wave numbers contribute insignificantly to the energy budget. This value of Kolmogorov length scale satisfies the assumption of particle-point approach for any choice of particle diameters of the configurations studied. The turbulence modulation phenomenon due to the presence of particles in the RVE can be also identified by comparing the spectra, as highlighted in Fig. 7e . In every case, the particles increase the energy of the scales in the order of the Kolmogorov scale. Even in the cases where a neutral ). Filled lines: n = 50 , dashed lines: n = 1.5 , dotted lines: particle-free flow effect is perceived evaluating the kinetic energy budget, the energy of large wavenumbers is increased by one or two orders versus unladen RVE. This reflects the interaction of particles with the flow at these scales, but without necessary influencing a change of the behavior at larger scales. In the case of turbulence augmentation, the increase is given for almost all scales due to the enhanced inertia of the whole fluid. Conversely, when global attenuation is computed, the mechanism seems to affect mainly the energy measured at medium and large scales. Standard variables for turbulence modulation The previous analysis of the simulation results discussed the turbulence modulation phenomena from the point of view of the type of interaction between the particles and turbulence and its relationship with the dimensionless parameters employed to construct the RVE database. However, in order to compare with experimental or other numerical references, it is mandatory to present the data using the standard variables used in the bibliography. The most widely known approach is perhaps the classification map of regimes of interaction between particles and turbulence proposed in [32] . There, the relationship between the particles concentration and a time ratio parameter is displayed considering the regions of turbulence augmentation and attenuation. In Fig. 8a every RVE simulation is plotted on the map, where a red circle denotes turbulence augmentation (turbulent intensity greater than 110% of the unladen case), while a blue square indicates turbulence attenuation (turbulent intensity lower than 80% of the unladen case). The remaining cases are considered neutral regarding turbulence modulation and marked with a black cross. The time ratio parameter relates the particle response time p = d 2 p ∕(18 ) and the turnover time of large eddies defined as e = L∕u rms . Regarding the experimental and numerical collections summarized in previous publications [13, 32] , the results of the simulations here performed agree with the general observation. There, the timescale ratio is defined as the ratio between the response time of particles ( p ) and the turnover time of large eddies ( e ). Flows with larger timescale ratios induce turbulence enhancement, while smaller timescale ratios lead to turbulence attenuation. References usually mention that the separation between them is defined by the timescale ratio p ∕ e = 1 for any mass loading, represented by the dashed green line in Fig. 8a . The results obtained here do not show such a fixed threshold. Instead, a dependence on the inter-particle distance is found. This value, which measures the average distance between particles, was already mentioned as a critical parameter for turbulence modulation [43] . It is easy to see that, to keep the same inter-particle distance, the concentration should be modified proportionally to the cube of the diameter. This could explain why the enhancement-attenuation threshold found on RVE simulations is not fixed to p ∕ e = 1 . Further studies are required to gain insight into these interactions and phenomena. Another reference diagram employed in the literature [12, 13, 44] displays the change in turbulence intensity as a function of the length scale ratio between the particle diameter and the integral length scale, as presented in Fig. 8b . The results obtained, where the particle diameter reported is the mean of the diameter's distribution, quantitatively agree with the variation in turbulent intensity. However, the transition between turbulence augmentation and attenuation is not as sharp as reported in the literature. The narrow range of flow conditions explored in the reported experiments could justify the mentioned differences. The results of equilibrium stress tensors stored in the database obtained by the RVE simulations can be considered as sample points of a continuous N-dimensional function F, where N is the number of inputs for the RVE, i.e. Therefore, it is necessary to define the multi-valuated and multi-dimensional function F to predict the fine scale response for non-simulated points. For this objective, the following actions are carried out to prepare the sample points: (a) The results obtained are extended using previous computed data corresponding to unladen RVEs [42] . This extension is necessary to take into account cases where the region analyzed in the coarse scale is unladen. In this regard, it was established that unladen results correspond to a disperse phase volume fraction of p = 10 −10 , i.e. a negligible concentration. (b) The mean diameter and the concentration inputs are re-scaled using a logarithmic function to favor a homogeneous sampling over these dimensions. (c) Every input parameter is normalized between the values 0 and 1 in order to have the same range of values for each of the inputs. This reduces the potential numerical issues or instabilities throughout the training and/or the searching processes. For the task of describing the F function, artificial neural networks were successfully employed in previous works [27] . However, the reliability of their predictions depends on the non-deterministic selection of the network structure and transference functions. Therefore, we turn to a straightforward strategy for data interpolation without decision parameters. The algorithm and the associated data structure employed here is the alternating digital tree (ADT) [45] . In this approach, an ADT is constructed by inserting the k N-dimensional sample points. Each sample point represents a simulated RVE configuration. Then, the interpolation of a non-simulated point can be reduced to the problem of determining the members of the set of k simulated points which Simulations with no significant changes of turbulence regarding particle-free case are marked with × lie within a prescribed subregion of a N-dimensional space, which is known as geometric searching. The ADT structure allows for the efficient solution of the geometric searching problem, leading to an algorithm whose computational burden is proportional to k log(k). Algorithm 2 briefly explains how to employ the synthetic fine scale model developed in the P-DNS framework. The iterative coupling between the coarse and the fine scales is established indicating which equations must be solved at each stage. Using this model, a transport equation for each component of the fine-scale stress tensor must be included in the coarsescale problem. In a time-step of the P-DNS, the computing time required to calculate the dimensionless parameters, evaluate the database, and apply the memory model is about 25% of the total cost. This overhead, where 75% corresponds to solving the six transport equations, is comparable to the extra cost when standard turbulence models are used. The turbulent flow in a channel confined between two parallel plates (Poiseuille flow) is studied in this subsection for the unladen and laden cases. Previous works [46, 47] have found that, under equivalent body forces applied, the mean The instability indexes, Id 1 and Id 2 , which depend on the local velocity gradient , the kinematic viscosity and cell reference size H c , are computed following the procedure explained in [26] . The local strain tensor is rotated to obtain a pure shear configuration as shown in Eq. (2.19) . The discrete phase volume fraction p is computed using Eq. (2.9), while the density ratio is fixed (as only air/water systems are considered here). Finally, to obtain the particle size distribution parameters, processing of the particles in the cell must be performed. First, the local complementary cumulative distribution function for the mass fraction ( Y c ) is computed. Second, the mean diameter d * is the value that satisfies Y c (d * ) = 1 − e −1 . This value must be non-dimensionalized with the cell size, i.e. d = d * ∕H c . Third, the spread parameter is obtained searching for the value of n that minimizes the error between the computed distribution Y c and the Rosin-Rammler distribution Y [Eq. (2.20) ]. To predict the actual fine-scale contribution at the current time, the memory model presented in [28] is employed. streamwise air velocity increases when inertial particles are introduced in the flow. In this section the target is first, to obtain reference solutions for the flow with and without particles using a direct numerical simulation (DNS) procedure. Then, the P-DNS methodology is employed to predict the flow behavior using a coarse mesh. Finally, these results are compared against those obtained using other numerical approaches for turbulence modeling. The case study consists of flow at Re = 57,000 in a domain mimicking a ventilation air duct of 1.5 × 0.5 × 0.5 m. Periodic boundary conditions are applied in longitudinal and transverse coordinate directions, while a pressure jump of p = 0.0151875 Pa is imposed between the inlet and the outlet driven the flow. Gravity is neglected to highlight the interaction between turbulence and particles without the additional effect of particle deposition. The initial solution for the air velocity field is obtained by performing first a simulation without particles under the same pressure conditions until a statistically stationary state of fully developed turbulence is reached. Then, solid particles of = 1000 kg/ m 3 and 200 μ m of diameter are inserted at random positions, homogeneously distributed over the entire computational domain and with initial velocity equal to the air velocity at the position of each particle. The volume fraction of the dispersed phase is = 0.0005 leading to a mass loading of = 0.5. For the DNS simulation, the domain is discretized with a mesh of 300 × 150 × 50 cells refined towards the walls. A DNS solution is obtained solving the equations for the carrier flow and the particles employing second-order time and spatial schemes, and without including any additional modeling for turbulence modeling or boundary layer treatment. On the other hand, a coarse mesh of 36 × 12 × 12 equal-size cells is employed to obtain the P-DNS solution. This uses the RVE synthetic model developed here, and the wall-RVE model [28] to improve the prediction of inertia stresses on wall-bounded cells. Also, solutions with a static Smagorinsky LES model and an unsteady RANS solver with a standard k-model are obtained in same coarse mesh for further comparisons. In the latter, standard wall functions are considered but without including any extra term for turbulence modulation. Statistically converged solutions, i.e. averages of the results computed by averaging over the two homogeneous directions and over a long enough simulation time, are presented in Fig. 9 . The mean streamwise velocity (on the left) and the mean kinetic energy (on the right) of the carrier phase along the wall-normal coordinate are shown. Due to the inertia introduced to the system by the particles, laden DNS solutions show the expected increase in kinetic energy. This inertia enhancement leads to a modification of the mean profile, which is sharpened, and the flow rate is increased. These mechanisms, regarding the energy value and its consequence on the velocity, are very well reproduced by the P-DNS solution, even on a very coarse mesh. The noticeable lack of accuracy of the other numerical approaches employed should be also noted. In the LES case, although some inertia augmentation is predicted, as the boundary layer is fairly unresolved the overall results are inaccurate. On the other hand, solving with the URANS model negligible differences are found between the laden and unladen results, which shows that ensemble averaged solutions fail to predict turbulence modulation. It is widely accepted that aerosol or respiratory droplet transmission is a critical factor for the rapid spread and continued circulation of SARS-CoV-2 virus in humans. In this context, computational techniques can elucidate the flow and propagation of viruses or other contaminants in order to mitigate or avoid infections. The target of this section is to perform a set of numerical experiments to evaluate the influence of the turbulence modeling when a saliva exhalation event is simulated. In particular, we are interested in evaluating the secure social distance of 1-2 m. This value, which was estimated from predictions of the distance traveled by expelled saliva droplets using standard numerical approaches among others [48] , could be over/underestimated due to their lack of modeling of the turbulence modulation phenomenon. The case studies presented in this section model respiratory events such as talking or coughing into a closed environment where the ambient temperature is assumed to be 20 • C. In order to simulate the saliva expulsion, a circular region of radius (r = 2.5 cm) roughly models a human's being the velocity and temperature defined as u(t) = 5f (t) m/s and T(t) = 20 + f (t)(37 − 20) • C respectively [1] . The saliva (water) droplets are injected into the domain from the same circular region with an initial direction that conforms a spray cone angle of 30 • . Particles are released during 0.2 s with an initial velocity equivalent to the surrounding flow and a temperature restricted to T p = 37 • C. Temperature differences between air and droplets expelled by the human body and the ambient temperature can not be neglected. Therefore the energy equation is added to the governing equations while buoyancy effects are modeled coupling the energy and momentum equations through the Boussinesq approximation [49] . Moreover, heat transfer from the surface of Lagrangian droplets to the surrounding fluid is considered using the Ranz-Marshall empirical correlation [1, 50] . Two conditions for the aerosols are modeled: (1) Exhalation with a large amount of water from the respiratory tract, and (2) the same droplets but after drying out under typical ambient conditions. The Rosin-Rammler distribution is employed to impose the number of aerosols as a function of the diameter. In (1) parameters d = 80 μ m and n = 8 are taken from [3] where a data fitting of the experimental measurements reported in [51] was performed. In that experimental study, exhaled droplet's mass and size due to talking and coughing were quantified. Parameters in (2) , d = 10μ m and n = 8 , were selected following the study in [52] where the transformation of particles' size distribution due to evaporation was predicted. With the aim of reducing the computational requirements and in agreement with the injection direction of the particles, the domain has the conic shape shown in Fig. 10a . At the outlet fixed static pressure is imposed, while on the remaining boundaries slip conditions are chosen. The domain is discretized using a mesh of 952K hexahedral cells, which is refined towards the inlet. Each flow condition is solved using both the P-DNS procedure and a static Smagorinsky LES approach. Because the main difference between these two methods is the treatment of the unresolved scales of turbulence, the comparison is useful to evaluate the influence of including the modeling of the turbulence modulation phenomenon. Note that for the cell length employed, it is assumed that the isotropic turbulence requirement of LES methods for non-simulated turbulent scales is satisfied. Table 10b summarizes the case studies analyzed. The quantitative comparison of the simulated cases is done by computing the time evolution of the accumulated saliva volume that reaches a certain distance from the inlet. For this purpose, six sample planes i at distances x k = 0.1k from the inlet in the streamwise direction are included. For a given plane k, the accumulated volume V k at time t is computed as: Figure 11a presents the predictions of P-DNS and LES for the penetration of droplets with the condition (1). Inertia allows more than 90% of particles to advance behind 1 , but most of them are quickly stopped by the air, and then sink slowly towards the floor in close proximity to the inlet. This is reflected by the results where less than 35% of saliva volume expelled can reach 0.3 m of distance, and less than 5% can reach 0.6 m. Predictions by the P-DNS and LES techniques are quite similar, and their agreement on computing the accumulated volume behind planes is a consequence of a very similar prediction of flow pattern evolution, which can be observed in Fig. 12 . There is a perceptible discrepancy between numerical predictions on the accumulated volume at plane 4 . The reason can be found observing the particles snapshot for t = 0.4 where a dense cluster of droplets has advanced more in the P-DNS solution than in the LES solution. Figure 13 shows a snapshot of the Id 1 field on the slice Z = 0 at time t = 0.4 s. This field actually is one of the parameters used by each cell to evaluate the local response of the fine scale. From the RVE database it can be seen that an absolute value of Id approximately below 50 does not produce instabilities at the fine scale. The figure indicates that most of the droplets are located outside the main flow stream where turbulence is more developed (higher Id numbers). Based on the RVE results obtained, as the concentration is low in these regions, the turbulence modulation phenomenon is not expected to be appreciable there. The comparison of the results obtained modeling this phenomenon (P-DNS) or not (LES) confirms this expectation. Figure 11b presents a comparison between the predictions of the accumulated volume using P-DNS and LES when condition (2) is simulated. Due to the smaller size of droplets, and therefore smaller Stokes number and settling velocities, most of the aerosols are totally dragged by the main flow. The inlet flow stops at t = 0.2 s, then the air velocity is slowly dissipated but its inertia carries a large number of droplets. For instance, at t = 1 s approximately 30% of the droplets have reached the last plane 6 . Penetration lengths are similar for both the P-DNS and LES predictions. However, more differences are seen than in condition (1) as the particle trajectories are more intricate since they follow the complex structures of the flow. As the aerosols follow the instantaneous airflow streamlines, the patterns described by their trajectories are more complex, as shown in Fig. 14 . A distinctive feature is the mushroom-like shape that the jet describes when a laden flow impacts a fluid at rest. As mentioned above, under this condition (2) more droplets are located on the main stream where larger gradients are found. This latter is a potential condition to find turbulence modulation. Figure 15 displays the volume fraction of the dispersed phase p on the slice Z = 0 at time t = 0.4 s. There, it can be seen that the volume fraction of the droplets is not larger than 10 ppm. As the dimensionless droplet size, i.e. the ratio between the particle diameter and cell size is approximately 10 − 3 ≤ d p ∕h ≤ 3 × 10 − 3 , the evaluation of the RVE database for these laden flow conditions returns results showing no turbulence modulation. In particular, we refer to Fig. 5a where laden RVE results with these flow conditions do not differ from the unladen case. The similarity between P-DNS and LES results reveals that the turbulence modulation could be neglected when the propagation of pathogens in air in the built environment is analyzed. The main reason is the low value of mass loading resulting after a saliva exhalation event in an open space. These results also validate the predictions, using standard turbulence models, reported in previous studies [1, 3, 7] . The complex nature of particulate systems leads to a failure of finding a universal indicator for turbulence modulation mechanisms, and therefore standard approaches for turbulence modeling could fail to predict the behavior of multiphase flows. A general computational methodology has been proposed to circumvent this shortcoming. It is based on the offline simulation of representative volume elements to quantify the change in the dynamic response of the carrier flow due to the presence of the particles under different flow configurations. The high fidelity solutions of pure shear RVEs were obtained by DNS simulations using the Lagrangian pointparticle approach. These results were validated by performing an analysis of their kinetic energy spectrum and by comparisons with published experimental and numerical results. An overview of the many findings mentioned claims that turbulence augmentation could be found when (2), solution with P-DNS. Concentration field p on slice Z = 0 at time t = 0.4 s large particles are employed, while turbulence attenuation could be detected when small particles are included. The conditional statements are due to the influence of other flow parameters involved, such as the dispersed phase volume fraction, the inter-particle distance, and the shear stress imposed. These parameters have shown to modify the turbulence modulation development and should be studied in more detail. In any case, the quantified data obtained, which is synthesized to favour its use on future multiscale computations, are, up to the authors' knowledge, the study with the widest range of water-air configurations analyzed. However, it is worth mentioning that the results obtained are restricted to physical phases with a density ratio of 1000 approximately, such as air-water systems. Also, the interaction level corresponds to a diluted solution (two-way coupling), where the dispersed phase is conformed by spherical and rigid particles. Most of these assumptions are directly related to having employed the particle-point approach to model the particulate phase. Several physical phenomena were neglected in the current RVE studies. For example, gravity, mass and energy transfers were not considered. The latter allows us to consider equilibrium conditions in the RVE and employing an existing memory model to compute the current inertia stress. This approach should be extended, for example, if evaporation processes were studied. Another open question is the reliability of using a large number of particles per computational parcel, as occurs in RVE configurations with a large concentration of small particles. Although this modeling makes the solution of these particular cases affordable, no numerical study is found in the literature about the sensitivity of the results to this assumption. A quite high number of RVE simulations have been solved (over 800) using high performance computing facilities. This offline procedure, which is performed only once for the air-water systems studied, took 40K core-hours to be completed (3 weeks on twenty computing nodes of four cores each). Even so, it is very likely that some dimensionless parameters could be sub-sampled, leading to an inaccurate prediction of the inertia stresses on configurations far from the sampling points. Nevertheless, the developed database was shown to be adequate for the purposes of the current work, i.e. to establish the methodology and evaluate the ability to reproduce turbulence modulation phenomena. The synthetic model developed could be then enhanced by performing new RVE simulations of unexplored configurations, improving in this way the accuracy of interpolations. The general aim of the P-DNS methodology is overcoming the computational burden associated with the direct numerical solutions of the governing equations in complex flow problems. Therefore, when they are affordable, DNS solutions could be employed as reference results to validate the P-DNS predictions on coarser meshes. In this online stage, the computing time required to evaluate the database and apply the memory model represents approximately 25% of the time-step cost. This extra computational burden of the P-DNS is comparable to the overhead when standard turbulence models are used. Evidence of the relevance of this methodology is found when the turbulent flow between two plates at p = 0.5 is studied. There, the modification of velocity and turbulent statistics of a base flow due to the inclusion of droplets were well predicted, on coarse meshes, only by the P-DNS method. However this influence on the results should not be generalized to the full range of numerical experiments. Whether the volume fraction or the mass loading of the particulate phase are below a certain threshold, turbulent statistics do not vary between laden or unladen cases. In such situations, a standard turbulence model can predict the main flow features. Aerosols spread from saliva exhalation to open spaces is a problem where turbulence modulation is negligible. This study has contributed to validate that using standard turbulence models can be accurate enough in these cases. Among the wide range of possibilities for future applications of the P-DNS approach, a next step in the study of multiphase flows is to avoid the treatment of individual particles at the coarse scale solution. This implies using a two-fluid like approach at that scale, which requires that dispersed phase velocity and average forces be also computed in the RVE. Homogenization procedures must also be performed to formulate the coarse scale equation for the dispersed phase. In this work, an extension of the multiscale Pseudo-Direct Numerical Simulation (P-DNS) methodology towards the study of multiphase flows was introduced. An RVE has been designed for the isolated modeling of the behavior of the finer scales of two-phase flows. The specific case of semidiluted water droplets in air was studied through a systematic off-line analysis of the RVEs. A database was built with mean statistics of the carrier flow dynamic response to a broad set of boundary conditions and features of the suspended particles. The database results confirmed that the onset and evolution of instabilities in multiphase flows can differ drastically from those in homogeneous flows. In agreement with previous numerical and experimental references, it was found that if mass loading is high enough, large particles can induce turbulence augmentation while small particles can induce attenuation. Moreover, the relationship between the dispersed phase volume fraction, the spread of diameters, and shear stresses applied to the system were quantified and reported. A synthetic model has been developed using an alternating digital tree to quickly interpolate high-dimensional data. This allowed us to generate high-fidelity macroscopic flow data to conduct multiscale simulations. In this regard, the P-DNS framework was evaluated in the laden turbulent flow between plates. Here it was shown that the P-DNS ability for turbulence modulation modeling improves the accuracy of the predictions for measured quantities. The approach presented is general, and may be seen as a first step in the development of a novel vademecum-based computational method to study multiphase flows. Any advance of computational techniques to help mitigate the effects of the current pandemic crisis should be welcome. In this work, advances towards the understanding and modeling of the turbulence modulation in multiphase flows have been presented. It was shown that this phenomenon rarely occurs when the pathogens are spread in air, which validates the use of standard numerical methodologies for turbulence modeling in these cases. However, the technologies here developed can be directly employed to give insight into other applications of interest in the battle against infectious diseases, such as the regional aerosol deposition in human airways, the design of dry powder inhalers, or the optimization of packaging for intravascular drug delivery, among others. Detailed simulation of viral propagation in the built environment Can indoor sports centers be allowed to re-open during the Covid-19 pandemic based on a certificate of equivalence? On coughing and airborne droplet transmission to humans Numerical modeling of the distribution of virus carrying saliva droplets during sneeze and cough Can CFD establish a connection to a milder Covid-19 disease in younger people? Aerosol deposition in lungs of different age groups based on Lagrangian particle tracking in turbulent flow Modelling aerosol transport and virus exposure with numerical simulations in relation to SARS-CoV-2 transmission by inhalation indoors High fidelity modeling of aerosol pathogen propagation in built environments with moving pedestrians Turbulent dispersed multiphase flow On models for turbulence modulation in fluidparticle flows A study of turbulence modulation models for gas-particle flows Classification of turbulence modification by dispersed spheres using a novel dimensionless number Turbulent modulation in particulate flow: a review of critical variables Particlesinduced turbulence: a critical review of physical concepts, numerical modelings and experimental investigations Fully resolved simulations of turbulence modulation by high-inertia particles in an isotropic turbulent flow Effects of turbulence modulation and gravity on particle collision statistics Direct numerical simulation of turbulent flows laden with droplets or bubbles Modeling the hydrodynamics of multiphase flow reactors: current status and challenges Multiscale modeling of multiphase-flow with complex interactions A multiple resolution approach using adaptive grids for fully resolved boundary layers on deformable gas-liquid interfaces at high schmidt numbers Multiscale computations of thin films in multiphase flows Particle-resolved direct numerical simulation for gas-solid flow model development Turbulence modulation in dilute particle-laden flow Modeling of collisional transport processes in spray dynamics CFD-dem simulation of turbulence modulation in horizontal pneumatic conveying Point-particle DNS and LES of particleladen turbulent flow-a state-of-the-art review Ryzhakov P (2020) A pseudo-DNS method for the simulation of incompressible fluid flows with instabilities at different scales The pseudodirect numerical simulation method for multi-scale problems in mechanics An improved version of the pseudo-direct numerical simulation method (P-DNS) for solving turbulent flows The variational multiscale method-a paradigm for computational mechanics Large eddy simulation and the variational multiscale method FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials On predicting particle-laden turbulent flows A local constitutive model for the discrete element method. Application to geomaterials and concrete Particle-laden turbulent flows: direct simulation and closure models The importance of the forces acting on particles in turbulent flows Integrable form of droplet drag coefficient Sobre la distribución del tamaño de las gotas en un aerosol (spray) Simulación de la inyección directa de combustible en motores de combustión interna A second-order in time and space particle-based method to solve flow problems on arbitrary meshes On mesh-particle techniques Cluster seshat (CIMEC UNL/CONICET) Pseudodirect numerical simulation (P-DNS) databases. Mendeley Data V2 Effect of inter-particle spacing on turbulence modulation by Lagrangian PIV Effect of particle size on modulating turbulent intensity An alternating digital tree (ADT) algorithm for 3D geometric searching and intersection problems Turbulence modification and heat transfer enhancement by inertial particles in turbulent channel flow Turbulence modulation for gas-particle flow in vertical tube and horizontal channel using four-way Eulerian-Lagrangian approach The ventilation of buildings and other mitigating measures for Covid-19: a focus on wintertime Computational methods for fluid dynamics Ranz and Marshall correlations limits on heat flow between a sphere and its surrounding gas at high temperature Exhaled droplets due to talking and coughing Mechanistic insights into the effect of humidity on airborne influenza virus survival, transmission and incidence Acknowledgements This research was partially funded by the projects PRECISE (BIA2017-83805-R) and PARAFLUIDS (PID2019-104528RB-I00) of the Agencia Española de Investigación (AEI) . The authors also acknowledge the financial support from the CERCA programme of the Generalitat de Catalunya, and from the Spanish Ministry of Economy and Competitiveness through the "Severo Ochoa Programme for Centres of Excellence in R&D" (CEX2018-000797-S). J. Gimenez thanks the Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT) from Argentina for funding this research via the project PICT-2018 N • 02464. On behalf of all authors, the corresponding author states that there is no conflict of interest.