Analysis of the flow field in cyclone separators

Analysis of the flow field in cyclone separators

Camp,,tees and Fluids, Vol. 2, pp. 249-260. PergamonPress, 1974.lh'inledin Great Britain ANALYSIS OF THE FLOW FIELD IN CYCLONE SEPARATORS C. T. CROWE...

596KB Sizes 0 Downloads 67 Views

Camp,,tees and Fluids, Vol. 2, pp. 249-260. PergamonPress, 1974.lh'inledin Great Britain

ANALYSIS OF THE FLOW FIELD IN CYCLONE SEPARATORS C. T. CROWE* a n d D. T. PRATt Department of Mechanical Engineering, Washington State University, Pullman, Washington, U.S.A. (Received 5 January 1973.)

/kl~lraet--A computer solution for subsonic, two-dimensional, gas-particle flows is applied to the study of the flow field in an axial.inlet peripheral-discharge cyclone separator. The analysis shows the effect of gas-particle momentum coupling on the gas flow field and fractional collection efficiency. The results explain the observed increase in overall collection efficiency with increased dust loading.

NOMENCLATURE a, b, c parameters in generalized differential equation (Equation 10) and d dp particle diameter ~,,/t= unit vector in r and z directions, respectively Stokes law correction factor mass flux M~ particle mass flow Mt shear stress moment mt mass of particle class " i " N, number of effective turns nj number density of particle class " i " P pressure r radius, radial coordinate T shear stress t time U particle velocity V gas velocity V,.,~ volume of a " t a n k " W~ inlet width Xj fraction of particle mass flow issuing from station " j " z axial coordinate = particle acceleration r/. efficiency ~Tt number flow of particle class " i " along trajectory viscosity /z,n effective viscosity p gas density p~ particle material density ~j bulk density stream function A1residence time in a tank to vorticity * Correspondence to Clayton T. Crow¢, Associate Professor, at Washington State University, Dept. of Mech. Eng.; Pullman, Washington, U.S.A., 99163. 249

250

C.T. CRow~and D. T. PRAT]

Subscripts

i o r t z

particle class " i " initial condition radial direction tangentialdirection axial direction INTRODUCTION

