Journal of Modern Physics
Vol. 2  No. 9 (2011) , Article ID: 7141 , 15 pages DOI:10.4236/jmp.2011.29120

A Hyperbolic Eulerian Model for Dilute Two-Phase Suspensions

Sarah Hank, Richard Saurel, Olivier Le Metayer

IUSTI, Aix Marseille University, Marseille Cedex, France Bastidon de la Caou, Roquevaire, France

E-mail: {Sarah.Hank, Richard.Saurel, Olivier.Lemetayer}@polytech.univ-mrs.fr

Received May 4, 2011; revised June 22, 2011; accepted July 1, 2011

Keywords: Turbulence, Riemann Solver, Pressureless Gas Dynamics

ABSTRACT

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

1. 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 NavierStokes 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 particles trajectories. Indeed, multivalued solutions are possible at a given point. These issues have been addressed by many authors [4-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.

2. 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

where 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,

(1)

and for the dispersed phase,

(2)

represents the particles apparent density, , where represents the condensed phase real density (considered here as incompressible). The total energies are defined by:

and.

In the carrier phase system, there is a coupling between the internal energy (), the density () and the pressure (P) with the help of a convex equation of state (EOS):

.

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. 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.

3. 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:

•      Reynolds average:

•      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:

Multiplying the evolution equation for K by,

.

We obtain,

.

It is now clear that K represents the turbulent pressure denoted by in the rest of the paper. We then define the “turbulent” entropy that obeys the following conservation law:

(with).

The flow model for the carrier phase thus reads in absence of interaction terms,

(3)

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:

(4)

, ,

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 0, the pressureless gas dynamics equations are recovered.

4. The Riemann Problem for the Dispersed Phase

The Riemann problem is considered first for the System (4) with the following notations:

,.

Figure 1. Wave diagram for the dispersed phase Riemann problem.

The Riemann problem resolution consists in determining states and as well as the various waves speeds starting from an initial discontinuity separating and (Figure 1). Let us first examine smooth solutions.

4.1. Simple Waves Velocities

System (4) can be written in quasi-linear form as:

with:

.

where represents the particles “turbulent sound speed”, such that:

.

The eigenvalues of this matrix correspond to the waves speeds. Unlike the pressureless model, 3 real and distinct eigenvalues appear. Consequently the system is strictly hyperbolic:

, ,.

4.2. Characteristic Relations and Riemann Invariants

With the help of left eigenvectors associated to each eigenvalue, the following characteristic relations are obtained:

.

These relations can be integrated to obtain:

,.

These relations apply across leftand right-facing compression and expansion waves, but not across discontinuities.

4.3. Discontinuities

Let us consider a discontinuity propagating at velocity σ. The Rankine-Hugoniot relations of System (4) read:

(5)

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

5. The Riemann Problem for the Gas Phase

The same analysis is carried out for System (3) yielding the following results:

•      Waves speeds:

, ,. (with, ,)

•      Riemann invariants

;.

•      Shock relations

(6)

Here direct integration of Riemann invariants is more difficult and an approximate Riemann solver is preferred, rather than the exact one. An excellent candidate is the HLLC solver [10], initially developed for the Euler equations but having a straightforward extension to the present turbulent model.

6. Dissipative Terms and Turbulence  Creation

The entropy creation for the gas and the particles has different origins with the present modeling. For the gas, the entropy creation results of drag effects while for the dispersed phase, it results of particles collisions, modeled 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 not for the thermodynamic entropy growth.

This separation principle was successfully used in [7] in another two phase-flow context.

6.1. Viscous Drag

For the sake of simplicity and in order to introduce as few parameters as possible in the model, we use the simplest form of viscous drag: the Stokes drag. This force reads:

.

where R is the particle radius (constant in the studied examples) and µ the gas viscosity. Let denote the particle Reynolds number and a drag coefficient such that [11]:

if, otherwise.

Under these notations the drag force for a cloud of droplet, reads per unit volume as:

(7)

The resulting gas-phase system now reads:

(8)

,

.

Here, a new parameter has been introduced for the turbulent temperature definition. However, it has no influence on turbulent pressure calculation. This can be demonstrated by calculating the turbulent pressure evolution equation from the turbulent entropy equation.

6.2. Turbulent Viscosity

As shown with System (1-2), the dissipation due to drag effects produces entropy in the gas phase only. In the following, we speculate that particles collisions are responsible for turbulence creation, resulting in particles cloud dispersion. These collisions 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 these collisions. A turbulent stress tensor involving a single parameter () is thus inserted in the particles subsystem (9).

(9)

with.

As the particles are assumed incompressible, their internal energy depends only of the temperature (or entropy). As the thermodynamic entropy is constant along trajectories in System (4) () its analogue in terms of internal energy, as the particles are incompressible, reads

resulting in the last equation of System (9), that can alternatively be written as,

In the present formulation, it has been assumed that all entropy production is stored as turbulence [7] as its origin is mechanical only. Combining the various equations of System (9), the turbulent entropy evolution equation reads:

(10)

It is worth to mention that numerical resolution of Equation (10) is not necessary. Indeed, difficulties are present to approximate the creation term. The turbulent entropy (or turbulent energy) is obtained from the particles total energy equation where the production terms appear in divergence form.

6.3. Model Summary

The new two-phase flow model is summarized hereafter:

(11)

with the following definitions and closure relations:

, , , , , , ,

7. Numerical Resolution

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

7.1. 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 subsystem with MUSCL type higher order extension (see again [10]). The numerical fluxes are obtained with the Riemann problem solution for each subsystem, solved along the normal vector of computational cells faces. The HLLC solver [12] is used for both subsystems.

Let us recall that with this formulation, the particles total energy equation is used to compute the turbulent particles energy. Their internal energy is determined from the last equation of System (11).

7.2. Drag Relation

At the end of each hyperbolic step, relaxation terms are integrated by considering the following ODE systems:

,

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

7.3. Diffusive Terms

The particle phase subsystem contains diffusive contributions:

, that are approximated by an explicit finite volume scheme. Details are given in Appendix A.

8. Validation and Comparison with the  Conventional Model

We now examine the new turbulent model behavior on relevant tests cases. In many situations, the conventional model (1-2) using pressureless gas dynamics equations is considered quite accurate and the predictions with the new model must be quite closed. To compare typical solutions, we first consider one-dimensional situations of two-phase flows in shock tubes. The conventional model (pressureless equations) is solved with the algorithm described in [4]. Obviously, the same modeling of drag effects is used in both models.

8.1. One-Dimensional Configurations

8.1.1. Shock Wave Interaction with a Particle Cloud

We consider the interaction of a shock wave with a particle cloud initially at rest. The particles apparent density is. For numerical reasons, we set on both sides of the cloud a weak but non zero particles apparent density (). The left part of the shock tube corresponds to the high pressure chamber, initially at P = 2 atm and with the initial gas density. The right part corresponds to the low pressure chamber with the initial pressure P = 1 atm and with the gas density. The initial situation is depicted in the Figure 2.

Reference results obtained with the pressureless 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 the turbulent 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 set to 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 Figure 4.

The weak particles volume fraction () 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 some differences, but these differences are located outside 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 velocity, pressure and density. The apparent particles density profiles also are closed. For a better comparison, we show 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 examined hereafter the influence of the turbulent viscosity parameter.

8.1.2. Influence of Turbulent Viscosity

We now consider the turbulent model only and use different values of the turbulent viscosity: kg/m/s, and. The apparent 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.

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.

8.1.3. Influence of the Initial Turbulent Pressure

We now consider three different initial turbulent pressures (, ,) for the same two-phase shock tube test. The turbulent viscosity is set to. Corresponding results are shown in the Figure 7 at time t = 8.6 ms.

Some influence appears, in particular a dissymmetry in the cloud, due to the fact that the turbulent pressure production 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 turbulent) yield very closed results. This is encouraging as the pressureless model is considered as quite accurate in such flow conditions. We now examine two-dimensional flow configurations.

