Mathematical Modeling of Gas-Solid Two-Phase Flows: Problems, Achievements and Perspectives (A Review)

: Mathematical modeling is the most important tool for constructing theories of different kinds of two-phase ﬂows. This review is devoted to the analysis of the introduction of mathematical modeling to two-phase ﬂows, where solid particles mainly serve as the dispersed phase. The main problems and features of the study of gas-solid two-phase ﬂows are included. The main characteristics of gas ﬂows with solid particles are discussed, and the classiﬁcation of two-phase ﬂows is developed based on these characteristics. The Lagrangian and Euler approaches to modeling the motion of a dispersed phase (particles) are described. A great deal of attention is paid to the consideration of numerical simulation methods that provide descriptions of turbulent gas ﬂow at different hierarchical levels (RANS, LES, and DNS), different levels of description of interphase interactions (one-way coupling (OWC), two-way coupling (TWC), and four-way coupling (FWC)), and different levels of interface resolution (partial-point (PP) and particle-resolved (PR)). Examples of studies carried out on the basis of the identiﬁed approaches are excluded, and they are also excluded for the mathematical modeling of various classes of gas-solid two-phase ﬂows.


Introduction
Continuum flows, which carry dispersed admixtures, include [1][2][3][4] sandstorms, tornadoes, volcanic eruptions, forest fires, and precipitation in the form of hail, snow, etc. Examples of technical devices that use two-phase currents include the paths of solid-fuel jet engines, pneumatic devices, and many others.
Today, we can state that there is continual growth of interest among researchers in the study of two-phase flows.This seems to be due to two factors.First, in recent years, there has been a tremendous growth in the possibilities of both the mathematical and physical (experimental) modeling of two-phase flows.Second, the range of problems under study for various types of two-phase flows is expanding.The second circumstance largely stems from the first circumstance.
This review differs from other reviews on this topic, which have largely been devoted to some narrower problems (for example, the problem of the influence of particles on gas turbulence, the problem of particle clustering, etc.).This review attempts to analyze the state of mathematical modeling in a broader sense.In the review, a classification of two-phase turbulent flows according to particle inertia is constructed.This classification covers almost the entire range of particle inertia.This classification is of great prognostic value since it offers new dimensionless criteria that allow one to analyze both the existing results at different qualitative levels and conduct new studies regarding various classes of two-phase flows determined by the concentration and inertia of particles.
The inertia of particles (which is primarily determined by their size and density) can vary in a colossal range (many orders of magnitude).One-phase flows are characterized by a number of space-time scales that are determined by the magnitude of their inherent flow velocity, flow regime (laminar, transient, and turbulent), flow geometry, etc.For the accurate modeling of particle motion, it is necessary to consider particles' interactions at different scales; these interactions are determined by (1) averaged movement, (2) widescale fluctuation motions, (3) fine-scale fluctuation motions, (4) different instabilities (for example, Tollmin-Schlichting instability in boundary layers, Taylor-Gertler instability in pipes, and Kelvin-Helmholtz instability in pure shear layers), etc.
It important to notice that many of these forces, in one form or another, contain the velocity of the carrier phase u(τ), which is a random variable in a turbulent flow.Therefore, a question often arises regarding the applicability of a particular expression that is obtained theoretically or empirically for other conditions to calculating the influence of forces (for example, for laminar flow or in the absence of velocity shift).

Multiplicity of Modeling Parameters
The main parameters are as follows: (1) three components of the average speed (U i , U j , and U k ), (2) three components of the fluctuation (rms) velocity ((u i 2 ) 1/2 , (u j 2 ) 1/2 , and (u k 2 ) ), (3) average temperature (T), ( 4) fluctuation (rms) in temperature ((t 2 ) ), (5) double correlations of various components of fluctuation velocities (components of the Reynolds stress tensor) (u i u j ), (6) double correlations of fluctuation velocity and fluctuation temperature (u i t , u j t ), etc.This multiplicity is explained by the fact that similar parameters for the dispersed phase are added to the parameters indicated above (for example, averaged (V i , V j , V k ) and fluctuation velocities ((v i 2 ) 1/2 , (v j 2 ) 1/2 , and (v k 2 ) 1/2 ) and particle temperatures (T p , (t p 2 ) ), etc.), including their size, size distribution, averaged and fluctuation (rms) concentrations (Φ, (φ 2 ) ), and the parameters for the carrier phase (in the presence of particles), the heat of phase transitions, and many others.

Multiplicity of Collision Processes
The main factors contributing to the occurrence of collisions between particles are listed as follows: (1) polydispersity, which leads to a difference in the averaged velocities; (2) the influence of the gradient of the averaged velocity of the carrier phase; (3) the influence of gravity (Archimedes); (4) the turbulent transport effect, which leads to the appearance of relative velocity between nearby particles; (5) the effect of clustering, i.e., an increase in the concentration of the dispersed phase in local regions of space; (6) electrostatic interactions; and (7) Brownian motion.

