A control-volume finite element method for dilute gas-solid particle flows

Computers Fluids Vol. 23, No. 8, pp. 1073--1096, 1994 Copyright © 1994 Elsevier Science Ltd Printed in Great Britain. All fights reserved 0045-7930(94)E0001-E

Pergamon

0045-7930/94 $7.00+ 0.00

A CONTROL-VOLUME DILUTE

FINITE

GAS-SOLID

ELEMENT PARTICLE

METHOD

FOR

FLOWS

C. MASSON and B. R. BALIGA Department of Mechanical Engineering, McGill University, 817 Sherbrooke Street West, Montreal, Quebec, Canada H3A 2K6 (Received 8 October 1992; in revised form 9 August 1993)

Abstract--The formulation of a co-located equal-order Control-Volume-based Finite Element Method (CVFEM) for the solution of two-fluid models of 2-D, planar or axisymmetric,incompressible,dilute gas-solid particle flows is presented. The proposed CVFEM is formulated by borrowing and extending ideas put forward in earlier CVFEMs for single-phase flows. In axisymmetricproblems, the calculation domain is discretized into torus-shaped elements and control volumes: in a longitudinal cross-sectional plane, or in planar problems,theseelementsare three-nodetriangles, and the control volumesare polygons obtained by joining the centroidsof the three-nodetriangles to the midpoints of the sides. In each element, mass-weighted skew upwind functions are used to interpolate the convectedscalar dependent variables and the volumeconcentrations.An iterative variable adjustment algorithm is used to solve the discretized equations. The capabilities of the proposedCVFEM are illustrated by its application to two test problems and one demonstration problem, using a simple two-fluidmodel for dilute gas-solid particle flows. The results are quite encouraging.

1. I N T R O D U C T I O N The objective of this paper is the adaptation of some available Control-Volume-based Finite Element Methods (CVFEM) for single-phase flows [1-8] in order to formulate a CVFEM for the solution of two-fluid models [9, 10] of 2-D, planar or axisymmetric, incompressible, dilute gas-solid particle flows in regular- or irregular-shaped geometries. The capabilities of the proposed CVFEM are illustrated in this paper by its application to two test problems and one demonstration problem. The results obtained for the first test problem are compared with a proposed analytical solution. In the second test problem, the CVFEM results are checked against those of another numerical investigation [11] and with results obtained using a well-established finite-volume staggered-grid method [12, 13]. In the demonstration problem, the proposed CVFEM is used to simulate laminar, 2-D, axisymmetdc, dilute gas-solid particle flows in a geometry that is characteristic of a split-flow inertial separator [14]. In dilute gas-solid particle flows, the motion of the particles is controlled by local aerodynamic forces that they experience, and the influence of particle-particle collisions is negligible. In contrast, the particle motion in dense gas-solid particle flows is controlled by particle-particle collisions, and aerodynamic forces on the particles have negligible influence on their trajectories. Further discussion of the distinguishing features of dilute gas-solid flows and quantitative qualification criteria are available in the works of Crowe [9, 10] and Ishii et al. [15]. Examples of gas-solid particle flows in the environment include dust and sand storms, and transport of particulate pollutants in the atmosphere. Industrial examples of such flows are found in pneumatic conveyor systems, pulverized coal fired furnaces, cyclone separators and classifiers, spray drying and cooling systems, and split-flow inertial separators for air intakes of helicopter gas-turbine engines. Crowe [9, 10] has provided comprehensive reviews of available mathematical models o f dilute gas-solid particle flows. Rigorous discussions of mathematical models of a variety of multiphase flows are available in the works of Ishii et al. [15], Anderson and Jackson [16], Whitaker [17], Ishii [18], Slattery [19], Drew [20], Crapiste et al. [21], and Jiang et al. [22]. A comparative evaluation of these various models is not the intent of this paper. Rather, as was stated earlier, the goal here is to propose and demonstrate the capabilities of a CVFEM for the solution of such models in complex 2-D geometries. 1073

1074

C. MASSONand B. R. BALIGA

CVFEMs have proved quite successful in the solution of single-phase fluid flow and heat transfer problems [1-8], but they have not yet been used for the prediction of gas-solid particle flows. The CVFEM proposed in this paper is based primarily on ideas contained in earlier CVFEMs put forward by Baliga and Patankar [1], Prakash and Patankar [5], Schneider and Raw [7, 8], and Saabas [6] in the context of single-phase flows in 2-D (planar) and 3-D geometries. The adaptation of these ideas for the solution of two-fluid models of the problems of interest is presented in this paper via concise descriptions of the various steps involved in the formulation of the proposed CVFEM. Two-fluid models of dilute gas-solid particle flows have been solved in the past mostly using finite-volume or finite-difference methods that use a primitive-variables formulation and staggered grids for the velocity components [9-11, 23, 24]. These methods are well established and are the basis of robust computer codes. In terms of geometrical flexibility, however, such methods are quite limited because they are typically formulated for problems in which the boundary surfaces lie along the axes of commonly used orthogonal coordinate systems. Furthermore, such staggered-grid primitive-variables formulations are not well-suited for extensions to complex irregular-shaped geometries via curvilinear orthogonal or nonorthogonal grids [25, 26]. The goal in this work is to overcome this limitation through the use of CVFEMs, which provide the geometrical flexibility of Finite Element Methods (FEMs) along with the ease of physical interpretation traditionally associated with Finite Volume Methods [2, 12]. Attention in this paper is limited to laminar gas-solid particle flows, but the proposed CVFEM is capable of solving some of the available models of turbulent gas-solid flows [27, 28]. The explicit demonstration of this capability, however, is not within the scope of this paper. 2. G O V E R N I N G EQUATIONS The most detailed model for gas-solid particle flows is the so-called complete local description [29]. This approach deals with the dynamics of each phase, and the interfaces between these phases, on the basis of first principles [29, 30]. Though the complete local description models the actual flow exactly, the simulation of practical problems using this model is beyond the capacity of available computers [10], and some simplifications are needed to obtain a tractable model. Such simplifications can be obtained by using some averaging processes [18, 20-22, 29]. These average descriptions are also called the continuum mechanical approach [20] and are applied to problems where complete details of the actual flow are not needed. In such cases, the average model describes each phase as a continuum, occupying simultaneously the same region in space: the gas phase is commonly referred to as the continuous phase, and the solid particles are denoted as the particulate phase. Numerous average descriptions have been proposed both for the continuous and particulate phases [11, 16-24,29-32]. The resulting models use an Eulerian formulation both for the continuous and the particulate phase. In such an Eulerian approach, in addition to the specification of how each phase interacts with itself through stresses, the mutual interaction of the two phases has to be included. These additional terms usually need further modelling, based on basic rules of physics and empirical data applicable for a specific flow topology (for example, stratified and dispersed flows): several rigorous, and fairly complex, models have been proposed [15, 19-22, 30], but some major sources of uncertainties remain. Jiang et al. [22] have proposed expressions for the stresses and the interaction terms based on an idealized model of the flow within the averaging region. The resulting model consists of 18 equations and 18 unknowns. Consideration is limited in this paper to 2-D, planar or axisymmetric, laminar, non-reacting, incompressible, dilute gas-solid particle flows in regular- or irregular-shaped geometries. It is also assumed that the particles are uniform solid spheres. A relatively simple mathematical model is used. It is a specialization to dilute concentrations of a general model proposed by Anderson and Jackson [16]. It is not the most rigorous or sophisticated mathematical model that we could have chosen, but it is considered to be adequate within the scope of the stated goals of this work. The chosen mathematical model consists of a set of six differential equations: a continuity equation and two momentum equations for each of the two phases. Two-dimensional axisymmetric versions of these equations are presented here with respect to the cylindrical (r - z ) coordinate system.

Control-Volume Finite Element Method

1075

2. I. Continuous phase z-momentum equation

3 t(£pGuG)..] O-'

-

3 3p _SDz.~_SGz ~Z (¢p GUGUG) + r 1~ r3 (r£pGvGuG)= --¢~Z

