A Hyperbolic Eulerian Model for Dilute Two-Phase Suspensions

Conventional modeling of two-phase dilute suspensions is achieved with the Euler equations for the gas phase and gas dynamics pressureless equations for the dispersed phase, the two systems being coupled by various relaxation terms. The gas phase equations form a hyperbolic system but the particle phase corresponds to a hyperbolic degenerated one. Numerical difficulties are thus present when dealing with the dilute phase system. In the present work, we consider the addition of turbulent effects in both phases in a thermodynamically consistent way. It results in two strictly hyperbolic systems describing phase’s dynamics. Another important feature is that the new model has improved physical capabilities. It is able, for example, to predict particle dispersion, while the conventional approach fails. These features are highlighted on several test problems involving particles jets dispersion and are compared against experimental data. With the help of a single parameter (a turbulent viscosity), excellent agreement is obtained for various experimental configurations studied by different authors


Introduction
Two phase diluted flows are present in many fundamental and industrial applications ranging from astrophysics, fluid mechanics, chemical engineering, combustion, nuclear engineering and so on.The Eulerian approach is widely used to deal with such flows, considering velocity non-equilibrium effects in a two-fluid description.It considers the gas phase governed by the Euler or Navier-Stokes equations where the volume occupied by the particles is neglected, and the dispersed phase (solid or liquid) governed by the pressureless gas dynamics equations [1,2].The coupling between those two sub-systems is achieved with relaxation terms expressing drag forces, heat and mass transfer effects.The dispersed phase system (pressureless gas dynamics equations) poses however serious difficulties.It is not possible to express the equations in characteristic variables as there is no set of independent eigenvectors.Shocks are unconventional as they don't respect the sub-characteristic criterion [3].These theoretical difficulties induce computational ones, as for example the treatment of reflective boundary conditions as well as any situation involving crossing parti-cles trajectories.Indeed, multivalued solutions are possible at a given point.These issues have been addressed by many authors [4][5][6] to cite a few.These references describe different approaches, sometimes based on slight modifications of the flow model, other times based on non-Eulerian formulations.We investigate in the present paper another flow model modification that consists in inserting turbulent effects in the particles cloud.It means that particles are evolving with the mean particle flow velocity and are subjected to velocity fluctuations, due for example to particles collisions.It results in the presence of a pressure term in the momentum equation that renders the particles system strictly hyperbolic.In the absence of particles collisions, here modeled with the help of a "turbulent viscosity", the turbulent pressure term vanishes and the pressureless gas dynamic equations are recovered.From a physical standpoint, turbulent pressure effects result in enhanced predictions of the corresponding flow model, compared to pressureless gas dynamic equations.Consider for example a free jet of particles evolving in a gas at rest.Modeling this jet with pressureless equations results in a straight jet, without any enlargement during its transport in the gas.When particles turbulence is present, dispersion occurs resulting in jet enlargement in perfect agreement with experimental measurements if proper turbulent viscosity is used.
The present paper is organized as follows.The conventional two-fluid model with pressureless equations for the particle phase is recalled in Section 2. Then turbulent effects are modeled in Section 3 with the thermodynamic method given in [7].The Riemann problem solution for the particles system is examined in Section 4, and those of the gas phase is considered in Section 5. Dissipative effects including relaxation ones and turbulent viscosity are inserted in Section 6.They result in turbulence production in both phases.Section 7 deals with some details about numerical implementation.Several test problems and comparison between the conventional pressureless model and the new turbulent model are considered in Section 8, first in one dimension, then in two dimensions, allowing comparisons with experimental data of free jet particle flows in air.Conclusions are given in Section 9.

Conventional Eulerian Flow Model for Dilute Suspensions
The following flow model considers each phase as a continuous media.In the carrier phase, molecular collisions are so intense that pressure appears as an external force to a given gas control volume.In the dispersed phase, the collisions are less intense and are here considered as negligible.This assumption is reasonable as the particle volume concentration is weak [8], say  represents the dispersed phase volume fraction.In addition, neglecting the volume occupied by the particles for the gas flow results in the following system [1,2].The evolution equations for carrier phase are the following ones, and for the dispersed phase, In the carrier phase system, there is a coupling between the internal energy ( ), the density ( e  ) and the pressure (P) with the help of a convex equation of state (EOS): ( , ) P P e   .
In the dispersed phase, such coupling is absent.The particles are incompressible and collisional effects are neglected.The only interaction force considered here is the drag force, modeled by the velocity relaxation parameter λ > 0. The power of the drag force is assumed to be transferred with the particles velocity, resulting in entropy production in the gas, and isentropic evolutions in the particles: It can be shown easily that the gas dynamic system is strictly hyperbolic with waves speed u, u + c and u -c.The sound speed c is defined by, where γ represents the specific heats ratio.The particles phase system is hyperbolic degenerated with a single characteristic speed p u .This poses both theoretical and numerical issues.The aim of the following section is precisely to correct the particles system in order to have better mathematical properties as well as improved physical meaning.To do this, turbulent effects are considered with the simplest thermodynamically consistent approach.