Multiplicity of Phase and Chemical Transformations
Phase transformations are not considered in this review because the subject of this review is gas flows with solid particles (gas-solid two-phase flow).However, for nonisothermal two-phase flows, particle melting may occur during interphase heat exchange.The melting of particles in the gas stream leads to the transition of a gas-solid two-phase flow into a gas-liquid two-phase flow.The subsequent process of the crystallization (solidification) of droplets can cause the flow to "return" to its initial state.

Multiplicity of Dimensionless Parameters
An example of such parameters is represented by the numerous Stokes numbers (see Section 3.3) that characterize the inertia of the dispersed phase with ratios of the various scales of the carrier gas, the Reynolds number of the particle, etc.

Main Characteristics of Two-Phase Flows
This section presents classifications of two-phase flows developed on their basis.

Particle Concentrations
The possible varieties of particle concentrations (classification) are given in [8,9].There are three classes of two-phase flows: (1) dilute two-phase flows without the reverse effect of the dispersed phase; (2) dilute two-phase flows with the reverse effect of the dispersed phase; and (3) dense two-phase flows with intense collisional interactions between particles.

One-Way Coupling
To model the motion of particles in dilute two-phase streams (dispersed phase Φ ≤ 10 −6 ), that is, streams with a small Φ, "one-way coupling" (OWC) is applied.

Four-Way Coupling
Further growth (Φ > 10 −3 ) requires the inclusion of the contribution of interparticle interactions to the process of the momentum and energy transfer of the dispersed phase [15][16][17][18].The chaotic motion of particles during their interaction is called "pseudoturbulence" to distinguish it from the actual turbulent fluctuations of particle velocities associated with their involvement in the turbulent motion of the carrier flow.
It should be noted that there is a clustering phenomenon in two-phase flows that consists of a sharp increase in the concentration of particles in local areas.This significant rise in Φ leads to an increase in the probability of particle collisions, even in a dilute two-phase flow.
In reference to what has been discussed above, it is clear that in flows with a small mass of dispersed phase content, in which particles do not undergo collisions and do not have a reverse effect on the flow of the carrier continuous medium, clustering phenomena can lead to flow restructuring.The formation of local areas of increased particle concentrations has been revealed experimentally or via calculation in various flows, including homogeneous isotropic turbulence [19,20], shear flows in pipes (channels) [21,22], flows in boundary layers [23], jet flows, traces behind streamlined bodies, flows around blunted bodies [24], and free, concentrated vortices [25][26][27].

Particles' Dynamic Relaxation Time
The inertia of particles is the time of dynamic relaxation τ p , which is represented in the following form where ρ p is the physical density; τ p0 is the time of dynamic relaxation of the Stokes particle; µ is the dynamic viscosity.

Stokes Numbers
There are three dimensionless criteria listed in [28,29]: Stk f , Stk L , and Stk K , representing Stokes numbers in averaged, large-scale, and small-scale fluctuation movements, respectively.
In the equation above, T i denotes a given characteristic time of the carrier phase.

Stokes Number in Time-Averaged Motion
We can apply the Stokes number to averaged motion, which we express as where T f is the characteristic time of the carrier phase in the averaged motion.

Stokes Number in Large-Scale Fluctuation Motions
In this case, the Stokes number assumes the following form where T L is the characteristic time of the carrier gas in a large-scale fluctuation motion (temporal Lagrangian integral turbulence scale).

Stokes Number in Small-Scale Fluctuation Motions
The inertia of particles in small-scale fluctuation motions can be characterized by the Stokes number we represent as where τ K is the Kolmogorov time scale of turbulence.

Classification of Turbulent Two-Phase Flows According to Particle Inertia
We will briefly describe the classification of turbulent two-phase flows according to particle inertia (see Figures 1 and 2) depending on Stokes numbers [28,29].

Classification of Turbulent Two-Phase Flows According to Particle Inertia
We will briefly describe the classification of turbulent two-phase flows accordin particle inertia (see Figures 1 and 2) depending on Stokes numbers [28,29].Flow around fixed "frozen" particles.In this case, the particles have an extrem large amount of inertia, which remains completely static and whose temperature d not change.An analogue of such a hypothetical class of two-phase flows is a one-ph

Classification of Turbulent Two-Phase Flows According to Particle Inertia
We will briefly describe the classification of turbulent two-phase flows according to particle inertia (see Figures 1 and 2) depending on Stokes numbers [28,29].Flow around fixed "frozen" particles.In this case, the particles have an extremely large amount of inertia, which remains completely static and whose temperature does not change.An analogue of such a hypothetical class of two-phase flows is a one-phase Flow around fixed "frozen" particles.In this case, the particles have an extremely large amount of inertia, which remains completely static and whose temperature does not change.An analogue of such a hypothetical class of two-phase flows is a one-phase flow in heat exchangers, where fixed pipes act as such particles, through which the working fluid moves.