3 [

3uG'~

10 (

t~uG~

tr"w)

(1)

r-momentum equation 3 3 1 3 3p vG 3t (Ep GvG) + ~Z (ep GUGVO) + r Or (rEp GVGvG) = -- E ~r -- So, + S ° - # - ~

(2) continuity equation

(ePG) +

(3)

(EpGuG)+ -rOr (rEpGv G) = 0

2.2. Particulate phase z-momentum equation 3 d d + 3 __ 13 3p ~t (~tp u ) 3z (otpOudud) + -r Or (rCtpdvdu d) = -- Ct-~Z + SDz + S~

(4)

r-momentum equation

3 (txpdvd) +

3t

13

(ozpdudvd)+ r Or (ro~pdvdvd)

3p

d

--~X~r "~ SDr "~-Sr

(5)

continuity equation

(0~Pd) "["

(ctpdud) + -rOt (r~tpdv d) = 0

(6)

where So~ = K(u C - Ud); SD, = K(v G -- v d)

3

/~

K=gcte-~CDRe

d pGIvd -- VGId ; Re d /z

(7) (8)

The six unknowns are u d, v d, u G, v G, p and 0c The superscripts G and d refer to the continuous (Gas) and particulate (dispersed) phases, respectively; pG and pd are the mass density of the continuous and particulate phases, respectively; p is the pressure;/~ is the dynamic viscosity of the continuous phase; u c and v G, and u d and v d, are the intrinsic volume-average [22] velocity components in the z and r directions of the continuous and particulate phases, respectively; p denotes the continuous-phase pressure; and S~, S0, S~ and S,d are the corresponding volumetric source terms, where body forces, such as gravitational force, can be included. The volume concentration of the particulate phase is denoted by ,c The volume concentration of the continuous phase, E, is related to ~ by: e + ~t = 1. The particles are considered to be of uniform size, with diameter d. Co is the drag coefficient, and Re d is a particle Reynolds number. The terms (SDz -- ~t(3p/3z)) and (SD, -- ~(3p/3r)) are the total mutual z- and r-direction forces, respectively, exerted by the continuous and particulate phases on one another, per unit volume of the mixture. Soz and SD, will be referred to as the momentum coupling terms. As has been discussed by Rudinger [33] and Maxey and Riley [34], when the density of the particulate phase is much greater than that of the continuous phase, the most important mutual force is the frictional drag: the expressions given in equations (7) and (8) for SD, and SD, include only the frictional drag, because it is assumed that for the gas-solid particle flows of interest here, pd/pG >1 103. Additional

1076

C. MASSONand B. R. BALIGA

details regarding these terms are available in the works and Riley [34]. Again, only a simple model has been not to undertake a detailed evaluation of the various propose and demonstrate a CVFEM for the solution particle flows.

of Ishii et al. [15], Rudinger [33], and Maxey used in this work because the goal here is available models; rather, the aim here is to of two-fluid models of 2-D dilute gas-solid

3. PROPOSED M E T H O D In the proposed CVFEM, all of the dependent variables are stored at the same nodes (co-located) and interpolated over the same elements (equal-order). As was stated earlier, this method was constructed by adapting and extending ideas from earlier CVFEMs for single-phase incompressible fluid flows proposed by Baliga and Patankar [1, 2], Prakash and Patankar [5], Schneider and Raw [7, 8], and Saabas [6]. The description of the proposed numerical method is given for 2-D axisymmetric problems. However, it is presented in a general manner, with respect to the 2-D cylindrical and Cartesian coordinate systems shown in Fig. 3(a). The discretized equations are written in a form that allows both the 2-D axisymmetric and planar (Cartesian) formulations to be obtained easily. The proposed expressions are appropriate for 2-D axisymmetric problems; and by setting r = 1 in the various discretized expressions, the 2-D planar (Cartesian) formulation is obtained. Concise descriptions of the various steps involved in the formulation of the proposed CVFEM are presented in this section, with suitable emphasis on features that are of particular interest in simulations of dilute gas-solid particle flows in complex geometries. 3. I. D o m a i n discretization

It is convenient to present the domain discretization procedure with respect to a longitudinal cross section of the axisymmetric domain of interest. This cross section is first divided into three-node triangular elements. Then the centroids of the element are joined to the midpoints of the corresponding sides. This creates polygonal control volumes around each node in the finite element mesh. The longitudinal cross section of a sample domain discretization is shown in Fig. 1: the solid lines denote the domain and element boundaries; the dashed lines represent the

t Fig. 1. Discretizationof the cross section of a typical axisymmetriecalculation domain.

Control-Volume Finite Element Method

(a)

1077

(b) Fig. 2. Typicalcontrol volumes.

control-volume faces; and the shaded areas show the control volumes associated within one internal node and one boundary node. The discretization of the longitudinal cross section is rotated through 2r~ radians about the axis of symmetry. The result is a discretization of the axisymmetric calculation domain into torus elements, of triangular cross section, and torus control volumes, of polygonal cross section. In the rest of the paper, for conciseness in the presentation, these torus elements and control volumes will be referred to as triangular (3-node) elements and polygonal control volumes, respectively. 3.2. Integral momentum conservation equations The discussion in this section is presented with respect to the continuous-phase z-momentum equation. The treatment of the other momentum equations is analogous. Consider a typical node i in the calculation domain: it could be an internal node, such as the one shown in Fig. 2(a), or a boundary node, similar to the one shown in Fig. 2(b). An integral formulation corresponding to equation (1), when applied to the polygonal control volume surrounding node i in Fig. 2, can be written as follows: d . n2nr ds +

d. n2nr ds +

e

+ SDz -- S~ + ~ (ep GuG) d~e~

aoc

+ [similar contributions from other elements surrounding node i] + [boundary contribution if applicable] = 0

(9)

where n is a unit vector normal to the differential length element ds and pointing outward with respect to the control volume, d is the combined convection~liffusion flux of uG: J = Jc + JD

(10)

Jc = (~P°V~) uG

(11)

Jo = - v V u°

(12)

With reference to equation (11), in expression (EpGV~), the subscript m is attached to the velocity vector Vm C in order to emphasize its connection to the mass flux: this velocity is interpolated in a special way, as is discussed in the next section. It should also be noted here that in the particulate-phase momentum equations, the diffusion flux, JD, is assumed to be negligible. The form of equation (9) emphasizes that it can be assembled by using an element-by-element procedure. 3.3. Interpolation functions The derivation of algebraic approximations to the integral conservation equations requires the specification of element-based interpolation functions for the dependent variables, continuousCAF 23/B--G

1078

c. MASSONand B. R. BALIGA

phase viscosity, source terms, mass densities, and volume concentrations. For conciseness in this presentation, these interpolation functions will be discussed here only for u c, V~, pC, /~, ct, p, So~. and S~. Other dependent variables are interpolated analogously. For convenience in the formulation of these functions, in each element, a local (x,y) coordinate system is defined such that the origin is at the centroid of the triangular element, the x axis is in the direction of z and the y axis is in the direction of r, as shown in Fig. 3(a). Some of the interpolation functions will be expressed with respect to this local coordinate system. 3.3. I. Viscosity, density and sources. In each triangular element, the centroidal values of # and pC are assumed to prevail over the corresponding element. The volumetric source term, S~, is stored at the vertices of the triangular elements and assumed to prevail over the corresponding portions of the control volumes within that element. The momentum coupling source term, SD.., is linearized, and expressed in the following general form [12]: (13)

Soz = (So=)c + (So:)pu G where (SDz)C :

-Ku °

(SDz)p =

K.

(14)

In each triangular element, the values of (Sm)c and (SDz)p are stored at the vertices, and are assumed to prevail over the corresponding portions of the control volumes within that element. This linearization of the momentum-coupling source term allows simultaneous (coupled), iterative solution of the continuous- and particulate-phase discretized momentum equations. This linearization and the iterative, coupled-solution algorithm has proven to be much less prone to convergence difficulties than other linearizations with a purely segregated iterative solution algorithm. Additional details are provided in Section 3.7. 3.3.2. Velocity in mass-flux terms. In the calculation of mass flow rates across the control-volume faces, the velocity is interpolated in a special way in each element, and this mass-flux related velocity is denoted by: Vm 6 = u~i + v~j

(15)

In the case of the particulate phase, the superscript, G, is replaced by the superscript, d. This special treatment, borrowed from the works of Prakash and Patankar [5] and Saabas [6], prevents the occurrence of spurious pressure oscillations in the proposed co-located equal-order CVFEM. The development of this interpolation is based on the discretized momentum conservation equations. Therefore, it will be presented later in this section. 3.3.3. uCin diffusion terms. In the derivation of algebraic approximations to the surface integrals of the diffusion flux in equation (9), u G is interpolated linearly in each element. Thus, with reference 3 F

xt-% 2

2 I

a

a

v

z

(b) Fig. 3. Typical triangular element.

Control-Volume Finite Element Method

1079

to the local (x, y) coordinate system in Fig. 3(a), the interpolation function for u G in element 123 is given by (16)

u G = A x + By + C

An equivalent, and perhaps more elegant, development of this linear interpolation on triangular elements could IX done using barycentric or area coordinates, traditionally employed in FEMs [35]. It should also be noted that with such linear interpolation functions, Delaunay triangulation is required to ensure that algebraic approximations of the diffusion transport terms contribute positively to the coefficients in the discretized equations. Barth [36] has presented a formal proof of this statement for 2-D planar problems. 3.3.4. u G in convection terms. In the derivation of algebraic approximations to the surface integrals of the convective flux in equation (9), a MAss-Weighted skew upwind interpolation scheme (MAW) for u G is used. The MAW scheme defines mass-weighted average values of u G at the integration points associated with each of the three control surfaces within a triangular element [Fig. 3(a)], namely, u °, usG, U~, in the following manner: Let

rh, =

j'o

£p°V~ " n,2nr ds

0

rh~ = f ) ePCV6m" n~2nr ds

m, =

fo

EpGV~• n,2nr ds

(17)

where n,, n~ and n, are the unit-normal vectors shown in Fig. 3(b). Then

[ (m')]

U~ + (1 --fro)UP where f m = min max - ~ ,

~fm

0 ,1 ,

if #l~lr > 0

u° = LfmU ~ + (1 --fm)UC2 where fro=rain max - ~ - ~ , 0 , 1 , i f r h , < 0

08)

\ m s ' 0) , 11 , if m , > 0 I frou° + (1 --fm)U ° where f m = minrmax(rh'/

u° = Lieu ~, + (1 -f~)u~

f,.u~ + (1 --fm)U~ w,oro

.,O =

I

L f"u~

+

(1 -f..lu °

']'

whereS'.:minPmax(-'(",O),L \ m:

where f,, = rain max

,

m,o), ,}, ,0

,

11, 3

if m~ < 0

(19)

if m, > 0 if m, <

0

(20)

These mass-weighted average values of u G at the integration points r, s and t are assumed to prevail over the corresponding control surface when the surface integrals of the convective-flux terms in equation (9) are evaluated. The proposed MAW scheme is designed to prevent negative coefficients in the algebraic discretization equations. This scheme is an adaptation of similar schemes proposed by Schneider and Raw [7] and Saabas [6]. It ensures, at the element level, that the extent to which the dependent variable at a node exterior to a control volume contributes to the convective outflow is less than or equal to its contribution to the inflow by convection. Thus, it is a sufficient condition to ensure that the algebraic approximations to the convective transport terms in equation (9) contribute positively to the coefficients in the discretization equation. It should also lx noted that the MAW scheme takes better account of the influence of the direction of the flow than the donor-cell scheme of Prakash [37], so it is less prone to false diffusion [6, 7, 37].

1080

C. MASSON and B. R. BAL1GA

The Flow-Oriented upwind interpolation scheme (FLO) of Baliga and Patankar [1, 2], and the streamline-upwind Petrov-Galerkin scheme (SUPG) of Brooks and Hughes [38] are more accurate than the MAW scheme. However, the FLO and SUPG schemes can produce negative coefficients in the discretization equations, especially when obtuse-angled elements and high element Peclet numbers are involved [6]. In the numerical solution of gas-solid particle flows, it is very useful to avoid negative coefficients in the discretization equations; otherwise, during the solution process, always-positive dependent variables, such as volume concentration of the continuous and particulate phases, could take on negative values and lead to a failure of the overall solution algorithm. In order to avoid this problem, only the MAW scheme was used in this work. A better approach may be to use the FLO scheme, whenever possible, and switch to the proposed MAW scheme when negative contributions to the coefficients in the discretization equations are encountered, at the element level [6]. However, such a hybrid approach is outside the scope of this work. The MAW scheme defined by equations (18)-(20) is highly implicit. This does not pose any special difficulties in the proposed derivation of the discretization equations, as presented in the next section, because it is based on successive-substitution, or Picard, linearization of the convective transport terms in the momentum equations. However, the MAW scheme would make Newtontype linearizations very difficult. It should also be noted that in this scheme, to obtain expressions for up, up, and u~ in terms of u~, up, and u~, a 3 x 3 matrix of element-interpolation coefficients must be inverted. Further details are available in the work of Saabas [6]. 3.3.5. Pressure p. Pressure is interpolated linearly in each element: (21)

p = dx + ey + f

3.3.6. Volume concentration a. In most of the available finite volume methods for dilute gas-solid particle flows, the function that is used to interpolate a is based on the upwind scheme [9, 11, 23, 24]. The donor-cell scheme of Prakash [37] is one way of implementing this idea in CVFEMs, but this scheme is prone to considerable false diffusion [6, 37]. In this work, an adaptation of the MAW scheme described previously has been implemented. The modified MAW scheme defines a material mass-weighted average of a at each of the three control surfaces of a triangular element [Fig. 3(b)], namely a,, ~s, a,, in a similar manner to that for the mass-weighted averages of u G. The values of a,, as, a, can be determined from equations (18) to (20) with the following modifications: replace up, up and up by a,, as and at, respectively; replace u~, u2c and u3c by a~, a2 and a3, respectively; and rh,, rh~. and rh, by M•,a, M~ • ~ and M• ,d, respectively, where flit d =

fo

p d V dm • nr2rtr

ds

• d _Ms

pdVdm"

M~ =

p W ~ " nt2nr ds

ns2nr ds

(22)

These material mass-weighted averages of a are assumed to prevail over the corresponding control surfaces within the element when the mass flow rates in the integral continuity and momentum equations are evaluated. E is computed using a + E = 1. 3.4. Discretized momentum equations

The discretization equations are obtained by first deriving algebraic approximations to the element contributions and the boundary contributions, if applicable, and then assembling these contributions appropriately. Again for conciseness, this derivation will be presented only for the continuous-phase z-momentum equation. The following discussion pertains to node 1 of the element shown in Fig. 3. In each element, there are diffusion, convection, source and unsteady contributions. The derivation of algebraic approximations to each of these contributions is presented separately.

Control-Volume Finite Element Method

1081

3.4.1. Element diffusion contribution. In each element, the diffusion flux Jo can be expressed in terms of its components in the z and r directions: ~uO\" { duO\.

<,3)

where i and j are unit vectors in the z and r directions, respectively. The interpolation function given in equation (16) is used to approximate JO" With reference to element 123 in Fig. 3(a) and the r-z and x - y coordinate systems, the diffusion contribution is approximated as follows:

f ~ JD" n2nr ds = 2~z# ~r° + ra [Aya -

Bx~]