Turbulent Dilute Two-Phase Flow Model
Our aim is to insert turbulent effects in System (1-2).To do this, we follow the method described in [9] that proceeds in several steps.Let us consider first the carrier phase, the extension to the dispersed phase being straightforward. ) Two different average definitions are used: Favre average is used for any variable weighted by the density: In the case of an isotropic flow, we define the quantity K by: After some calculations, we obtain for the carrier phase: The gas total energy is now defined by: where n represents the number of degrees of freedom in which velocity fluctuations develop (usually, n = 3).To close the system, an evolution equation for the quantity K is needed.Combination of the energy, momentum and mass equations results in: The Gibbs identity reads: Imposing ds = 0 (for a flow evolving without shock waves and dissipation) and combining the two last relations, an evolution equation for K is obtained: Let us define the "turbulent" polytropic coefficient: We obtain, It is now clear that K represents the turbulent pressure denoted by t in the rest of the paper.We then define the "turbulent" entropy that obeys the following conservation law: ).
The flow model for the carrier phase thus reads in absence of interaction terms, ( ) with the following definitions and equations of state, Applying the same method to the pressureless dispersed phase Subsystem (2), the following "turbulent" particles flow model is obtained:  It is shown in the following sections, that both subsystems are hyperbolic and well posed for the Riemann problem resolution.It is clear from System (4) and thermodynamic definitions that when pt 0, the pressureless gas dynamics equations are recovered.s 

The Riemann Problem for the Dispersed Phase
The Riemann problem is considered first for the System (4) with the following notations: 0  W and r W (Figure 1).Let us first examine smooth solutions

Si le Waves Velocities
where .
The eigenvalues of this matrix correspond to th waves speeds.Unlike the pressureless model, 3 real and distinct eigenvalues appear.Consequently the system is st e rictly hyperbolic:

Characteristic Relations and Riemann Invariants
From the Riemann problem invariants and shock relations, it is clear that the dispersed phase Riemann problem and associated solvers are straightforward extension of for System (3) yielding ideal gas Riemann solvers for the Euler equations.Details may be found in [10].