Lagrangian and Eulerian Modeling of Two-Phase Flows
Mathematical modeling methods play an important role in the study of the processes of the motion of solid particles.Detailing a large number of processes in which information about each individual is not always indisputable can lead to a decrease in the reliability of the created model.

Reasons for Considering of the Two-Phase Nature of Tornados
It is clear that the attempts to describe all the varieties of two-phase flows using all models can hardly be justified.As a consequence, for certain classes of flows (see Section 3), which are characterized primarily by the concentration of the dispersed phase and Stokes numbers, specific models should be preferred.

Lagrangian Modeling
The system of equations is as follows (example taken from [30]) where x p is the position vector (radius vector) of the particle, v is the instantaneous velocity vector of the particle, ω p is the angular velocity vector of the particle, m p is the mass of the particle, and I is the moment of inertia of the particle.Equation ( 8) describes the change in angular velocity of the particle due to viscous interactions with the surrounding gas.
Due to the viscosity of the liquid, a moment of rotation T acts on a rotating particle.

Eulerian Modeling
Let us briefly consider the current approaches to constructing continuum equations for the motion of dispersed impurity and analyze the features of describing their behavior for different classes of two-phase flows.
Algebraic and differential models.There are two main approaches to determining the velocity correlations of dispersed phases.One of them is presented in [31,32] where A is the function of a particle's involvement in gas fluctuation movement.
The other approach consists of applying gradient relations like the Boussinesq relations for a single-phase flow [33] or in the form presented in [34,35] where ν p is the turbulent viscosity coefficient of the dispersed phase.Different methods for determining the value of ν p have been described in the literature [34,35].
According to this first method, to transform stochastic equations like the Langevin equation into a kinetic equation for a group of particles, a probability density function (PDF) that describes the coordinate x, velocity v, and temperature distribution of particles t p is introduced: where averaging is not carried out over time but over realizations of the random fields of the carrying gas flow.
The equation for the PDF is then used, as presented in [36]: where Here, T pL and T pLt are the times of interaction between particles (droplets) and energyintensive fluctuations of velocity and temperature, respectively, for the non-inertial impurities T pL = T L and T pLt = T Lt .
The system of Equations ( 13)-( 15) is not closed, as the equations contain information related to particle involvement in fluctuating motion turbulent stresses v i v j and turbulent heat flux v j t p in the dispersed phase as well as the turbulent diffusion of momentum and heat arising from non-uniform particle concentrations.
A mathematical description of momentum and heat transfer processes in the dispersed phase of varying complexity was developed in [36]: where is the energy of the fluctuations of the dispersed phase velocity.
The second method.The methods presented in [48,49] allow for the acquisition of equations for joint PDF distributions of the velocity and temperature of the dispersed impurity [50].
The third method.The third method constitutes the construction of a closed kinetic equation based on the expansion of the characteristic functional into a series of cumulants [51,52].

Advantages and Limitations of Lagrangian and Eulerian Modeling
Let us consider Euler-Lagrange and Euler-Euler models with respect to describing the motion of flows of continuous media with solid particles, droplets, and bubbles [36].
The advantage of Euler-Lagrange (trajectory and stochastic) models is their ability to obtain detailed statistical information about the motion of individual particles via integrating the equations of motion.
It should be noted that with an increase in the concentration of the dispersed phase, there are also difficulties in using Euler-Lagrange models [53].

Description of the Gas Flow Carrying the Particles
An increase in the volume fraction of the dispersed phase can affect the carrier medium (see Section 3.1.2).Let us consider the motion of a continuous medium (gas) in the presence of particles when the particles begin to have a reverse influence on the medium's characteristics.

Actual Equations
The relevant equations are as follows: The continuity equation (Equation ( 17)) has a form similar to the single-phase flow equation.

Time-Averaged Equations
The resulting averaged continuity, motion, and energy equations are expressed as follows: Let us assume that the distributions of the averaged velocities and particle concentrations are known.We need to determine the turbulent gas stresses u i u j and the turbulent heat flux u j t as well as the correlations between particle concentration fluctuations and gas velocity and temperature fluctuations φ u i and φ t , which can be represented as follows [54][55][56]:

Equations for the Reynolds Stresses
Subtracting Equations (20)- (22) from Equations ( 17)-( 19) yields the following equations: Equation ( 25) differs from the corresponding equation of a single-phase flow via the presence of the last group of terms on the right-hand side, which takes into account the dynamic influence of particles on the carrier flow.
The system of Equations ( 20), ( 21), (23), and ( 25) is unclosed because Equation ( 25) contains unknown triple correlations of velocity fluctuations of the carrier phase, serving as correlations related to fluctuations in particle concentration and velocity.
One-parameter models.The model based on the turbulence energy equation is the most common (as in the case of single-phase flow).It is necessary to multiply the equation of fluctuation motion by u i , sum over i, and then average the result.The equation will then have the following form: The Equation ( 26) can be rewritten concisely: where the additional dissipation ε p, caused by the presence of a dispersed phase, assumes the following form: There have been several studies (e.g., [60][61][62]) in which the authors attempted to estimate the magnitude of the terms on the right-hand side (28) for different classes of two-phase flows.This means the second and third terms on the right-hand side (28) are small compared to the first term.Thus, in the implementation of quasi-equilibrium and non-equilibrium flows (see Figures 1 and 2), the first term on the right-hand side (28) will play a determining role in the process of dissipation: Considering this mechanism [63-65] leads to Equation ( 27) assuming the following form where P p is the additional generation caused by the presence of particles.Two-parameter models.As in the study of single-phase turbulent flows, the two-parameter k − ε-turbulence model has become the most widespread.
By analogy with the equation for a single-phase flow in the case of two-phase flow, we obtain where ε εp is the reduction in dissipation due to the presence of particles.The equation for ε εp is most often represented as follows [66,67]: where the constant C ε3 can assume the following values: C ε3 = 1.0 [68], C ε3 = 1.2 [66], and C ε3 = 1.9 [69].

Methods of Numerically Modeling Two-Phase Flows
The main methods of numerically modeling two-phase flows are presented below.

Particle-Resolved DNS
Particle-resolved direct numerical simulation (PR DNS) is the method that most fully describes the physics of two-phase flows.In this method, flow around each particle is allowed to occur.In this case, the behavior of each particle is determined by both external acting forces and the aerodynamic drag force from the carrier gas (determined in the calculation process).This method is also applicable to the calculation of more complex two-phase flows carrying droplets or bubbles, where the interfacial surface may deform.This deformation is calculated using the aerodynamic force determined in the calculation process.
A well-known limitation of this method is the following circumstance.It is possible to calculate the movement of gas around each particle when the step of the computational grid is small compared to the particle size.The application of this method is complicated when the particle size exceeds the size of the smallest turbulent vortices (Kolmogorov microscale) and the number of particles is large.
To date, various numerical methods and algorithms have been developed to implement PR DNS.In [70], this method was used to calculate the force acting on a single stationary particle in decaying homogeneous isotropic turbulence (DHIT).Effective methods of implementing PR DNS include the immersed boundary method [71], which uses a Cartesian grid throughout the computational domain, and the lattice Boltzmann method [72], which also uses a Cartesian grid that is not aligned in terms of particle shape.Another method is Physalis [73], which uses a local analytical solution to determine the flow around each particle.

Particle Point Methods
Lagrange's methods of description are the oldest methods of describing the motion of particles.These methods can be used to calculate the motion of millions of particles.The condition for the applicability of Lagrange's approaches is the smallness of the particle size compared to the Kolmogorov spatial scale.In this case, the particles can be considered as point particles.
The most important characteristic of particle inertia is dynamic relaxation time τ p .In the case of small values of τ p , the instantaneous velocity of the particle is close to the corresponding velocity of the carrier gas, and the particles are tracers.In this case, an equilibrium flow is realized (Section 3).With an increase in τ p , the particles cannot fully track the turbulent fluctuations of the gas, and a quasi-equilibrium flow is realized.In this case, to describe the motion of particles, it is necessary to integrate the equations of their motion.
Lagrange's models can have different levels of description of the turbulence of the carrier gas, ranging from Reynolds-averaged Navier-Stokes equations (RANS), wherein only fields of averaged turbulence characteristics are calculated, to large-eddy simulation (LES) and direct numerical simulation (DNS), wherein only large vortices and vortices of all scales are resolved (up to Kolmogorov), respectively (see Figure 3).The particle concentration determines the required level of description of the interfacial interaction (see Figures 1 and 2): (1) the mode of movement of single particles, where their presence does not have a reverse effect on the characteristics of a non-existent gas (one-way coupling-OWC); (2) the mode of weakly dusty flow (dilute two-phase flows), with a reverse effect of particles (two-way coupling-TWC), and the mode of highly dusty flow (dense twophase flows), where collisions of particles with each other play a significant role (four-way coupling-FWC).
centration determines the required level of description of the interfacial interaction (see Figures 1 and 2): (1) the mode of movement of single particles, where their presence does not have a reverse effect on the characteristics of a non-existent gas (one-way coupling-OWC); (2) the mode of weakly dusty flow (dilute two-phase flows), with a reverse effect of particles (two-way coupling-TWC), and the mode of highly dusty flow (dense two-phase flows), where collisions of particles with each other play a significant role (four-way coupling-FWC).