_ r o + rc 3D" n2nr ds = 2nl~ ----f-- [ - A y c + Bxc]

fo :

(24) (25)

3.4.2. Element convection contribution. In each element, the convection contribution is approximated as: J c n2r~r ds = (1 - ~,)M?u~;

J c . n2~r ds = (1 - ~,)~I?u?

(26)

In these equations, up and u,G are obtained using the MAW interpolation of u C, given by equations (18)-(20); a, and a, are computed using the modified MAW scheme; and ~:/~ and ~:/~ are given by expressions similar to equation (22). It should also be noted that (1 -e,)~:/G and (1 -et).~:/o are the continuous-phase mass flow rates mr and m,, respectively, as given in equation (18). 3.4.3. Element source contribution. The volume integral involving the source terms is approximated as follows:

Ilaoc(E~'~-SDz--SG) d~//':<£lC~Z)ele'~gl(U?--ud)--(SG)l>3V'laoc

(27'

where (S~)~ and K~ are the stored nodal values associated with the element 123, and (~ao, is the volume defined by the points 1, a, o, and c: (tao, = 2n [DET[ (2rj + 2ro + .30

ra + rc)

(28)

where DET = xl Y2 + x2Y3 + x3yl

-- Yl x2 -

y2x3 - y3xt

(29)

3.4.4. Element unsteady contribution. A fully-implicit unsteady formulation [12] is used in this work. It should be noted here that in the proposed method, the solution of steady-state problems is also obtained through the use of a fully-implicit unsteady formulation. This approach is related to the solution of the coupled, non-linear, steady-state equations using iterative methods with under-relaxation [12]. In the two-fluid model considered in this work, there are two sets of governing equations, one related to the continuous phase and one to the particulate phase. In the iterative solution of the steady-state equations, it is often necessary to prescribe different relaxation factors for each of the two sets of discretization equations, in order to ensure similar evolution of the corresponding solutions. In the unsteady formulation, however, only the time step is needed to be given. This time step will naturally ensure a common evolution of the solutions for each set of equations. In the choice of a suitable time step, guidance may be obtained from physical quantities such as the value of the particle relaxation time [33]. The volume integral involving the unsteady term is approximated as follows: f,

c~ PG~l:la°' (CU~--E*U~*) aoc~l (EPGUC)d~/" = At

(30)

where At is the time step, and E* and u~* are the values of E and u G, respectively, at the previous time step. 3.4.5. Discretized equations. Finally, the discretized momentum equations are derived and assembled using element-by-element procedures akin to those used in FEMs. The resulting u G and

1082

c. MASSONand B. R. BALIGA

v G discretization equations for the continuous phase at a node i can be cast in the following general forms:

(acT+Ki~¢~+

~,

,at

( ac ~'+ Ki le'c~+ eiOG~_~

+K~ud~v+

~ 6V.b ] V? = ~., aC,b

b" + K i v ~ v

(~-P)(31)

A-----f-

+ e*P°"F'~VV?*A~ + e i C " ( - ~ y )

(32)

The bar over the pressure gradients denotes volume-averaged value associated with the control volume. The contribution of the momentum-interaction term is stated explicitly in these equations in order to clearly represent the coupling between the momentum equations of the two phases. Further manipulations of the unsteady contributions are needed in the calculation of mass flow rates [see equations (36)-(37)]; therefore, they are also stated explicitly in the previous equations.

3.5. Discretized equations for p and The continuous- and particulate-phase discretized equations, which are employed volume concentration, 0¢, respectively; (` is The integral mass conservation equation volume surrounding node i in Fig. 2, can

(`pGVmG"n2~r ds +

continuity equations are used to obtain two sets of to compute the continuous-phase pressure, p, and the computed using c¢ + (` = 1. for the continuous phase, when applied to the control be written as follows:

(`pGV~" n2~r ds +

o,.-~(ep°) d ~

+ [similar contributions from other elements surrounding node i] + [boundary contributions if applicable] = 0

(33)

The volume integral involving the unsteady term in equation (33) is approximated using equation (30) with the following specialization: u~ = 1 and utc* = 1. In each element, the mass flux, (`p GVm, G can be expressed in terms of its components in the z and r directions:

ep°V~ = (`pGu~i+ ep%~

(34)

Interpolation functions for u~ and v~ have to be first prescribed in order to approximate the mass flux contribution at the integration points on the control-volume faces within the elements. To facilitate this, equations (31) and (32) are rewritten in the following manner: u? = c,i +

(~)

v?=#i+d~(-~--y)

(35,

where , c b" ,Y___,ac.bu,,b+ ui = .b

P °"F'c" ((`*u?* - (`,u?)

ac~ + K/U~v

;

d~'-

v G p G,,F-cv ~.acnbVnb + b v + K~Vd3V'¢.~+ ~ ((`* V?* -- (.iV?) ~,=.b ; d~=

ac~ + KiC¢v

(`///'c"

ac~ + K i ~

(36)

¢i.,t/.c.~

ac~'+ K / F .

(37)

Within each element, for the evaluation of the mass fluxes at the integration points r, s, and t (Fig. 3) on the faces a-o, b-o, and o-c, respectively, the velocity components are interpolated using the recommendation of Prakash and Patankar [5] and Saabas [6]:

Control-Volume Finite Element Method

1083

Integration-point values of ~, 3, d u and d ~ are calculated using linear interpolation of the corresponding values at the vertices of the element. Similar mass-flux interpolation schemes have been used by Peric et al. [39] and Rice and Schnipke [40] in the formulation of co-located finite volume and finite element methods, respectively. As was stated earlier, this special interpolation of u~ and v~ prevents the occurrence of spurious checkerboard-type distributions [5, 6, 39, 40]. It should be noted, however, that this interpolation procedure may not be well-suited for Newtontype linearizations of the convective terms in the momentum equations. In the derivation of algebraic approximations to integrals of mass flow rates in equation (33), u~ and v~ are interpolated in each element by the functions given in equation (38). The same functions are also used to approximate the terms that represent mass fluxes in the momentum equations. Using these interpolation functions, the contributions of element 123 (Fig. 3) to the mass flow rates that appear in equation (33) can be approximated as follows:

fa

'epCV~.n2nr

ds - 1 -~ ~, {(xa)[pGVCo+p c v~]-(ya)[p G CUaC + p Gua]} G

f o CP CVCm " n2nr ds = 1 -2 ~t, {(Yc)[P° CUo+ P G

GU¢]_(X¢)[p C GVo+pGVcC]}

(39) (40)

where Uo, c vo, G u~, G va, G UcG and Vc~ are obtained using equation (38). Note that the subscript, m, has been dropped on the right-hand side of equations (39) and (40), in order to simplify the notation. ct, and 0it are obtained using the modified MAW scheme. Addition of the contributions of other elements surrounding the point i yields the complete discretized pressure equation. The complete discretized equation for the particulate-phase volume concentration, at, is derived from the particulate-phase continuity equation in a very similar manner: the modified MAW is used to interpolate ~t, and the velocity components, u~ and vm, d are interpolated linearly in each element. Compact representations of the discretized p- and ~t-equations for a typical node i are the following: acfpi = ~ ac~bPnb + b p

(41)

nb a c ~ x i = ~ aC~b~Xnb+ b ~ nb

(42)

Carver [13] suggests subtraction of the continuous-phase continuity equation from the particulate-phase continuity equation to derive a discretization equation for ~, and an addition of these equations in the derivation of the discretization equation for p, so as to explicitly account for the coupling between the phases in the calculation of ~t and p. This treatment is appropriate only when local mass conservation is ensured over each control-volume for each phase, individually, as in the finite volume method used by Carver [13]. In the proposed co-located equal-order CVFEM, for problems that involve inflows and outflows, ~ is prescribed at all nodes located on the inflow boundaries, and p is prescribed at one (or more) node(s) located at the outflow boundaries: thus, for the control volumes surrounding the nodes on the inflow boundaries, local mass conservation of the particulate phase is not necessarily satisfied; and local mass conservation of the continuous phase is not necessarily respected for the control volumes associated with the nodes on the outflow boundaries. Thus, at nodes on the inflow and outflow boundaries, the treatment proposed by Carver [13] could not be incorporated into the proposed CVFEM. Instead, as was described earlier in this section, the discretization equations for ~t is obtained from the continuity equation for the particulate phase; and only the continuous-phase continuity equation is used to derive the discretization equations for p. Therefore, the coupling between the two phases is not directly accounted for in the calculation of ~t and p, but this did not lead to any major difficulties. It should also be noted that the proposed discretization equations for ~t do not need any special treatment to ensure physically realistic solutions, while the linear combination approach suggested by Carver [13] does: the b ~ of the proposed discretized concentration equation [equation (42)] is always equal to or greater than zero and so are all the coefficients in this equation. This feature ensures that

1084

c. MASSONand B. R. BALIGA

/> 0; and since ~ ,~ 1 in the dilute gas-solid particle flows of interest in this work, there is no need to incorporate any special procedures to ensure that ~ < 1 during the iterative solution procedure. Another important feature of the formation proposed in this paper is that in equations (36) and (37), the unsteady contribution is all included in the numerator of fi and t~: thus, when the steady-state solution is reached, the unsteady contribution is zero, and the mass flux interpolation is independent o f the time step. Each of these variables can be defined in such a way that parts of the unsteady terms appear in the numerators of fi and b, and the other parts are included in the denominators of fi, d u, f, and d ~'. This would be inappropriate since the mass flux interpolation would then depend on the time step even under steady-state conditions. 3.6. Boundary conditions In this work, it is assumed that the domain boundaries remain at fixed spatial locations, and they could coincide with rigid walls, symmetry surfaces, or fluid inlet and outlet regions. At rigid walls, the intrinsic volume-average velocity of the continuous phase approaches the velocity of the wall, and, therefore, the no-slip condition can be applied. The solid particles will, in general, move and roll relative to the wall. In the present analysis, it is assumed that the particles collide with the wall through specular elastic collisions: this is equivalent to the assumption of the impermeable and slip (zero shear) conditions at the wall for the particulate phase [41]. There is no outflow or inflow of ~ at rigid walls. Again, these assumptions are adequate in the context of the stated goals of this work. More sophisticated wall boundary conditions, such as those discussed by Jiang et al. [22], can be incorporated into the proposed method, but this is outside the scope of this paper. At symmetry surfaces, the normal velocity component and the normal gradient of the tangential velocity component are both zero, for both the continuous and the particulate phases. There is no outflow or inflow of ~ at symmetry surfaces. At inflow boundaries, u ~, v c, u a, v a, and ~ are all prescribed. Special treatments are needed for ri, 3, d u and d ~ on boundaries with prescribed velocities, such as walls and inflow boundaries. At nodes which lie on such boundaries, d u and d ~' are set to zero, and, therefore, fi=Usp

t;=Vsp

(43)

At outflow boundaries, the treatment recommended by Patankar [12] is used: it is assumed that convection is the dominant transport process and diffusion is negligible. This is handled by setting the diffusion flux equal to zero at the outflow boundaries. As has been discussed by Patankar [12], in incompressible flows, pressure is a relative variable: its absolute value is not significant; only differences in pressure are meaningful. Thus it is necessary to specify only one arbitrary value of pressure at any convenient node in the calculation domain. In the proposed method, the value of pressure was fixed to a convenient, or suitable, value at one (or more, if appropriate) node(s) on the outflow boundary. 3. 7. Solution o f the discret&ation equations The discretization equations, equations (31), (32), (41), and (42), form sets of coupled nonlinear algebraic equations. In this work, at each time step, a modified version of an iterative variable adjustment procedure proposed by Saabas [6] was used to solve these equations: 1. Start with guessed or available velocity and concentration fields. 2. Calculate coefficients in the discretized unsteady momentum equations without including contributions of the pressure-gradient terms. 3. Calculate ric, ~c, du~, d,~, rid, t3d, d ~ and d v~. 4. Calculate coefficients in the discretized equations for pressure, and solve these equations to obtain updated values of p. 5. Add contributions of the pressure-gradient terms to the appropriate coefficients of the discretized z-momentum equations calculated in Step 2, and solve for u G and u d, simultaneously. 6. Add contributions of the pressure-gradient terms to the appropriate coefficients of the discretized r-momentum equations calculated in Step 2, and solve for v a and v d, simultaneously. 7. Calculate coefficients in the discretized equations for ~, and solve these equations to obtain updated values of ~.

Control-Volume Finite Element Method

1085

Table I. Examples of Co Red relations

Expression

Type

Ca Red

Validity

Stoke's Solution [ 4 3 ] Oseen's Solution [ 4 3 ] Wallis' Expression [ 4 4 ]

Theoretical Theoretical Empirical

24 24(1 + ~ Red) 24(I +-~ Red~7)

Red < I Red < 5 Red < 1000

Ref. [II]

24(I + ~ Re ez'3)

8. Calculate E ( = 1 - ~). 9. Return to Step 2, and repeat until appropriate convergence criteria are satisfied. If the unsteady formulation is used only to facilitate the solution of steady-state problems, then it is not necessary to do Step 9 in this procedure. Rather, Steps 1-8 are repeated until steady-state conditions prevail. In this work, the solution was considered to be converged when the non-dimensional average, absolute, residue for each set of discretization equations was less than 10 -1° . Depending on the problem, global values, such as the separator efficiency, were also monitored, and it was stipulated that the absolute value of the relative change from one iteration to the next (or from one time step to the next in steady-state problems) should be less than 1 0 - 6 : in most of the calculations, however, it was found that the convergence criterion based on the non-dimensional average residue is the more demanding one. In this work, in order to faciliate the implementation and testing of the proposed CVFEM, structured grids were used: the nodes in the finite element mesh lie along nonorthogonal lines that allow (/, J) indexing. In Steps 4 and 7, a line Gauss-Seidel algorithm based on the tri-diagonal matrix algorithm [12] was used to solve the discretization equations for p and ct, respectively. In Steps 5 and 6, a line Gauss-Seidel method based on a coup&d-equation line solver [12, 42] was used. The coupled solution of u ~ and u d, and v C and v d, is crucial to obtain a robust algorithm. These ideas could be extended to unstructured grid implementations by employing point, rather than line, Gauss-Seidel methods. More sophisticated approaches, such as direct, or block-decomposition, techniques based on sparse matrix algorithms may also prove viable for use with unstructured grids. To improve the rate of convergence, block-correction procedures and multigrid techniques, for example, can be included. In this work, however, these options were not considered. 4. A P P L I C A T I O N S Some of the capabilities of the proposed CVFEM are illustrated in this section by its application to two test problems and one demonstration problem. The results obtained for the first test problem are compared with a proposed analytical solution. In the second test problem, the CVFEM results are compared with those of another numerical investigation [11] and with results obtained using a well-established finite-volume staggered grid method [12, 13]. In the demonstration problem, the proposed CVFEM is used to simulate laminar, 2-D, axisymmetric, dilute gas-solid particle flows in a geometry that is characteristic of a split-flow inertial separator [14]. Before the presentation and discussion of the results, it is instructive to identify some of the parameters typically involved in the test and demonstration problems. One parameter is the Reynolds number, Re = p Cu~ L//~: L is a characteristic length of interest; and u~ is a characteristic velocity of the continuous phase. The ratio of the densities of the particulate and continuous phases, 7 = P d/pC, is also a parameter. Another important parameter is the ratio of the characteristic times, z G and z d, for the continuous and particulate phases, respectively. This parameter is called the Stokes number, Sk = zo/T c. The continuous-phase characteristic time is given by z C = L / u ~ . The particulate-phase characteristic time is taken to be the particle relaxation time given by [33]: zd --

4 pdd2 3/~CD Re d

(44)

The expression for Co Re d is typically obtained from theoretical or experimental studies on a single spherical particle. Table 1 gives some well-known expressions. 4.1. Specified solution in a square domain 4. I.I. Problem statement. The first test case has no real physical significance. It was used only to verify the formulation and implementation of the proposed CVFEM. The procedure consists

1086

C. MASSON and B. R. BALIGA

first to propose a concentration distribution, a mass conserving velocity field for each phase, and a pressure field. This ensures that the continuity equations are satisfied. In this test, the calculation domain is a square enclosure of side L. The following steady-state solution was proposed: U d ~ U d°

V d ~ U d°

uG = u C°y

vG

O~ ~ O~°

0 =

(45)

=p °x + y

--i--

p

where u d°, ct°, u °°, and pO are prescribed constants and pd/pC = 103. The solution is written in terms of a Cartesian (x, y) coordinate system. This proposed solution satisfies the momentum equations only for the following non-zero volumetric source terms: Sxo _ - p ° ( 1 - ~ ° ) + S yC = p ° ( 1

K

(y) u G°

- u d°

(46)

- 0t °) - K u do

(47)

f

S d = p°ot° + K ~u do - u G o Y-L' ~)

(48)

Say = p°ot° + Ku d°

(49)

This dilute gas-solid particle flow is governed by equations (1)-(8), with the volumetric source terms given by equations (46)-(49). The continuous-phase boundary conditions are: (i) Couette velocity profile at the inlet plane x = 0; (ii) outflow treatment [12] at the outlet plane x = L; and (iii) given x- and y-components of velocity at planes y = 0 and y = L. The particulate-phase boundary conditions are: (i) uniform inlet velocity and concentration profiles at the inlet planes x = 0 and y -- 0; and (ii) outflow treatment at the outlet planes x = L and y = L. 4.1.2. Results. The simulations of this problem using the proposed CVFEM, with coarse and fine grids, give the exact solutions [equation (45)]. This behaviour is expected since with this specified solution, the interpolation functions used in the proposed numerical methods give exact values of the various fluxes and sources. Nevertheless, this test problem and these successful simulations were very useful: they clearly indicated the validity of the implementation and the capability of the proposed CVFEM to solve the mathematical model of dilute gas-solid particle flows given by equations (1)-(8). 4.2. Flow in a channel with a restriction 4.2.1. Problem statement. In the second test problem, steady, laminar, dilute gas-solid particle

flow in a channel with a restriction is investigated. The influence of gravity is considered negligible. This problem is similar to that proposed and analyzed by Di Giacinto et al. [11]. A schematic

I_

~L/2_~

6L

~--

_1

L/6

-I

(

2L/9

L

5, Fig. 4. Geometry of the restriction.

Control-Volume Finite Element Method

1087

Table2. Valuesofparametersforflowin a channel with a restriction Re

Sk

=i.

100 100 100

10 -2 10 -2 10 -I

10 -3 5 x 10 -3 5 x 10 -3

1000 1000 1000

illustration of the problem is given in Fig. 4. The boundary conditions are: (i) uniform inlet profile for =; (ii) identical Poiseuille inlet velocity profiles for both phases; (iii) outflow treatment at the outlet plane; (iv) no-slip condition at the walls for the continuous phase; and (v) slip condition at the walls for the particulate phase. The nondimensional parameters considered in this problem are given in Table 2. In the CVFEM simulations, only one-half channel was modelled, assuming that the symmetry condition applies at the centerline. 4.2.2. Results. Three simulations are presented for this test problem (see Table 2). They have been selected to illustrate the effects of the Stokes number, Sk, and the inlet concentration, =i,, on the flow behaviour. All of the presented results correspond to a Reynolds number of 100, with the channel height (Fig. 4) as the characteristic length L, and continuous-phase inlet centerline velocity, G ud, as the characteristic velocity u~. The grid used had 73 x 37 nodes; the results were considered essentially grid independent since calculation on finer grids produced negligible changes in the centerline velocities. Figures 5-8 present comparisons of the centerline variations of u C, u d, p, and computed using the proposed CVFEM and a finite-volume staggered-grid method [12, 13]. Where applicable, the results of Di Giacinto et al. [1 I] are also shown in these figures. As can be observed in Figs 5-8, the proposed CVFEM and finite-volume/staggered-grid solutions are in good agreement. Quantitative comparisons with the results of Di Giacinto et al. [11] are not very good, but there is agreement in the qualitative behavior of these solutions as a function of Sk and 0tin. For Sk = 1 0 - 2 and ¢in = 1 0 - 3 , the finite-volume/staggered-grid solution using a 73 x 37 grid is in close agreement with the solution proposed by Di Giacinto et al. [11], while, in the other cases, these solutions are only qualitatively similar. However, for Sk = 10-~ and ~i, = 5 x 10 -3, our finite-volume/staggered-grid solution on a coarser grid (37 x 10) was found to be in good agreement with the solution of Di Giacinto et aL [11]. This indicates that some of the results of Di Giacinto et al. [11] may not be grid independent. For Sk = 10 -2, the particulate phase is almost in equilibrium with the continuous phase: the results can therefore be analyzed using a simple homogeneous two-phase flow theory [44], where suitable average velocity and density are determined to treat the two phases as a pseudofluid that obeys the usual equations of single-phase flow with the average properties. The average viscosity

Sk

1.6

CVFEM o o o o o Stag. Grid ,~_ ***** Ref. [ 1 1 ] /~-~E"n - - - CVFEM ~"~x"a ooooo Stag. Grid . T~"~-.'a .... CVFEM • 0 \ ~-L'~ Aa,Aa,a, Stag. Grid ~ ='~.00000 Ref. [ 1 1 ]

-

1.51.4-

v~

m

'-1

*~,.o 1.1 1.0

10-

o o~'<

¢t~

~

~

0

0

0

T

I

I

0

1

2

Fig. 5. Centerline

I

.3

z/L

5xl

0-3 0 --~

~'~:'~.

0

=

~tn 10 -s

10 -1 5 x l

-~.~.<....

1.31.2-

10 -2

I

I

I

4-

5

6

continuous-phasevelocity--73 x 37 grid,

1088

C. MASSONand B. R. BALIGA Sk

1.6 1.5 -o'g

n m

.".n "A- .-'~ "-'f].

1 .4

~-,,,1

- CVFEM o o o o o S t a g . Grid --CVFEM o o o o o S t a g . Grid .... CVFEM Z~Z~AAAS t a g . Grid

eqn

10 -z 10 -3 2 105x10-" 10 -I

5 x 1 0 -~

.3

3

1 .2 1.1 1.0

-j m

TO

I

I

1

2

3

I

I

I

4

5

6

z/L Fig. 6. Centerlineparticulate-phase velocity--73-37 grid. is approximately equal to the gas viscosity for very low solid concentration, and the average density is given by pm=

EpG

+ ~pd

(50)

Therefore, the effects of increasing the solid concentration, with all other conditions unchanged, is similar to the increase of the Reynolds number in single-phase flow when ~ > 1. This is exactly the behaviour of the proposed solutions: as the Reynolds number increases, the recirculating zone behind the restriction becomes larger; recovery of the continuous-phase velocity profile takes place over a larger axial length; and the velocity profiles in the core region become flatter, which results in lower centerline velocities. One of the effects of the inlet volume concentration, c(i,, can be seen by analyzing the results for Sk = l0 -2. As ~i. is increased, the entrance and total pressure drops increase to compensate the augmentation of the drag due to the presence of a larger number of particles.

As the Stokes number increases, the particulate phase is no longer in equilibrium with the Sk - CVFEM o o o o o S t a g . _ Grid

10 -2

¢rtt*,#~. Ref. [ 1 1 J -- - - CVFEM 1.0

-'~

0.5-

10 -z 5 x 1 0 -3 0 -1 1 5 x 10 -3

,-,or, u,-, S t a a . Grid .... CVFIZM Z~AA,,Z~ S t a g . Grid

"

00000Ref.

a.11.()_3

[11J

0.0-0.5

-

=

-1.0

-

C~

-1.5

-

-2.0

-

-2.5 -3.0

-

(3_

0

~,~<>

-.3.5 0

B""

er""

_&.-

.a.~ ja¢'- "- -

I

I

I

I

I

I

1

2

3

4

5

6

Fig. 7. Centerline

z/L

pressure---73

x 37 grid.

Control-Volume Finite Element Method

.E

1.10 1.09

-q -

1.08

-

1.07

-

1.06

-

1.05 1.04 1.05

-

1.02

-

1.01

i/

iI

A

A

A

&

/~

A

1089

A

A

A

A

I

~n

Sk - ooooo

I

J

CVFEM S t a g . Grid

10 -z

~ln 10 -x

CVFEM

1 0 -2

5 x 1 0 -3

oooooStag. ....

Grid 10 -1 5 x 1 0 -5 Grid

CVFEM

z~nz~n S t a g g e r e d

-~H I

1.00 0

I

I

I

I

I

I

1

2

5

4

5

6

Fig. 8. Variation of

Or/Grin

z/L

along the centerline---73 x 37 grid.

continuous phase, which can be seen by a larger difference in u G and u d in the region of the restriction. The centerline value of ct at the outlet reaches values bigger than the inlet concentration, ~i,, which indicates an accumulation of particles along the centerline of the channel downstream of the restriction. 4.3. Split-flow inertial separator

This test problem is presented to illustrate the capability of the proposed numerical method to solve problems that involve flows in complex geometries. The proposed CVFEM is used in this section to simulate steady flows in an idealized split-flow inertial particle separator. Such separators are usually installed at the inlet of helicopter gas-turbine engines in order to prevent sand and foreign-object ingestion. Such ingestion is responsible for a large proportion of early damage and unscheduled maintenance. The main advantage of such separators is their low-maintenance requirements [14]. In this problem, a dilute gas-solid particle flow in an idealized inertial separator is investigated. The idealization comes from the assumption of zero swirl in the flow and zero gravity which allows an axisymmetric analysis. In a real separator, swirl is introduced by inlet blades in order to increase the separation efficiency. Furthermore, only laminar flows are considered in this work. Thus, these simulations are not intended to model a practical separator. Rather, they are used mainly to demonstrate the capabilities of the proposed CVFEM. It should be noted, however, that this idealized problem is also well-suited for a demonstration of some of the underlying physics of dilute gas-solid particle flows in this geometry. Thus, in the discussion of the results of the CVFEM simulations, an assessment of the effects of the various non-dimensional parameters on overall pressure drop and separation efficiency is included.

Re

Fig. 9. Geometry of the separator.

1090

C. MASSONand B. R. BALIGA

4.3.1. Problem statement. A schematic illustration of the idealized problem is given in Fig. 9. The equivalent axisymmetric geometry with swirling flow has been analyzed in the past [45], using a finite element method to solve a one-way coupling model, in which the fluid flow was computed using Euler equations. Viscous flow analysis of a similar Cartesian separator has been realized [46] using one-way coupling and the Lagrangian formulation for the particulate phase. In this work, the two-fluid model is used. The inlet is at the left end, and there are two outlets at the right end. The goal is to deviate the particles in the bypass duct, while ensuring only particle-free gas enters the main duct. The bypass ratio, b, defined as the ratio of the continuous-phase mass flow rate through the bypass duct to the continuous-phase inlet mass flow rate, is controlled by the pressure difference between the main-duct outlet plane and the bypass outlet plane. The separator efficiency, r/er, indicates, for a given bypass ratio b, the effectiveness of a given separator. This efficiency is defined as the ratio of the particulate-phase mass flow rate through the bypass duct to the particulate-phase inlet mass flow rate. The boundary conditions are: (i) uniform inlet volume concentration, ~in; uniform inlet velocity profiles for u c, u d, and v G = u d = 0 ; (ii) outflow treatment [12] at the outlet planes, with given pressure difference between the bypass and main outlets; and (iii) no-slip condition at the walls for the continuous phase, and slip condition at the walls for the particulate phase. 4.3.2. Non-dimensional parameters. The parameters involved in the problem of interest are the Reynolds number, Re, the Stokes number, Sk, the ratio of the densities of the particulate and continuous phases, 7, the inlet volume concentration of the particulate phase, ~in, and the bypass ratio, b. There are many geometric parameters in this problem. The set chosen for this study is illustrated in Fig. 9: the ratio of the inlet external radius, Re, to the inlet internal radius, Ri, is equal to 2; the other geometric parameters may be obtained from Fig. 9, since this figure is drawn to scale. Re is based on DH, the inlet hydraulic diameter ( = 2Re - 2Ri). The particulate phase enters the separator at the same uniform velocity as that of the continuous phase. All simulations in this study were done with Re = 200. One combination of y and ~i. was considered: y = 1000 and ctm = 10 -3. Three values of the bypass ratio were studied: b = 10%, b = 20%, and b = 30%. The influence of the Stokes number was investigated by conducting simulations in the range 10 -3 ~ 0.2, the value of ~ef asymptotes to 100% with increasing values of Sk. 4.3.3.2. Static-Pressure Drop in the Main Duct The variation of nondimensionalized static-pressure drop in the main duct, Ap*, with Stokes number, Sk, is presented in Fig. 11. In this figure, the bypass ratio, b, is a parameter, and Ap* is defined as follows: A p *g = Pin,c - - Pout,M G G 2 0 . 5 p (Uin)

where Pi.,c is the static pressure at the central node in the inlet plane of the separator; Pout.u is the static pressure at the outlet plane of the main duct; pc is the density of the continuous phase; and

1091

Control-Volume Finite Element Method

1009O 80 c:c:ob 7O 6O 5t 50 ,4-. ~D 4O 30-' 20" 10 I I 0 0.001

= 10~

I

I IIII

p

I

I

I

I

0.01

! I[11

I

i

I

l lli~

i

I

0.1

Sk Fig. 10. Variation of separator efficiency with Stokes number. I

G is the prescribed uniform velocity of the continuous phase at the inlet plane of the separator. In all simulations, Pin,c was essentially equal to the area-weighted average of the static pressure at the inlet plane. Also presented in Fig. 11 are results pertaining to flow, in the same separator, of a homogeneous mixture (line with long dashes) and a single-phase fluid (line with short dashes). The density and viscosity of the homogeneous mixture were obtained using Pmix = ~inP d Jr £inPG;

~/mix= ~ ( l Jr 2.50~in ) ~,~

These equations are quite appropriate for the dilute gas-solid particle flows considered in this study. For the single-phase fluid flow, the density and viscosity were set equal to those of the continuous phase in the gas-solid particle flow. The results in Fig. 11 show that Ap* decreases as b increases, for a fixed Sk. This is to be expected because as b increases, the amount of fluid flow in the bypass duct increases and that through the main duct decreases. For a fixed value of b, Ap* asymptotes to the homogeneous-mixture solution

50-

40-

10~ 20~

o_ 30 <::l 20

10

0.001

20~

i

i

I

II,lll

I

]

I

0.01

llllll

I

0.1

I

I

llllll

1

Sk Fig. I 1. Variation of static-pressure drop in the main duct with Stokes number.

1092

C. MASSONand B. R. BAL1GA

10-

. . . . .

41'

>m tD..

~

~ ~

--30~

........................................

30~

6

-2o~


---~- . . . . . . . . . . . . . . . . . .

2

i

0.001

i i i itll[

i

I

i I lllJ I

0.01

l

l

10~

i ~II11 i

0.1

1

Sk Fig. 12. Variation of static-pressure drop in the bypass duct with Stokes number. as Sk decreases below a value of 0.01. This is also an expected result: as Sk decreases, the particulate and continuous phases move towards equilibrium conditions with one another; and at very small values of Sk, two-fluid models of dilute gas-particle flows become equivalent to the homogeneousmixture model. At large values of Sk (~>0.06), Ap* values given by the two-fluid model fall below those produced by the homogeneous flow model. This is because the separator efficiency increases with increasing Sk; and for Sk/> 0.06, the concentration of particles in the main duct beyond the split-off point is significantly below the uniform concentration (=~n) that is assumed to prevail in the homogeneous flow model. The results in Fig. 11 also show that for a fixed value of b, starting from Sk = 0.001, as Sk increases, Ap* first increases until it reaches a maximum, and then it decreases monotonically. This is due to the opposite variations with Sk of the size of the particles (or particle inertia) and the particle number density (or particle-gas contact surface area). For a fixed 0ti,, at low values of Sk, there is a large amount of small particles; high values of Sk, on the other hand, correspond to a relatively smaller amount of large particles. High values of Sk are also accompanied by high values of r/el, and this further reduces the particle number density in the main duct beyond the split-off point. Increases in size and number density of the particles are both accompanied by increases in the drag force exerted by the particles on the gas, or vice versa. The results in Fig. 11 show that for 0.001 ~< Sk ~<0.02, the influence of increasing particle size dominates the opposite influence of decreasing particle number density. The reverse is true for Sk/> 0.02. In the range of parameters considered, the Ap* values for the gas-solid particle flow are all larger than the corresponding values for the single-phase fluid flow, but this difference decreases at high values of Sk. This is also expected, because as Sk increases, for fixed ~m and b, r/or increases, so less particles flow through the main duct beyond the split-off point.

4.3.3.3. Static-Pressure Drop in the Bypass Duct The variation of nondimensionalized static-pressure drop in the bypass duct, Ap*v, with Stokes number, Sk, is presented in Fig. 12. In this figure, the bypass ratio, b, is a parameter, and Ap*v is defined as follows Pin,c Pout,BY Ap*v = 0.5pG(u~)2 -

-

where Po.t,Bv is the static pressure at the outlet plane of the bypass duct; and the physical meanings ofp~,.,, pC, and uiG,are the same as those given earlier. Also presented in Fig. 12 are results obtained

Control-Volume Finite Element Method

1093

with a homogeneous-mixture model (line with long dashes) and results pertaining to a single-phase fluid flow (line with short dashes) in the same separator. The results in Fig. 12 show that Ap*y increases as b increases, for a fixed Sk. This is expected because with an increase in b, the flow through the bypass duct increases. For a fixed value of b, Ap*y asymptotes to the homogeneous-mixture solution as Sk decreases below a value of 0.01. The explanation for this is the same as that provided earlier to explain this feature of Ap* results. For Sk/> 0.01, the variation of Ap~'y with Sk is different for different values orb. These differences are caused by the combined influences of the relative speed of the particles with respect to the gas at the inlet-plane of the bypass duct (Fig. 9); the particle concentration distribution in the bypass duct; the particle number density in the bypass duct; and the particle size (or inertia). At b = 10%, Ap*v decreases with increases in Sk, indicating that the dominant influence is that of the corresponding decrease in particle number density: indeed, for Sk/> 0.04, the two-fluid-model Ap*y values are lower than those obtained for single-phase fluid flow. For these particular cases, detailed examinations of the continuous- and particulate-phase speeds, vG=JuG2"[-U G2 and Vd = ~ v d2, (see Fig. 13) showed that the former are significantly larger than the latter at the inlet-plane of the bypass duct, just after the split-off point: so the particles facilitate the fluid flow in the bypass duct, rather than impose a drag on it. For b = 20%, Ap~v decreases with increases in Sk, but reaches a minimum value in the vicinity of Sk = 0.08: this is because the decrease in particle number density caused by increasing Sk (at a fixed 0ti,) is balanced at this point by the increase in this variable produced by increases in nef. At b = 30%, the influences of the (~/O(In

0

1 = I

r

,,

2 [ I

I

(X/~ln

3 I

~"~'i " ~ t

4 = I

5

0

1

I

2

3

4

5

f/"

r " ~

Bypass

Bypass Outlet

I'XX~

Inlet

~* /

",/,

.

-. _

I

i

I

• ,/UGil VTu,.

-.y-- o o,.

I

I

I

=

I

I

/

=

I

I

'

I

'

T

I

'

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

v/u .

v/#,.

(o) 0 r

10

~/QIn

O{/C(In

20

30

0

I I I I I I I I I I I I I I

, ~L...

"~,.~

Bypass Inlet

10

20

30

I I I I I I I I i I t I I

r ~'~,,~

lyuPfl~S

',,. / i

a

........ ]yo,.

0.0

J l l ~ l l J l l 0.2 0.4 0.6 0.8

i

I

i

I

I

I

'

I

I

1.0 O.O 0.2 0.4 0.6 0.8 1.0

(b) Fig. 13. Profiles of V°/u9m, vd#d G and =/=in at the bypass inlet and outlet planes, b = 10%: (a) S k -- 2 x l 0 -2, (b) S k -- 6 x 10-2. CAF 23/8--H

1094

C. MASSONand B. R. BALIGA 0 i

O(/O(In

O{/O(in

1

2

I

I

3

i

I

4 i

r ~

I

5

1

i

I

Bypass Inlet

' i ~

0

2

I

r ~

i

3

I

i

I

I

i

I

i

0.8

I

I

1.2

I

i

1.6 2.0 0.0

v/u?.

i

~

I

I

_

I

5 i

Bypass Outlet

"""'--

~

/~ 0.0 0.4

4

I

_

vTu,.I

I

0.4 0.8

1.2 1.6 2.0

v/u?.

(<:,) O{JOiln 0 r

25

(X//O~ln

50

75 1O0 llllllJlllla IILLII --~'Z'~-~ _ "~- - - Bypass

0

20

10

r ~ ,,'

\\

Va/u?. \ \

i--

VT,u.;"

t ..... i

I

i

I

I

/

i

I

i

40

Bypass Outlet

" ~

I; 0.0 0.4 0.8

30

iaIIIIIIilIililllll

w

I

<
I

1.2 1.6 2.0 0.0 0.4 0.8

i

I

1.2

i-I-i

1.6 2.0

v/u?.

v/u?.

(b) Fig. 14. Profiles of V°lu~, Va/u~ and e/ein at the bypass inlet and outlet planes, b = 30%: (a) Sk=2x 10-2;(b) Sk=2x 10-I. aforementioned factors on Ap*v balance each other out very well for 10-3~< Sk ~<0.04, so Ap*v is essentially constant in this range of Sk. However, for Sk > 0.04, increases in Sk lead to substantial increases in Ap*v. For these particular cases, detailed examinations of the particulate- and continuous-phase speeds, and volume concentrations (see Fig. 14) showed the following: (1) at the inlet-plane of the bypass duct, just after the split-offpoint, the concentration of the particulate-phase is very high in the low-velocity region adjacent to the upper wall; (2) downstream of this point, these low-velocity particles move to the central regions of the bypass duct, mainly because of the geometry of this particular separator; (3) this, in turn, causes a considerable drag on the continuous phase, and requires high values of Ap*v. 5. C O N C L U S I O N In this paper, some of the ideas contained in available CVFEMs for single-phase flows [1-8] have been amalgamated and extended to formulate a co-located, equal-order, CVFEM for the solution of two-fluid models of 2-D, planar or axisymmetric, incompressible, dilute gas-solid particle flows in regular- or irregular-shaped geometries. The various steps in the formulation of this CVFEM have been presented concisely, with emphasis and discussion of features that are particularly relevant, or necessary, in the simulation of dilute gas-solid particle flows. In order to focus attention on the formulation of the proposed CVFEM, the discussions in this paper are presented in the context of a simple two-fluid model that is a specialization, to dilute gas-solid particle flows, of

Control-Volume Finite Element Method

1095

a general m o d e l p r o p o s e d b y A n d e r s o n a n d J a c k s o n [16]. Reference has been m a d e to m o r e s o p h i s t i c a t e d a n d r i g o r o u s m o d e l s in the literature [15-22], b u t a c o m p a r a t i v e e v a l u a t i o n o f these m o d e l s is n o t within the scope o f this paper. N o v e l features o f the p r o p o s e d C V F E M vis-a-vis available C V F E M s [1-8] include the following: a f o r m u l a t i o n t h a t can h a n d l e p l a n a r o r a x i s y m m e t r i c 2-D flows; a modified M A s s - W e i g h t e d skew u p w i n d scheme ( M A W ) for the i n t e r p o l a t i o n o f the v o l u m e c o n c e n t r a t i o n , ~, o f the p a r t i c u l a t e phase; a p p r o p r i a t e i n c o r p o r a t i o n o f the discretized u n s t e a d y terms in the definitions o f ~G a n d ~3°; a n d an iterative variable a d j u s t m e n t a l g o r i t h m in which linearized discretization e q u a t i o n s for u G a n d u d, a n d v G a n d v d, are solved simultaneously, using a line G a u s s - S e i d e l a l g o r i t h m b a s e d on a c o u p l e d - e q u a t i o n line solver. In a d d i t i o n , the sets o f discretized e q u a t i o n s for p a n d ct are derived f r o m the integral mass c o n s e r v a t i o n e q u a t i o n s for the c o n t i n u o u s a n d p a r t i c u l a t e phases, respectively, r a t h e r t h a n a linear c o m b i n a t i o n o f these equations, as is d o n e in some available finite v o l u m e m e t h o d s [13, 24]. T h e reasons for a d o p t i n g this a p p r o a c h , a n d the features o f the p r o p o s e d C V F E M t h a t allow this a p p r o a c h w i t h o u t a d d i t i o n a l special treatments, have been discussed in this p a p e r . T h e capabilities o f the p r o p o s e d C V F E M have been illustrated in this p a p e r by its a p p l i c a t i o n to two test p r o b l e m s a n d one d e m o n s t r a t i o n p r o b l e m . This is the first a p p l i c a t i o n o f C V F E M s to p r o b l e m s involving dilute g a s - s o l i d particle flows. T h e results are quite encouraging. T h u s the p r o p o s e d C V F E M c o u l d serve as a useful tool in the c o m p a r a t i v e e v a l u a t i o n s o f available twofluid m o d e l s o f g a s - s o l i d particle flows [9, 10, 15-24, 27-32], a n d a p p l i c a t i o n s o f these m o d e l s to practical p r o b l e m s .

Acknowledgements--The support of the Natural Sciences and Engineering Research Council of Canada (NSERC), in the

form of a Postgraduate Scholarship to the first author and an operating grant to Professor Baliga, is gratefully acknowledged. The first author is also grateful to Hydro-Qurbec for financial support through a fellowship for postgraduate studies at McGill University. The research activities of Professor Baliga are also supported by the Centre de Recherche en Calcul Appliqu6 (CERCA) in Montrral, Qurbec.

REFERENCES 1. B. R. Baliga and S. V. Patankar, A control-volume finite element method for two-dimensional fluid flow and heat transfer. Numer. Heat Transfer 6, 245 (1983). 2. B. R. Baliga and S. V. Patankar, Elliptic systems: finite element method II, In Handbook of Numerical Heat Transfer (Edited by Minkowycz et al.), Chap. II, pp. 421-461. Wiley, New York (1988). 3. N. A. Hookey and B. R. Baliga, Evaluation and enhancements of some control volume finite-element methods: Part II--Fluid flow problems. Numer. Heat Transfer 14, 273 0988). 4. B. LeDain-Muir and B. R. Baliga, Solution of three-dimensional convection-diffusion problems using tetrahedral elements and flow-oriented upwind interpolation functions. Numer. Heat Transfer 9, 143 0986). 5. C. Prakash and S. V. Patankar, A control volume-based finite-element method for solving the Navier Stokes equations using equal-order velocity-pressure interpolation. Numer. Heat Transfer 8, 259 (1985). 6. H. J. Saabas, A control volume finite element method for three-dimensional, incompressible, viscous fluid flow. Ph.D. Thesis, Dept. of Mechanical Engineering, McGilt University, Montreal (1991). 7. G. E. Schneider and M. J. Raw, A skewed positive-influence coefficient upwinding procedure for control-volume-based finite element convection diffusion computation. Numer. Heat Transfer 9, l 0986). 8. G. E. Schneider and M. J. Raw, Control-volume finite-element method for heat transfer and fluid flow using colocated variables--1. Computational procedure. Numer. Heat Transfer 11, 363 (1987). 9. C. T. Crowe, REVIEW--Numerical models for dilute gas-particle flows. A S M E J. Fluid Eng. 104, 297 0982). 10. C. T. Crowe, The state-of-the-art in the development of numerical models for dispersed phase flows. Proc. Int. Conf. on Multiphase Flows '91, Tsukuba, Japan, 49-60 (1991). I I. M. Di Giacinto, F. Sabetta and R. Piva, Two-way coupling in dilute gas-particle flows. A S M E J. Fluids Eng. 104, 305 (1982). 12. S. V. Patankar, Numerical Heat Transfer and Fluid Flow. McGraw-Hill, New York (1980). 13. M. B. Carver, Numerical computation of phase separation in two-fluid flow. A S M E J. Fluids Eng. 106, 147 (1984). 14. M. L. Yaffe, T700 aims at low combat maintenance. Aviation Week Space Technol., 28 January, 45 (1974). 15. R. Ishii, Y. Umeda and M. Yuhi, Numerical analysis of gas-particle two-phase flows. J. Fluid Mech. 203, 475 0989). 16. T. B. Anderson and R. Jackson, A fluid mechanical description of fluidized beds. I&EC Fundamentals 6, 527 (1967). 17. S. Whitaker, The transport equations for multiphase systems. Chem. Eng. Sci. 2~, 139 (1973). 18. M. Ishii, Thermo-Fluid Dynamic Theory of Two Phase Flow. Eyrolles, Paris. 19. J. C. Slattery, Momentum, Energy, and Mass Transfer in Continua. Krieger, Malabar (1981). 20. D. A. Drew, Mathematical modeling of two-phase flow. Ann. Rev. Fluid Mech. 15, 261 0983). 2 I. G. H. Crapiste, E. Rotstein and S. Whitaker, General closure scheme for the method of volume averaging. Chem. Eng. Sci. 41, 227 (1986). 22. T. S. Jiang, M. H. Kim, V. J. Krcmes~ Jr. and J. C. Slattery, The local volume-average equations of motion for a suspension of non-neutrally buoyant spheres. Chem. Eng. Commun. 50, 1 0987).

1096

C. MASSON and B. R. BALIGA

23. F. Durst, D. Milojevic and B. Schnung, Eulerian and Lagrangian predictions of particulate two-phase flows: a numerical study. Appl. Math. Model. 8, 101 (1984). 24. D. B. Spalding, Numerical computation of multiphase flows. Workshop Lecture Notes, Purdue University (1978). 25. W. Shyy and T. C. Vu, On the adoption of velocity variable and grid system for fluid flow computation in curvilinear coordinates. J. Comp. Phys. 92, 82 (1991). 26. F. Moukalled and S. Acharya, Application of an adaptive grid procedure for the calculation of turbulent separated flows. Computers Fluids 21, 455 (1992). 27. S. E. Elghobashi and T. W. Abou-Arab, A two-equation turbulence model for two-phase flows. Phys. Fluids 26, 931 (1983). 28. C. P. Chen and P. E. Wood, A turbulence closure model for dilute gas-particle flows. Can. J. Chem Eng. 63, 349 (1985). 29. H. B. Stewart and B. Wendroff, Two-phase flow: models and methods. J. Comput. Phys. 56, 363 (1984). 30. G. Hetsroni, Handbook of Multiphase Systems. McGraw-Hill, New York (1982). 31. J. X. Bouillard, R. W. Lyczkowski and D. Gidaspow, Porosity distributions in a fluidized bed with an immersed obstacle. AIChE J. 35, 908 (1989). 32. B. T. Chao, W. T. Sha and S. L. Soo, On inertial coupling in dynamic equations of components in a mixture. Int. J. Multiphase Flow 4, 219 (1978). 33. G. Rudinger, Relaxation in gas-particle flow. Nonequilibrium Flows (Edited by Wegener) Chap. 3 (1969). 34. M. R. Maxey and J. J. Riley, Equation of motion for a small rigid sphere in a nonuniform flow. Phys. Fluids 26, 883 (1983). 35. O. C. Zienkiewicz, The Finite Element Method. McGraw-Hill, London (1977). 36. T. J. Barth, On unstructured grids and solvers. Computational Fluid Dynamics, Von Karman Institute Lecture Series, 1990-03 (1990). 37. C. Prakash, Examination of the upwind (donor-cell) formulation in control volume finite-element methods for fluid flow and heat transfer. Numer. Heat Transfer 11, 401 (1987). 38. A. Brooks and T. J. R. Hughes, Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comp. Meth. Appl. Mech. Eng. 32, 199 (1982). 39. M. Peric, R. Kessler and G. Scheuerer, Comparison of finite-volume numerical methods with staggered and co-located grids. Computers Fluids 16, 389 (1988). 40. J. G. Rice and R. J. Schnipke, An equal-order velocity-pressure formulation that does not exhibit spurious pressure modes. Comp. Meth. Appl. Mech. Eng 58, 135 (1986). 41. J. L. Sinclair and R. Jackson, Gas-particle flow in a vertical pipe with particle-particle interactions. A1ChE J. 35, 1473 (1989). 42. P. F. Galpin, J. P. Van Doormaal and G. D. Raithby, Solution of the incompressible mass and momentum equations by applications of a coupled equation line solver. Int. J. Numer. Meths. Fluids 5, 615 (1985). 43. H. Schlichting, Boundary-layer Theory. McGraw-Hill, New York (1979). 44. G. B. Wallis, One-dimensional Two-phase Flow. McGraw Hill, New York (1969). 45. D. S. Breitman, E. G. Dueck and W. G. Habashi, Analysis of a split-flow inertial particle separator by finite elements. J. Aircraft 22, 135 (1985). 46. M. Zedan, P. Hartman, A. Mostafa and A. Sehra, Viscous flow analysis for advanced inlet particle separators. AIAA 90-2136, 8pp. (1990).