8.2. Two-Dimensional Configurations

We now consider a more sophisticated situation where extra effects appear. It consists in the injection of a parti-

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. 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.

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 injector (C1, C2, C3). The configurations used by these authors

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.

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.

8.2.1. Boundary Conditions

The computational domain contains 6 different boundaries 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-

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

Table 1. Experimental data for solid particle jet injection in air from [13,14].

responds to a supersonic particles injection with respect to the turbulent speed of sound. Thus all data of the Table 1 are imposed at this boundary. Boundary 4 is treated as a tank for both gas and particle phases. Numerical 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. System (11), in the stationary case is integrated in a control volume which is represented to the left of Figure 10.

8.2.2. Gas Phase

The steady gas phase equations reduce to:

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.

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):

(12)

This relation is undetermined, indeed and. Nevertheless it can be used to integrate the energy and turbulent entropy equations. For the total energy, the integration gives, with the help of (12),

, (13)

where, and.

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. Schematic representation of the experimental facility of [13,14] for solid particles jets injection in air.

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:

(14)

It is the second available relation, but an additional

Figure 9. Boundary conditions.

Figure 10. Schematic representation of the tank boundary condition.

equation is required to close the full system. This one corresponds to thermodynamic entropy conservation between the tank and the inlet section, expressed as:

(15)

Then a Riemann solver is built where the left facing wave jump conditions (Riemann invariants or shock relations) are replaced by System (13-15). The right facing wave obeys to conventional relations, as well as the contact discontinuity.

8.2.3. Particle Phase

In the same way, the system associated to the particle phase is integrated in the control volume. The stationary system reads:

Following the same methodology, the following relations are obtained:

(16)

(17)

(18)

(19)

with.

A Riemann solver is built where the left facing wave jump conditions are replaced by System (17-19).

8.2.4. 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 volume fraction in the entire domain, corresponding to the apparent density of , associated to a volume fraction of the order of. Initial turbulent pressure of the particles is set to the same value as in gas phase (= 10 Pa).

8.2.5. Model Parameters

The same drag force correlation (7) with the same gas viscosity is used as in the preceding 1D example (kg/m/s). In all computations, the turbulent viscosity is constant and equal to kg/m/s.

8.2.6. Qualitative Comparisons

The particle jet dynamics computed alternatively with the conventional pressureless model and the new turbulent one are shown in the Figure 11.

Large differences between the two computed jets are clearly visible. With the pressureless conventional model, no jet enlargement occurs. The particles trajectories correspond to straight lines. With the turbulent model, jet spreading occurs.

 

8.2.7. Comparison with the Experimental Data of [13]

A computational domain involving 280 × 100 cells is used on the geometrical configuration shown in the Figure 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.

8.2.8. Comparison with the Experimental Data of [14]

The same computations are done for the configuration 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

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. Computed results in lines versus experimental data of [13]. Normalized apparent density cross cut at abscissa C1 = 10 D, C2 = 20 D and C3 = 30 D. The jet width increases during its propagation, due to collisional turbulent effects that are perfectly reproduced by the model. The apparent density along the centerline is denoted by.

agreement. For the last graph of Figure 13 the agreement is not as good as for the other graphs. But it can be noticed on the experimental curve that some noise is present. The experimental cross cut of the apparent particle density contains oscillations. It means that an error bar is

Figure 13. Computed results in lines versus experimental data of [14]. Normalized apparent density cross cut at abscissa C1 = 4 D, C2 = 10 D and C3 = 15 D. The jet width increases during its propagation, due to collisional turbulent effects that are perfectly reproduced by the model. The apparent density along the centerline is denoted by.

present in the experiments. Also, the variable plot on this graph corresponds to the normalized density. At the normalized distance of 15D, both and

Figure 14. Results comparison with the different models on the configuration studied in [14].

become small. Thus, the ratio becomes inaccurate.

The results of Figures 13 and 14 are obtained with the same turbulent viscosity of kg/m/s. This agreement is particularly interesting as the 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 pressureless 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.

9. 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 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 been successfully reproduced by the model with the same turbulent viscosity coefficient, of the order of kg/m/s, while the experimental conditions (particles diameter, particles density and injector diameter) are varying considerably.