THE CYCLONE separator is used extensively in industry to separate dust and condensedphase pollutants from gas streams. A vortex motion created by the flow entering the separator causes a centrifugal force to act on the particles, which moves them toward the v,all and ultimately into a collection bin. The advantages of the cyclone separator are its simplicity, inexpensiveness and the fact that it can be fabricated from almost any material. The disadvantages are its characteristically low collection efficiency and high pressure loss. Although cyclone separators have been widely used for many years, very little is known about the nature of the internal flow field. The design of a cyclone separator is based primarily on experience, previous designs or empirical relationships. One such relationship relates the " c u t size", or the particle size for which 50 per cent collection is achieved, to separa tor operating variables[I ]: [2~tN¢ Vi(pp - P) ~

(1)

where W~ is the inlet width, Vi the inlet velocity and Ne is the effective number of turns in the cyclone. Of course, the number of effective turns which the gas makes during its traverse through the separator is not known a p r i o r i , and must be estimated from experience or from what has been evaluated for a similar design. In fact Ne can be assigned a quantitative value only after dp,o has been measured and equation (I) used to calculate N~. Besides the empirical expression for particle " c u t size", other expressions have been proposed for pressure loss and for the effects of flow rate, gas viscosity and the particle mass density on collection efficiency. In each expression, some empirical factor must be chosen before the expression can be used. The choice of the factor may be difficult if the equation is to be applied to an unconventional design. One parameter which affects the performance of cyclone separators in an interesting fashion is "dust loading" or the weight of dust per unit volume of gas. Increasing the dust loading increases collection efficiency and reduces the pressure drop. This trend suggests that separators should be operated at high loadings to be most effective. However, high loadings also lead to serious erosion and plugging problems. The increase in collection efficiency with increased dust loading is attributed to the large particles impacting and sweeping the smaller particles toward the wall. Baxter[l], after analyzing the data from several different studies on dust loading versus collection efficiency, suggested the following empirical equation : 100 -- rh = (pn~/pn2).l~ 2

100

--

(2)

r/2

where r/ is the efficiency and P8 is the bulk density of the dust. The purpose of the present paper is to analyze the flow field within a cyclone separator in order to provide a better understanding of the influence of the various engineering design

Analysis of the flow field in cyclone separators

251

parameters on separator performance. In particular, the effect of dust loading on the flow field and, in turn, its effect on collection efficiency will be considered. ANALYSTS OF GAS-PARTICLE FLOWS The flow field within a cyclone separator can be regarded as a subsonic, two-dimensional gas-particle flow. The analysis of such a flow field is complicated by the fact that not only must the conservation equations account for the mass, momentum and energy of each phase, but additional equations are required to relate the momentum and energy coupling between phases. The aerodynamic drag on a particle due to particle-gas velocity differences gives rise to momentum coupling between phases. Correspondingly, convective and radiative heat transfer between the gas and particle results in energy coupling. In this analysis, however, it will be assumed that the flow is isothermal and the dust particles and gas are in thermal equilibrium. This assumption, which is valid for cyclone separators operating under typical conditions, corresponds to the condition for which the energy equation is already satisfied and need not be considered. A calculational scheme to obtain stable solutions for internal, two-dimensional, subsonic flows of a single phase has been developed by Spalding[2] and his co-workers. Solutions based on iterative methods had been obtained earlier for unsteady two-dimensional inviscid[3] and laminar[4] flows, but instabilities were encountered at high Reynolds numbers. Incorporating certain features in the analysis of steady flows which improved the stability of the unsteady flow solutions, Spalding et al. developed a solution procedure stable for all Reynolds numbers and applicable to two-dimensional flows with non-uniform properties. A finite difference formulation of the differential equations known as the "tank-and-tube" formulation provides systems of simultaneous algebraic equations for the flow properties which are solved on a digital computer. Crowe and Pratt[5] have shown how this "tankand-tube" formulation can be extended to the analysis of gas-particle flows, and it is this scheme which is applied to the present problem. The extension of the "tank-and-tube" approach to dust-laden flows is reviewed briefly in the following section. TANK AND TUBE APPROACH Basically there are two ways the differential equations for two-dimensional fluid flow can be formulated; either the component velocities and pressure, or vorticity and stream function can be chosen as the dependent variables. The latter was chosen for reasons of computational ei~ciency, and the corresponding differential equations were derived. There are two methods by which finite difference equations can be derived from differential equations; Taylor series expansion or integration over finite areas assuming certain distributions of the variables between nodes. Spalding et a/.[2] employed the latter method and refers to the resulting finite difference equation as the "tank-and-tube" formulation. That is, the area around each node can be likened to a tank connected by tubes to adjacent tanks. In this way, the equation for the dependent variables is valid for a cell of finite size as opposed to a point for the differential formulation. It is through this concept that particle effects can be incorporated into the solution procedure. In deriving the governing differential equations for fluid flow, it is assumed that the fluid is a continuum so that the elemental volume can be shrunk to an arbitrary point. The continuum assumption is not valid for the discrete distribution of particles found in gasparticle flows. In deriving the governing differential equations for particle-laden flows, the

252

C.T. CROWEand D. l'.

PRA.Fr

elemental volume must be small enough to be homogeneous, yet large enough to contain a large number of particles. These constraints are often incompatible. The tank-and-tube formulation circumvents the "continuum '" problem associated with the differential formulation of the equations. The tank itself is the elemental volume through which the particles pass and for which the governing equations in finite difference form must be satisfied. Consider the particle which passes through the tank illustrated in Fig. I. If the particle is accelerated in the z-direction by the aerodynamic drag of the fluid, then there is a corresponding deficit of z-momentum of the fluid in the tank and, according to the tank-and-tube assumptions, is distributed uniformly throughout the tank.

(I, J ÷ I) Tank boundary

\

I . . . . . .

7

i

__L_

(t-l,d}

(I,'J)

B i

/

L ....

/

trajectory i

(t, J -IJ ---

Z

Fig. 1. Tank-and-tube configuration enclosing node (I, J). In the analysis of dust-laden flow in a cyclone separator, it is valid to assume that the particles occupy negligible volume. That is, no effort is made to account for the change in gas-mass in a tank because of the volume occupied by the particles. Because the particles occupy negligible volume and do not contribute mass to the flow, the differential equation for the stream function of a pure gas flow is unchanged by the presence of the particles. The equation for the stream function in cylindrical polar coordinates is

#~

~

+ ~r \~r ~ f + ~ = °

(3)

where co is the vorticity. T h e momentum equations for gas-particle flow can be written as follows*: I 8 8 : (DrVr Vz) + - - ( p V z r qr 8Z

3P 2) 4- E r l i m i ~ . + ~ : T: i "'

(4)

*Although equations (4) and (5) are differential formulations of the two-phase flow equations, they arc to be integrated over a tank area. That is, no difficultiesarise since n~can be regarded as the number of particles of size class " i " per unit tank volume.

Analysis of the flow field in cycloneseparators

253

for the z-direction, and gP

rl Ord (prV 2) + OOz O (p It, Vz) - p TVt2 + ~'i nimi~t" + -&r = T,

(5)

for the r-direction where ~z, and ~,, are the accelerations of particle class " i " in the z and r directions respectively. 7", and Tz are the component shear stresses and are identical to those presented in Ref.[2] The particle accelerations are given by

ctz, = dUz,/dt

(6)

~t,, = d U,,/dt- U,,2 /r

(7)

and Taking the partial derivative of equation (5) with respect to z and subtracting the partial derivative of equation (4) with respect to r gives

div

+ - grad

r

2

+ r- -:-cz

• isop nimia"

(p Vt 2)

--~-~z -- r ~r

nimict,,

=

div grad

+ S,~ (8)

where e, the mass flux vector and the operator " i s o " is defined as O

O

iso = e, ° ---. Oz --ez Or S,~ is a term which derives from the shear stress relations and depends on the distribution of velocity and effective viscosity. Because it is often zero, and because the laws relating to turbulent effective viscosity are not precise, S,o is omitted from the equation. It is argued by Spalding[2] et al. that the error in so doing is small. However, provision is made to include S~, when its evaluation is possible and when it is considered necessary. Since the cyclone separator has a swirl component of velocity, a third equation is needed which relates the change in the moment of tangential momentum to the moment due to viscous stresses; namely,

div(drV,) + ~

m,n,

d(rU,,)/dt =

r~m.

(V,/r) +

r

/2ef f ~ r

(Vt/r)

(9)

The governing differential equations can be expressed as special cases of the more general equation,

(I0) where q~ is the dependent variable in question and a, b, c and dare the appropriate variables. For example, if equation (10) is used for the vorticity equation, equation (8), then = (o/r

CAF Vol. 2 No. 3/4--B

a:.r

2

b=r

2

C =

/Aeff

254

c . T . CROWE and D. T. PRATT

and d = -c- (p V t ) - r cz ~



2

8

8r

Or 0

2



"~z

'

The first two terms in equation (10) are referred to as the " c o n v e c t i o n " terms, the third and fourth as the "diffusion" terms and the fifth as the " s o u r c e " term. The only effect of particles in the flow field is to add the particle acceleration terms to the " s o u r c e " term. That is, the particles represent a " s o u r c e " (or sink) of momentum to the gas flow field. A typical node in the flow field is illustrated in Fig. i where the dotted line identifies the tank boundary. Equation (10) is integrated over the tank area assuming that a, d, n l , ~,,, ez, and the dependent variable are constants and that b and c vary linearly with r and z. Evaluating the integral of the convection terms requires values of the dependent variable on the tank boundaries (in the " t u b e s ' ) . It is assumed that the value of the dependent variable in a tube is equal to that in the tank from which the flow comes. This "upwind difference" scheme is responsible for the unconditional stability of the solution procedure. However, the truncation error associated with the upstream difference scheme introduces a false diffusion of the dependent variable which may be regarded as an artificial viscosity. When the flow is turbulent, the artificial viscosity is usually much smaller than the effective viscosity and, therefore, not significant. Integration of equation (10) around each node results in a system of simultaneous algebraic equations for the dependent variable. The boundary conditions are then established which, for elliptic equations, requires values of the dependent variable or its derivative along the entire flow boundary. The equations are then solved using the Gauss-Seidel successive substitution method. Pressure can be recovered by integrating the component momentum equations after the solution has been obtained. PARTICLE TRAJECTORIES In order to determine the magnitude of the particles' contribution to the source term in the vorticity equation, equation (8), it is necessary to know the particles' trajectories, at what locations they cross tank boundaries and their average acceleration in the tank. The equations of motion in the axial, radial, and tangential direction for particle class " i " can be written as dU. rni dt "' - 3 n f p dv,(V: - U:,) m, \ dt

mi\

2'

(11)

= 3~fpdp,( V, - U,,)

dV,, v, 'rv,.'/t dt +

= 3~f~dr,(V t-

U,,)

(12) (13)

The body force term due to gravity has been neglected because, in a cyclone separator, the centrifugal and aerodynamic drag forces are relatively much larger than the gravitational forces. Also it can be shown that the particle motion due to the centrifugal force is an order

Analysis of the flow field in cyclone separators

255

of magnitude larger than that due to turbulent dispersion. The factorfis defined as CD R e / 2 4 and, for spherical particles, can be evaluated by the empirical relationship[6] f=

1 + 0-15 R e °'687

(14)

up to a Reynolds number of 1,003, which encompasses the range for dust particles in a cyclone separator. Particle trajectories can be obtained by integrating the equations of motion twice with respect to time. However, economy and accuracy are improved if they are integrated once analytically. Assuming the f a c t o r f a n d the gas velocity are constant* over the time interval At, one obtains by integration U~, = Vz - ( 11"=- Uz,o)eXp(- El At)

(15)

U,, = V, - (Vr - Ur,o)eXp(-EiAt) + [1 - e x p ( - . E ~ A t ) ] U , , 2 / r E ~ U,, = V, - ( V , - U,,o)exp( - E , At) - [1 - e x p ( - E ,

At)]U,, Ur,/rE,

(16) (17)

where the subscript " o " corresponds to the initial condition and E i = 18pf/ppdp, 2

The increment in distance is then obtained by multiplying the average velocity by the time interval. SOLUTION PROCEDURE The program must first be supplied with the necessary input data such as the computational grid pattern (distribution of tanks), boundary conditions and other data which relate to the gaseous and particulate phases. Specifically, the stream function at the wall of the separator is set equal to zero and the stream function along the axis is determined by the gas mass-flow rate through the separator. The vorticity variation normal to the wall is approximated as linear, the magnitude of the gradient depending on the effective viscosity. On the centerline the vorticity is zero and varies directly with the radius with departure from the centerline. It is also necessary to choose a discrete size distribution for the particles and locations of particle entry. For example, uniform entry of the particles with the flow must be approximated by discrete entry locations with fraction Xj of the particle mass entering at station " j " . Thus, at station " j " , the number of particles of size class " i " which issue into the flow per unit time is tii = .~j l(4 p/mi

(18)

where ~/p is the total mass flux of particles per unit time. Particle mass conservation requires that n~ be constant along the entire trajectory for particle size class " i " starting at station " j " . Of course, the solution improves with an increasing number of particle sizes and entry locations, but the calculations become more expensive, so that a compromise must be established. The solution is started by establishing the flow field in the absence of particles (ni = 0). Using this flow field, the particle trajectory equations [Equations (15), (16) and (17)] are * It is consistent with the accuracy of the tank-and-tube formulation to assume a uniform velocityin each tank.

256

C . T . CROWE and D. T. PRAa-r

integrated to locate the particle trajectories in the flow. The number density of particle size class " i " in each tank is determined from (19)

nl = i l i A r / V t a . k

where Az is the residence time of a particle in the tank, and Vta,k is the tank volume. The accelerations ct,, and ~tz,are average values for the particle travel through the tank. The flow field is recalculated incorporating the two-phase source terms presented in equation (10). Then new particle trajectories are calculated, the source terms corrected and the process is repeated until further repetition fails to improve the solution. APPLICATION

TO A CYCLONE

SEPARATOR

There are several types of cyclone separators classified by the manner in which the flow enters the separator and the dust is discharged. The separator to be analyzed here is the axial-inlet, peripheral-discharge type sketched in Fig. 2. Swirl vanes at the inlet produce a Dust laden gas

Swirl vanes

Gas discharge tube

Dust Peripheral dust collector

Fig. 2. Schematic diagram of an axial-inlet peripheral-discharge cyclone separator.

rotational motion of the dust cloud and the centrifugal forces accelerate the dust particles toward the wall. The gas streamlines converge and exit through the gas-discharge pipe. The particles move to and subsequently along the separator wall toward the peripheral collector and then are removed. Only a small portion of the total gas flow through the separator is bled offwith the dust particles. Of course, some of the smaller dust particles do not reach the wall and are carried out by the primary gas flow. The collection et~ciency of the separator is defined as the mass of dust collected divided by the total mass of dust which enters he separator.

Analysis of the flow field in cycloneseparators

257

The cyclone separator studied in this paper was 2.4 in. dia and 9"6 in. long. The i. dia. of the inlet port and the dia of the gas discharge tube were 1"4 in. The mixture entered the separator at 44 ft/sec through swirl vanes oriented at 45 °. The computational grid pattern used for the calculations is illustrated in Fig. 3. The 11 by 11 grid was selected to give a reasonably fine description of the flow field without requiring ~/

-Ztt



/

9 8

/

/

/

/

_I

9"6 In. /

/

/

, /I

/

Pori'icle entry locotions

I///

6-"

--

T

1.2 in.

0.7 iri.

4 /

3/ 2 J

/

__.. 1 2

L.

I

_1_

1

3

4

5

6

.L ~/ 7

I 8

Fig. 3. Configuration, computational grid pattern and particle entry locations used for computer solution.

excessive computer time. Indeed, the solution can be improved by choosing a finer grid although Spalding[2] and his coworkers found that the solution corresponding to eleven grid lines can give quite accurate predictions of vorticity distributions. The purpose of this paper, however, is to illustrate how the computer model can be applied to the analysis of the flow field in a cyclone separator and not to accurately predict the efficiency of a particular separator. The gas was assumed to be air at standard conditions. The gas mass bled off through the particle collector was neglected although it could be included if a more detailed solution were desired. The inlet flow of particles was approximated by three entry locations as shown in Fig. 3. It was also assumed that the particles and gas were in velocity equilibrium at entry. The particle trajectory calculations were initiated at the entry ports and executed using a time interval of 0.0001 sec. It took approximately 30 time intervals for a particle to traverse a tank. The discrete particle-size distribution used for the calculations is shown in Fig. 4 together with the continuous distribution which it approximates. This somewhat refined distribution of smaller size particles was needed to properly evaluate the collection efficiency since it is only these smaller particles which escape. It was also assumed that the particles were spherical and had a specific gravity of 4. The eddy viscosity model suggested by Spalding[2] and his coworkers for the effective viscosity in cylindrical combustion chambers was employed in the analysis of the cyclone separator because of the geometric similarity between the two. The important feature of the gas-particle flow analysis used in this study is the fact that it accounts for the particle-gas momentum coupling. Thus, calculations were done at several

258

C.T. CROWEand D. T. PRATT --

60

40--

Actual

distr,bution

20--

=L

I0

/

/

/

10/~// //

--

8 --

E

o

/ //21/.t

/

6.3/.LJ " /

4'7/.t/F

6 --

5.7 F.

4

4-g

/

/ E

/ /

,

I

I

I0

Weight

3,C

per cent

i

I

I

!

50

70

90

99

less

then

size

steered

Fig. 4. Discrete particle-size distribution used for dust.

dust loadings to illustrate the capability of the analytic approach. If the effect of the particles on the flow field were neglected, as often done in gas-particle flow analyses, then dust loading would have no effect on the flow field. The dust loadings considered in this study were 0, i 5, 30, and 60 and 120 grains/ft 3. A typical run required 6 min on an IBM 360--67. RESULTS

Typical gas streamlines and particle trajectories predicted by the calculations are shown in Fig. 5. As expected, the gas streamlines diverge after leaving the annular entry region and then converge again as the gas exits through the gas-discharge tube. The particles move I

r"

I

2IFI

I ,/

I

/

l

/lO~

. . . .

]

- -

Gas

---

Por~icle

/

/

I

l

l

. . . . .

I

I

/

6.3p. \

1 .

streamline trojectory

Fig. 5. Predicted gas streamlines and particle trajectories at a dust loading of 30 grains/ft 3.

toward the wall due to the centrifugal forces; the larger, more massive, particles impact the wall in a short distance while the smaller particles diverge little from the gas streamlines. Although it is not shown in Fig. 5, the calculations predict a small on-axis recirculation zone adjacent to the upstream face of the separator. Few particles would reach and become trapped in this zone, however.

Analysis of the flow field in cyclone separators

259

The effect of dust loading on the radial and axial gas velocities at a station 3.6 in downstream of the inlet port is illustrated in Fig. 6. Comparing the axial velocity distribution of the flow with no dust particles and that for a dust loading of 60 grains/ft 3, one notes that the %\\%\

\ \\

\ \ \ \ \ \

...//

Dust-laden

I

',

[

\

~,, ~ ",.

\

\ \\

\ \ \\

i

i -0 75

\ \\\\\\

oo6

qas

....

-/

D u s t - laden gas

0-04-

Dust,.eqa -I'0

\ \ \ \

/

/

/

/,. J

-0'50

Rodi01 v e l o c i t y ,

j -0"25

ft/sec

W!

",x, 0

5

I0

Axial v e l o c i t y ,

15

J

J

20

25

ft/sec

Fig. 6. Effect o f dust l o a d i n g on the radial and axial gas velocities at 3-6 in f r o m inlet.

effect of" dust loading is to cause a slight increase in axial velocity near the periphery of" the separator and a reduction in velocity along the centerline. The increase in velocity near the wall is the result of" the larger particles "dragging" the gas along in the axial direction. The decreased velocity near the centerline occurs because the mass flow rate is the same in both cases and mass conservation must be satisfied. The effect of" dust loading on the radial component of" gas velocity is more striking. As one notes in Fig. 6, the radial velocity is negative which means the flow is moving toward the centerline. When dust particles are introduced into the flow the radial speed of" the gas toward the centerline is reduced. In fact, the reduction in radial speed varies monatomically with increased dust loading. This reduction in radial speed renders the gas flow field less effective in moving the dust particles towards the centerline and allowing their escape through the gas-discharge tube. Thus, the effect of" increased dust loading is to increase collection efficiency which qualitatively agrees with and explains the experimental results. The effect of" dust loading on the collection efficiency as predicted by the computer solution for the particle sizes used in the analysis is shown in the Table. The column for zero dust loading corresponds to the efficiency one would predict by neglecting gas-particle momentum coupling. One notes that as dust loading is increased the collection efficiency of` the smaller particles increases accordingly. In fact, all the 4"7/~ particles are collected at a dust loading of" 120 grains/f`t 3. One also notes that the cut size is 3 microns at a dust loading of" 60 grains/f`t 3 and that it decreases with increased dust loading. Applying I.~pple's equation [Equation (I)] predicts a cut size of` 2.9 # which, considering the approximations made in its derivation, shows remarkably good agreement with the more detailed computer solution. There is no provision in l_~pple's equation, however, to predict the effect of" dust loading on cut size. The effect of dust loading on the overall collection efficiency is illustrated in Fig. 7. One notes as before a monatomic increase in efficiency with dust loading. It is difficult to corn-

C . T . CROWE a n d D . T. PRA'Vr

260

pare the predicted overall collection efficiency with Baxter's empirical equation, equation (2), because the type of separator and range of dust loading for which it is applicable is unknown. In fact, Baxter's equation suggests that as the dust loading is decreased toward zero, the collection efficiency approaches 100% which, of course, is incorrect.

96

"g =

94

g u 92 9O

0

15

I

30

4'5 Dust

=

I

60 I0oding,

75

90

105

~20

gift 3

Fig. 7. Predicted dependence of overall collection efficiency on dust loading.

A series of controlled experiments are currently being developed at Washington State University to measure the collection efficiency of axial-inlet peripheral-discharge separators to assess the accuracy of the model. CONCLUSIONS

The model developed to analyze the flow field in a cyclone separator shows the effect'of dust loading on the flow field and explains the observed increase in collection efficiency with increased dust loading. The model can easily be used in its present form to predict the effect on efficiency of other designs and operating parameters such as separator length, diameter, swirl vane angle, inlet velocity and so on. In addition, the effect of particle size distribution and material density on collection efficiency is easily within the capabilities of the present program. This analytic model should serve as a most useful tool to the designer of dustcontrol equipment alleviating the need to rely solely on empiricism. Acknowledgements--The authors wish to acknowledge the support of the W,S.U. Computer Center in carrying out this program. REFERENCES 1. Stern A. C., Air Pollution, 3, Academic Press, New York 0968). 2. Gosman A. D., Pun W. M., Runchal A. K., Spalding D. B. and Wolfshtein M., Heat and Mass Transfer in Reeireulati~ Flows, Academic Press, New York (1969). 3. Blair A., Metropolis N., Taub A. H. and Tsingou M., A study of a numerical solution to a two-dimensional hydrodynamic problem, Physics and Mathematics, LA-2165 (1957). 4. Baraket H. Z. and Clark J. A., Analytic and experimental study of transient laminar natural convection flows in partially filled containers, Proe. 3rd Int. Heat Transfer Conf., Chicago, 11, (1966). 5. Crowe C. T. and Pratt D. T., Two-dimensional gas particle flows, Proc. 1972 Heat Transfer and Fluid Mechanics Institute, Stanford University Press, p. 386, (1972). 6. Wallis G. B., One-Dimensional Two-Phase Flow, McGraw-Hill New York (1969).