Direct Numerical Simulation
To date, there has been a significant amount of work in which researchers have studied various problems regarding the physics of two-phase flows using DNS and by describing interphase interactions and interphase boundary at various levels.
One of the first papers in which the behavior of point particles (PP DNS) under damped homogeneous isotropic turbulence decaying homogeneous isotropic turbulence (DHIT) was studied was [74].In this study, the motion of 432 particles was studied at a very small Reynolds number ( 35 Re < λ ).Only linear aerodynamic drag was taken into account in the equations of particle motion.
In later studies [75,76] devoted to the study of particle motion, both regarding forced homogeneous isotropic turbulence forced homogeneous isotropic turbulence (FHIT) and with respect to damped homogeneous isotropic turbulence (DHIT), emphasis was placed on the study of the possibilities of various methods of interpolation (linear interpolation, high-order Lagrangian interpolation, and high-order Hermite interpolation) of the gas velocity at the location of the particle.
A more complex case of turbulent two-phase flow turbulent flow in a channel is considered in [77,78].In [77], in addition to the aerodynamic drag force, the Safman force was also taken into account, and in [78], a more advanced Fourier-Chebyshev pseudo-spectral method was used to interpolate the gas velocity at the particle's location.To date, there have been numerous studies using the PP OWC DNS method of two-phase flows in the channel [79], in pipes [80][81][82], under FHIT [19,20], and under DHIT [83].
With an increase in the concentration of particles, the particles begin to have the opposite effect on the characteristics of the carrier gas flow (see Section 3), so TWC DNS is necessary.This introduces additional difficulties in mathematical modeling.Firstly, in the equation of the motion of a particle, it is not the initial (inherent in a single-phase flow) velocity of the gas that should be present but the "new" velocity of the flow caused by the presence of particles.In [84], it was suggested that the difference between these velocities is small if the diameter of the particles is smaller than the size of the numerical grid, L d p < .This condition is almost always satisfied in the case of PP DNS.Secondly, it is necessary to introduce a source term in the equations of gas motion [85].If the particle

Direct Numerical Simulation
To date, there has been a significant amount of work in which researchers have studied various problems regarding the physics of two-phase flows using DNS and by describing interphase interactions and interphase boundary at various levels.
One of the first papers in which the behavior of point particles (PP DNS) under damped homogeneous isotropic turbulence decaying homogeneous isotropic turbulence (DHIT) was studied was [74].In this study, the motion of 432 particles was studied at a very small Reynolds number (Re λ < 35).Only linear aerodynamic drag was taken into account in the equations of particle motion.
In later studies [75,76] devoted to the study of particle motion, both regarding forced homogeneous isotropic turbulence forced homogeneous isotropic turbulence (FHIT) and with respect to damped homogeneous isotropic turbulence (DHIT), emphasis was placed on the study of the possibilities of various methods of interpolation (linear interpolation, high-order Lagrangian interpolation, and high-order Hermite interpolation) of the gas velocity at the location of the particle.
A more complex case of turbulent two-phase flow turbulent flow in a channel is considered in [77,78].In [77], in addition to the aerodynamic drag force, the Safman force was also taken into account, and in [78], a more advanced Fourier-Chebyshev pseudospectral method was used to interpolate the gas velocity at the particle's location.To date, there have been numerous studies using the PP OWC DNS method of two-phase flows in the channel [79], in pipes [80][81][82], under FHIT [19,20], and under DHIT [83].
With an increase in the concentration of particles, the particles begin to have the opposite effect on the characteristics of the carrier gas flow (see Section 3), so TWC DNS is necessary.This introduces additional difficulties in mathematical modeling.Firstly, in the equation of the motion of a particle, it is not the initial (inherent in a single-phase flow) velocity of the gas that should be present but the "new" velocity of the flow caused by the presence of particles.In [84], it was suggested that the difference between these velocities is small if the diameter of the particles is smaller than the size of the numerical grid, d p < L. This condition is almost always satisfied in the case of PP DNS.Secondly, it is necessary to introduce a source term in the equations of gas motion [85].If the particle is smaller than the Kolmogorov scale (d p < η K ), then there are no special problems.Otherwise, (d p > η K ) raises the question of the relevance of the assumption of point particles.In [86,87], calculations of a two-phase flow containing a lot of very small particles at Φ = O(10 −4 ) were carried out, and the number of particles was commensurate with the number of cells of the computational grid.
Examples of studies in which PP TWC DNS modeling was performed include [88][89][90][91].In [90], turbulent flow in a channel was studied.The volume concentration of particles was equal to Φ ≈ 10 −4 .It was assumed that the particles were of the Stokes type (adhering to the linear law of resistance).It was found that in the case of small particles (d p < η K ), their presence suppressed turbulence, and, on the contrary, the presence of relatively large particles (d p > η K ) caused turbulence intensification.In [89,90], a two-phase flow in a channel at Re τ = 180, determined from the half-height of the channel, was simulated.It was revealed that the presence of particles reduces resistance and leads to an increase in longitudinal fluctuations of the gas velocity.At the same time, the presence of particles caused a decrease in gas velocity fluctuations in the other two directions and significantly reduced Reynolds stresses.In [91], a two-phase turbulent flow in a channel was simulated for the same Reynolds number (Re τ = 180), taking into account the nonlinearity in the particle drag law (regarding non-Stokes particles).It was found that particles with small Stokes numbers increased the intensity of turbulence, Reynolds stresses, and the level of viscous dissipation.At the same time, particles with large Stokes numbers led to a decrease in the intensity of turbulence.
The numerical concentration of particles N 0 and the number of particles in the Kolmogorov vortex N η are related, as N η = N 0 η 3 K .The calculations [92] carried out allowed for the clear identification of two regimes.At Stk K < 1, the presence of particles results in a decrease in the decay of turbulent energy (first mode).At Stk K > 1, particles accelerate the decay of turbulence (second mode).In [93], the results regarding PR TWC DNS with respect to the direction of turbulent two-phase upward flow in a vertical channel are presented.
A further increase in particle concentration necessitates the consideration of interparticle collisions (see Section 3), which requires conducting FWC DNS.Intense interparticle collisions influence particle motion statistics and, consequently, their backreaction with respect to gas flow.This greatly complicates mathematical modeling.Currently, several stochastic approaches have been developed in order to depart from simple deterministic calculations of pairwise particle collisions, which require immense computational time.
Examples of studies in which PP FWC DNS modeling was performed include [94,95].In [94], the mathematical modeling of turbulent two-phase flow in a vertical pipe in the presence of small heavy particles was carried out over a wide range of variations in particle mass concentration (M = 0.1 − 30).Various modeling techniques for real wall roughness were used to better match the results with experimental data.It was found that the results strongly depend on the model of wall roughness used rather than on the variation of the parameters characterizing the inter-particle collision process.The calculations also revealed a decrease in turbulence intensity with an increase in particle mass concentration.In [95], the modeling of turbulent two-phase downward flow in a channel was performed at Re τ = 642 and particle mass concentration M = 0.8.The calculations were carried out for smooth and rough walls, where roughness was modeled by placing fixed tiny particles on the wall.It was discovered that rough walls enhance the suppression of turbulence caused by the presence of particles in the flow.
In [96], the interaction between a stationary homogeneous isotropic turbulent (HIT) flow and inertial particles while accounting for inter-particle collisions (PP FWC DNS) was studied via direct numerical simulation (DNS).The calculations were performed for a periodic cubic box with a size of 128 3 for two values of the Reynolds Taylor number (Re λ = 35.4 and Re λ = 58) while varying the volume concentration of particles (from Φ = 1.37 • 10 −5 to Φ = 8.22 • 10 −5 ) and the Stokes number (Stk K = 0.19 − 12.7).Elastic spherical particles with a diameter of d p = 67.6 µm, corresponding to d p /η K = 0.1, served as the dispersed phase.The Stokes number was varied by changing the particle density over a wide range, namely, ρ p = 150 − 18, 000 kg/m 3 .The results [96] showed that dissipation decreases up to 32% with an increase in the Stokes number and the volume concentration of particles.It was shown that this maximum reduction in dissipation is overestimated by 7% when accounting for inter-particle collisions.The spectral analysis revealed a transfer of energy from large to small scales due to particle flow, which explains the difference in dissipation.

