Ind. Eng. Chem. Res. 1999, 38, 1177-1182
1177
Unsteady Compressible Flow of Granular Materials† Tommaso Astarita Department of Energetics, Thermofluodynamics and Environmental Control, Universita` di Napoli “Federico II”, Naples, Italy
Raffaella Ocone*,‡ School of Chemical, Environmental and Mining Engineering, University of Nottingham, Nottingham, U.K.
Elementary compressible flow problems of granular materials have been analyzed in the past few years (Ocone, R.; Astarita, G. J. Rheol. 1993, 37, 727-742), and even only moderately complex problems present subtle difficulties, because of inelasticity of particle-particle collisions. Some analogies with classical gas dynamics exist, and those make it possible to approach basic problems using tools used there. Some results can be obtained from such an approach: standing compression and rarefaction waves, thermal relaxation tails, and the like. Introduction The study of pressure waves in granular materials (GM) has recently attracted much attention because of its connection to a large amount of practical problems. It is well-known, for instance, that wave propagation in fluidized beds is related to their (or lack thereof) stability.2,3 Most of the theories, which can be found in the literature, are based on an empirical or semiempirical approach, being the complete equations of motion very difficult to solve. Even in systems where the interstitial fluid does not play a major role, the mathematical complete problem is extremely difficult to solve, both analytically and numerically. In this work, following an analogous route developed in gas dynamics (GD), and using thermodynamic tools developed by Ocone and Astarita,1 we approach the problem of pressure waves stability. The results can furnish important insights in handling practical problems arising in silos, discharge from silos, chute flow, hoppers, and other industrial solid handling units. As gases, granular materials are intrinsically compressible. Although the intrinsic particle density is constant, the solid volume fraction, hence the local averaged density of the particulate phase, may well change in space and time. One of the fundamental tools in classical gas dynamics is the thermodynamic theory, which makes it possible to define the concept of entropy and from this the speed of sound. Even if compressible materials always undergo compressible flow, in solving practical problems, one assumes incompressible flow; this is due to the fact that if the flow is only slightly compressible, the assumption of compressible flow can produce numerical errors in the solution procedure.4 In classical GD the difference between an incompressible and a compressible flow may be easily defined and quantitatively verified, introducing the concept of steady and isentropic flow. In this way, it is possible to reduce the independent thermodynamic variables to only one; indeed, a unique relation between † This paper was submitted for the Roy Jackson special issue. ‡ Present address: Department of Mechanical and Chemical Engineering at Heriot-Watt University, Edinburgh EH14 4AS, U.K. Phone: (0044)-131-451 3777. Fax: (0044)-131-451 3129.
density and temperature and density and pressure, and between pressure and temperature exists. Furthermore, after having defined the Mach number, M, as the ratio of the local velocity and the speed of sound in the same point, it is possible to consider the ratio of the static and stagnation values of each thermodynamic variable which is a function only of the Mach number. The concept of the stagnation condition is very useful in compressible flow: if we adiabatically and isentropically decrease the velocity, the thermodynamic variables will change in values until they reach a new value which corresponds to zero velocity. Such a value is called the stagnation value. The ratio between the static and stagnation density furnishes an insight of “how” compressible the flow is. Since the density ratio is a monotonically decreasing function of the Mach number, by fixing the accuracy, it is possible to find M such that the flow can be assumed incompressible. In granular flow the situation is more complex; even if a theory of granular thermodynamics1 has been recently developed, because of the inelasticity of particleparticle and/or particle-wall collisions, there is always a dissipation of internal pseudoenergy with a consequent change of the entropy level (even for a situation of rest). A way to overcome this problem is to assume that there is a positive influx of energy from the walls, considering, for instance, vibrating walls. Since we mostly wish to discuss the logical status of the unsteady compressible granular flow, we use what is probably the simplest possible example of the problem: a basic flow that it is maintained in a thermalized state by vibration of a constant section conduit. Furthermore, we have chosen to use, for the constitutive equations for granular flow rheology, the simplest ones that still make sense at least in the limit of a rather dilute system. An in-depth analysis of the assumptions made in this paper may be found in Astarita et al.5 The simplest problem in classical GD, viz. the velocity of sound (i.e. the celerity of propagation of an infinitesimal discontinuity of pressure or density) can be extended easily to granular flow.6 An extension of the classical Hugoniot7,8 analysis of the propagation of finite discontinuities is also possible.9 In such an analysis, the discontinuities are assumed to exist; however, their logical status is unclear, since contrary to the case of
10.1021/ie980393z CCC: $18.00 © 1999 American Chemical Society Published on Web 02/25/1999
1178 Ind. Eng. Chem. Res., Vol. 38, No. 4, 1999
classical gas dynamics, in granular flow theory the Hugoniot equations do not enjoy the status of a proper asymptote of the finite standing-wave case.10 The complete problem (i.e., the stability of such discontinuity) may not be solved in an exact way and it is necessary to use the full unsteady-state balances of mass, momentum, and fluctuating energy. Governing Equations The constitutive equations for granular flow rheology that we will use in this work have been extensively described by Astarita et al.5,11 In analogy with GD, we neglect the effects of viscosity and thermal conductivity so that there are only three constitutive relations that need to be specified:
U)T
(1)
p ) νFpT
(2)
I ) pxT/D
(3)
D ) dp/(1 - e)
(4)
with
where T is the granular pseudotemperature, U the internal energy, Fp the particle density, ν the solid volume fraction, p the particulate phase pressure, D an intrinsic length scale, I the dissipation rate of pseudoenergy, e the coefficient of restitution for particleparticle collisions, and dp the particle diameter. Since the particle density is constant, it is possible to write eq 1 as
p ) FT
(5)
with F ) Fpν the particulate phase density. Equations 1 and 5 are a particular case of constitutive equations for a perfect gas. Indeed, if one assumes that the gas constant is equal to 1, and that the specific heat at constant volume is also 1, the equations of state of a perfect gas reduce to eqs (1) and (5). Under these hypotheses, the ratio of the two specific heats, γ, comes out to be 2. However, the third constitutive equation has not any counterpart in GD (in this case the coefficient of restitution is exactly equal to 1), and introduces an intrinsic scale length, which while significantly larger than dp (e is likely to be close to unity), is nevertheless likely to be small with respect to the duct transverse dimension. A fundamental parameter in the analysis of compressible flow is the Mach number. Its definition is based upon the concept of speed of sound, namely the rate of propagation of small disturbances. Despite that these concepts are quite trivial in GD, the extension to GM poses interesting conceptual problems. Indeed, the concept of “granular” entropy, an essential concept in the definition of speed of sound, has been introduced only a few years ago by Ocone and Astarita.1 In a successive paper Ocone and Astarita6 found a relation for the particulate speed of sound, a, which under our hypotheses reduces to
a ) x(2T)
(6)
The Mach number, defined as M ) v/a (with v the average particulate velocity), is a fundamental criterion
Figure 1. M2 versus M1 for classical gas dynamics (dotted line) and granular flow (solid line).
for the flow definition. A flow for which M < 1 is called subsonic: every disturbance can travel upstream of the main flow and it may influence every point of the flow field. Indeed, it is a well-known result in GD that the fluid velocity in a tube is a function of both the upstream and the downstream pressures (the equations are elliptic in the space coordinates). The opposite situation occurs in a supersonic flow (M > 1): the outlet pressure does not play any role in determining the fluid velocity distribution. The “reversible” unsteady-state mass, momentum, and energy balances for a one-dimensional flow in the x direction are in conservative form, with v the velocity, t time, and E the total energy U + v2/2 (by reversible we mean that thermal conductivity and bulk viscosity have been neglected; in GD, but not in GM, this implies that the flow, in the absence of shock waves, is actually reversible in a thermodynamic sense):
Ft + (Fv)x ) 0
(7)
(Fv)t + (Fv2 + p)x ) 0
(8)
Et + (FvE + vp)x ) (Q - I)
(9)
where the fundamental difference with GD is the presence of the lower-order term (Q - I) which makes the system (7)-(9) inhomogeneous.12 Q represents the rate of energy input per unit volume. Equations (7)-(9) are coupled with the constitutive equations (1), (3), and (5). In GD, the (homogeneous) system corresponding to eqs (7)-(9) is such that the energy balance reduces to the result that, in the absence of shocks, if a system starts from an isentropic state, it stays so at all times (the substantial derivative of entropy, DS/Dt, is zero). Hence, the system of partial differential equations (PDE) reduces to a 2 × 2 one for which the Riemann invariants can be found; from that point on, the problem reduces to the integration of ordinary differential equations (ODE) along the easily determined characteristic lines.13 In the GM theory one does not obtain the result DS/Dt ) 0 because of the presence of the (Q - I) term in the energy balance: hence, the classical shortcut of GD is not available, since the complete 3 × 3 system of PDEs does not admit Riemann invariants.10
Ind. Eng. Chem. Res., Vol. 38, No. 4, 1999 1179
Figure 2. Weak shock wave.
Figure 3. Strong shock wave.
Analysis The steady-state form of the equations in GD implies that the mass flow rate, the impulse, and the total enthalpy are constant; these quantities are also equal on both sides of a shock in the Rankine-Hugoniot classical analysis of shock waves:
Fv ) c1
(10)
p + Fv2 ) c2
(11)
U + p/F + v2/2 ) c3
(12)
From a physical point of view eq (10) is trivial; eq (11) states that, in the absence of dissipative effects, momentum is conserved; eq (12) states the conservation of total enthalpy. The latter is similar to the Bernoulli equation for simple fluids.
In the analysis of shock waves in GM theory, one cannot impose that all three quantities are equal on the two sides of the shock, as well as that Q ) I on both sides: the problem appears to be overdetermined. In practice, Ocone and Astarita9 let go of the condition Q ) I behind the shock, and they obtained modified Rankine-Hugoniot shock equations which, however, do not satisfy the condition Q ) I behind the shock. Hence, the Ocone and Astarita9 analysis yields the short (in the D scale) range solution of the standing wave problem, with a thermal relaxation expected to take place behind the shock.14 The correct way to approach the analysis of GM standing waves is to let go of the condition that U + p/F + v2/2 has the same value on both sides of the wave, rather than of the condition that Q ) I on both sides. In fact, both ahead and behind any standing wave, the total enthalpy has to be constant along the flow direction, and Q must balance I. Neither one of these
1180 Ind. Eng. Chem. Res., Vol. 38, No. 4, 1999
Figure 4. Mach number versus y for an expansion wave as a function of time.
conditions holds within the wave, andswhich is far more important-the value of U + p/F + v2/2 behind the wave is certainly not going to be equal to the one ahead of the wave (if Q ) I has to hold behind the wave). In classical GD, eqs (10)-(12) determine a unique relation between the Mach number ahead, M1, and behind, M2, the wave:
(M12 - M22)(1 - γM12M22 + (γ - 1)/2(M12 + M22)) ) 0 (13) It is possible to find a similar relation also if we substitute eq (12) with Q ) I:
subsonic wave (both the Mach number on the left and on the right side is less than unity) have been analyzed by Astarita et al.5; in this paper only the stability of supersonic-subsonic waves is considered. By supersonicsubsonic wave we indicate the situation where the Mach number either on the left or on the right is greater than 1 (the other one being less than 1). It has been shown5 that, in a subsonic wave, points in which the flow is supersonic cannot exist, and the flow is everywhere subsonic. Equations (10)-(12), taking into account eqs (1), (3), and (5), can be written in an equivalent form:
Ft + vFx + Fvx ) 0
(15)
vt + vvx + (1/F)px ) 0
(16)
Tt + (p/F)vx + vTx ) (Q - I)/F
(17)
(M1 - M2)(2γM1M2 - 1 + γ2M1M2(M12 + M1M2 + M22)) ) 0 (14) As expected, both eqs (13) and (14) admit the trivial solution M1 ) M2. The relevant solution may be found considering only the second part of the left-hand side of these equations. The difference between GD (dotted line) and GM (solid line) is shown in Figure 1, where M2 is plotted as a function of M1. In classical GD, at steady state, a shock wave exists only if the Mach number ahead the wave is greater than 1, and the shock will produce a compression. Points to the left of M1 ) 1 do not have a physical meaning. Both graphs are symmetrical with respect to the line of equation M2 ) M1; however, while in GD, as may be expected from the previous discussion, the curve intersects this line for M1 ) 1, in GM the point of intersection is M1 ) 1/x(3γ): there could be a completely subsonic wave. The stability and the difficulties that arise for a
In GM theory one has two more external scales than one has in GD: the intrinsic length scale D appearing in eq 3 and the rate of energy input per unit volume Q. This in turn implies that the set of eqs (2), (3), and (15)(17) can be made dimensionless by assigning only the mass flow rate, G:
τ ) tx(Q/GD);
b ) Fx(QD/G3);
φ ) vx(G/QD) (18)
θ ) TG/QD;
i ) I/Q;
β ) p/x(QGD);
which reduces the set of equations to
y ) x/D
Ind. Eng. Chem. Res., Vol. 38, No. 4, 1999 1181
Figure 5. Temperature versus y for an expansion wave as a function of time.
bτ + φby + bφy ) 0
(19)
φτ + φφy + (1/b)βy ) 0
(20)
θτ + (β/b)φy + φθy ) 1 - i
(21)
β ) bθ
(22)
i ) β xθ
(23)
In GD it is impossible to find a length scale, and when the diffusive effects are neglected, the width of a shock wave is zero. In GM the situation is completely different and the length scale plays an important role. Results and Discussion As already said, we are only interested in waves that start, or bring to a supersonic Mach number. Figure 2 represents what, in GD, may be called a shock wave, showing the solution reached at steady state. A step between the left and right fixed (in time) values is taken for all variables as the initial condition. Moreover, the two values are chosen such that they satisfy eq (14). Given the above dimensionless form of the equations, both the mass flow rate and the dissipation rate at the two boundary points are equal to 1. All the profiles start with a zero derivative, and then the curves present a first region (A) where a very sharp slope is observed. Thus, a second region (B) can be envisaged where the variations become smoother, and in the limit the derivative goes to zero. This behavior is more clear in the temperature profile, where the first derivative changes sign and the temperature, after having reached
a maximum, decreases to an absolute minimum. The presence of the two regions is related to two different length scales: the first, very small, is due to artificial viscosity effects introduced by the numerical integration; the second is related to relaxation phenomena occurring after the compression wave.5 It is interesting to note that the maximum temperature is practically equal to the value that one would calculate from the Ocone and Astarita9 analysis of the Rankine-Hugoniot shock wave, and the width of the relaxation region is of the same order of magnitude found by Astarita et al.5 for completely subsonic waves. This last result is of practical importance; indeed, one can get all the information needed on the wave, solving algebraic equations as those in Ocone and Astarita,9 and then separately consider the relaxation zone for completely subsonic waves, starting from a nonequilibrium state (for nonequilibrium state we mean a situation where the condition Q ) I is not satisfied). In other words, the global problem can be solved by solving two simpler problems. The difference between the two regions is even more evident in Figure 3 where a “shock wave” with a higher Mach number on the left side is reported. A well-known result in GD is that for a diffusive fluid the width of a shock wave decreases as M1 increases; analogously, here the length scale of region A drastically reduces with respect to the one reported in Figure 2, while the relaxation scale decreases only slightly. Much more complex is the behavior of an expansion wave, and this is the reason numerical results are presented in a different way. In Figures 4 and 5, the distribution of the Mach number and temperature along y is reported at different times. Starting from the left-
1182 Ind. Eng. Chem. Res., Vol. 38, No. 4, 1999
hand side, each curve is taken at 1000 time unit differences. The phenomenon corresponds to a relatively low time and by no means does it represent a steady solution; as the time increases the curves become smoother (the phenomenon tends to die out). First of all, it is worth noticing that compared to the previous figures (Figures 2 and 3) the length scale is of various orders of magnitude greater: this proved that expansion waves are unstable. From an analysis of the curves three different regions can be envisaged: a first one (A) where both M and θ increase; a second one (B) where the Mach number decreases and the temperature increases; a third one (C) where M increases and θ decreases. It is quite difficult to make final conclusions from the results discussed; nevertheless, some conjectures can be made. Region (C) appears to be quite stable, and it is possible to assimilate this region to a wave which propagates at a higher speed within the undisturbed fluid on the right-hand side. The speed of such a wave can be easily obtained from the figures; indeed, within the region, the distance between two successive curves remains constant, and practically equal to 1000. The local maximum of the Mach number (end of region (A)) increases with time, and the reverse is true for the local minimum (end of region (B)). A similar behavior is observed for the temperature, but in this case the absolute maximum decreases with time. Conclusions The two possible types of supersonic waves (compression and rarefaction) in granular flow have been analyzed through a numerical approach. The constitutive equations used are the ones which correspond to dilute systems (low values of the solid volume fraction). Compression waves are found to be stable. Two different physical phenomena, with different length scales, govern the problem: the first one is similar to what is normally observed in classical gas dynamics (i.e., a sudden compression, which immediately brings the thermodynamic variable to a different state) is observed. A relaxation zone, characterized by a much larger value of the length scale, follows the compression region. The occurrence of the two zones makes it possible to separate the initial problem into two simpler ones. Supersonic compression waves cause, even for values of the Mach number ahead of the wave very close to unity, a very strong compression and they should therefore be avoided. Expansion waves are unstable; such a result, of very practical importance, is in accordance with classical gas dynamics. The instability makes it impossible to have an internal mechanism able to accelerate the flow (and then to reach the final
velocity). Again, this is in accordance with classical gas dynamics; analogously, we believe that a nozzle is needed in this case as well. We are currently working on the latter problem and we are confident that techniques used for gases would work in granular materials as well. Dedication We are pleased to dedicate this paper to Professor Roy Jackson. One of us (R.O.) has had the honor of being introduced to the field while being a graduate student of his. Literature Cited (1) Ocone, R.; Astarita, G. A Pseudo-Thermodynamic Theory of Granular Flow Rheology. J. Rheol. 1993, 37, 727-742. (2) Jean, R. H.; Fan, L. S. On the Model Equation of Gibilaro and Foscolo with Corrected Bouyancy Forces. Powder Technol. 1982, 72, 201-205. (3) Foscolo, P. U.; Di Felice, R.; Gibilaro, L. G. The Pressure Field in an Unsteady-State Fluidized Bed. AIchE J. 1989, 35, 1921-1926. (4) Zucrow, M. J.; Hoffman, D. J. Gas Dynamics; J. Wiley: New York 1976; Vol. 1. (5) Astarita, T.; Ocone, R.; Astarita, G. The Rayleigh Approach to the Rheology of Compressible Granular Flow. J. Rheol. 1997, 41, 513-527. (6) Ocone, R.; Astarita, G. On Waves of Particulate Phase Pressure in Granular Materials. J. Rheol. 1994, 38, 129-139. (7) Hugoniot, H. Me´moire sur la Propagation du Mouvement dans un Fluid Indefini. J. Math. Pur. Appl. 1887, 3, 477-491; 1887, 4, 153-167. (8) Hugoniot, H. Me´moire sur la Propagation du Mouvement dans les Corps et spe´cialement dans les Gaz Parfaits. J. Ecole Polyt. 1888, 58, 1-17. (9) Ocone, R.; Astarita, G. Compression and Rarefaction Waves in Granular Flow. Powder Technol. 1995, 82, 231-237. (10) Astarita, T.; Ocone, R.; Astarita, G. Compressible Flow of Granular Materials. In Particulate Flows: Processing and Rheology; IMA Volumes in Mathematics and Its Applications; Drew, D. A., Joseph, D. D., Passman, S. L., Eds.; Springer-Verlag: New York, 1997; Vol. 98, 1-22. (11) Astarita, T.; Ocone, R.; Astarita, G. Extension of Dynamics of Polymeric Liquid Methodology to New Areas. J. Non-Newtonian Fluid Mech. 1998, 76, 5-25. (12) Morse, P. M.; Feshback, H. Methods of Theoretical Physics; McGraw-Hill: New York, 1953. (13) Courant, R.; Hilbert, D. Methods of Mathematical Physics; Interscience: New York, 1952. (14) Astarita, T.; Ocone, R.; Astarita, G. Thermal Relaxation Behind Pressure Waves in Granular Flow. Fluidization VIII; Tours, France, 1995; Vol. I, pp 539-547.
Received for review June 17, 1998 Revised manuscript received September 14, 1998 Accepted September 15, 1998 IE980393Z
