Powder Technology 150 (2005) 42 – 55 www.elsevier.com/locate/powtec
Numerical modeling of gas–particle flow using a comprehensive kinetic theory with turbulence modulation C.K. Chana,*, Y.C. Guob, K.S. Lauc a
Department of Applied Mathematics, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong b Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China c Department of Applied Physics, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong Received 5 March 2003; received in revised form 15 September 2004; accepted 1 November 2004 Available online 25 January 2005
Abstract In most fluidized beds, both solids flux and gas phase Reynolds number are high and the flow are usually turbulent. It is therefore necessary to consider both the effects of particle–particle collisions and particle phase turbulence in any mathematical model for simulating gas–particle flows. A comprehensive model is developed in the present work in which a two-equation (k–e) turbulence model is used for calculating the gas phase. In addition, a transport equation of particle phase turbulent kinetic energy is proposed and used for modeling the particle phase turbulence (k p model). Similar to that of the single gas phase, effective viscosity of the particle phase is the sum of the laminar viscosity caused by particle–particle collisions described by kinetic theory and the turbulent viscosity caused by collections of particles described by the k p model. The proposed model is used to predict gas–particle flows in a vertical pipe. Results obtained using this model compare well with experimental data. D 2004 Elsevier B.V. All rights reserved. Keywords: Gas–particle flows; Kinetic theory; Turbulent kinetic energy; Circulating fluidized beds
1. Introduction Eulerian approach is a common method for describing gas particle flows in circulating fluidized beds. This approach has been reviewed comprehensively by Gidaspow [1] and Enwald et al. [2]. In Eulerian approach, numerical models can be classified into three categories according to how the effective viscosity of the particle phase is determined. The first is the particle phase constant viscosity model used by several authors such as Tsuo and Gidaspow [3], Kuipers et al. [4,5], Theologos and Markatos [6,7] and Nieuwland et al. [8]. The advantage of this model is its simplicity, but it does not take into account the effect of nonuniform characteristics of the particle phase viscosity. The second model is an extension of the gas phase k–e turbulence model to the gas–particle flows. Chen [9]
* Corresponding author. Tel.: +852 2766 6919; fax: +852 2362 9045. E-mail address:
[email protected] (C.K. Chan). 0032-5910/$ - see front matter D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.powtec.2004.11.011
proposed an algebraic model of particle turbulent viscosity, derived from the concept of particle tracking fluid. A similar model was also proposed by Elghobashi and Abou-Arab [10] in which the effect of particle phase on gas phase turbulence was considered by adding additional terms to the transport equations of gas phase turbulent kinetic energy and its dissipation rate. Zhou and Huang [11] derived a particle phase turbulent kinetic energy (k p) equation that was determined by convection, diffusion, production and interaction with gas phase turbulence. By coupling with the gas phase k–e turbulence model, Zhou and Huang developed a k–e–k p model for simulating gas–particle flows. As the above models do not include the mechanism of particle– particle collisions, they are generally restricted to dilute gas– particle flows. The third model is the particle kinetic theory model, which is based on kinetic theories of non-uniform dense gases as described by Chapman and Cowling [12]. It assumes that the random motion of particles is analogous to the thermal motion of molecules in a gas and the starting
C.K. Chan et al. / Powder Technology 150 (2005) 42–55
Nomenclature C d drag coefficient C pp empirical modeling constant C 1, C 2, C l empirical constants of gas phase k–e turbulence model d p particle diameter e restitution coefficient g i acceleration due to gravity in the ith direction g 0 radial distribution function G k production term of the gas phase turbulent kinetic energy equation G p source term of the particle effect on gas phase turbulent kinetic energy G s modulus of particle phase or mass flux k turbulent kinetic energy L p particle mean free path l turbulence scale M g molecular weight of gas p pressure P p particle phase pressure due to particle–particle collision r radial coordinates R radius of riser R g universal gas constant Re Reynolds number T gas phase temperature u, v axial and transverse velocities Y v velocity vector
point is the Boltzmann integral–differential equation for a velocity distribution of particles. Jenkins and Savage [13] derived a conservation equation for the so-called granular temperature based on a Maxwellian velocity distribution for a single particle with the granular temperature proportional to the mean square of the random particle velocity. For the gas phase, the thermal motion and molecular collisions give rise to gas phase laminar viscosity and gas phase pressure. The random motion and interactions between individual particles give rise to the particle phase viscosity and particle phase pressure, which can be regarded as laminar characteristics of the particle phase and depend strongly on the granular temperature. Sinclair and Jackson [14] developed a model for fully developed gas–solid flow in a vertical pipe based on the particle kinetic theory model in which the interaction of the particles and the gas is restricted to a mutual drag force. Their computational results demonstrated that the interaction between individual particles produces lateral segregation of solids in the radial direction. Pita and Sundaresan [15] extended this work and investigated the fully developed flow of gas–particle suspensions in vertical pipes of different diameters. It was shown that this model captures
43
Greek letters b gas–particle flow drag coefficient c dissipation due to inelastic particle collision C H transport coefficient of particle pseudo-temperature e gas phase turbulent kinetic energy dissipation rate eg gas phase volume fraction ep particle phase volume fraction e p,max maximum value of particle phase volume fraction H particle pseudo-temperature g diffusion coefficient l viscosity m kinematic viscosity n bulk viscosity u generalized dependent variable P i,k particle phase stress tensor q density re turbulent Schmidt number for turbulent kinetic energy dissipation rate s i,j gas phase stress tensor s rp particle mean motion relaxation time xp particle phase turbulent kinetic energy dissipation rate Subscripts g gas phase i, j, k coordinate directions in inlet l laminar p particle phase T turbulence
the existence of steady state multiplicity when different pressure gradients can be obtained for the same gas and solid fluxes. However, the model manifests an unsatisfactory degree of sensitivity to the elasticity of the particle– particle collisions. Pita and Sundaresan [16] also studied numerically the steady developing flow of gas–particle suspensions in a vertical riser. Results show that the inlet granular temperature has remarkable influence on the flow pattern. The core–annulus flow can be predicted only when the inlet granular temperature is high and the rate of segregation is directly proportional to the granular temperature. The above model developed by Sinclair and Jackson [14] gives reasonable results when elastic particle–particle collisions are considered and reveals remarkable rich varieties of behavior in a vertical pipe over a wide range of flow conditions. Presently, the kinetic theory approach seems to be the most promising mathematical approach to describe the behaviour of particle phase flow in gas–particle flow systems [17]. Based on particle phase kinetic theory model, several investigators such as Ding and Gidaspow [18], Nieuwland et al. [19], Samuelsberg and Hjertager [20], Benyahia et al. [21], Detamore et al. [22] and Louge
44
C.K. Chan et al. / Powder Technology 150 (2005) 42–55
et al. [23] have carried out numerical simulations on bubbling fluidization, circulating fluidization and pneumatic transport. The kinetic theory model has also been modified by many researchers such as Lun [24], Grace and Sun [25], Kim and Arastoopour [26] to account for the slightly elastic and rough particles, non-uniformly sized particles and cohesive particles. Gidaspow [1] extended the model to both dense and dilute gas–particle flows starting with the Boltzmann equation when a nonMaxwellian velocity distribution was employed. As all the above formulations were based on particle–particle collisions, they can be classified as laminar particle phase model. Dasgupta et al. [27] extended the laminar particle phase model of Sinclair and Jackson [14] to a turbulent particle phase model. They proposed a particle phase turbulent kinetic energy equation for a fully developed gas–particle flow in a vertical pipe and a k–e model was used for gas phase turbulence. However, the effects of particle phase turbulence on the granular temperature, assumed to be constant, were neglected. Thus, these models were focused on the turbulent mechanism of particle fluctuations rather than the mechanism of particle–particle collisions as described by the kinetic theory. Hrenya and Sinclair [28] extended the model proposed by Dasgupta et al. [27] to account for two mechanisms that give rise to the lateral segregation of particles, i.e. interactions associated with individual particles based on kinetic theory and interactions associated with collections of particles based on the analogy with gas phase turbulent flows. The predicted results showed that if the effects of particle phase turbulence on the governing equations alone are considered, the model still manifests an excessive sensitivity to elastic particle–particle collisions. By considering the effects of particle phase turbulence on the kinetic theory, the expected segregation patterns can be predicted and does not have any unsatisfactory degree of sensitivity to elastic particle–particle collisions. In addition, Zheng et al. [29] used a k–e–k p–e p–H model to study gas–particle flow in a riser. However, the turbulent kinetic energy dissipation rate for the particle phase was not provided. Zhang and Reese [30] incorporated the kinetic theory into a single equation model to study the turbulence modulation of particle phase. In their paper, two kinds of turbulent kinetic energy equations of the gas phase, as developed by Crowe and Gillandt [31] and Chen and Wood [32], were investigated and shown to predict the turbulence modulation due to solid particles well. However, these studies were on dilute gas–solid flow with low mass loading ratio of below 2%. In this paper, a comprehensive mathematical model for developing gas–particle flows is proposed based on the particle phase turbulence model described by particle turbulent kinetic energy. The granular temperature (pseudotemperature) and particle phase turbulent kinetic energy are solved according to their respective transport equations. The
governing equations are time-averaged for both the gas phase and the particle phase. Effective viscosity of the particle phase is determined by the laminar and turbulent viscosities. The laminar viscosity due to particle–particle collisions is determined by kinetic theory and the turbulent viscosity due to collections of particles is described by the k p model. Flows in risers are very turbulent with high Reynolds numbers and the traditional kinetic theory generally does not take into account of time-averaged turbulence behaviour. The present model fills in this gap and is capable of taking full account of the different forces arising from the interactions between the mean and fluctuating velocity of both the gas phase and particle phase. These turbulence modulation include 1.
drag force from interactions between mean velocities of gas phase and particle phase; 2. gas phase Reynolds stress from the interactions between gas phase velocity fluctuations; 3. particle phase Reynolds stress due to interactions between particles velocity fluctuations; 4. change of turbulent kinetic energy of both the gas phase and particle phase due to interactions between the gas phase and particle phase velocity fluctuations; 5. particle phase pressure and shear stress from interactions between individual particles. In order to verify the present model, experimental results by Bader et al. [33] are used to compare with our predictions.
2. Governing equations for gas phase As given by Gidaspow [1] and Samuelsberg and Hjertager [20] and using a k–e turbulence model, the time-averaged governing equations together with the constitutive relations for gas phase are given as follows. Gas phase continuity equation B ag qg mj ¼ 0; Bxj
ð1Þ
where a and q represent the volume fraction and density of gas phase and subscripts g refers to the solid phase. Gas phase momentum equation Bsi;j B Bp ag qg mi mj ¼ ag þ þ b mpi mi þ ag qg gi ; Bxj Bxi Bxj ð2Þ where s i,j is the gas phase shear stress and is related to the Bmj i gradients of velocity components by si;j ¼ lg;e ð Bm Bxj þ Bxi Þ Bm 2 k 3 di;j Bxk ; lg;e ¼ lg;l þ lg;T , is the effective viscosity, l g,T = C A a g q g k 2/e is the gas phase turbulence viscosity determined by the k–e turbulence model and l g,l is the gas phase laminar viscosity.
C.K. Chan et al. / Powder Technology 150 (2005) 42–55 Table 1 Empirical constants for gas phase turbulence model
Particle phase momentum equation
C1
C2
Cl
rk
re
1.44
1.92
0.09
1.0
1.3
Gas phase turbulent kinetic energy equation B B lg;e Bk a q mj k ¼ þ Gk ag qg e þ Gp ; Bxj g g Bxj rk Bxj
B B ap qp mpi þ ap qp mpi mpj Bt Bxj ¼ ap
ð3Þ
i where G k is the production term given by Gk ¼ lg;T ð Bm Bxj þ Bmj Bmi Bxi Þ Bxj : The last term on the right-hand side of Eq. (3) is the source term due to interaction between the gas phase and the particle phase. This source term is caused by fluctuation correlation of the gas phase velocity and particle phase velocity through momentum transfer between the gas and particle phases. As described by Zhou [34], this term in dilute gas– pffiffiffiffiffiffiffi 2a q particle flow is model as Gp ¼ sprp p Cpg kkp k , where k p is the particle turbulent kinetic energy and s rp is the particle phase relaxation time. However, our present study considers the single particle fluctuation described by the kinetic theory and the source term dueq toffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the presence of 2ap qp particles is modeled as Gp ¼ srp Cpg k kp þ H k ; where H is the particle phase pseudo-temperature, C pg is an empirical modeling constant, and s rp is given dp2 qp Re0:6667 1 1 þ p6 by srp ¼ 18l . Hence, with increasing parg;l ticle diameter d p, the particle relaxation time increases and leads to a decrease in G p, which is also due to the fact that larger particle size leads to increase in the velocity slip between the gas phase and particle phase. The particle phase turbulent kinetic energy and pseudo temperature are considered as turbulence modulation of the particle phase on the gas phase. Gas phase turbulent kinetic energy dissipation rate equation lg;e Be B B a q mj e ¼ Bxj g g Bxj re Bxj
e þ ð4Þ C1 ðGk þ Gp Þ C2 ag qg e ; k
where the empirical constants associated with the k–e turbulence model as given in Table 1 are based on Reynolds and Cebeci [35].
3. Particle phase kinetic theory Using the kinetic theory as developed by Gidaspow [1], the governing equations for the laminar particle phase are given as follows. Particle phase continuity equation B a ap qp þ ap qp mpj ¼ 0: Bt Bxj
45
ð5Þ
BPi;j Bp þ b mpi mi þ ap qp gi ; Bxi Bxj
ð6Þ
Bm Bm Bm where Pi;j ¼ Pp di;j þ fp di;j Bxpjj þ lp;l ð Bxpij þ Bxpji Þ 23 Bmpj di;j Bxj is the particle phase shear stress, P p=a pq p [1+2(1+e)a p g 0]H is the pressure due to particle–particle collision, e is the coefficient of restitution and g 0 is a radial distribution 1function which can be expressed as g0 ¼ ap 3 3 1 5 ½1 ð ap;max Þ : In the above expression of P i,j , l p,l and f p are particle phase laminar viscosity and bulk viscosity respectively such that pffiffiffiffiffiffiffi
2 10qp dp pH 4 lp;l ¼ 1 þ ð1 þ eÞg0 ap 5 96ð1 þ eÞg0 rffiffiffiffiffi 4 2 H þ ap qp dp g0 ð1 þ eÞ ; 5 p and rffiffiffiffiffi 4 2 H : fp ¼ ap qp dp g0 ð1 þ eÞ 3 p In Eq. (2) and (5), b is the drag coefficient of gas–particle flow, according to the volume fraction of particle phase and is defined as 8 a2p lg;e ap qg jY mg Y mp j > > < 150 ag ðdp /p Þ2 þ 1:75 dp /p ; ap b0:8 b¼ 3 ap qg jY mp j 2:65 mg Y > > ag ; ap z0:8 : Cd 4 dp / p C d is the drag coefficient between the gas and a single particle, defined as (
0:687 24 1 þ 0:15Re ; Rep V 1000 p Cd ¼ Rep 0:44; Rep N1000 with the particle Reynolds number given by Rep ¼ j mYg mYp j ag qg dp /p lg;l
.
Particle phase pseudo-temperature equation B B ap qp H þ ap qp mpj H Bt Bxj Bmpj 2 B 2 BH 2 þ CH c; ¼ Pi;j 3 3 Bxj Bxj 3 Bxi
ð7Þ
qffiffiffiffi Bm where c ¼ 3ð1 e2 Þa2p qp g0 H d4p Hp Bxpkk is the dissipation term of pseudo-temperature pffiffiffiffiffidue to elastic collisions 2 150q dp pH 6 between particles. CH ¼ 384ðp1þe Þg0 1 þ 5 ð1 þ eÞg0 ap þ
46
C.K. Chan et al. / Powder Technology 150 (2005) 42–55
qffiffiffiffi
2a2p qp dp g0 ð1 þ eÞ Hp is the transport coefficient of pseudo-temperature.
4. Particle phase turbulence model Similar to the method suggested by Khalil [36] and based on the particle phase momentum equation obtained from the kinetic theory, taking ap ¼ ap þ a Vp ; mpi ¼ mpi þ mVpi ; mpj ¼ mpj þ mVpj and substituting them into the instantaneous equations for m pi and m pj , the following timeaveraged momentum equations for the particle phase can be obtained, by adopting the symbol bPQ for mean quantities, as
B B ap qp mpi þ ap qp mpk mpi Bt Bxk Bp BPi;k ¼ ap þ þ bðmi mpi Þ þ ap qp gi Bxi Bxk B qp mpk aVp mVpi þ qp mpi aVp mVpk þ qp ap mVpi mVpk Bxk B BpV ; ð8Þ q aVp mVpi aVp þ qp aVp mVpi mVpk Bt p Bxi
B B ap qp mpj þ ap qp mpk mpj Bt Bxk BPj;k Bp ¼ ap þ þ b mj mpj þ ap qp gj Bxj Bxk B qp mpk aVp mVpj þ qp mpj aVp mVpk þ qp ap mVpj mVpk Bxk B BpV qp aVp mVpj aVp þ qp aVp mVpj mVpk : ð9Þ Bt Bxj
In order to close Eq. (8) and (9), it is necessary to simulate the unknown second-order and third-order correlations in terms of lower-order correlations or mean-flow properties. As suggested by Zhou [34], one approach is to use particle phase turbulent viscosity model similar to that of the gas phase and to derive the transport equation of particle phase turbulent kinetic energy k p. To obtain the transport equation for mVpi mVpj , it is necessary to multiply the transport equation for mpi by mpj and the transport equation for mpj by mpi . The conservation equation for mpi mpj may be obtained. Similarly, multiplying Eq. (8) by mpj and multiplying Eq. (9) by mpi , the conservation equation for m——– pi m pj may be obtained. Lastly, the transport equation of mVpi mVpj may be obtained by taking the difference of the time-averaged conservation equation for m pi m pj and the conservation
mpj . The exact k p equation can be written as equation for P mpi P B B ap q p k p þ ap qp mpk kp Bt Bxk qp ap mVpk mVpi mVpi Bkp B 4 lp þ fp ¼ mVpi pVdi;k Bxk 3 Bxk 2 BmVpi BmVpi 4 Bp l þ fp eVp mVpi 3 p Bxi Bxk Bxk BpV BqV ap mVpi þ aVp mVpi þ aVp mVpi qp gi Bxi Bxi þ b mVi mVpi kp qp mpk aVp mVpi þ qp mpi aVp mVpk B þ qp ap mVpi mVpk þ qp aVp mVpi mVpk mpi Bxk Bqp P mpk aVp mVpi mVpi þ aVp mVpk mVpi mVpi 2Bxk B P mpi qp aVp mVpi Bt aq Bqp
P P p þ aVp mVpk mpi mpj aVp mVpi mVpi ð10Þ 2Bxk 2at It is then necessary to close the correlation terms such as aVp mVpi and aVp mVpj . The simplest method is to use gradient modeling. The expressions for aVp mVpi and aVp mVpj can be expressed as aVp mVpi ¼
mp;T Bap ; rp Bxi
ð11Þ
aVp mVpj ¼
mp;T Bap ; rp Bxj
ð12Þ
where m p,T is the turbulent kinematic viscosity. The turbulent stress production term can be modeled as qp ap mVpi mVpk
Bmpi Bmpk Bmpk Bmpi ¼ lp;T þ Bxk Bxi Bxk Bxk
ð13Þ
The interaction between gas phase and particle phase caused by velocity slip can be modeled as pffiffiffiffiffiffiffi b mVi mVpi kp ¼ b Cpg kkp kp :
ð14Þ
The diffusion term can be modeled as qp ap mVpk mVpi mVpi Bkp 4 lp;l þ fp mVpi pVdi;k 3 Bxk 2 lp;T þ gp Bkp B ¼ ; Bxk rp Bxk ð15Þ
B Bxk
where g p=(4/3)l p,l+f p is the diffusion coefficient.
C.K. Chan et al. / Powder Technology 150 (2005) 42–55
Similar to the gas phase, the dissipation term can be modeled using the concept of Kolmogorov’ dimensional analysis as shown by Zhou and Huang [11] such that BmVpi BmVpi 4 lp;l þ fp ¼ CDp qp ap xp ; ð16Þ 3 Bxk Bxk where x p is the turbulent kinetic energy dissipation rate. Since q pa px p has the dimension of l p,Tk p/l 2, the dissipation term can be modeled as 2 BmVpi BmVpi 4 lp;l þ fp ¼ CDp qp ap kp3 =l; ð17Þ 3 Bxk Bxk where l is the particle phase turbulence scale. Assuming that this scale is proportional to the turbulence scale of gas phase which is determined by turbulent kinetic energy and its dissipation rate, 3
k2 lg ¼ ; e
3
k2 l ¼ Clp ; e
the dissipation rate can then be modeled as 32 BmVpi BmVpi kp 4 lp;l þ fp ¼ Cpe qp ap e : 3 Bxk Bxk k
ð18Þ
ð19Þ
The particle phase turbulent viscosity l p,T and effective viscosity l p,e are defined as 1
3
lp;T ¼ Clp qp ap kp2 k 2 =e;
ð20Þ
and lp;e ¼ lp;T þ lp;l :
ð21Þ
Neglecting the third- and higher-order correlation terms corresponding to the fluctuation of particle phase volume fraction aVp, and substituting Eqs. (11)–(21) into Eq. (10), the transport equation for particle phase turbulent kinetic energy can be obtained as B B ap qp kp þ ap qp mpk kp Bt Bxk gp þ lp;T Bkp Bmpk Bmpi Bmpi B ¼ þ þ lp;T Bxk rp Bxk Bxi Bxk Bxk 32 kp mp;T Bap Cpe qp ap e q gi k rp Bxi p
pffiffiffiffiffiffiffi mp;T Bap Bmpi þ b Cpg kkk kp þ qp rp Bxi Bt Bqp mp;T Bap mp;T Bap Bmpi þ qp mpk þ qp mpi þ rp Bxi rp Bxk Bxk 2Bxk mp;T Bap mp;T Bap Bp
mpi mpi þ : ð22Þ rp Bxk rp Bxi Bxi While there is no definite value for the empirical constants r p, a range from 0.35 to 0.7 has been suggested
47
by Zhou and Huang [11]. In this paper, an intermediate value of r p=0.5 was chosen. C pe is also an empirical constant used in the particle phase turbulence model. It was found that the choice of this constant affects the dissipation rate in the governing equation of particle phase turbulent kinetic energy. Numerical experiment was conducted on C pe and it was found that C pe =1.5 provided the particle velocities closest to experimental data. Therefore, in this paper, the empirical constants associated with the particle phase turbulence model are taken as r p=0.5 and C pe =1.5.
5. Steady and time-averaged governing equations for particle phase The basic equations obtained from kinetic theory for particle phase are applicable to laminar particle phase flow and they can be considered as the instantaneous equations of particle phase flow. For turbulent gas–particle flows, it is necessary to derive the time-averaged equations for both gas phase and particle phase. The time-averaged equations for gas phase are given by Eqs. (1)–(4). Using the above particle phase k p turbulence model and the time-averaged equations of the kinetic theory, the steady governing equations for particle phase are given as follows. Particle phase continuity equation B B ap qp mpk ¼ Bxk Bxk
qp mp;T Bap : rp Bxk
ð23Þ
Particle phase momentum equation BPi;k B Bp ap qp mpk mpi ¼ ap þ þ b mi mpi Bxk Bxi Bxk Bmpi Bmpk B þ þ ap qp gi þ lp;T Bxk Bxk Bxi qp mp;T Bap B B þ mpk þ Bx Bxk rp Bx k i qp mp;T Bap
mpi ; ð24Þ rp Bxk
Table 2 Conditions used in the simulation for comparison with experiments Bader et al. Particle type Particle diameter Particle density Diameter of the riser Height of the riser Volume fraction at inlet Particle mass flux Operating pressure Superficial gas velocity
FCC 76 Am 1714 kg/m3 0.305 m 11.505 m 0.25 98 kg/(m2s) 1 atm 3.7 m/s
48
C.K. Chan et al. / Powder Technology 150 (2005) 42–55
where Pi;k
Bmpk ¼ Pp di;k þ fp di;k Bxk
Bmpi Bmpk 2 Bmpk di;k þ lp;e þ : Bxk Bxi 3 Bxk
Comparing with the particle phase shear stress P i,j described by kinetic theory, the particle phase laminar viscosity l p,l is replaced by effective viscosity l p,e in the P expression of Pi,j . Particle phase pseudo-temperature equation 2 Bmpk B 2 ap qp mpk H ¼ Pi;k þ Bxk 3 3 Bxi B CH þ lp;T BH 2 c;
Bxk Bxk 3 rp
this paper, restitution coefficient in the pipe is taken as e=0.99 and value for particle wall collision is taken as e w=0.8. The particles are considered as spherical with a diameter of 76 Am. A grid system of 150 30 as shown in Fig. 1 is used. The grid size is fixed in the axial direction and varies in the radial direction from 0.0088 m to 0.0028 m. Both gas-phase and particle-phase conservation equations in Eulerian coordinates are integrated in the computational cell to obtain the finite-difference equations. The gas phase equations are solved by means of the SIMPLE algorithm by Patankar [37] with pressure–velocity corrections and a similar procedure is used for the particle phase, but without pressure–velocity corrections. In addition, multiple iterations between the gas phase and the particle phase are adopted based on a two-way coupling.
ð25Þ
where the dissipation term c is similar to that in the kinetic theory.
6. Numerical procedure In order to validate the present model, simulations are made to compare with results of a circulating fluidized bed by Bader et al. [33]. Their experimental conditions as shown in Table 2 are chosen for the numerical simulation. As the exact value of restitution coefficient is difficult to be established, values used in other studies generally range from 0.8 to 1.0. In
7. Boundary conditions The boundary conditions for gas phase and particle phase applied in the simulation are given as follows. Gas phase boundary conditions ! ! ! !
Inlet condition: u g,in(r)=u g,max[1(r/R 2)], v g,in=0. Axial condition: v g=0, Bu g/Br=0 (u g=u g,k,e). Wall condition: u g=v g=k=e=0. Exit condition: Bu g/Bx=0 (u g=u g,v g,k,e). Particle phase boundary conditions
! ! !
!
Inlet condition: u p,in=G s/(a p,inq p), v p,in =0. k p,in=u p,in2, H=10.0. Axial condition: v p=0, Bu p/Br=0 (u p=a p,u p,k p,H). Wall condition as suggested by Soo [38]: u p=L p(Bu p/Br)jr=R , L p=d p/a p1/3, v p=0, Bu p/Br=0 (u p=a p,k p,H). Exit condition: Bu p/Bx=0 (u p=u p,v p,a pk p,H).
8. Results and discussion 8.1. Comparison of different particle phase mathematical model
Fig. 1. Computational grid.
The proposed particle phase mathematical model consists of particle phase kinetic theory and particle phase turbulent kinetic energy model. If the value of C lp is set as zero, and the value of particle phase turbulent kinetic energy (k p) is also kept as zero, the present comprehensive mathematical model becomes the ordinary particle phase kinetic theory. On the other hand, if the value of particle phase pseudotemperature (H) is set as zero, without solving the governing equation of H, the present comprehensive mathematical model changes to simple particle phase k p turbulence model, which is much similar to those dilute particle phase k p turbulence models [11,34,39].
C.K. Chan et al. / Powder Technology 150 (2005) 42–55
Fig. 2 shows the distributions of computed particle phase volume fraction and axial velocity at an axial location of 4.1 m using three different particle phase mathematical models together with the experimental data of Bader et al. [33]. The ordinary particle phase kinetic theory can predict the main feature of core–annulus flow in the vertical pipe, but the predicted particle phase volume fraction just has its high values in annulus region very closed to the wall of the pipe. The particle phase volume fraction rises abruptly near the wall, and has a nearly uniform value within most of the inner region. This phenomenon can also be found in the predictions of particle phase volume fraction by other researchers’ studies using particle phase kinetic theory [20]. However, the experimental data shows that the particle phase volume fraction has a relatively high value in the large annulus region from the nondimensionalized radial location of 0.8 to 1.0.
Fig. 2. Distribution of particle phase volume fraction and axial velocity: (a) particle phase volume fraction; (b) particle phase axial velocity.
49
Comparing with the experimental data, the computed maximum particle phase axial velocity at the axis of the pipe is 14.3 m/s when using particle phase kinetic theory, which is much higher than the experimental data. As the main mechanism of particle–particle collision effect in dense gas–particle flows is not accounted for, the predictions of particle phase k p turbulence model show that k p model cannot predict the main feature of core–annulus flow pattern in the vertical pipe. The computed particle phase volume fraction in the large annulus region is much lower than the experimental values as well as that predicted by the kinetic theory. Furthermore, the predicted particle phase downward flow annulus region is too small and the predicted particle phase axial velocity at the central core is also slower than the experimental values and those predicted by the kinetic theory. Based on the particle phase kinetic theory, the present comprehensive model takes into account the particle phase turbulence effect. It can be seen that the present comprehensive model is capable to predict the main feature of core–annulus flow pattern in the vertical pipe. In addition, the particle phase turbulence reduces the high axial velocity at the central region and also increases the particle phase volume fraction within most part of the pipe annulus. Although compared with experimental data, there are some discrepancies in the near wall region, the present model gives more reasonable particle phase axial velocity distribution in the core region and more reasonable particle phase volume fraction distribution in the annulus region of the vertical pipe than those predicted by the other models. Particle phase pseudo-temperature obtained from kinetic theory and particle phase turbulent kinetic energy obtained from k p turbulence model affect the flow characteristics in a way of particle phase viscosity that exists in the particle phase continuity equation and momentum equation. Fig. 3(a) shows the distribution of computed particle phase effective viscosity, turbulent viscosity and laminar viscosity using three different models respectively. The particle phase turbulent viscosity calculated using k p turbulence model is much larger than the particle phase laminar viscosity calculated using kinetic theory, and the present comprehensive model gives the largest value of particle phase effective viscosity of the three models. Fig. 3(b) shows the predicted results of particle phase effective viscosity and its component using comprehensive model. The two viscosity components of computed particle phase turbulent viscosity and laminar viscosity are similar to the predictions using k p turbulence model and kinetic theory independently as shown in Fig. 3(a). The particle phase turbulent viscosity has a high value in most parts of the pipe and decreases rapidly to almost zero at the wall. Gidaspow and Lu [40] provided experimental data of particle phase viscosity obtained using different methods [1,41]. They were in the range of 0.06 to 0.1 kg/(m s). In this paper, based on a predicted value of particle volume fraction of 0.02, the predicted particle phase viscosity is around 0.06 to 0.07 kg/(m s), which is in good agreement
50
C.K. Chan et al. / Powder Technology 150 (2005) 42–55
uniform distribution of velocity and mass fraction. That is the reason why the particle phase k p turbulence model is not capable of predicting reasonable results of core–annulus flow pattern in the vertical pipe under dense gas–particle flow condition. Similarly, as the present comprehensive model gives the largest value of particle phase effective viscosity, it does not give rise to a more uniform distribution of particle phase velocity and volume fraction. In fact the predictions of the present comprehensive model give reasonable results of core–annulus flow pattern and much higher particle phase volume fraction than the value calculated using the particle phase k p turbulence model. The main reason is that the present comprehensive model considers the particle–particle collision effect described by the kinetic theory. The main difference between the present comprehensive model and the k p turbulence model is that there is a collision pressure term in the momentum equation of the present model. Fig. 4 shows the distribution of computed particle phase pressure due to particle–particle collisions. Predicted results show that the particle–particle collision pressure is nearly zero at the wall and reaches a maximum at the centre of the vertical pipe. Previous experimental results of Polashenski and Chen [43] on packed bed and bubbling bed showed higher solid pressures. However, their superficial gas velocities, ranging from 0.03 m/s to 0.9 m/s were much lower than in the present paper for a riser. The negative gradient of the particle–particle collision pressure is one of the reasons for the accumulation of particles near the annulus wall region. 8.2. Hydrodynamic of gas phase flow
Fig. 3. Distribution of particle phase viscosities: (a) particle phase viscosities using different models; (b) particle phase effective viscosity and its component using the kinetic theory with the k p model.
with data of Gidaspow and Lu [40]. In our study, slip conditions are used for the particle phase wall boundary conditions and results of particle phase pseudo temperature and turbulence kinetic energy are found to be very small. Although experimental results of Polashenski and Chen [42] showed that solid viscosity at the wall was not zero, they were found to be much smaller than in the central region of the pipe. This phenomenon is due to the weak particle–particle interaction near the wall of the pipe. In addition, the particle phase turbulent viscosity is larger than the particle phase laminar viscosity. Therefore, the particle phase effective viscosity is dominated by particle phase turbulence. It is recognized that large turbulent viscosity enhances the momentum and mass transfer and gives rise to more
Fig. 5 shows the distributions of computed gas phase axial velocity at an axial location of 4.1 m using three different particle phase mathematical models together with the experimental data of Bader et al. [33]. In general, the
Fig. 4. Distribution of particle phase pressure due to particle–particle collision.
C.K. Chan et al. / Powder Technology 150 (2005) 42–55
51
Ahmadi. Furthermore, the distribution turbulent kinetic energy for both the gas phase and the particle phase in this paper are similar to those predicted by Cao and Ahmadi. When the modified kinetic theory is used with the k p turbulence model, the particle phase pseudo-temperature as well as the sum of particle phase turbulent kinetic energy and pseudo-temperature are larger thanthe p gas ffiffiffiffiffiffiphase ffi turbulent 2ap qp kinetic The term G kH k or Gp ¼ ¼ C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
energy. p pg srp 2ap qp C þ H k is the production term in Eq. k k pg p srp (3), to modulate the gas phase turbulent kinetic energy for gas particle flow. It can be seen that the gas phase turbulent kinetic energy calculated in the gas–particle flow is larger than that in single gas flow. In addition, the gas phase turbulent kinetic energy in single gas flow attains a peak value near the wall and drops to zero with a steep gradient at Fig. 5. Distribution of gas phase axial velocity.
distribution of predicted gas phase axial velocity is similar to that of the particle phase as shown in Fig. 2(b). As 76 Am FCC particles are used in the simulation, the particle terminal velocity of 0.25 m/s gives rise to very small velocity slip between the gas phase and particle phase. When the modified kinetic theory is used with the k p turbulence model for the particle phase, the predicted gas phase flow patterns are also core–annulus flows. On the other hand, using the k p turbulence model alone for the particle phase, the predicted gas phase is not in a core– annulus flow pattern. Compared with the experimental data, the modified kinetic theory with k p turbulence model predicts more reasonable velocity distribution in the central core region. The gas and particle phases mainly interact in two ways. One is the drag force between the gas phase and particle phase described in the momentum equations, and another is the turbulence interaction between the gas phase and particle phase. In this paper, a k–e turbulence model is used for modeling the gas phase turbulent flow, and the interaction between the two phases is described by the gas phase turbulent kinetic energy and its dissipation rate equations. The computed gas phase turbulent kinetic energy using different particle phase mathematical models in gas–particle flow as well as the single gas phase flow in a vertical pipe are shown in Fig. 6(a). Cao and Ahmadi [44] developed a thermodynamics model to predict relatively dense gas–particle flow in a vertical channel. The predicted maximum value of particle phase turbulent kinetic energy was about 1.5 m2/s2. Their results, based on particle phase mass flux of 20.4 kg/(m2 s) were in good agreement with experimental data of Miller and Gidaspow [41]. In our paper, the particle phase mass flux is 98 kg/(m2 s), which is much larger that that of Cao and Ahmadi. As higher mass flux leads to lower fluctuating kinetic energy of both the gas phase and the particle phase. The maximum predicted value of the particle phase kinetic energy of 0.3 m2/s2 is comparable with that of Cao and
Fig. 6. Distribution of gas phase turbulent kinetic energy and its dissipation rate: (a) gas phase turbulent kinetic energy; (b) gas phase turbulent kinetic energy dissipation rate.
52
C.K. Chan et al. / Powder Technology 150 (2005) 42–55
the wall, whereas the gradient of gas phase turbulent kinetic energy in gas–particle flow is relatively smooth due to the effect of the particle phase. On the contrary, when only the k p turbulence model is used, as the particle phase turbulent kinetic energy is smaller than the gas phase kinetic pffiffiffiffiffiffi ffi turbulent 2a q energy, the term of Gp ¼ sprp p Cpg kkp k is a reduction term in Eq. (3). This reduces the gas phase turbulent kinetic energy for gas particle flow. In fact, experimental studies [45] on hydrodynamics of gas–solid flows in similar vertical pipe show that the gas phase turbulent kinetic energy is not as small as predicted by the k p turbulence model. This indicates again that reasonable results cannot be obtained by using the k p turbulence model for the particle phase alone. Using three different mathematical models for the particle phase, similar phenomena can be found in the distributions of gas phase turbulent kinetic energy dissipation rate in gas–
particle flows and in single phase gas flow as shown in Fig. 6(b). 8.3. Sensitivity analysis of particle phase turbulence model constants Model constants associated with k p turbulence model for particle phase are chosen according to the values used in k p turbulence model for dilute gas–particle flows [39]. The value of C lp is an important factor in the expression of particle phase turbulent viscosity. In the k p turbulence model for dilute gas–particle flows [39], the value of C lp is chosen as 0.0064. In this paper, the parameter sensitivity for C lp is investigated. Fig. 7 shows how a p, U p, k p and H are affected by changing C lp from 0.003 to 0.007. It can be seen that C lp clearly influences the radial distribution of
Fig. 7. Sensitivity of model parameter C lp: (a) particle phase volume fraction; (b) particle phase axial velocity; (c) particle phase turbulent kinetic energy; (d) particle phase pseudo-temperature.
C.K. Chan et al. / Powder Technology 150 (2005) 42–55
particle phase hydrodynamics. With an increase in C lp, the particle phase turbulent viscosity becomes larger, the particle volume fraction becomes smaller in the near wall region, and the distribution of particle phase axial velocity becomes flatter. As the gradient of particle phase axial velocity reduces, the turbulent stress production term and particle phase shear stress in Eq. (22) and (25) respectively become small, reducing the particle phase turbulent kinetic energy and particle phase pseudo-temperature as shown in Fig. 7(c) and (d). When the value of C lp is chosen as 0.005, the distribution of computed particle phase axial velocity is closest to the experimental data among the five cases. Fig. 8 shows the parametric sensitivity of the model predictions with respect to C pg. This constant affects the source term in the gas phase turbulent kinetic energy and its dissipation rate equations, and also affects the source term in
53
the particle phase turbulent kinetic energy equation. It can be seen that C pg has a great effect on the radial distribution of particle phase hydrodynamics. With an increase in C pg, the particle volume fraction becomes smaller, the distribution of particle phase axial velocity becomes flatter, and the particle phase pseudo-temperature also becomes smaller. This phenomenon is similar to the case of increase in the value of C pg. However, particle phase turbulent kinetic energy increases with increasing C pg. This is different from the phenomenon when the value of C pg is increased. As gas phase turbulent kinetic energy is larger than theffiffiffiffiffiffiparticle p ffi phase kinetic energy, the source term of b Cpg kkp kp caused by the fluctuation velocity difference in the governing equation of particle phase kinetic energy is a production term. With increased value of C pg, the particle phase turbulent kinetic energy also increases. When the value of
Fig. 8. Sensitivity of model parameter C pg: (a) particle phase volume fraction; (b) particle phase axial velocity; (c) particle phase turbulent kinetic energy; (d) particle phase pseudo-temperature.
54
C.K. Chan et al. / Powder Technology 150 (2005) 42–55
C pg is chosen as 0.75, which is the same value as in modeling of dilute gas–particle flows [39], the distribution of computed particle phase axial velocity is the closest among the five cases when compared with experimental data.
9. Conclusions A comprehensive mathematical model with turbulence modulation is developed for simulating gas–particle flows in a vertical pipe. By adding a source term into the gas phase turbulent kinetic energy and dissipation rate equations to account for the particle phase, a modified two-equation (k– e) turbulence model is developed for calculating the gas phase flows. Based on particle kinetic theory, a transport equation for particle phase turbulent kinetic energy is proposed and used for modeling the particle phase turbulence. Time-averaged governing equations for particle phase are derived. Computed results show that the particle phase turbulent viscosity gives significant contribution to particle phase effective viscosity. The predicted particle phase volume fraction and velocities of both gas phase and particle phase are in good agreement with experimental data and the present model is capable of predicting core–annulus flow pattern of both gas and particle phases. Compared with the usual particle phase kinetic theory and particle phase k p turbulence model, although there are still some discrepancies in the near wall annulus region of the vertical pipe, the present model gives more reasonable particle phase axial velocity distribution in the core region and more reasonable particle phase volume fraction distribution in the large annulus region. The sensitivity of the model predictions with respect to C pg and C pg indicates that turbulence is an important factor in gas–particle flow behaviour. In the modeling of particle phase flow in gas–particle flows, the particle phase turbulence should be taken into account carefully.
Acknowledgements This work was supported by the Research Committee of The Hong Kong Polytechnic University under Grant No. A-PA81 and G-YC18.
References [1] D. Gidaspow, Multiphase Flow and Fluidization—Continuum and Kinetic Theory Descriptions, Academic Press, New York, 1994. [2] H. Enwald, E. Peirano, A.E. Almstedt, Eulerian two-phase flow theory applied to fluidization, Int. J. Multiph. Flow 22 (1996) 21 – 66. [3] Y.P. Tsuo, D. Gidaspow, Computation of flow patterns in circulating fluidized beds, AIChE J. 36 (6) (1990) 885 – 896. [4] J.A.M. Kuipers, K.J. van Duin, F.P.H. van Beckum, W.P.M. van Swaaij, A numerical model of gas fluidized beds, Chem. Eng. Sci. 47 (1992) 1913 – 1924.
[5] J.A.M. Kuipers, K.J. van Duin, F.P.H. van Beckum, W.P.M. van Swaaij, Computer simulation of the hydrodynamics of a twodimension gas-fluidized bed, Comput. Chem. Eng. 17 (1993) 839 – 858. [6] K.N. Theologos, N.C. Markatos, Modeling of flow and heat transfer in fluidized catalytic cracking riser-type reactors, Trans. IchemE 70 (3) (1992) 239 – 245 (Part A). [7] K.N. Theologos, N.C. Markatos, Advanced modeling of fluid catalytic cracking riser-type reactors, AIChE J. 39 (6) (1993) 1007 – 1017. [8] J.J. Nieuwland, M.S. van Annaland, J.A.M. Kuipers, W.P.M. van Swaaij, Hydrodynamic modeling of gas/particle flows in riser reactors, AIChE J. 42 (6) (1996) 1569 – 1582. [9] C.P. Chen, Studies in Two-Phase Turbulence Closure Modeling, PhD thesis, Michigan State University, 1983. [10] S.E. Elghobashi, T.W. Abou-Arab, A two-equation turbulence model for two-phase flows, Phys. Fluids 26 (4) (1983) 931 – 938. [11] L.X. Zhou, X.Q. Huang, Prediction of confined turbulent gas–particle jet by an energy equation model of particle turbulence, Sci. China, Ser. A 33 (1) (1990) 52 – 59. [12] S. Chapman, T.G. Cowling, The Mathematical Theory of NonUniform Gases, 3rd ed., Cambridge University Press, 1970. [13] J.T. Jenkins, S.B. Savage, A theory for the rapid flow of identical, smooth nearly elastic, spherical particles, J. Fluid Mech. 130 (1983) 187 – 202. [14] J.L. Sinclair, R. Jackson, Gas–particle flow in a vertical pipe with particle–particle interactions, AIChE J. 35 (9) (1989) 1473 – 1486. [15] J.A. Pita, S. Sundaresan, Gas–solid flow in vertical tubes, AIChE J. 37 (7) (1991) 1009 – 1018. [16] J.A. Pita, S. Sundaresan, Developing flow of a gas–particle mixture in a vertical riser, AIChE J. 39 (4) (1993) 541 – 552. [17] H. Arastoopour, Numerical simulation and experimental analysis of gas/solid flow systems: 1999 Fluor–Daniel plenary lecture, Powder Technol. 119 (2–3) (2001) 59 – 67. [18] J. Ding, D. Gidaspow, Bubbling fluidization model using kinetic theory of granular flow, AIChE J. 36 (4) (1990) 523 – 538. [19] J.J. Nieuwland, P. Huizenga, J.A.M. Kuipers, W.P.M van Swaaij, Hydrodynamic modelling of circulating fluidised beds, Chem. Eng. Sci. 49 (1994) 5803 – 5811. [20] A. Samuelsberg, B.H. Hjertager, Computational modeling of gas/ particle flow in a riser, AIChE J. 42 (6) (1996) 1536 – 1546. [21] S. Benyahiaa, H. Arastoopour, T.M. Knowlton, H. Massah, Simulation of particles and gas flow behavior in the riser section of a circulating fluidized bed using the kinetic theory approach for the particulate phase, Powder Technol. 112 (1–2) (2000) 24 – 33. [22] M.S. Detamore, M.A. Swanson, K.R. Frender, C.M. Hrenya, A kinetic-theory analysis of the scale-up of circulating fluidized beds, Powder Technol. 116 (2–3) (2001) 190 – 203. [23] M.Y. Louge, E. Mastorakos, J.Y. Jenkins, The role of particle collisions in pneumatic transport, J. Fluid Mech. 231 (1991) 345 – 359. [24] C.K.K. Lun, Kinetic theory for granular flow of dense, slightly inelastic, slightly rough spheres, J. Fluid Mech. 233 (1991) 539 – 559. [25] J.R. Grace, G. Sun, Influence of particle size distribution on the performance of fluidized beds reactors, Can. J. Chem. Eng. 69 (1991) 1126 – 1134. [26] H.S. Kim, H. Arastoopour, Simulation of FCC particles flow behavior in a CFB using modified kinetic theory, Can. J. Chem. Eng. 73 (5) (1995) 603 – 611. [27] S. Dasgupta, R. Jackson, S. Sundaresan, Turbulent gas–particle flow in vertical risers, AIChE J. 40 (2) (1994) 215 – 228. [28] C.M. Hrenya, J.L. Sinclair, Effects of particle–phase turbulence in gas–solid flows, AIChE J. 43 (4) (1997) 853 – 869. [29] Y. Zheng, X. Wan, A. Qian, F. Wei, Y. Jin, Numerical simulation of the gas–particle turbulent flow in riser reactor based on k–e–k p–e p–H two-fluid model, Chem. Eng. Sci. 56 (2001) 6813 – 6822. [30] Y. Zhang, J.M. Reese, Particle–gas turbulence interactions in a kinetic theory approach to granular flows, Int. J. Multiph. Flow 27 (2001) 1945 – 1964.
C.K. Chan et al. / Powder Technology 150 (2005) 42–55 [31] C.T. Crowe, I. Gillandt, Turbulence Modulation of Fluid–Particle Flows—A Basic Approach, 3rd International Conference on Multiphase Flows, ICMF’98, Lyon, France, 1998. [32] P.E. Chen, C.E. Wood, A turbulence closure model for dilute gas– particle flows, Can. J. Chem. Eng. 63 (1985) 349 – 360. [33] R. Bader, J. Findlay, T.M. Knowlton, Gas/solids flow patterns in a 30.5-cm-dia. circulating fluidized bed, Proc. Int. Circulating Fluidized Bed Conf., Compiegne, France, 1988, pp. 123 – 137. [34] L.X. Zhou, Theory and Numerical Modeling of Turbulent Gas– Particle Flows and Combustion, Science Press and CRC Press, 1993. [35] W.C. Reynolds, T. Cebeci, Calculation of Turbulent Flows, Springer Verlag, 1976. [36] E.E. Khalil, Modeling of Furnace and Combustors, Abacus Press, 1982. [37] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, 1980. [38] S.L. Soo, Fluid Dynamics of Multiphase System, Blaisdell Waltham Press, 1967.
55
[39] C.M. Liao, W.Y. Lin, L.X. Zhou, Simulation of particle–fluid turbulence interaction in sudden-expansion flows, Powder Technol. 90 (1) (1997) 29 – 38. [40] D. Gidaspow, H.L. Lu, Collisional viscosity of FCC particles in a CFB, AIChE J. 42 (9) (1996) 2503 – 2510. [41] A. Miller, D. Gidaspow, Dense, vertical gas–solid flow in a pipe, AIChE J. 38 (11) (1992) 1801 – 1815. [42] W. Polashenski, J.C. Chen, Measurement of particle phase stresses in fast fluidized beds, Ind. Eng. Chem. Res. 38 (1999) 705 – 713. [43] W. Polashenski, J.C. Chen, Normal solid stress in fluidized beds, Powder Technol. 90 (1) (1997) 13 – 23. [44] J. Cao, G. Ahmadi, Gas–particle two-phase turbulent flow in a vertical duct, Int. J. Multiph. Flow 21 (6) (1995) 1203 – 1228. [45] Y.L. Yang, Experimental and Theoretical Studies on Hydrodynamics in Concurrent Upward and Downward Circulating Fluidized Beds, PhD dissertation, Tsinghua University, Beijing, China, 1991.