Large Eddy Simulation
The use of the Reynolds-averaged Navier-Stokes (RANS) equations requires far fewer computational resources, which is its undeniable advantage.This approach has been successfully used in practical calculations.However, the Reynolds equations and turbulence models used to solve the equations do not have acceptable universality and, therefore, cannot be used to solve a wide range of applied problems.
Large Eddy Modeling (LES) is a compromise between DNS and RANS.This approach is limited to the study of flows only on scales exceeding some given value.In the LES model, the Navier-Stokes equations, which are filtered in space, are solved, and only large eddies are allowed to move.Small eddies have a more versatile structure and are modeled using subgrid scale models.
The LES-based solution contains richer information than the RANS-based solution.It contains not only the average flow characteristics (velocity, temperature, pressure, and concentration fields) and Reynolds stress distributions but also spectral characteristics (velocity, temperature, and pressure fluctuation spectra), two-point moments, and temporal and spatial scales of turbulence.Many of these characteristics are important for engineering applications.For example, temperature fluctuations play a fundamental role in the calculation of chemically reactive flows.
LES is similar to DNS, but the grid used in the process is much larger.Small vortices are approximated using a subgrid-scale (subgrid-scale) model of turbulence.The most commonly used model is the dynamic Smagorinsky model of vortex viscosity [97].Other well-known models are based on scale similarity assumption [98], Taylor series expansion [99], or approximate deconvolution [100].
One of the early works that used the PP OWC LES method was the study presented in [101].In this work, particle dispersion was investigated regarding the case of homogeneous shear flow.The authors did not use the term LES, but they considered the spatially averaged Navier-Stokes equation for the gas and used time-and space-varying coefficients for the small-scale vortices.The calculations were carried out for only 48 passive particles, and the influence of subgrid scales on their motion was not considered.
The work presented in [102] investigated particle dispersion in a turbulent pipe flow using PP OWC LES and DNS methods for different Reynolds numbers.The equation of particle motion considered drag force, lift force, and buoyancy force.Due to very low particle volume concentrations, their back-reactions on the gas and interparticle collisions were not considered.Moreover, the influence of subgrid scales of the gas velocity was also not considered.The main conclusion of this work was that the dynamic relaxation time of particles plays an important role in their sedimentation.
In [103] studied particle motion in a vertical channel with a very low particle volume concentration using the PP OWC LES method.The dynamic Smagorinsky approach, previously developed in [104], was used as a subgrid-scale model.A comparison of the results obtained with those of DNS-based modeling showed good agreement.It should be noted that this work investigated the influence of subgrid-scale velocities on particle settling.For this purpose, an additional equation for the transport of kinetic energy of subgrid-scale turbulence was used, revealing only a minor effect on the calculation results.
In [105], the authors performed calculations of a two-phase flow for the case of forced homogeneous isotropic turbulence (FHIT), for which the reverse influence of particles on gas was taken into account, i.e., using the PP TWC LES method.The authors applied various subgrid-scale models to the equations of motion of the carrying gas.A very important conclusion was drawn: an increase in particle mass concentration leads to a decrease in the weighting coefficients in the dynamic model of vortex viscosity.As a consequence, the calculation error due to the use of subgrid-scale models for the two-phase flow was reduced compared to the single-phase flow.
The PP FWC LES method was used to account for particle collisions in [106] in the study of two-phase flow in a vertical channel at Re τ = 644 and a volume concentration of up to Φ = 1.4 × 10 −4 .The impact of drag force, gravitational force, and lift forces (due to the presence of gas velocity shear and particle rotation) on particle behavior was considered in the work.A deterministic model was used to account for particle collisions.Conclusions were drawn about the significant influence of inter-particle collisions on the statistical characteristics of particle motion, including the concentration magnitude.
In [107], two-phase flow calculations were performed using the PP FWC LES method in a channel with a very high particle volume concentration Φ = 1.3 × 10 −2 .Among all the forces, only the drag force and gravitational force were considered.The calculations showed that particles have a colossal effect on turbulence, leading to a thinning of the boundary layer, an increase in gas velocity fluctuations in the longitudinal direction, and, conversely, a reduction in gas fluctuations in the two other directions.
In [108], the parameters of a two-phase flow in a channel were calculated at a particle volume concentration of Φ = 4.8 × 10 −4 and a Reynolds flow rate of Re = 42, 000, both of which were set based on the height of the channel.The authors separately considered the effects of particle back-influence on gas and inter-particle collisions (PP TWC LES and PP FWC LES).They also emphasized the use of various particle collision models (hard-sphere and soft-sphere), different wall conditions (smooth and rough), and different subgrid viscosity models (the Smagorinsky model and a dynamic model).The calculation results showed that the differences when using different collision and subgrid models were insignificant.At the same time, the consideration of particle collisions and wall roughness leads to better agreement with the available experimental data.
In [109], PP FWC LES was performed for a two-phase flow with particles at a volume concentration of Φ = 7.3 × 10 −5 and a Reynolds flow rate of Re = 11, 900, both of which were set based on half of the height of the channel.The authors used a subgrid model developed earlier in [110] for the particle motion equation as well as a deterministic model to calculate inter-particle collisions.It was shown that with such a small volume concentration of particles, their influence on gas turbulence is negligible.At the same time, it was found that the influence of particle collisions plays a significant role.A good agreement was found between the results and the DNS data described in [86] as well as with the experimental data.
The authors of [109] later performed PP FWC LES simulations of a two-phase flow [111] in a horizontal pipe at a Reynolds number of Re = 120, 000, which was set based on the pipe diameter.The peculiarity of this study was the consideration of particle polydispersity and particle rotation as well as the inclusion of not only the drag force but also the lift force of Saffman and the Magnus force.Wall roughness was modeled by introducing coefficients of normal and tangential velocity restitution that differ from unity and by taking into account the so-called shadow effect at small wall collision angles.
In [112], PP FWC LES of a two-phase flow in a channel was performed with the presence of particle agglomeration effects.The main technique that allowed for the consideration of the appearance of particle agglomerates in the flow after their collision was the introduction of van der Waals forces, which are responsible for the phenomenon of cohesion.Various aerodynamic and energy systems can serve as examples of the use of two-phase flows in the future [113][114][115][116][117][118][119].It should be mentioned that mixing and chemical reactions can occur in a two-phase flow.The coupling of CFD with chemicals can be used to evaluate the performance of devices [120][121][122].
The following conclusions can be drawn from the above description and analysis of works devoted to the mathematical modeling of two-phase flows.
As is known in this field, the Reynolds number is the most important criterion for single-phase flows.It is known that high Reynolds number values limit the use of the DNS method, as the requirements for computing power increase sharply.The Reynolds number of a particle Re p determines the mode of flow around a particle (from the Stokes mode to the mode of formation of turbulent wakes behind moving particles) and is the most important criterion for two-phase flows.It is important to note that the Reynolds number of a particle can be determined not only from the difference in average velocities but also from the difference in fluctuating velocities between the carrier gas and the particles.
In the overwhelming majority of works, the ratio of a particle's diameter to the Kolmogorov space scale of turbulence is used as the only "two-phase" criterion.This is obviously not sufficient, especially in the case of studying the flow in channels (pipes), where there may be a difference in the average velocities of the gas and particles.It seems appropriate to use other criteria more widely, such as the Reynolds number of a particle Re p , the Stokes number in average motion Stk f , and the Stokes number in large-scale fluctuating motion Stk L .This will allow for the mathematical modeling of various classes of flows in accordance with the developed classification of two-phase flows (see Figures 1 and 2).