REFERENCES

  1. F. E. Marble, “Dynamics of Dusty Gases,” Annual Review of Fluid Mechanics, Vol. 2, No. 1, 1970, pp. 397-446. doi:10.1146/annurev.fl.02.010170.002145
  2. Y. B. Zeldovich, “Gravitational Instability: An Approximate Theory for Large Density Perturbations,” Astron & Astrophys, Vol. 5, No. 84, 1970, pp. 168.
  3. P. D. Lax, “Weak Solutions of Nonlinear Hyperbolic Equations and Their Numerical Computation,” Communications on Pure and Applied Mathematics, Vol. 7, No. 1, 1954, pp. 159-193.
  4. R. Saurel, E. Daniel and C. Loraud, “Two-Phase Flows: Second-Order Schemes and Boundary Conditions,” AIAA Journal, Vol. 32, No. 6, 1994, pp. 1214-1221. doi:10.2514/3.12122
  5. Y. Brenier and E Grenier, “Sticky Particles and Scalar Conservation Laws,” SIAM Journal on Numerical Analysis, Vol. 35, No. 6, 1998, pp. 2317-2328. doi:10.1137/S0036142997317353
  6. A. Chertock, A. Kurganov and Y. Rykov, “A New Sticky Particle Method for Pressureless Gas Dynamics,” SIAM Journal on Numerical Analysis, Vol. 45, No. 6, 2007, pp. 2408-2441. doi:10.1137/050644124
  7. R. Saurel, S. Gavrilyuk and F. Renaud, “A Multiphase Model with Internal Degrees of Freedom: Application to Shock-Bubble Interaction,” Journal of Fluid Mechanics, Vol. 495, No. 1, 2003, pp. 283-321. doi:10.1017/S002211200300630X
  8. G. Rudinger, “Some Effects of Finite Particle Volume on the Dynamics of Gas-Particle Mixtures (Gas Particle Mixture with Finite Particle Volume Affecting Frozen and Equilibrium Flows Behind Shock Wave),” AIAA Journal, Vol. 3, 1965, pp. 1217-1222.
  9. R. Saurel, A. Chinnayya and F. Renaud, “Thermodynamic Analysis and Numerical Resolution of a Turbulent-Fully Ionized Plasma Flow Model,” Shock Waves, Vol. 13, No. 4, 2004, pp. 283-297. doi:10.1007/s00193-003-0216-z
  10. E. F. Toro, “Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction,” Springer Verlag, Berlin, 1997.
  11. L. Schiller and Z. Naumann, “A Drag Coefficient Correlation,” VDI Zeitung, Vol. 77, 1935, pp. 318-320.
  12. E. F. Toro, M. Spruce and W. Speares, “Restoration of the Contact Surface in the HLL-Riemann Solver,” Shock Waves, Vol. 4, No. 1, 1994, pp. 25-34. doi:10.1007/BF01414629
  13. K. Hishida, K. Kaneko and M. Maeda, “Turbulence Structure of a Gas-Solid Two-Phase Circular Jet,” Trans. JSME, Vol. 51, 1985, pp. 2330-2337. doi:10.1299/kikaib.51.2330
  14. Y. Tsuji, Y. Morikawa, T. Tanaka, K. Karimine and S. Nishida, “Measurement of an Axisymmetric Jet Laden with Coarse Particles,” International Journal of Multiphase Flow, Vol. 14, No. 5, 1988, pp. 565-574. doi:10.1016/0301-9322(88)90058-4

A. Numerical method for diffusive terms

Time explicit integration of corresponding terms is considered:

Space integration is detailed in 1D, multi-D extension being similar. The space integration is transformed successively as,

, where.

Two neighboring computational cells i and i + 1, separated by cell boundary i + 1/2 are schematized in the Figure 15. The velocity derivative has precisely to be expressed at these cell boundaries.

The following intercell conditions are used:

where the superscripts (–) and (+) denotes respectively the left and the right sides immediately closed to the boundary. Thus,

,

The cell boundary velocity derivative reads consequently,

,

In the one-dimensional case, viscous derivative terms are thus approximated by:

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