The Riemann Problem for the Gas Phase
he same analysis is carried out T the following results: Waves speeds: (with Here direct integration of Riemann invari difficult and an approximate Riemann solver is preferred, rather than the exact one.An excellent candidate is the H tions ants is more LLC solver [10], initially developed for the Euler equabut having a straightforward extension to the pre-he entropy creation for the gas and the particles has e gas, the on results of drag effects while for the ispersed phase, it results of particles collisions, modeled t for th rder to introduce as w parameters as possible in the model, we use the simdrag: the Stokes drag.This force ads: sent turbulent model.

Dissipative Terms and Turbulence Creation
T different origins with the present modeling.For th entropy creati d hereafter with the help of a "turbulent viscosity".Also, it is important to note that the only dissipative effects considered in the following have a mechanical origin and not a thermal one.Thus the entropy production will be stored as turbulence creation only and no e thermodynamic entropy growth.This separation principle was successfully used in [7] in another two phase-flow context.

Viscous Drag
For the sake of simplicity and in o fe plest form of viscous re 6π ( ) where R is the particle radius (constant in the studied examples) and µ the gas viscosity.Let Re p denote the particle Reynolds number and coefficient such that [11]: Under these notations the drag force for a cloud droplet, reads per unit volume as: The resulting gas-phase system now reads: Here, a new parameter has been introduced for the turbulent temperature nition.However, it influence on turbulent pressure calculation.This ca demonstrated by calculating the turbulent pressure evolution equation from the turbulent entropy equation.

Turbulent Viscosity
As shown with System (1-2), the dissipation due to dr effects produces entropy in the gas phase only.In the following, we speculate that particles collisions are responsible for turbulence creation, resulting in particles cl isio a vt C defi has no n be ag oud dispersion.These coll ns are not intense enough compared to molecular collisions that occur in the gas phase.Thus a thermodynamic pressure is improper to model particles collisions.A transport coefficient, such as a "turbulent viscosity" is more adequate to model hese collisions.A turbulent stress tensor involving t single parameter ( t  ) is thus ubsystem (9).
inserted in the particles s  entropy production is stored as tu ce [7] as its origin is mechanical only.Combinin various equations of System (9), the turbulent entropy evolution equation reads: e creation term.The turbulent entropy (or turbulent energy) is obtained from the particles total energy equation where the production terms appear in divergence form.


It is worth to mention that numerical resolution of Equation ( 10) is not necessary.Indeed, difficulties are present to approximate th

Model Summary
The new two-phase flow model is summarized hereafter: with the following definitions and closure relations:

Numerical Resolution
Operator splitting is used to solve the various hyperbolic, diffusion and relaxation terms.

Hyperbolic Step
System (11), in absence of relaxation and diffusion corresponds to a system of conservation laws, split in two subsystems. ),( A Godunov type method is used to solve each subsysm with MUSCL type higher order extension (see again vector of computational cells faces.The HLLC solver [12] is used for both subsystems. Let us recall that with formulation, th rticles total energy equation is used to compute the turbulent particles energy.Their internal energy is determined from the last equation of System (11).

Drag Relation
At the end of each hyperbolic step, relaxation terms are in te [10]).The numerical fluxes are obtained with the Riemann problem solution for each subsystem, solved along the normal this e pa tegrated by considering the following ODE systems:

   
A high-order Runge-Kutta method is used in this aim.

Diffusive Terms
The particle phase subsystem contains diffusive contributions:

Validation and Comparison with the Conventional Model
W levant tests cases.In many situations, the conventional reless gas dynamics equations is nsidered quite accurate and the predictions with the that are approximated by an explicit finite volume scheme.Details are given in Appendix A. e now examine the new turbulent model behavior on re model (1-2) using pressu co we first consider one-dimensional situations of two-phase flows in shock tubes.The conventional model (pressureless equations) is solved with the alg thm de- n with a Particle Cloud essureless e turbu se new model must be quite closed.To compare typical solutions, ori scribed in [4].Obviously, the same modeling of drag effects is used in both models.We consider the interaction of a shock wave with a particle cloud initially at rest.The particles apparent density is ).The left part of the shock tube corresponds to the high pressure chamber, initially at P = 2 atm and with the initial gas density  ) on both sides of the cloud implies a fast acceleration of the particles outside the cloud (see the particles velocity graph of Figure 4), but the small particles proportion has no influence on the gas motion.
The particle velocity profiles of Figures 3 and 4 have me so differences, but these diff sid erences are located out-comparison, we sh ined hereafter the influen e the cloud, where the particles are absent.So, these differences have no importance.For the same reason turbulent pressure profiles must be observed within the cloud.
Figures 3 and 4 show the same behavior for gas ve-locity, pressure and density.The apparent particles density profiles also are closed.For a better ow the particles variables in the particles cloud only and on the same graphs in the Figure 5.
Excellent agreement is observed between the two computations, showing that in one dimension, the two models are equivalent.We exam ce of the turbulent viscosity parameter.

Influence of Turbulent Viscosity
We now consider the turbulent model only and use different values of the turbulent viscosity: For this 1D test problem, it appears clearly that the turbulent viscosity has no influence on the cloud dynamics.In all these computations, the initial turbulent pressure was set to 10 Pa.We now investigate the influence of this initial condition to the results.Some influence appears, in particular a dissymmetry in the cloud, due to the fact that the turbulent pressure pr g res ncouraging as the in figuration oduction is different on the left and the right sides of the cloud.These differences are not significant enough to detect unphysical behavior.
In the one-dimensional case, whatever the turbulent parameters are, the two models (pressureless and turbunt) yield very closed results.This is e le pressureless model is considered as quite accurate such flow conditions.We now examine two-dimensional flow con s.

Two-Dimensional Configurations
We now consider a more sophisticated situation where extra effects appear.It consists in the injection of a parti-  cle jet in air at rest.In order to exclude the effects of dynamics fragmentation with liquid jets, that contain extra physics, we consider solid particles jets.Careful experimental studies were done independently [13,14] with particles made of different materials.Corresponding experimental data will be used for model's validation.Both experimental configurations can be summarized as shown in the Figure 8.A particle injector of a few millimeters diameter is used to inject a free particle jet into an open cavity initially filled with air at atmospheric conditions.
In [13,14] experiments, cross cut of apparent density profiles are given at various locations from the injecto (C1, C2, C3).The configurations used by these authors r