Conclusions
Two-phase flows have a colossal distribution in nature and are widely used in practice.The extreme complexity of the physics of two-phase turbulent flows is determined by the factors described in detail in Section 2.2.In addition, in two-phase systems, it is necessary to consider the processes of mixing and chemical reactions [120][121][122], which are of great importance for the operation of a wide range of technical devices.All of the above information complicates the mathematical modeling of such flows.
In the last 20-30 years, there has been a tremendous growth in interest among numerous researchers in the numerical modeling of two-phase flows with particles.As a result, there has been significant progress in improving the methods and approaches for the mathematical modeling of such flows.Currently, there are advanced methods such as particle-resolved direct numerical simulation (PR DNS), which allows for the determination of local gas velocities influenced by the presence of particles and interphase interaction forces.This method has well-known limitations associated with its small number of particles and their "coarseness".However, DNS is a very computationally intensive method for solving practical problems.Therefore, in the near future, RANS and LES methods and the modeling of particle motion based on the Euler approach are likely to be more in demand and require further improvements.
In conclusion, we have formulated, in our opinion, the following promising directions for further progress in the field of the mathematical modeling of two-phase flows with particles: (1) The development of mathematical modeling methods for two-phase flows with relatively large particles (non-equilibrium flows) that only interact with large energycarrying vortices and are characterized by dynamic slippage (velocity difference) in relation to average motion.(2) The development of mathematical modeling methods for two-phase flows with large particles, which form turbulent wakes behind them.With the increase in particle concentration, these turbulent wakes will interfere with each other, and the particles will undergo collisions.(3) The development of mathematical modeling methods for two-phase flows containing particles of different sizes (polydisperse particles).Such flows are of interest to practicing engineers.Particles of different sizes will have different velocities and different effects on gas flow and tend to collide with each other at lower concentrations.
The development of mathematical modeling methods for two-phase flows with particles complicated by phase transitions (melting and subsequent evaporation) and chemical reactions (primarily combustion reactions).

