- Email: [email protected]
Particle phase boundary layer theory in vertical two-phase gas–solid flow
Physica A 255 (1998) 26–47
Particle phase boundary layer theory in vertical two-phase gas–solid ow Y.A. Sergeev a; ∗ , P.L. Iske b , V.N. Kurdyumov c , J.H.J. Moors d a Department
of Engineering Mathematics, University of Newcastle, Stephenson Building, Newcastle upon Tyne NE1 7RU, UK b Knowledge Management, ABN-AMRO, Global Transaction Services, P.O. Box 263, 1000 EA Amsterdam, The Netherlands c Universidad Politecnica de Madrid, E.T.S.I. Aeronauticos, Plaza Cardenal Cisneros, 3, 28040, Madrid, Spain d Faculty of Electrical Engineering, Eindhoven University of Technology, P.O. Box 513, NL 5600 MB Eindhoven, The Netherlands Received 23 October 1997; received in revised form 5 January 1998
Abstract To gain understanding of the two-phase vertical gas–solid ow behaviour, the equations of gas– solid suspension ow are used to develop a model to predict the radial pro les of the volumetric fraction and the velocity of the solid particles. We derived two coupled di erential equations for the volume fraction of the particles and the uctuations of the particle velocity. These could be solved approximately, assuming the existence of a lean phase in the centre (the core) and a dense particle phase near the wall (the annulus). In this boundary layer theory, the method of matched asymptotic expansions could be employed. As a result, the obtained expressions have a generic character and the observed qualitative trends are insensitive with respect to the detailed c modelling of the particle–turbulence and particle–particle interactions. 1998 Elsevier Science B.V. All rights reserved. Keywords: Kinetic theory of gas–particulate suspensions; Multiphase ow; Two-phase gas– particulate ow
1. Introduction Dispersed ow phenomena are encountered in a wide variety of industrial applications. In two-phase dispersed gas–solid as well as gas–liquid pipe ow a so-called core–annulus behaviour can be observed. This means that the high-density component (i.e. the solid particles or the liquid drops) is found preferentially near the wall. This ∗
Corresponding author. Fax: 0191 222 5498; e-mail: [email protected].
c 1998 Elsevier Science B.V. All rights reserved 0378-4371/98/$19.00 Copyright PII S 0 3 7 8 - 4 3 7 1 ( 9 8 ) 0 0 0 6 2 - 4
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
27
a ects the ow pattern: the central lean phase (the core) owing faster than the annulus dense phase. Under certain conditions down ow near the wall can even occur. An accurate computational model, based on a correct description of the physical phenomena occurring in a vertical two-phase ow, is required to interpret and extrapolate the experimental data. However, to predict qualitatively the impact of a change in the geometrical and physical parameters on the ow behaviour, a less detailed model can be used. Therefore, we developed a model, which only contains the essential physics. The model focuses on the steady state fully developed ow. This case represents a one-dimensional problem. This problem has been addressed in many other publications. Sinclair and Jackson [1] discussed a model for fully developed ow in which the particles interact with the gas. The velocity pro le of the gas is modi ed by the presence of the particles, but remains parabolic, i.e. high velocity in the centre of the pipe and low velocity at the wall. The velocity pro le of the particles is induced by the drag forces on the particles. Because of the velocity gradients, uctuations in the particle velocity will be generated as particles in adjacent layers collide with each other. The random motion of the particles is associated with a quasi-temperature, in analogy with the kinetic theory of gases. It was shown by Pita and Sundaresan [2] that in the model of Sinclair and Jackson the inelasticity of the collisions causes an unrealistically strong e ect on the
ow phenomena. Other who used either the collisional approach or a combination of turbulence modelling and the collisional approach include Simonin [3], Louge et al. [4] and Ding and Gidaspow [6]. The work of the rst two authors so far, is restricted to (very) low particle volume fractions (i.e. smaller than 0.5%), whereas the latter authors concentrated on applying kinetic theory which was derived to describe rapid granular
ow. Based on the collisional theory of gas–particulate suspensions, the asymptotic methods have been successfully applied in the recent work by Sergeev and Zhurov [5] to develop the approximate analytical theory of vertical gas–particulate ow at low and moderate pressure gradients. Because of the over-sensitivity of these collisional models towards the restitution coecient of the particle–particle collisions, Dasgupta et al. [7] recently extended the model of Sinclair and Jackson in order to incorporate the interaction between the turbulent gas velocity eld and turbulent motion of the particles as source for the
uctuations. However, some aspects of the model remain speculative, e.g. the closure relations that are used are the same as used in the granular kinetic theory with some adjustable parameters. This means that in the modelling of the particle uctuations still collisions are held responsible, whereas it is very unlikely that the (direct) particle collisions are important. Based on the particle diameter, the Reynolds number is of the order one, and in this regime particles hardly collide. In order to account for the turbulence, sometimes a conventional Reynolds decomposition and time-averaging procedure is applied to the two-continua transport equations for mass, momentum and energy. This introduces various unknown correlations between the uctuations in the velocities, pressure and volume fractions of both phases, for which closure relations have to be speci ed. A generalisation of the one-phase k–
28
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
model to turbulent two-phase ows is often used to provide the required closure relations [ 7 –10]. In this approach, the turbulent properties of the uid phase follow from the transport equations for the turbulence intensity k1 and its dissipation rate 1 , which are solved in the model. The corresponding properties of the particulate phase are determined from suitably chosen correlation functions. The nature of the equations to be solved is such that usually severe numerical problems are encountered when solving them. The system of non-linear partial di erential equations form two-point boundaryvalue problem for which standard numerical recipes not always work. Furthermore, problems occur due to multiplicity of the solutions. It seems that the majority of the models predict the same qualitative behaviour of the gas–solid ow, including the segregation phenomenon. But the question is how sensitive these qualitative results are with respect to the detailed microscopic description, such as particle–particle interactions and particle–turbulence interactions. In this article we show that under very general assumptions the qualitative behaviour of two-phase gas– solid ow can be predicted. Existing models, such as the ones based on kinetic theory for granular ow or those based on detailed turbulence modelling, can all be shown to be realisations of the present model. Because of its relative simplicity the model can be used to predict trends in the two-phase ow behaviour. On the other hand, because there is no detailed microscopic description of the dispersed ow, the model yields qualitative results only. Some model parameters appear which have to be determined by a more fundamental study or by tting experimental data. To study the behaviour of relatively dense gas–solid systems we take the “Eulerian– Eulerian” approach in which the gas and the particle phase are treated as two interpenetrating uids, which are governed by their own separate transport equations valid at each point of the ow. These transport equations are derived from the microscopic conservation laws for mass, momentum and energy by an averaging procedure, for which di erent schemes have been proposed (see, e.g., [11]). The advantage of the continuum description is that it contains a limited number of independent variables. However, since the details of the microscopic ow-behaviour have been “averaged-out”, one has to postulate additional constitutive equations for the e ective stress tensors and interfacial transfer terms for mass, momentum and energy. To nd an approximate solution of the model equations, we base ourselves on the observation that in developed ow a lean region in the centre of the pipe (the core) and a dense region near the wall (the annulus) are found. Because the annulus is rather thin it seems appropriate to study the problem by applying techniques used in boundary-layer theory.
2. Description of the model In this section we will rst review the transport equations for the mean quantities used in the two-continua model. Then we will present the model assumptions and nally the equations will be reduced to their one dimensional steady-state form for the mixture.
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
29
2.1. Basic equations The basic transport equations solved in the two-continua model represent the conservation of mass and momentum for each of the individual phases separately. The variables solved from the mass and momentum conservation equations are the volume fractions g and s of the gas and solids phase, respectively and the velocities ug and us . We will summarise the relevant equations below: Mixture continuity: g + s = 1:
(1)
Gas phase continuity: @ g + ∇ · ( g ug ) = 0: @t
(2)
Solids phase continuity: @ s + ∇ · ( s us ) = 0: @t
(3)
Gas phase conservation of momentum: g
@ ( g ug ) + g ∇ · ( g ug ug ) = g ∇ · g + g g g + sg (us − ug ): @t
(4)
Solid phase conservation of momentum: @ ( s us ) + s ∇ · ( s us us ) = s ∇ · g + ∇ · s + s s g + sg (ug − us ): (5) @t In these equations, g and s are the gas and solid densities, respectively, g and s are the e ective stress tensors of the gas and solids phase, and sg is the viscous drag function between the two phases. Now we consider the steady plane-parallel vertical ow of gas–solid mixture in a vertical cylindrical pipe (i.e. no time dependence and no dependence on the z- and -coordinates). With the above assumptions the basic system of equations (1) – (5) reduces to the equations of conservation of the gas and solids axial momentum and the equation of conservation of radial momentum of solid particles. From here on the volumetric fraction of the solid particles will be denoted by . These equations are, respectively, s
1 d dP − (1 − ) (rg ) − (1 − )g g + sg (us − ug ) = 0; dz r dr rz 1 d 1 d dP g − (rrz (rs ) = 0; ) − s g + sg (ug − us ) − − dz r dr r dr rz d g d s rr = rr = 0; dr dr − (1 − )
(6) (7) (8)
where z is the vertical axis and P is the gas pressure. A further simpli cation can be obtained by the addition of Eqs. (6) and (7). The shear stresses in the gas phase are
30
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
neglected with respect to the shear stresses in the solids phase (which has been veri ed to be a good approximation), with as result the following axial momentum equation for the mixture: dP 1 d + [(1 − )g + s ]g + (rs ) = 0; dz r dr rz
(9)
which does not contain the viscous drag function sg . Inside the mixture the drag forces cancel. The radial momentum equation for the mixture implies d s = 0: dr rr
(10)
2.2. Constitutive equations Now we consider the closure for the system of equations (9) and (10). For that s and the bulk purpose we will write down constitutive equations for the shear stress rz s stress rr . We now introduce a uctuation energy (or “quasi-temperature”), , which is equal to one third of the mean square of the velocity uctuations of the particles: 1 = us0 us0 : 3
(11)
In this expression, us0 denotes the uctuation in the particle velocity. A remark should be made here. For the values of parameters relevant for these ows the value of the viscous drag coecient sg is very large, which means that the gas velocity ug di ers very little from the solids velocity us in the bulk ow, as well as in the ow close to the wall (we will take us = ug = u). The set of equations (9)–(10) is complemented with a conservation equation for the
uctuation energy in the solids phase. This equation is written as 1 d s du (rq) + rz + = 0; r dr dr
(12)
where q is the pseudo-heat ux and the rate of dissipation of the energy of random particulate motion which will be discussed separately. The interpretation of as a thermodynamic property of the system suggests that results, obtained in the kinetic theory for gases, can be employed to nd analogous expressions for the closures for our model. In this work, expressions for the stresses will be derived which can be related to the kinetic theory for dense gases [12]. However, we will show that in order to pursue our analysis, we do not necessarily need the detailed expressions for e.g. the viscosity and the heat conduction coecient, but that much more simple relations, showing the right asymptotic behaviour, are sucient. The constitutive relations for s and q are written as s rz = −
du ; dr
(13)
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
31
s rr = − C;
(14)
d : dr
(15)
q= −
In these equations denotes the apparent viscosity of the solids phase, C is a proportionality constant (similar to the ideal gas law), and is the quasi-temperature–heat conduction coecient in the solids phase. Inspired by the beforementioned analogy with the kinetic theory for gases we propose the following expressions for the transport coecients, appearing in the previous equations: √ (16) = f1 ( )s ( ) ; C = f2 ( )s ;
(17)
√
= f3 ( )s ( ) :
(18)
Eqs. (9), (10), (12) with the constitutive relations (13) – (15) represent the closed system for the unknown parameters , u and . Of course, explicit expressions for the functions f1 , f2 and f3 as well as for the “mean free path” (i.e. the distance over which the particles exchange their momentum and uctuation energy) have to be known, before getting quantitative results. For the analysis presented here, however, we only need the following properties of the functions f1 , f2 and f3 : lim f1 ( )( ) = o1 ;
(19)
lim
f2 ( ) = o2 ; →0
(20)
lim f3 ( )( ) = o3 ;
(21)
→0
→0
where o1 , o2 and o3 are positive and nite constants. It is interesting to note that, although there is no physical basis for the use of the collision theory in these systems (in which the gas plays an important role and the particles hardly collide), the expressions derived in the collision theory do satisfy the constraints (19)–(21). In particular, for f2 ( ) one derived f2 ( ) = (1 + 4 g0 );
g0 = (1 − ( = max )1=3 )−1 ;
(22)
in which max is the volume fraction corresponding to random close packing (which is estimated at max ' 0:65). Expression (22) can be observed to give correct limiting behaviour described by Eq. (20), with o2 = 1. In the kinetic theory for granular ow, the limiting values of 1 and 3 are nite: the mean free path diverges as −1 for small , but the functions f1 ( ) and f3 ( ) become linear in the volume fraction. We will now discuss the rate of dissipation of the uctuation energy , which appears in Eq. (12). We write
= !3=2
(23)
32
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
which is consistent with expressions for dissipation caused by rotations. It is interesting to notice that the same form for ! is obtained in the framework of the collisional theory by Lun et al. [13] where ! turns out to be a function of . Now we introduce dimensionless variables. The scale of length will be R, the radius of the pipe and velocities will be scaled with V, the super cial gas velocity, and the quasi-temperature will be scaled with V 2 . As a scale for the mass density we use s , the mass density of the solids phase. Expressed in the dimensionless variables (denoted by a tilde) the set of equations (9), (10), (12) reads 1 1d dP˜ + [(1 − ) + ] − d˜z F r˜ r˜
p du˜ ˜ rf ˜ 1 ( )( ) ˜ = 0; d˜r
d ˜ = 0; (f2 ( )) d˜r 1 d r˜ d˜r
(24) (25)
p d˜ ˜ rf ˜ 3 ( )( ) ˜ d˜r
!
p du˜ 2 ˜ + f1 ( )( ) ˜ + ˜ = 0: d˜r
(26)
In Eq. (24) the Froude number F appears, given by F = gR=V 2 . The density ratio = g =s is usually very small (about 10−3 ). Since ' 10−2 –10−1 it can be concluded that (1 − ) is small compared to and can be neglected. From now on we will skip the tilde in the notation for dimensionless variables; where confusion is possible we will specify whether we use dimensional or dimensionless variables. Now we consider the boundary conditions for the problem. The closed system of equations (24)– (26) will be considered in combination with the boundary conditions: d du = =0 dr dr u=0 +
at r = 0;
(27)
at r = 1; d =0 dr
(28) at r = 1:
(29)
The last boundary condition gives a balance between the ux (second term) and the rate of dissipation of uctuation energy at the wall ( rst term). Eqs. (24)–(26) can be conveniently rewritten in a form containing only two unknown functions, (r) and (r). Integrating Eq. (24) and using the boundary condition (27) yields for the velocity gradient: 1 1 1 dP du √ r = + dr f1 ( )( ) 2 dz rF
Zr 0
r 0 dr 0 :
(30)
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
33
Substituting Eq. (30) in Eq. (26) we obtain the system of equations for the functions (r) and (r) √ d 1 d rf3 ( )( ) r dr dr 2 Zr 1 dP 1 1 √ r + r 0 dr 0 + = 0; (31) + rF f1 ( )( ) 2 dz 0
d f2 ( ) = 0: dr
(32)
One more boundary condition is required. Such a condition can be, for example, an integral condition for the total ux of the solids phase. The dimensionless form of this condition is as follows: Z1 Qs = 2
us r dr;
(33)
0
where Qs is assumed to be given. Since the condition (33) is of the integral form, it causes certain inconveniences to nd a solution of the problem. Instead of (33) we shall use the supplementary boundary condition for in the form = ∗
at r = 1;
(34)
where ∗ is supposed to be known. The problem (21), (22), with the boundary conditions (27)–(29) and (34) is closed. As soon as the solution of this problem is obtained, Qs can be found from Eqs. (33) and (34) as a function of ∗ : Qs = Qs ( ∗ ).
3. Particle phase boundary layer theory 3.1. Introduction As mentioned before, the two-phase dispersed ow often shows a core-annulus behaviour. The rather thin annulus can be interpreted as an ‘extended’ boundary layer, therefore we are inspired to apply a boundary layer theory. To make this point clear, we consider the boundary condition (29). We will consider as the small parameter that is required to employ a boundary layer approach. (In the collision theory would have a value of the order of several tens (dimensionless) particle diameters, which indeed is rather small (about 10−4 –10−3 )). Due to the fact that the small parameter is in front of the derivative in the boundary condition, the formation of a boundary layer of the particle phase in the vicinity of the wall (i.e. as r → 1) must be expected.
34
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
An estimate of the thickness of the boundary layer as a function of ∗ = |r = 1 can already be given: from Eq. (32) we nd for r → 1: 0 ∗ d ∗ ∗ d = 0; (35) f2 ( ) + f2 ( ) dr r=1 dr r=1 1 f2 ( ∗ ) d ; (36) = dr r=1 f20 ( ∗ ) ≡
∗ f0 ( ∗ ) = ∗ 2 ∗ = O( ) (d =dr)r=1 f2 ( )
(37)
where (0 ) ≡ d=d . So, the thickness of the boundary layer is of the same order as . For reasons of convenience we once more give the set of equations with the boundary conditions from which the radial pro les of the particle uctuations and volume fraction can be obtained: √ d 1 d rf3 ( )( ) r dr dr 2 Zr 1 dP 1 1 √ r + r 0 dr 0 + = 0; (38) + rF f1 ( )( ) 2 dz 0
d f2 ( ) = 0 ; dr d d = =0 dr dr = ∗ +
(39) at r = 0;
at r = 1;
d =0 dr
at r = 1:
(40) (41) (42)
Due to the non-linearity, it turns out that it is dicult to analyse this system of equations numerically. In particular, the in uence of physical and geometrical parameters is hard to determine. Therefore, we prefer to formulate an approach that gives an approximate solution, but which allows us to predict trends in the ow behaviour. We will use the method of singular perturbations (see, e.g., [14]) in the form of matched asymptotic expansions for the problem (38) – (42). In accordance with this method we consider the following two domains: The boundary layer, i.e. the vicinity of the wall, 1 − r¡ = O( ). Below we denote the particle phase concentration and quasi-temperature in this region as l and l , respectively. The bulk, i.e. the central part of the ow, 06r¡1 − = 1 − O( ); the particle concentration and the quasi-temperature are written as b and b , respectively. In each of these regions the solution is sought in the form of asymptotic expansions of the functions (r) and (r) in terms of the parameter . We will seek the principal
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
35
terms of such expansions, only keeping the terms in lowest order in . In our analysis, rst we derive an equation for the quasi-temperature . In this paper we assume that the pressure drop is larger than required to compensate the gravitational forces on the solids in the mixture. This is the situation that is most studied in the literature [1,2,6]. We assume dP = O( 0 ): (43) dz 3.2. Quasi-temperature in the boundary layer We start our analysis by considering ow in the boundary layer. Here we make the following a priori assumption ∗ = O( m );
(44)
where m¿0. The proof of this assumption as well as the value of m will be given in the process of solving the problem (38)–(42). Since 061 − r6O( ) within the boundary layer, we introduce a new variable: y=
1−r :
(45)
This variable is found in the interval [0; −1 ]; within the boundary layer we have 06y6O( 0 ). Now we consider Eq. (38). Since ¡1, the integral is of the order of unity (i.e. O( 0 )) or smaller. The dissipation is supposed not to be larger than the second term. Keeping only terms of the lowest order in , only the rst term in Eq. (38) remains, so we can write p dl d f3 ( l )( l ) l = 0: (46) dy dy So, we obtain ( l )f3 ( l )
d3=2 l = C1 : dy
(47)
The boundary conditions at the wall read l = ∗ l =
dl dy
at y = 0; at y = 0:
(48) (49)
From d f2 ( l )l = 0 dy
(50)
it obviously follows f23=2 ( l )3=2 l = C2 :
(51)
36
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
Using Eqs. (47) and (51) and the boundary conditions (48) and (49) we obtain: C2 =
2 f23=2 ( ∗ ) C1 : 3 f3 ( ∗ )( ∗ )
(52)
It follows that the radial distribution of the quasi-temperature in the annulus can be written in the form: 2=3 f2 ( ∗ ) 2 C12=3 : (53) l = 2=3 ∗ 2=3 ∗ 3 f2 ( l )f3 ( ) ( ) The constant C1 is still unknown and must be obtained by asymptotic matching of the outer solution, b , and the boundary layer solution, l : b (r)|r→1 ∝ l (y)|y→∞ ;
(54)
where the symbol ∝ is used to indicate the coincidence of leading asymptotic terms with respect to . In order to give an asymptotic expression for the quasi-temperature we rst study the distribution of the solids fraction in the boundary layer. Combining Eqs. (47) and (51) leaves us with ( l )
d l + q( ∗ ) = 0; dy
(55)
where ( ) =
3 f3 ( )f20 ( )( ) 2 f25=2 ( )
(56)
q( ) =
3 f3 ( )( ) : 2 f23=2 ( )
(57)
and
From Eq. (55) it follows directly that the solids volume fraction in the boundary layer does not depend on any parameter except for the solids concentration at the wall ∗ ! From the properties of the function ( ) it also follows that limy→∞ l = 0. Therefore, to nd the asymptotic behaviour of the particle concentration at the outer boundary of the boundary layer (as y → ∞) and, consequently, the order of magnitude of the solids volume fraction in the bulk, the properties (20) and (21) can be used so that in Eq. (56) f3 ( )( ) and f2 ( ) can be replaced by f3 (0)(0) = 3
and
f2 ( ) = 2 ;
(58)
respectively, where 3 and 2 are the dimensionless quantities corresponding to the parameters o3 and o2 (note that 2 = o2 ). Then the boundary layer equation (55) takes the form 3=2 2 d l = − q( ∗ ) 2 5=2 : dy 3 3 l
(59)
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
37
From Eq. (59) the asymptotic behaviour of the solids volume fraction at the conventional boundary between the bulk and the boundary layer follows in the form l ∼
2=3
3 1 1 2=3 ∗ q ( ) 2 y2=3
as y → ∞:
(60)
Here we arrive at the important conclusion that the bulk is indeed much more dilute than the boundary layer; the order of magnitude of solids volume fraction in the bulk: b = O( 2=3 ):
(61)
We can also give the asymptotic expression for the quasi-temperature in the boundary layer. From Eq. (47) it becomes clear that 2=3 C1 y2=3 for y → ∞: (62) l (y) = 3 3.3. Quasi-temperature in the bulk Now we consider the solution in the bulk region. From the estimate (61) and assumption (43) it follows that the integral term in Eq. (38) can be neglected compared to the gas pressure gradient term. Since the solid volume fraction in the bulk is low the combinations f3 ( b )( b ) and f1 ( b )( b ) in Eq. (38) can be replaced by the (non-zero!) zero order terms in their Taylor expansions, i.e. by 3 and 1 respectively, where 1 = f1 (0)(0) is the dimensionless quantity corresponding to the parameter can now be written as p (dP=dz)2 r 2
1 d db √ r b + + = 0: r dr dr 4 1 3 3 b
(63) o1
in Eq. (19). Eq. (38) (64)
Eq. (64) has to be solved subject to the condition of symmetry (40) and the condition of asymptotic matching with the boundary layer solution (54). To nd the unknown constant C1 , that appeared in the solution for the quasi-temperature in the boundary layer, we require that the solution b for the bulk region of Eq. (64) must have the asymptotic behaviour of the form (62) as r → 1. This means that b → 0 as r → 1, so that instead of Eq. (64), with the boundary conditions (40) and (42) formulated within the interval 06r61 − O( ), we can consider Eq. (64) within the interval 06r61, together with the boundary condition (40) and b = 0
at r = 1:
(65)
We can write Eq. (64) in a convenient form by introducing a new variable , de ned by −1 dP √ (66) = 2 1 3 b : dz
38
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
The resulting equation reads: r2 1 d p d r + √ + 0 = 0 r dr dr
(67)
with 0 the modi ed uctuation energy dissipation rate in terms of . The boundary conditions for Eq. (67) are d = 0; (68) dr r=0 |r=1 = 0:
(69)
Eq. (67) can be slightly rewritten: 2 p d 1 d d r + + r 2 + 0 = 0: r dr dr 2 dr
(70)
Note that if we neglect the dissipation of uctuation energy, the solution of Eq. (70) yields a universal radial pro le of , which does not depend on any parameter! The asymptotic analysis of Eq. (70) without dissipation with the boundary conditions (68) and (69) yields: = A(1 − r)2=3
as r → 1:
(71)
This can be validated by substituting Eq. (71) into Eq. (70), observing that the divergent terms cancel. The conclusion is that the behaviour of the particle uctuations in the bulk near the wall is consistent with the behaviour of the particle uctuations in the boundary layer, given by Eq. (62). The method of matched asymptotic expansions can now be applied, since A is a xed constant that can be determined numerically. In case of dissipation, the constant A will depend on the rate of dissipation. If we write (compare to Eq. (23)):
0 = !3=2
(72)
we can obtain the function A(!) (A(0) ' 0:969) using a transformation of the radius: r → 1 − rˆ3=2 . The form (72) is consistent with formulae used in the kinetic theory of gases. In Fig. 1 we plotted the function A(!). In Fig. 2 we plotted the solution of Eq. (70) for several values of !. As follows from the results represented in Fig. 2, in absence of the dissipation the radial quasi-temperature pro le is monotonic. The dissipation leads to the appearance of maximum of the particle uctuations at a certain distance from the wall. 3.4. Matched asymptotic solutions We are now able to nd the quasi-temperature in the boundary layer, because the constant C1 in Eq. (60) can now be determined. Equating (62) and (71) and using the
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
39
Fig. 1. The matching factor A(!).
Fig. 2. Quasi-temperature pro les for di erent values of the dissipation factor !.
scaling (66) yields: C1 =
A3=2 3=4 1=4 3 (4 1 )
3=2 dP : dz
(73)
For the uctuation energy in the wall region we obtain from Eq. (53) l =
2=3 2=3 f2 ( ∗ ) 1 3 2 A(!) 2=3 ∗ ∗ 3 f3 ( )( ) f2 ( l ) (4 1 3 )1=2
dP : dz
(74)
Now we come to the a priori assumption (44): we recall that ∗ = O( m ). From Eq. (74) we nd: ∗ =
2=3 2=3 1 3 2 A(!) 2=3 3 f3 ( ∗ )( ∗ ) (4 1 3 )1=2
dP : dz
(75)
It follows then that m = 2=3:
(76)
40
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
It is now possible to derive a relation between the volume fraction of the particles in the bulk, b , and the volume fraction at the wall, ∗ . For the bulk we use the solution of Eq. (70), with the result |dP=dz| (r; !): b (r) = √ 4 1 3 Now we use Eq. (39) yielding: 2=3 2=3 3 2 2=3 A(!)f2 ( ∗ ) = f2 ( b (r))(r; !): 3 f3 ( ∗ )( ∗ )
(77)
(78)
Since the leading expansion term of the solids concentration is of the order 2=3 , using Eq. (78) and properties (20) and (21) to nd the radial distribution of in this region we obtain 2=3 2=3 A(!) 3 f2 ( ∗ ) 2 2=3 : (79) b (r) = 3 (r; !) 2 (f3 ( ∗ )( ∗ ))2=3 It might be useful to estimate the solids volume fraction at the centre of the pipe, b (0), as a function of the dissipation rate, !. From Eq. (79) it follows: 2=3 2=3 A(!) 3 f2 ( ∗ ) 2 2=3 : (80) b (0; !) = 3 (0; !) 2 (f3 ( ∗ )( ∗ ))2=3 To conclude this section we note that the following approximate form for the solids quasi-temperature, tot , valid within the whole interval 06r61 can be suggested. The function tot should satisfy tot = b
as r → 0;
(81)
tot = l
as r → 1:
(82)
We choose:
"
(1−r)= # (1−r)= 1 1 l : b + tot = 1 − 2 2
(83)
It should be noted that tot is not a solution to the equations. However, from this function, together with the boundary condition (41) for the particle volume fraction at the wall, we can calculate tot , the radial distribution of the solids volume fraction, using Eq. (39). With the functions tot and tot , we are able to determine the radial velocity pro le by integration of Eq. (30). 4. Velocities Below, the radial distribution of the velocity of the solid particles in the annulus and the bulk region will be obtained; since ub ' ug the formulae give simultaneously the
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
41
pro le of the gas velocity. The (solids) velocity can be obtained from Eq. (30), with the boundary condition at the wall (28). Since = O( 2=3 ) in the bulk region, and the thickness of the boundary layer in which = O( 0 ) is of the order , the integral term in Eq. (30) can be neglected for all r. Then we get the simpli ed expression: r dP=dz du = √ : dr 2 f1 ( )( )
(84)
By integrating this formula we obtain the velocity. In the boundary layer we obtain 1 dP ul (r) = − 2 dz
Z1 √ r
r0 dr 0 ; l f1 ( l )( l )
(85)
where we used u(1) = 0. The length of the integration interval is of the order and the integrand is O( −1=3 ). The conclusion is that the velocity in the boundary layer is O( 2=3 ). In the bulk the velocity is given by 1 dP ub (r) = − 2 dz 1 dP − 2 dz
1− Z
√
r0 dr 0 b f1 ( b )( b )
√
r0 dr 0 : l f1 ( l )( l )
r
Z1 1−
(86)
The rst integral in the RHS is O( 0 ). We have already seen that the second integral is O( 2=3 ), so it can be neglected in Eq. (86). Because = O( ) we can replace the upper limit 1 − in the integral by 1. Using the variable , introduced in Eq. (66), and replacing b in the combinations f3 ( b )( b ) and f1 ( b )( b ) by 3 and 1 respectively, we obtain 1=4 ub (r) = √ 3 3=4 2 1
s Z1 0 dP p r dr 0 : dz (r 0 ; !)
(87)
r
From this formula it can be concluded that the velocity in the bulk is proportional to the square root of the pressure gradient. Using the obtained radial distributions of solids velocity (87) and volume fraction (79) the total dimensionless solids ux given by Eq. (33) can be calculated as a function of the particle concentration at the wall ∗ . This gives ∗ = ∗ (Qs ) and, therefore, solves the problem in its original formulation.
(88)
42
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
To calculate Qs we have a contribution from the bulk and from the boundary layer: 1− Z Z1 ub b r dr + 2 ul l r dr: Qs = Qs; b + Qs; l = 2 0
(89)
1−
It is not dicult to demonstrate that the second integral can be neglected with respect to the rst one, which means that the solids ux is determined by the ux in the bulk only Z1 Qs ' Qs; b ' 2
ub b r dr:
(90)
0
From Eqs. (79) and (87) we obtain: 1 s Z1 Z 11=12 7=6 2=3 ∗ 0 2 A 3 f2 ( ) r r dP p dr 0 dr: Qs = dz (r) 3=4 0 2=3 ∗ ∗ 2=3 (r ) 3 1 (f3 ( )( )) 0
(91)
r
The double integral can be numerically evaluated (it depends on the dissipation factor !). We can write Qs in a convenient form s dP f2 ( ∗ ) 2=3 : (92) Qs = s (!) dz ∗ ∗ 2=3 (f3 ( )( )) Finally, we can calculate the dimensionless gas ux, which is given by Z1 Qg = 2
Z1 u(1 − )r dr ' 2
0
Z1 ub (1 − b )r dr ' 2
0
ur dr;
(93)
0
where we only took the contribution of the lowest order in ; in Eq. (93) is the density ratio; the velocity is given by Eq. (87) and we nd s s Z1 Z1 0 dP 2 1=4 dP r 3 r p dr 0 dr = g (!) : (94) Qg = dz 3=4 0 dz (r ) 1 0
r
We note that the particle volume fraction at the wall can be explicitly expressed through the total solids ux in case of relatively dilute two-phase suspension. Assuming low particle concentration throughout the entire pipe including the wall region we may write in Eq. (91) f2 ( ∗ ) ' 2 ∗ ;
f3 ( ∗ )( ∗ ) ' 3 :
(95)
Neglecting also the dissipation, the numerical value of A ' 0:969 and the integral in Eq. (91) can be estimated as 1 Z1 Z r r0 p (96) dr 0 dr ' 1:219 (r; 0) (r 0 ; 0) 0
r
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
43
so that for the particle concentration at the wall we nd: ∗ ' 0:25
3=4 1
2=3 2 1=4 3
Qs :
(97)
We complete this section with a remark on the shear stress and stability of the fully developed vertical ow of suspension of solid particles. The shear stress in the solids phase is given by the constitutive relations (13) and (16). Using the solution for the bulk (77), (79) and (87) it can be easily shown that dr z ¿0: (98) dr r=0 Using the boundary layer solution of Eq. (55) and the radial distributions of the solids quasi-temperature and velocity in the boundary layer, given by Eqs. (74) and (85), respectively, together with the value of the quasi-temperature at the wall (75) and the boundary condition for the solids velocity (28) the leading asymptotic term of the radial gradient of the e ective shear stress in the solids phase can be found in the form 1 dr z = −B( ∗ ) 2=3 ¡0; (99) dr r=1 where B = O(1) is the positive function of the solids volume fraction at the wall. From inequalities (98) and (99) it follows that the e ective shear stress in the solids phase reaches the maximum within the ow region. As shown by Joseph et al. [15] this might be evidence of the stability of the annular structure of the considered gas-particle ow in the vertical pipe. We also calculate here the shear stress at the wall. We have seen that for the velocity gradient at the wall one can write dP=dz du : (100) = √ dr r=1 2 ∗ f1 ( ∗ )( ∗ ) From the constitutive relations (13) and (16) we obtain the following remarkably simple relation between the wall stress and the pressure drop: r z |r=1 = −
1 dP : 2 dz
(101)
In the next section we will discuss results obtained by applying the approximate method, described in the previous sections, and compare these results with results obtained by numerical solution of the full set of Eqs. (38) – (42). 5. Accuracy, qualitative trends and conclusions A boundary layer theory was developed in order to obtain more insight in the in uence of physical and geometrical parameters on the model results. We will rst
44
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
make a comparison between results of the numerical scheme and results obtained by the method of matched asymptotic expansions. Though the theory does not depend on the explicit form of the closure relations, a choice for speci c closure relations has to be made in order to perform the calculations. Because the closure relations from the kinetic theory for granular ow are most commonly used, the values for o1 , o2 and o3 , de ned in properties (19)–(21), are based on the expressions given by Lun et al. [13]. For the suspension of ne solid particles (the diameter dp ∼ 10−4 m) these values read o1 ' 9:23 × 10−6 m;
o2 = 1;
o3 ' 3:46 × 10−5 m:
(102)
We take these values for our calculations. It should be noted, however, that the predicted trends are not sensitive towards the actual choice of the closures. To calculate the dimensionless parameters we consider the fully developed ow of solid particles of the diameter given above and the density s = 1:4 × 103 kg=m3 in the vertical pipe of the radius R = 0:5 m. The scale velocity is assumed to be V = 10 m=s, the gas density g = 1:2 kg=m3 , and the particle volume fraction at the wall ∗ = 0:2. First, the accuracy of the approximate solution described in Sections 3 and 4 will be tested. The solution of the set of equations (31) – (32) for the particle uctuation energy, the particle volume fraction and Eq. (30) for the velocity of the particles are calculated with both the approximate and numerical methods. The rst comparison is made at a relatively high pressure drop, dP=dz = −3000 Pa=m (the corresponding ˜ z˜ = −0:0214). In this case the particle volume fracdimensionless pressure drop dP=d tion obtained from the approximate solution di ers at most 3% from the complete numerical solution. The uctuation energy and the velocity di er a few percent. In Fig. 3 the radial distributions of the quasi-temperature, the particle volume fraction and the velocity are given. The particle volume fraction at the wall is found to be much larger than the particle volume fraction in the middle of the pipe, which is a result that is characteristic for all calculations. This is related to the radial distribution of the uctuation energy. In the middle of the pipe the uctuations are larger than at the wall. In fully developed ow particles prefer being in regions where the uctuations are smaller (cf. the ideal gas) and so they segregate towards the wall. Radial distributions for the particle uctuation energy, the particle volume fraction and the velocity for decreasing pressure drops show the same characteristics. The velocity and the quasi-temperature are decreasing with decreasing pressure drop, as can be concluded by comparing the results for a pressure drop of 3000 Pa=m (Fig. 3a) and 500 Pa=m (b) ˜ z˜ = −0:00357). It must be noted that the accuracy of the boundary (dimensionless dP=d layer theory decreases with decreasing pressure drop; the relative error between the numerical calculations and the approximate theory is slightly higher (about 5%) in the case (b). The important trend to be studied is the in uence of the particle diameter. In Fig. 4 the radial distributions of the velocity and the particle uctuation energy (quasitemperature) are given for various particle diameters. The in uence is studied by changing the particle diameter and adjusting ∗ and the dimensional pressure gradient dP=dz
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
45
Fig. 3. Comparison between the approximate and the numerical solution. The dimensionless quasitemperature, ˜ = =V 2 , the particle volume fraction, , and the dimensionless velocity, u˜ = u=V . Dashed line – numerical solution of the complete set of equations; solid line – approximate asymptotic solution. = 0:0068, ∗ = 0:2. (a) dP=dz = −3000 Pa=m; (b) dP=dz = −500 Pa=m.
in such a way that the super cial velocities of the gas and solids, i.e. Usup; g =
s Qg ; R 2
Usup; s =
s Qs ; R 2
(103)
respectively, where Qg and Qs are the dimensionless gas and solid uxes given by Eqs. (93) and (33), respectively, remain the same (the values of the parameters are given in Table 1 for all the calculations). We assume that the densities of the particles and the gas are those given above, and the pipe radius is R = 0:5 m. The super cial velocity of the gas is used as a scale velocity V in the solution described in Sections 3 and 4, V = Usup; g . We can conclude that the velocity in the centre of the pipe is larger for small particles. The velocity pro le is thus more peaked for smaller particles. Smaller particles uctuate more than larger particles (Fig. 4b). The calculated pro le of the quasi-temperature can now be used to nd the particle volume fraction pro le
46
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
Fig. 4. The in uence of the particle diameter obtained with the boundary layer theory on the radial pro les of (a) velocity, (b) uctuation energy (quasi-temperature). (···) dp = 150 m; (- - -) dp = 100 m; (– – –) dp = 70 m; (——) dp = 30 m. Table 1 Values for dimensional parameters used in calculations to study the in uence of particle diameter dp (m)
∗
dP=dz (Pa=m)
Usup; g (m=s)
Usup; s (m=s)
30 70 100 150
0:236 0:136 0:100 0:0635
−141:4 −158:9 −176:45 −208:6
20:00 20:00 20:00 19:99
0:1465 0:1466 0:1468 0:1467
in the bulk and the boundary layer from the approximate solution given in Section 3. It is dicult to illustrate graphically the in uence of the particle diameter on the radial distribution of the volume fraction, so we rather describe the found trends qualitatively. The particle volume fraction in the centre of the pipe decreases when the particle diameter is increased from 30 to 70 to 100 m and increases again when the diameter is further increased to 150 m. The particle volume fraction at the wall has to be larger for small particles to ensure that the super cial velocity of the solids phase is a constant. The ratio of the particle volume fraction at the wall to the particle volume fraction in the centre of the pipe is larger for small particles. From the obtained results it follows that the segregation is stronger for smaller particles, and the velocity pro le
attens with decreasing particle diameter.
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
47
In this paper we discussed a model in which the two-phase conservation equations for fully developed ow are solved by introducing a quasi-temperature, which is a measure for the particle uctuations. Closure relations for the viscosity, heat conductivity and the pressure were expressed in terms of this quasi-temperature. The equations are generic, i.e. the form of the equations does not depend on the microscopic model for the transport coecients such as the viscosity and the heat conductivity. Based on the observation that a large solids volume fraction exists in the wall region and a small solids volume fraction exists in the centre of the pipe a boundary layer theory for the solids phase has been developed. In this theory asymptotic solutions for the quasi-temperature and the volume fraction of the particles could be derived. This enables us to perform a qualitative analysis of the in uence of physical and geometrical parameters on a character of distribution of particle concentration and velocity in a rapid vertical two-phase pipe ow. A comparison was made between the results from the approximate asymptotic theory and results obtained by a direct numerical solution of the full set of equations. It was concluded that for studying trends qualitatively, the approximate method gives reasonable results, none of which are contradictory to experimental observations. References [1] J.L. Sinclair, R. Jackson, Gas–particle ow in a vertical pipe with particle–particle interactions, A.I.Ch.E. J. 35 (1989) 1473. [2] J.A. Pita, S. Sundaresan, Gas–solid ow in vertical tubes, A.I.Ch.E. J. 37 (1991) 1009. [3] O. Simonin, Prediction of the dispersed phase turbulence in particle-laden jets, Proc. ASME 1991 Conf. on FED, vol. 121: Gas–Solid Flows, 1991, p. 197. [4] M. Louge, E. Masorakos, J.T. Jenkins, The role of particle collisions in pneumatic transport, J. Fluid Mech. 231 (1991) 345. [5] Y.A. Sergeev, A.I. Zhurov, Asymptotic theory of two-phase gas–solid ow through a vertical tube at moderate pressure gradient, Physica A 236 (1997) 243. [6] J. Ding, D. Gidaspow, A bubbling uidization model using kinetic theory of granular ow, A.I.Ch.E. J. 36 (1990) 523. [7] S. Dasgupta, R. Jackson, S. Sundaresan, Turbulent gas–particle ow in vertical risers, A.I.Ch.E. J. 40 (1994) 215. [8] S.E. Elgobashi, T.W. Abou-Arab, A two-equation turbulence model for two-phase ow, Phys. Fluids 26 (1983) 931. [9] M.A. Rizk, S.E. Elgobashi, A two-equation turbulence model for dispersed dilute con ned two-phase
ows, Int. J. Multiphase Flow 15 (1989) 119. [10] C.P. Chen, P.E. Wood, A turbulence closure model for dilute gas–particle ows, Can. J. Chem. Eng. 63 (1985) 349. [11] D.A. Drew, Mathematical modelling of two-phase ow, Ann. Rev. Fluid Mech. 15 (1983) 261. [12] S. Chapman, T.G. Cowling, The Mathematical Theory of Non-uniform Gases, 3rd ed., Cambridge University Press, Cambridge, 1970. [13] C.K.K. Lun, S.B. Savage, D.J. Je rey, N. Chepurniy, Kinetic theories for granular ow: inelastic particles in Couette ows and slightly inelastic particles in a general ow eld, J. Fluid Mech. 140 (1984) 223. [14] M.D. Van Dyke, Perturbation Methods in Fluid Mechanics, Parabolic Press, Stanford, 1975. [15] D.D. Joseph, A.C. Bannwart, Y.J. Liu, Stability of annular ow and slugging, Int. J. Multiphase Flow, in press.
Particle phase boundary layer theory in vertical two-phase gas–solid ow Y.A. Sergeev a; ∗ , P.L. Iske b , V.N. Kurdyumov c , J.H.J. Moors d a Department
of Engineering Mathematics, University of Newcastle, Stephenson Building, Newcastle upon Tyne NE1 7RU, UK b Knowledge Management, ABN-AMRO, Global Transaction Services, P.O. Box 263, 1000 EA Amsterdam, The Netherlands c Universidad Politecnica de Madrid, E.T.S.I. Aeronauticos, Plaza Cardenal Cisneros, 3, 28040, Madrid, Spain d Faculty of Electrical Engineering, Eindhoven University of Technology, P.O. Box 513, NL 5600 MB Eindhoven, The Netherlands Received 23 October 1997; received in revised form 5 January 1998
Abstract To gain understanding of the two-phase vertical gas–solid ow behaviour, the equations of gas– solid suspension ow are used to develop a model to predict the radial pro les of the volumetric fraction and the velocity of the solid particles. We derived two coupled di erential equations for the volume fraction of the particles and the uctuations of the particle velocity. These could be solved approximately, assuming the existence of a lean phase in the centre (the core) and a dense particle phase near the wall (the annulus). In this boundary layer theory, the method of matched asymptotic expansions could be employed. As a result, the obtained expressions have a generic character and the observed qualitative trends are insensitive with respect to the detailed c modelling of the particle–turbulence and particle–particle interactions. 1998 Elsevier Science B.V. All rights reserved. Keywords: Kinetic theory of gas–particulate suspensions; Multiphase ow; Two-phase gas– particulate ow
1. Introduction Dispersed ow phenomena are encountered in a wide variety of industrial applications. In two-phase dispersed gas–solid as well as gas–liquid pipe ow a so-called core–annulus behaviour can be observed. This means that the high-density component (i.e. the solid particles or the liquid drops) is found preferentially near the wall. This ∗
Corresponding author. Fax: 0191 222 5498; e-mail: [email protected].
c 1998 Elsevier Science B.V. All rights reserved 0378-4371/98/$19.00 Copyright PII S 0 3 7 8 - 4 3 7 1 ( 9 8 ) 0 0 0 6 2 - 4
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
27
a ects the ow pattern: the central lean phase (the core) owing faster than the annulus dense phase. Under certain conditions down ow near the wall can even occur. An accurate computational model, based on a correct description of the physical phenomena occurring in a vertical two-phase ow, is required to interpret and extrapolate the experimental data. However, to predict qualitatively the impact of a change in the geometrical and physical parameters on the ow behaviour, a less detailed model can be used. Therefore, we developed a model, which only contains the essential physics. The model focuses on the steady state fully developed ow. This case represents a one-dimensional problem. This problem has been addressed in many other publications. Sinclair and Jackson [1] discussed a model for fully developed ow in which the particles interact with the gas. The velocity pro le of the gas is modi ed by the presence of the particles, but remains parabolic, i.e. high velocity in the centre of the pipe and low velocity at the wall. The velocity pro le of the particles is induced by the drag forces on the particles. Because of the velocity gradients, uctuations in the particle velocity will be generated as particles in adjacent layers collide with each other. The random motion of the particles is associated with a quasi-temperature, in analogy with the kinetic theory of gases. It was shown by Pita and Sundaresan [2] that in the model of Sinclair and Jackson the inelasticity of the collisions causes an unrealistically strong e ect on the
ow phenomena. Other who used either the collisional approach or a combination of turbulence modelling and the collisional approach include Simonin [3], Louge et al. [4] and Ding and Gidaspow [6]. The work of the rst two authors so far, is restricted to (very) low particle volume fractions (i.e. smaller than 0.5%), whereas the latter authors concentrated on applying kinetic theory which was derived to describe rapid granular
ow. Based on the collisional theory of gas–particulate suspensions, the asymptotic methods have been successfully applied in the recent work by Sergeev and Zhurov [5] to develop the approximate analytical theory of vertical gas–particulate ow at low and moderate pressure gradients. Because of the over-sensitivity of these collisional models towards the restitution coecient of the particle–particle collisions, Dasgupta et al. [7] recently extended the model of Sinclair and Jackson in order to incorporate the interaction between the turbulent gas velocity eld and turbulent motion of the particles as source for the
uctuations. However, some aspects of the model remain speculative, e.g. the closure relations that are used are the same as used in the granular kinetic theory with some adjustable parameters. This means that in the modelling of the particle uctuations still collisions are held responsible, whereas it is very unlikely that the (direct) particle collisions are important. Based on the particle diameter, the Reynolds number is of the order one, and in this regime particles hardly collide. In order to account for the turbulence, sometimes a conventional Reynolds decomposition and time-averaging procedure is applied to the two-continua transport equations for mass, momentum and energy. This introduces various unknown correlations between the uctuations in the velocities, pressure and volume fractions of both phases, for which closure relations have to be speci ed. A generalisation of the one-phase k–
28
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
model to turbulent two-phase ows is often used to provide the required closure relations [ 7 –10]. In this approach, the turbulent properties of the uid phase follow from the transport equations for the turbulence intensity k1 and its dissipation rate 1 , which are solved in the model. The corresponding properties of the particulate phase are determined from suitably chosen correlation functions. The nature of the equations to be solved is such that usually severe numerical problems are encountered when solving them. The system of non-linear partial di erential equations form two-point boundaryvalue problem for which standard numerical recipes not always work. Furthermore, problems occur due to multiplicity of the solutions. It seems that the majority of the models predict the same qualitative behaviour of the gas–solid ow, including the segregation phenomenon. But the question is how sensitive these qualitative results are with respect to the detailed microscopic description, such as particle–particle interactions and particle–turbulence interactions. In this article we show that under very general assumptions the qualitative behaviour of two-phase gas– solid ow can be predicted. Existing models, such as the ones based on kinetic theory for granular ow or those based on detailed turbulence modelling, can all be shown to be realisations of the present model. Because of its relative simplicity the model can be used to predict trends in the two-phase ow behaviour. On the other hand, because there is no detailed microscopic description of the dispersed ow, the model yields qualitative results only. Some model parameters appear which have to be determined by a more fundamental study or by tting experimental data. To study the behaviour of relatively dense gas–solid systems we take the “Eulerian– Eulerian” approach in which the gas and the particle phase are treated as two interpenetrating uids, which are governed by their own separate transport equations valid at each point of the ow. These transport equations are derived from the microscopic conservation laws for mass, momentum and energy by an averaging procedure, for which di erent schemes have been proposed (see, e.g., [11]). The advantage of the continuum description is that it contains a limited number of independent variables. However, since the details of the microscopic ow-behaviour have been “averaged-out”, one has to postulate additional constitutive equations for the e ective stress tensors and interfacial transfer terms for mass, momentum and energy. To nd an approximate solution of the model equations, we base ourselves on the observation that in developed ow a lean region in the centre of the pipe (the core) and a dense region near the wall (the annulus) are found. Because the annulus is rather thin it seems appropriate to study the problem by applying techniques used in boundary-layer theory.
2. Description of the model In this section we will rst review the transport equations for the mean quantities used in the two-continua model. Then we will present the model assumptions and nally the equations will be reduced to their one dimensional steady-state form for the mixture.
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
29
2.1. Basic equations The basic transport equations solved in the two-continua model represent the conservation of mass and momentum for each of the individual phases separately. The variables solved from the mass and momentum conservation equations are the volume fractions g and s of the gas and solids phase, respectively and the velocities ug and us . We will summarise the relevant equations below: Mixture continuity: g + s = 1:
(1)
Gas phase continuity: @ g + ∇ · ( g ug ) = 0: @t
(2)
Solids phase continuity: @ s + ∇ · ( s us ) = 0: @t
(3)
Gas phase conservation of momentum: g
@ ( g ug ) + g ∇ · ( g ug ug ) = g ∇ · g + g g g + sg (us − ug ): @t
(4)
Solid phase conservation of momentum: @ ( s us ) + s ∇ · ( s us us ) = s ∇ · g + ∇ · s + s s g + sg (ug − us ): (5) @t In these equations, g and s are the gas and solid densities, respectively, g and s are the e ective stress tensors of the gas and solids phase, and sg is the viscous drag function between the two phases. Now we consider the steady plane-parallel vertical ow of gas–solid mixture in a vertical cylindrical pipe (i.e. no time dependence and no dependence on the z- and -coordinates). With the above assumptions the basic system of equations (1) – (5) reduces to the equations of conservation of the gas and solids axial momentum and the equation of conservation of radial momentum of solid particles. From here on the volumetric fraction of the solid particles will be denoted by . These equations are, respectively, s
1 d dP − (1 − ) (rg ) − (1 − )g g + sg (us − ug ) = 0; dz r dr rz 1 d 1 d dP g − (rrz (rs ) = 0; ) − s g + sg (ug − us ) − − dz r dr r dr rz d g d s rr = rr = 0; dr dr − (1 − )
(6) (7) (8)
where z is the vertical axis and P is the gas pressure. A further simpli cation can be obtained by the addition of Eqs. (6) and (7). The shear stresses in the gas phase are
30
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
neglected with respect to the shear stresses in the solids phase (which has been veri ed to be a good approximation), with as result the following axial momentum equation for the mixture: dP 1 d + [(1 − )g + s ]g + (rs ) = 0; dz r dr rz
(9)
which does not contain the viscous drag function sg . Inside the mixture the drag forces cancel. The radial momentum equation for the mixture implies d s = 0: dr rr
(10)
2.2. Constitutive equations Now we consider the closure for the system of equations (9) and (10). For that s and the bulk purpose we will write down constitutive equations for the shear stress rz s stress rr . We now introduce a uctuation energy (or “quasi-temperature”), , which is equal to one third of the mean square of the velocity uctuations of the particles: 1 = us0 us0 : 3
(11)
In this expression, us0 denotes the uctuation in the particle velocity. A remark should be made here. For the values of parameters relevant for these ows the value of the viscous drag coecient sg is very large, which means that the gas velocity ug di ers very little from the solids velocity us in the bulk ow, as well as in the ow close to the wall (we will take us = ug = u). The set of equations (9)–(10) is complemented with a conservation equation for the
uctuation energy in the solids phase. This equation is written as 1 d s du (rq) + rz + = 0; r dr dr
(12)
where q is the pseudo-heat ux and the rate of dissipation of the energy of random particulate motion which will be discussed separately. The interpretation of as a thermodynamic property of the system suggests that results, obtained in the kinetic theory for gases, can be employed to nd analogous expressions for the closures for our model. In this work, expressions for the stresses will be derived which can be related to the kinetic theory for dense gases [12]. However, we will show that in order to pursue our analysis, we do not necessarily need the detailed expressions for e.g. the viscosity and the heat conduction coecient, but that much more simple relations, showing the right asymptotic behaviour, are sucient. The constitutive relations for s and q are written as s rz = −
du ; dr
(13)
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
31
s rr = − C;
(14)
d : dr
(15)
q= −
In these equations denotes the apparent viscosity of the solids phase, C is a proportionality constant (similar to the ideal gas law), and is the quasi-temperature–heat conduction coecient in the solids phase. Inspired by the beforementioned analogy with the kinetic theory for gases we propose the following expressions for the transport coecients, appearing in the previous equations: √ (16) = f1 ( )s ( ) ; C = f2 ( )s ;
(17)
√
= f3 ( )s ( ) :
(18)
Eqs. (9), (10), (12) with the constitutive relations (13) – (15) represent the closed system for the unknown parameters , u and . Of course, explicit expressions for the functions f1 , f2 and f3 as well as for the “mean free path” (i.e. the distance over which the particles exchange their momentum and uctuation energy) have to be known, before getting quantitative results. For the analysis presented here, however, we only need the following properties of the functions f1 , f2 and f3 : lim f1 ( )( ) = o1 ;
(19)
lim
f2 ( ) = o2 ; →0
(20)
lim f3 ( )( ) = o3 ;
(21)
→0
→0
where o1 , o2 and o3 are positive and nite constants. It is interesting to note that, although there is no physical basis for the use of the collision theory in these systems (in which the gas plays an important role and the particles hardly collide), the expressions derived in the collision theory do satisfy the constraints (19)–(21). In particular, for f2 ( ) one derived f2 ( ) = (1 + 4 g0 );
g0 = (1 − ( = max )1=3 )−1 ;
(22)
in which max is the volume fraction corresponding to random close packing (which is estimated at max ' 0:65). Expression (22) can be observed to give correct limiting behaviour described by Eq. (20), with o2 = 1. In the kinetic theory for granular ow, the limiting values of 1 and 3 are nite: the mean free path diverges as −1 for small , but the functions f1 ( ) and f3 ( ) become linear in the volume fraction. We will now discuss the rate of dissipation of the uctuation energy , which appears in Eq. (12). We write
= !3=2
(23)
32
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
which is consistent with expressions for dissipation caused by rotations. It is interesting to notice that the same form for ! is obtained in the framework of the collisional theory by Lun et al. [13] where ! turns out to be a function of . Now we introduce dimensionless variables. The scale of length will be R, the radius of the pipe and velocities will be scaled with V, the super cial gas velocity, and the quasi-temperature will be scaled with V 2 . As a scale for the mass density we use s , the mass density of the solids phase. Expressed in the dimensionless variables (denoted by a tilde) the set of equations (9), (10), (12) reads 1 1d dP˜ + [(1 − ) + ] − d˜z F r˜ r˜
p du˜ ˜ rf ˜ 1 ( )( ) ˜ = 0; d˜r
d ˜ = 0; (f2 ( )) d˜r 1 d r˜ d˜r
(24) (25)
p d˜ ˜ rf ˜ 3 ( )( ) ˜ d˜r
!
p du˜ 2 ˜ + f1 ( )( ) ˜ + ˜ = 0: d˜r
(26)
In Eq. (24) the Froude number F appears, given by F = gR=V 2 . The density ratio = g =s is usually very small (about 10−3 ). Since ' 10−2 –10−1 it can be concluded that (1 − ) is small compared to and can be neglected. From now on we will skip the tilde in the notation for dimensionless variables; where confusion is possible we will specify whether we use dimensional or dimensionless variables. Now we consider the boundary conditions for the problem. The closed system of equations (24)– (26) will be considered in combination with the boundary conditions: d du = =0 dr dr u=0 +
at r = 0;
(27)
at r = 1; d =0 dr
(28) at r = 1:
(29)
The last boundary condition gives a balance between the ux (second term) and the rate of dissipation of uctuation energy at the wall ( rst term). Eqs. (24)–(26) can be conveniently rewritten in a form containing only two unknown functions, (r) and (r). Integrating Eq. (24) and using the boundary condition (27) yields for the velocity gradient: 1 1 1 dP du √ r = + dr f1 ( )( ) 2 dz rF
Zr 0
r 0 dr 0 :
(30)
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
33
Substituting Eq. (30) in Eq. (26) we obtain the system of equations for the functions (r) and (r) √ d 1 d rf3 ( )( ) r dr dr 2 Zr 1 dP 1 1 √ r + r 0 dr 0 + = 0; (31) + rF f1 ( )( ) 2 dz 0
d f2 ( ) = 0: dr
(32)
One more boundary condition is required. Such a condition can be, for example, an integral condition for the total ux of the solids phase. The dimensionless form of this condition is as follows: Z1 Qs = 2
us r dr;
(33)
0
where Qs is assumed to be given. Since the condition (33) is of the integral form, it causes certain inconveniences to nd a solution of the problem. Instead of (33) we shall use the supplementary boundary condition for in the form = ∗
at r = 1;
(34)
where ∗ is supposed to be known. The problem (21), (22), with the boundary conditions (27)–(29) and (34) is closed. As soon as the solution of this problem is obtained, Qs can be found from Eqs. (33) and (34) as a function of ∗ : Qs = Qs ( ∗ ).
3. Particle phase boundary layer theory 3.1. Introduction As mentioned before, the two-phase dispersed ow often shows a core-annulus behaviour. The rather thin annulus can be interpreted as an ‘extended’ boundary layer, therefore we are inspired to apply a boundary layer theory. To make this point clear, we consider the boundary condition (29). We will consider as the small parameter that is required to employ a boundary layer approach. (In the collision theory would have a value of the order of several tens (dimensionless) particle diameters, which indeed is rather small (about 10−4 –10−3 )). Due to the fact that the small parameter is in front of the derivative in the boundary condition, the formation of a boundary layer of the particle phase in the vicinity of the wall (i.e. as r → 1) must be expected.
34
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
An estimate of the thickness of the boundary layer as a function of ∗ = |r = 1 can already be given: from Eq. (32) we nd for r → 1: 0 ∗ d ∗ ∗ d = 0; (35) f2 ( ) + f2 ( ) dr r=1 dr r=1 1 f2 ( ∗ ) d ; (36) = dr r=1 f20 ( ∗ ) ≡
∗ f0 ( ∗ ) = ∗ 2 ∗ = O( ) (d =dr)r=1 f2 ( )
(37)
where (0 ) ≡ d=d . So, the thickness of the boundary layer is of the same order as . For reasons of convenience we once more give the set of equations with the boundary conditions from which the radial pro les of the particle uctuations and volume fraction can be obtained: √ d 1 d rf3 ( )( ) r dr dr 2 Zr 1 dP 1 1 √ r + r 0 dr 0 + = 0; (38) + rF f1 ( )( ) 2 dz 0
d f2 ( ) = 0 ; dr d d = =0 dr dr = ∗ +
(39) at r = 0;
at r = 1;
d =0 dr
at r = 1:
(40) (41) (42)
Due to the non-linearity, it turns out that it is dicult to analyse this system of equations numerically. In particular, the in uence of physical and geometrical parameters is hard to determine. Therefore, we prefer to formulate an approach that gives an approximate solution, but which allows us to predict trends in the ow behaviour. We will use the method of singular perturbations (see, e.g., [14]) in the form of matched asymptotic expansions for the problem (38) – (42). In accordance with this method we consider the following two domains: The boundary layer, i.e. the vicinity of the wall, 1 − r¡ = O( ). Below we denote the particle phase concentration and quasi-temperature in this region as l and l , respectively. The bulk, i.e. the central part of the ow, 06r¡1 − = 1 − O( ); the particle concentration and the quasi-temperature are written as b and b , respectively. In each of these regions the solution is sought in the form of asymptotic expansions of the functions (r) and (r) in terms of the parameter . We will seek the principal
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
35
terms of such expansions, only keeping the terms in lowest order in . In our analysis, rst we derive an equation for the quasi-temperature . In this paper we assume that the pressure drop is larger than required to compensate the gravitational forces on the solids in the mixture. This is the situation that is most studied in the literature [1,2,6]. We assume dP = O( 0 ): (43) dz 3.2. Quasi-temperature in the boundary layer We start our analysis by considering ow in the boundary layer. Here we make the following a priori assumption ∗ = O( m );
(44)
where m¿0. The proof of this assumption as well as the value of m will be given in the process of solving the problem (38)–(42). Since 061 − r6O( ) within the boundary layer, we introduce a new variable: y=
1−r :
(45)
This variable is found in the interval [0; −1 ]; within the boundary layer we have 06y6O( 0 ). Now we consider Eq. (38). Since ¡1, the integral is of the order of unity (i.e. O( 0 )) or smaller. The dissipation is supposed not to be larger than the second term. Keeping only terms of the lowest order in , only the rst term in Eq. (38) remains, so we can write p dl d f3 ( l )( l ) l = 0: (46) dy dy So, we obtain ( l )f3 ( l )
d3=2 l = C1 : dy
(47)
The boundary conditions at the wall read l = ∗ l =
dl dy
at y = 0; at y = 0:
(48) (49)
From d f2 ( l )l = 0 dy
(50)
it obviously follows f23=2 ( l )3=2 l = C2 :
(51)
36
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
Using Eqs. (47) and (51) and the boundary conditions (48) and (49) we obtain: C2 =
2 f23=2 ( ∗ ) C1 : 3 f3 ( ∗ )( ∗ )
(52)
It follows that the radial distribution of the quasi-temperature in the annulus can be written in the form: 2=3 f2 ( ∗ ) 2 C12=3 : (53) l = 2=3 ∗ 2=3 ∗ 3 f2 ( l )f3 ( ) ( ) The constant C1 is still unknown and must be obtained by asymptotic matching of the outer solution, b , and the boundary layer solution, l : b (r)|r→1 ∝ l (y)|y→∞ ;
(54)
where the symbol ∝ is used to indicate the coincidence of leading asymptotic terms with respect to . In order to give an asymptotic expression for the quasi-temperature we rst study the distribution of the solids fraction in the boundary layer. Combining Eqs. (47) and (51) leaves us with ( l )
d l + q( ∗ ) = 0; dy
(55)
where ( ) =
3 f3 ( )f20 ( )( ) 2 f25=2 ( )
(56)
q( ) =
3 f3 ( )( ) : 2 f23=2 ( )
(57)
and
From Eq. (55) it follows directly that the solids volume fraction in the boundary layer does not depend on any parameter except for the solids concentration at the wall ∗ ! From the properties of the function ( ) it also follows that limy→∞ l = 0. Therefore, to nd the asymptotic behaviour of the particle concentration at the outer boundary of the boundary layer (as y → ∞) and, consequently, the order of magnitude of the solids volume fraction in the bulk, the properties (20) and (21) can be used so that in Eq. (56) f3 ( )( ) and f2 ( ) can be replaced by f3 (0)(0) = 3
and
f2 ( ) = 2 ;
(58)
respectively, where 3 and 2 are the dimensionless quantities corresponding to the parameters o3 and o2 (note that 2 = o2 ). Then the boundary layer equation (55) takes the form 3=2 2 d l = − q( ∗ ) 2 5=2 : dy 3 3 l
(59)
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
37
From Eq. (59) the asymptotic behaviour of the solids volume fraction at the conventional boundary between the bulk and the boundary layer follows in the form l ∼
2=3
3 1 1 2=3 ∗ q ( ) 2 y2=3
as y → ∞:
(60)
Here we arrive at the important conclusion that the bulk is indeed much more dilute than the boundary layer; the order of magnitude of solids volume fraction in the bulk: b = O( 2=3 ):
(61)
We can also give the asymptotic expression for the quasi-temperature in the boundary layer. From Eq. (47) it becomes clear that 2=3 C1 y2=3 for y → ∞: (62) l (y) = 3 3.3. Quasi-temperature in the bulk Now we consider the solution in the bulk region. From the estimate (61) and assumption (43) it follows that the integral term in Eq. (38) can be neglected compared to the gas pressure gradient term. Since the solid volume fraction in the bulk is low the combinations f3 ( b )( b ) and f1 ( b )( b ) in Eq. (38) can be replaced by the (non-zero!) zero order terms in their Taylor expansions, i.e. by 3 and 1 respectively, where 1 = f1 (0)(0) is the dimensionless quantity corresponding to the parameter can now be written as p (dP=dz)2 r 2
1 d db √ r b + + = 0: r dr dr 4 1 3 3 b
(63) o1
in Eq. (19). Eq. (38) (64)
Eq. (64) has to be solved subject to the condition of symmetry (40) and the condition of asymptotic matching with the boundary layer solution (54). To nd the unknown constant C1 , that appeared in the solution for the quasi-temperature in the boundary layer, we require that the solution b for the bulk region of Eq. (64) must have the asymptotic behaviour of the form (62) as r → 1. This means that b → 0 as r → 1, so that instead of Eq. (64), with the boundary conditions (40) and (42) formulated within the interval 06r61 − O( ), we can consider Eq. (64) within the interval 06r61, together with the boundary condition (40) and b = 0
at r = 1:
(65)
We can write Eq. (64) in a convenient form by introducing a new variable , de ned by −1 dP √ (66) = 2 1 3 b : dz
38
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
The resulting equation reads: r2 1 d p d r + √ + 0 = 0 r dr dr
(67)
with 0 the modi ed uctuation energy dissipation rate in terms of . The boundary conditions for Eq. (67) are d = 0; (68) dr r=0 |r=1 = 0:
(69)
Eq. (67) can be slightly rewritten: 2 p d 1 d d r + + r 2 + 0 = 0: r dr dr 2 dr
(70)
Note that if we neglect the dissipation of uctuation energy, the solution of Eq. (70) yields a universal radial pro le of , which does not depend on any parameter! The asymptotic analysis of Eq. (70) without dissipation with the boundary conditions (68) and (69) yields: = A(1 − r)2=3
as r → 1:
(71)
This can be validated by substituting Eq. (71) into Eq. (70), observing that the divergent terms cancel. The conclusion is that the behaviour of the particle uctuations in the bulk near the wall is consistent with the behaviour of the particle uctuations in the boundary layer, given by Eq. (62). The method of matched asymptotic expansions can now be applied, since A is a xed constant that can be determined numerically. In case of dissipation, the constant A will depend on the rate of dissipation. If we write (compare to Eq. (23)):
0 = !3=2
(72)
we can obtain the function A(!) (A(0) ' 0:969) using a transformation of the radius: r → 1 − rˆ3=2 . The form (72) is consistent with formulae used in the kinetic theory of gases. In Fig. 1 we plotted the function A(!). In Fig. 2 we plotted the solution of Eq. (70) for several values of !. As follows from the results represented in Fig. 2, in absence of the dissipation the radial quasi-temperature pro le is monotonic. The dissipation leads to the appearance of maximum of the particle uctuations at a certain distance from the wall. 3.4. Matched asymptotic solutions We are now able to nd the quasi-temperature in the boundary layer, because the constant C1 in Eq. (60) can now be determined. Equating (62) and (71) and using the
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
39
Fig. 1. The matching factor A(!).
Fig. 2. Quasi-temperature pro les for di erent values of the dissipation factor !.
scaling (66) yields: C1 =
A3=2 3=4 1=4 3 (4 1 )
3=2 dP : dz
(73)
For the uctuation energy in the wall region we obtain from Eq. (53) l =
2=3 2=3 f2 ( ∗ ) 1 3 2 A(!) 2=3 ∗ ∗ 3 f3 ( )( ) f2 ( l ) (4 1 3 )1=2
dP : dz
(74)
Now we come to the a priori assumption (44): we recall that ∗ = O( m ). From Eq. (74) we nd: ∗ =
2=3 2=3 1 3 2 A(!) 2=3 3 f3 ( ∗ )( ∗ ) (4 1 3 )1=2
dP : dz
(75)
It follows then that m = 2=3:
(76)
40
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
It is now possible to derive a relation between the volume fraction of the particles in the bulk, b , and the volume fraction at the wall, ∗ . For the bulk we use the solution of Eq. (70), with the result |dP=dz| (r; !): b (r) = √ 4 1 3 Now we use Eq. (39) yielding: 2=3 2=3 3 2 2=3 A(!)f2 ( ∗ ) = f2 ( b (r))(r; !): 3 f3 ( ∗ )( ∗ )
(77)
(78)
Since the leading expansion term of the solids concentration is of the order 2=3 , using Eq. (78) and properties (20) and (21) to nd the radial distribution of in this region we obtain 2=3 2=3 A(!) 3 f2 ( ∗ ) 2 2=3 : (79) b (r) = 3 (r; !) 2 (f3 ( ∗ )( ∗ ))2=3 It might be useful to estimate the solids volume fraction at the centre of the pipe, b (0), as a function of the dissipation rate, !. From Eq. (79) it follows: 2=3 2=3 A(!) 3 f2 ( ∗ ) 2 2=3 : (80) b (0; !) = 3 (0; !) 2 (f3 ( ∗ )( ∗ ))2=3 To conclude this section we note that the following approximate form for the solids quasi-temperature, tot , valid within the whole interval 06r61 can be suggested. The function tot should satisfy tot = b
as r → 0;
(81)
tot = l
as r → 1:
(82)
We choose:
"
(1−r)= # (1−r)= 1 1 l : b + tot = 1 − 2 2
(83)
It should be noted that tot is not a solution to the equations. However, from this function, together with the boundary condition (41) for the particle volume fraction at the wall, we can calculate tot , the radial distribution of the solids volume fraction, using Eq. (39). With the functions tot and tot , we are able to determine the radial velocity pro le by integration of Eq. (30). 4. Velocities Below, the radial distribution of the velocity of the solid particles in the annulus and the bulk region will be obtained; since ub ' ug the formulae give simultaneously the
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
41
pro le of the gas velocity. The (solids) velocity can be obtained from Eq. (30), with the boundary condition at the wall (28). Since = O( 2=3 ) in the bulk region, and the thickness of the boundary layer in which = O( 0 ) is of the order , the integral term in Eq. (30) can be neglected for all r. Then we get the simpli ed expression: r dP=dz du = √ : dr 2 f1 ( )( )
(84)
By integrating this formula we obtain the velocity. In the boundary layer we obtain 1 dP ul (r) = − 2 dz
Z1 √ r
r0 dr 0 ; l f1 ( l )( l )
(85)
where we used u(1) = 0. The length of the integration interval is of the order and the integrand is O( −1=3 ). The conclusion is that the velocity in the boundary layer is O( 2=3 ). In the bulk the velocity is given by 1 dP ub (r) = − 2 dz 1 dP − 2 dz
1− Z
√
r0 dr 0 b f1 ( b )( b )
√
r0 dr 0 : l f1 ( l )( l )
r
Z1 1−
(86)
The rst integral in the RHS is O( 0 ). We have already seen that the second integral is O( 2=3 ), so it can be neglected in Eq. (86). Because = O( ) we can replace the upper limit 1 − in the integral by 1. Using the variable , introduced in Eq. (66), and replacing b in the combinations f3 ( b )( b ) and f1 ( b )( b ) by 3 and 1 respectively, we obtain 1=4 ub (r) = √ 3 3=4 2 1
s Z1 0 dP p r dr 0 : dz (r 0 ; !)
(87)
r
From this formula it can be concluded that the velocity in the bulk is proportional to the square root of the pressure gradient. Using the obtained radial distributions of solids velocity (87) and volume fraction (79) the total dimensionless solids ux given by Eq. (33) can be calculated as a function of the particle concentration at the wall ∗ . This gives ∗ = ∗ (Qs ) and, therefore, solves the problem in its original formulation.
(88)
42
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
To calculate Qs we have a contribution from the bulk and from the boundary layer: 1− Z Z1 ub b r dr + 2 ul l r dr: Qs = Qs; b + Qs; l = 2 0
(89)
1−
It is not dicult to demonstrate that the second integral can be neglected with respect to the rst one, which means that the solids ux is determined by the ux in the bulk only Z1 Qs ' Qs; b ' 2
ub b r dr:
(90)
0
From Eqs. (79) and (87) we obtain: 1 s Z1 Z 11=12 7=6 2=3 ∗ 0 2 A 3 f2 ( ) r r dP p dr 0 dr: Qs = dz (r) 3=4 0 2=3 ∗ ∗ 2=3 (r ) 3 1 (f3 ( )( )) 0
(91)
r
The double integral can be numerically evaluated (it depends on the dissipation factor !). We can write Qs in a convenient form s dP f2 ( ∗ ) 2=3 : (92) Qs = s (!) dz ∗ ∗ 2=3 (f3 ( )( )) Finally, we can calculate the dimensionless gas ux, which is given by Z1 Qg = 2
Z1 u(1 − )r dr ' 2
0
Z1 ub (1 − b )r dr ' 2
0
ur dr;
(93)
0
where we only took the contribution of the lowest order in ; in Eq. (93) is the density ratio; the velocity is given by Eq. (87) and we nd s s Z1 Z1 0 dP 2 1=4 dP r 3 r p dr 0 dr = g (!) : (94) Qg = dz 3=4 0 dz (r ) 1 0
r
We note that the particle volume fraction at the wall can be explicitly expressed through the total solids ux in case of relatively dilute two-phase suspension. Assuming low particle concentration throughout the entire pipe including the wall region we may write in Eq. (91) f2 ( ∗ ) ' 2 ∗ ;
f3 ( ∗ )( ∗ ) ' 3 :
(95)
Neglecting also the dissipation, the numerical value of A ' 0:969 and the integral in Eq. (91) can be estimated as 1 Z1 Z r r0 p (96) dr 0 dr ' 1:219 (r; 0) (r 0 ; 0) 0
r
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
43
so that for the particle concentration at the wall we nd: ∗ ' 0:25
3=4 1
2=3 2 1=4 3
Qs :
(97)
We complete this section with a remark on the shear stress and stability of the fully developed vertical ow of suspension of solid particles. The shear stress in the solids phase is given by the constitutive relations (13) and (16). Using the solution for the bulk (77), (79) and (87) it can be easily shown that dr z ¿0: (98) dr r=0 Using the boundary layer solution of Eq. (55) and the radial distributions of the solids quasi-temperature and velocity in the boundary layer, given by Eqs. (74) and (85), respectively, together with the value of the quasi-temperature at the wall (75) and the boundary condition for the solids velocity (28) the leading asymptotic term of the radial gradient of the e ective shear stress in the solids phase can be found in the form 1 dr z = −B( ∗ ) 2=3 ¡0; (99) dr r=1 where B = O(1) is the positive function of the solids volume fraction at the wall. From inequalities (98) and (99) it follows that the e ective shear stress in the solids phase reaches the maximum within the ow region. As shown by Joseph et al. [15] this might be evidence of the stability of the annular structure of the considered gas-particle ow in the vertical pipe. We also calculate here the shear stress at the wall. We have seen that for the velocity gradient at the wall one can write dP=dz du : (100) = √ dr r=1 2 ∗ f1 ( ∗ )( ∗ ) From the constitutive relations (13) and (16) we obtain the following remarkably simple relation between the wall stress and the pressure drop: r z |r=1 = −
1 dP : 2 dz
(101)
In the next section we will discuss results obtained by applying the approximate method, described in the previous sections, and compare these results with results obtained by numerical solution of the full set of Eqs. (38) – (42). 5. Accuracy, qualitative trends and conclusions A boundary layer theory was developed in order to obtain more insight in the in uence of physical and geometrical parameters on the model results. We will rst
44
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
make a comparison between results of the numerical scheme and results obtained by the method of matched asymptotic expansions. Though the theory does not depend on the explicit form of the closure relations, a choice for speci c closure relations has to be made in order to perform the calculations. Because the closure relations from the kinetic theory for granular ow are most commonly used, the values for o1 , o2 and o3 , de ned in properties (19)–(21), are based on the expressions given by Lun et al. [13]. For the suspension of ne solid particles (the diameter dp ∼ 10−4 m) these values read o1 ' 9:23 × 10−6 m;
o2 = 1;
o3 ' 3:46 × 10−5 m:
(102)
We take these values for our calculations. It should be noted, however, that the predicted trends are not sensitive towards the actual choice of the closures. To calculate the dimensionless parameters we consider the fully developed ow of solid particles of the diameter given above and the density s = 1:4 × 103 kg=m3 in the vertical pipe of the radius R = 0:5 m. The scale velocity is assumed to be V = 10 m=s, the gas density g = 1:2 kg=m3 , and the particle volume fraction at the wall ∗ = 0:2. First, the accuracy of the approximate solution described in Sections 3 and 4 will be tested. The solution of the set of equations (31) – (32) for the particle uctuation energy, the particle volume fraction and Eq. (30) for the velocity of the particles are calculated with both the approximate and numerical methods. The rst comparison is made at a relatively high pressure drop, dP=dz = −3000 Pa=m (the corresponding ˜ z˜ = −0:0214). In this case the particle volume fracdimensionless pressure drop dP=d tion obtained from the approximate solution di ers at most 3% from the complete numerical solution. The uctuation energy and the velocity di er a few percent. In Fig. 3 the radial distributions of the quasi-temperature, the particle volume fraction and the velocity are given. The particle volume fraction at the wall is found to be much larger than the particle volume fraction in the middle of the pipe, which is a result that is characteristic for all calculations. This is related to the radial distribution of the uctuation energy. In the middle of the pipe the uctuations are larger than at the wall. In fully developed ow particles prefer being in regions where the uctuations are smaller (cf. the ideal gas) and so they segregate towards the wall. Radial distributions for the particle uctuation energy, the particle volume fraction and the velocity for decreasing pressure drops show the same characteristics. The velocity and the quasi-temperature are decreasing with decreasing pressure drop, as can be concluded by comparing the results for a pressure drop of 3000 Pa=m (Fig. 3a) and 500 Pa=m (b) ˜ z˜ = −0:00357). It must be noted that the accuracy of the boundary (dimensionless dP=d layer theory decreases with decreasing pressure drop; the relative error between the numerical calculations and the approximate theory is slightly higher (about 5%) in the case (b). The important trend to be studied is the in uence of the particle diameter. In Fig. 4 the radial distributions of the velocity and the particle uctuation energy (quasitemperature) are given for various particle diameters. The in uence is studied by changing the particle diameter and adjusting ∗ and the dimensional pressure gradient dP=dz
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
45
Fig. 3. Comparison between the approximate and the numerical solution. The dimensionless quasitemperature, ˜ = =V 2 , the particle volume fraction, , and the dimensionless velocity, u˜ = u=V . Dashed line – numerical solution of the complete set of equations; solid line – approximate asymptotic solution. = 0:0068, ∗ = 0:2. (a) dP=dz = −3000 Pa=m; (b) dP=dz = −500 Pa=m.
in such a way that the super cial velocities of the gas and solids, i.e. Usup; g =
s Qg ; R 2
Usup; s =
s Qs ; R 2
(103)
respectively, where Qg and Qs are the dimensionless gas and solid uxes given by Eqs. (93) and (33), respectively, remain the same (the values of the parameters are given in Table 1 for all the calculations). We assume that the densities of the particles and the gas are those given above, and the pipe radius is R = 0:5 m. The super cial velocity of the gas is used as a scale velocity V in the solution described in Sections 3 and 4, V = Usup; g . We can conclude that the velocity in the centre of the pipe is larger for small particles. The velocity pro le is thus more peaked for smaller particles. Smaller particles uctuate more than larger particles (Fig. 4b). The calculated pro le of the quasi-temperature can now be used to nd the particle volume fraction pro le
46
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
Fig. 4. The in uence of the particle diameter obtained with the boundary layer theory on the radial pro les of (a) velocity, (b) uctuation energy (quasi-temperature). (···) dp = 150 m; (- - -) dp = 100 m; (– – –) dp = 70 m; (——) dp = 30 m. Table 1 Values for dimensional parameters used in calculations to study the in uence of particle diameter dp (m)
∗
dP=dz (Pa=m)
Usup; g (m=s)
Usup; s (m=s)
30 70 100 150
0:236 0:136 0:100 0:0635
−141:4 −158:9 −176:45 −208:6
20:00 20:00 20:00 19:99
0:1465 0:1466 0:1468 0:1467
in the bulk and the boundary layer from the approximate solution given in Section 3. It is dicult to illustrate graphically the in uence of the particle diameter on the radial distribution of the volume fraction, so we rather describe the found trends qualitatively. The particle volume fraction in the centre of the pipe decreases when the particle diameter is increased from 30 to 70 to 100 m and increases again when the diameter is further increased to 150 m. The particle volume fraction at the wall has to be larger for small particles to ensure that the super cial velocity of the solids phase is a constant. The ratio of the particle volume fraction at the wall to the particle volume fraction in the centre of the pipe is larger for small particles. From the obtained results it follows that the segregation is stronger for smaller particles, and the velocity pro le
attens with decreasing particle diameter.
Y.A. Sergeev et al. / Physica A 255 (1998) 26–47
47
In this paper we discussed a model in which the two-phase conservation equations for fully developed ow are solved by introducing a quasi-temperature, which is a measure for the particle uctuations. Closure relations for the viscosity, heat conductivity and the pressure were expressed in terms of this quasi-temperature. The equations are generic, i.e. the form of the equations does not depend on the microscopic model for the transport coecients such as the viscosity and the heat conductivity. Based on the observation that a large solids volume fraction exists in the wall region and a small solids volume fraction exists in the centre of the pipe a boundary layer theory for the solids phase has been developed. In this theory asymptotic solutions for the quasi-temperature and the volume fraction of the particles could be derived. This enables us to perform a qualitative analysis of the in uence of physical and geometrical parameters on a character of distribution of particle concentration and velocity in a rapid vertical two-phase pipe ow. A comparison was made between the results from the approximate asymptotic theory and results obtained by a direct numerical solution of the full set of equations. It was concluded that for studying trends qualitatively, the approximate method gives reasonable results, none of which are contradictory to experimental observations. References [1] J.L. Sinclair, R. Jackson, Gas–particle ow in a vertical pipe with particle–particle interactions, A.I.Ch.E. J. 35 (1989) 1473. [2] J.A. Pita, S. Sundaresan, Gas–solid ow in vertical tubes, A.I.Ch.E. J. 37 (1991) 1009. [3] O. Simonin, Prediction of the dispersed phase turbulence in particle-laden jets, Proc. ASME 1991 Conf. on FED, vol. 121: Gas–Solid Flows, 1991, p. 197. [4] M. Louge, E. Masorakos, J.T. Jenkins, The role of particle collisions in pneumatic transport, J. Fluid Mech. 231 (1991) 345. [5] Y.A. Sergeev, A.I. Zhurov, Asymptotic theory of two-phase gas–solid ow through a vertical tube at moderate pressure gradient, Physica A 236 (1997) 243. [6] J. Ding, D. Gidaspow, A bubbling uidization model using kinetic theory of granular ow, A.I.Ch.E. J. 36 (1990) 523. [7] S. Dasgupta, R. Jackson, S. Sundaresan, Turbulent gas–particle ow in vertical risers, A.I.Ch.E. J. 40 (1994) 215. [8] S.E. Elgobashi, T.W. Abou-Arab, A two-equation turbulence model for two-phase ow, Phys. Fluids 26 (1983) 931. [9] M.A. Rizk, S.E. Elgobashi, A two-equation turbulence model for dispersed dilute con ned two-phase
ows, Int. J. Multiphase Flow 15 (1989) 119. [10] C.P. Chen, P.E. Wood, A turbulence closure model for dilute gas–particle ows, Can. J. Chem. Eng. 63 (1985) 349. [11] D.A. Drew, Mathematical modelling of two-phase ow, Ann. Rev. Fluid Mech. 15 (1983) 261. [12] S. Chapman, T.G. Cowling, The Mathematical Theory of Non-uniform Gases, 3rd ed., Cambridge University Press, Cambridge, 1970. [13] C.K.K. Lun, S.B. Savage, D.J. Je rey, N. Chepurniy, Kinetic theories for granular ow: inelastic particles in Couette ows and slightly inelastic particles in a general ow eld, J. Fluid Mech. 140 (1984) 223. [14] M.D. Van Dyke, Perturbation Methods in Fluid Mechanics, Parabolic Press, Stanford, 1975. [15] D.D. Joseph, A.C. Bannwart, Y.J. Liu, Stability of annular ow and slugging, Int. J. Multiphase Flow, in press.