are reported in the Table 1.
To deal with numerical simulations of corresponding test problems, initial and boundary conditions as well as flow model parameters have to be specified.

Boundary Conditions
The computational domain contains 6 different bounda-ries shown in Figure 9.
All the computations are done in 2D axisymmetric configuration.The boundary condition 1 is treated by symmetry condition.Boundaries 2 and 3 are treated with non reflective conditions.Boundary 5 is treated with wall condition.Boundaries 4 and 6 correspond respectively to gas and particles inflows.In particular, Boundary 6 cor - Table 1.Experimental data for solid particle jet injection in air from [13,14].[14] [13] Particle diameter ( m) 500 64 Material density (kg/m 3 ) 1020 2590 Injection speed (m/s) 24 30 Injection apparent density (kg/m 3 ) 0.2 1.0 Injector diameter D (mm) 20 13 responds to a supersonic particles injection with respect to the turbulent speed of sound.Thus all data of the Ta- ble 1 are imposed at this boundary.Boundary 4 is treated as a tank for both gas and particle phases.Numerica ystem (11),

Gas Phase
The ste ons redu : l treatment of this specific boundary needs particular care.Part of the solution is obtained by assuming a steady flow between the tank and the inlet section.S in the stationary case is integrated in a control volume which is represented to the left of Figure 10.

ady gas phase equati ce to
Using the slip condition on the lateral surfaces ( = 0), the integration of the mass equation gives (variables with subscript i, corresponds to variables at the inlet surface, those with subscript 0 corresponds to tank variables): 0 This relation is undetermined, indeed and .Nevertheless it can be used turbulent entropy equations.Fo en 0 S   to integrate the enr the total 0 0 u  ergy and ergy, the integration gives, with the help of (12),   A first relation is thus available to solve the boundary Riemann problem.The momentum conservation equation cannot be integrated as the pressure integral is unknown on the lateral surface.Integration of turbulent entropy equation using (12) results in: It is the second available equation is required to close the full system.This one corresponds to thermodynamic tween the tank and the inlet section, e entropy conservation bexpressed as: Then a Riemann solver is built where the left facing wave jump conditions (Riemann invariants or shock relations) are replaced by System (13)(14)(15).The right facing wave obeys t l as the contact discontinui

Particle Phase
In the same way, the system associated to the particle phase is integrated in the control volume.The stationary system reads:   A Riemann solver is built where the left facing wave relation, but an additional jump conditions are replaced by System (17-19).

Initial Conditions
The gas is initially at rest under atmospheric conditions with a low turbulent pressure = 10 Pa).The particles have an initial non-zero vol fraction in the entire do ( t P ume main, corresponding to the apparent density of

Model Parameters
The same drag force correlation ( viscosity is used as in the preceding 1D example kg/m/s).In all compu viscosity is constant and equal to 7) with the same gas Large differences between the two computed jets are clearly visible.With the pressureless conventional model, no jet enlargement occurs.The particles traj respond s

Comparison with the Experimental Data of [13]
A computational domain involving 280 × 100 cells is used on the geometrical configuration sh e 8 with parameters given in the Table 1.The normalized apparent density profiles are shown in the Figure 12 and compared with experimental data at steady state.

Comparison
T reported in [14].The normalized apparent density profiles at steady sate are shown in the Figure 13 and compared to experimental.
In this case too, the results obtained with the turbulent model and the experimental data are in a very good e particle jet dynamics computed alternatively with the conventional pressureless model and the new turbu-  The results of Figures 13 and 14 are obtained with the same turbulent viscosity of kg/m/s.This agreement is particularly interesting as t e experiments of [13] and [14] deal with very different materials density, different particles diameter and injector diameters.
When the same computations are done with the pres-sureless model (1-2), no jet enlargement appears, as shown in the Figure 14.The same remark holds when the turbulent viscosity is set to zero in the new model.Consequently the new model admits the same solution as the pressureless equations in the limit of vanishing turbulent viscosity and is able to improve considerably the predictions and capabilities with non-zero turbulent viscosity.

Conclusions
A new two-phase flow model for dilute particles suspensions has been developed.Compared to conventional pressureless particle dynamic models the new one is n s cessfully reproduced by the model with the same turbulent viscosity coefficient, of the order of The cell boundary velocity derivative reads consequently, In the one-dimensional case, viscous derivative terms are thus approximated by:

Figure 1 .
Figure 1.Wave diagram for the dispersed phase Riemann problem.
across left-and right-facing compression and expansion ves, but t across discontinuities.er a discontinuity propagating at velocity σ.ations of System (4) read:


are assumed incompressible, their internal energy depends only of the temperature (or entropy).As the thermodynamic entropy is constant along trajectories in System (4last equation of Syste(9), that can alternatively be written as, In the present formulation, it has been assumed that all rbulen g the S. HANK ET AL.

.
For numerical reasons, we set on both sides of the cloud a weak but non zero particles apparent density corresponds to the low pressure chamber with the initial pressure P = 1 atm and with the gas density 3 1 kg/m   .The initial situation is depicted in the Figure 2.  Reference results obtained with the pr model (1-2) A shock wave propagates to the right in the gas phase while a rarefaction wave propagates to the left.The shock then interacts with the particle cloud and a diffraction phenomenon occurs.A shock is transmitted to the cloud and a weak reflected shock is emitted to the left from the cloud interface.The particle cloud compresses and moves to the right.The various flow variables profiles are shown in the Figure 3 at time 2.22 ms.We now examine the turbulent model numerical solution in the same configuration. Turbulent model results With th lent model the initial turbulent pressures have to be set.We set for both phases a weak initial pressure of 10 Pa.The turbulent viscosity is arbitrarily 3 2.10 kg/m/s.In the one-dimensional case, the value of this parameter has no influence as will be shown later.The flow variables profiles are shown at the same time (t = 2.22 ms) in the Figure4.The weak particles volume fraction ( particles density and the particles velocity profiles are plot on the same graph for each value of the turbulent viscosity.Corresponding results are shown in the Figure 6.

Figure 7
at time t = 8.6 ms.

Figure 2 .
Figure 2. Shock tube initial situation with a pressure ratio of 2 and a particle cloud located between abscissa x = 1 m and x = 1.2 m.

Figure 3 .
Figure 3. Gas and particles variables profiles at times t = 0 and t = 2.22 ms for the two phase shock tube test case with the conventional pressureless model.

Figure 4 .
Figure 4. Gas and particles variables profiles at times t = 0 and t = 2.22 ms for the two phase shock tube test case with the new turbulent model.

Figure 5 .
Figure 5.Comparison of pressureless and turbulent models solutions on the shock tube test problem.Apparent density and velocity profiles are in perfect agreement.

Figure 6 .
Figure 6.Influence of the turbulent viscosity on the shock tube test case results.The apparent particles density as well as the particles velocity in the particle cloud shows no influence.The results are compared at time t = 8.6 ms.

Figure 7 .
Figure 7. Influence of the turbulent initial pressure in the dispersed phase.The particles apparent density and velocity profiles in the particles cloud are shown at time t = 8.6 ms for three initial pressures.

Figure 8 .
Figure 8. Schematic representation of the experimental facility of [13,14] for solid particles jets injection in air.

Figure 9 .Figure 10 .
Figure 9. Boundary conditions. with to straight lines.With the turbulent model, jet preading occurs.own in the Fig- ur with the Experimental Data of [14] he same computations are done for the configuration lent one are shown in the Figure 11.

Figure 11 .
Figure 11.Apparent density contours for the test problem of [13] (see Table 1 for corresponding data).The pressureless model results are shown on top while the turbulent one are shown on bottom of the figure.Particles jet behavior is shown at times t1 = 2.8 ms, t2 = 9.6 ms, t3 = 18.2 ms and t4 = 30 ms.This last time corresponds to steady state.

Figure 12 .Figure 14 .
Figure 12.Computed results in lines versus experimental data of [13].Normalized apparent density 0 /   p p 30 D. to collisi y the mod noted by

3 2. 10
 h strictly hyperbolic.It also has enhanced physical capabilities, thanks to a simple modeling of particles collisions.Indeed, turbulence production in the new model is linked to a turbulent viscosity, aimed to mimic particles collisions.Particles jets dispersion experiments from different authors and under different configurations have bee uc s A. Numerical method for diffusive term the left and the right sides immediately closed to the boundary.Thus, Time explicit integration of corresponding terms is considered: cells i and i + 1, separated by cell boundary i + 1/2 are schematized in the Figure 15.The velocity derivative * 1/ 2 i q has precisely to be expressed at these cell boundaries.followingintercell conditions are used:

Figure 15 .
Figure 15.Schematic representation of cell boundary i + 1/2 separating left and right computational cells i and i + 1. )