Nomenclature d p
particle diameter, m η K Kolmogorov length scale, m x p particle radius vector, m u vector of actual velocity of gas, m/s v vector of actual velocity of particle, m/s u i , u j , u k actual velocity components of gas, m/s v i , v j , v k actual velocity components of particle, m/s U i , U j , U k time-averaged velocity components of gas, m/s V i , V j , V k time-averaged velocity components of particle, m/s u i , u j , u k fluctuation velocity components of gas, m/s v i , v j , v k fluctuation velocity components of particle, m/s

τ
is the Kolmogorov time scale of turbulence.

Figure 1 .
Figure 1.Classification of turbulent two-phase flows' dependence on particle inertia.

Figure 2 .
Figure 2. Simplified schemes of turbulent two-phase flows of different classes depending on p cle inertia: (a) equilibrium flow, (b) quasi-equilibrium flow, (c) nonequilibrium flow, (d) flow large particles, (e) and flow around fixed "frozen" particles.

Figure 1 .
Figure 1.Classification of turbulent two-phase flows' dependence on particle inertia.

Figure 1 .
Figure 1.Classification of turbulent two-phase flows' dependence on particle inertia.

Figure 2 .
Figure 2. Simplified schemes of turbulent two-phase flows of different classes depending on particle inertia: (a) equilibrium flow, (b) quasi-equilibrium flow, (c) nonequilibrium flow, (d) flow with large particles, (e) and flow around fixed "frozen" particles.

Figure 2 .
Figure 2. Simplified schemes of turbulent two-phase flows of different classes depending on particle inertia: (a) equilibrium flow, (b) quasi-equilibrium flow, (c) nonequilibrium flow, (d) flow with large particles, (e) and flow around fixed "frozen" particles.

Figure 3 .
Figure 3. Classification of approaches to numerical simulation of two-phase flows depending on different levels of turbulence description and interphase coupling.

Figure 3 .
Figure 3. Classification of approaches to numerical simulation of two-phase flows depending on different levels of turbulence description and interphase coupling.