Comparison of Eulerian and Lagrangian Monte Carlo PDF methods for turbulent diffusion flames

Comparison of Eulerian and Lagrangian Monte Carlo PDF Methods for Turbulent Diffusion Flames ¨ BUS* H. MO

Institut fu ¨r Thermodynamik der Luft- und Raumfahrt, Universita ¨t Stuttgart, Pfaffenwaldring 31, 70569 Stuttgart, Germany

and ¨ GGEMANN P. GERLINGER and D. BRU

Lehrstuhl fu ¨r Technische Thermodynamik und Transportprozesse, Universita ¨t Bayreuth, Universita ¨tsstr. 30, 95440 Bayreuth, Germany Methods based on Monte Carlo solutions of the transport equation for the joint probability density function (pdf) of scalars are considered to be among the most promising approaches for modeling turbulent reactive flows. For the treatment of convection and diffusion two approaches are well established, either from Eulerian or Lagrangian point of view. Simply speaking, the former is computationally less expensive but also less accurate than the latter one. In this paper we provide an extension of a previous analysis of numerical diffusion of Eulerian schemes to two-dimensional flows in order to shed light on the shortcomings of this approach. For convergence acceleration, local time stepping is extended to Lagrangian pdf methods. Furthermore, the adaption of the Continuous Mixing Model to unequal particle weights is presented. The differences between the Eulerian and Lagrangian approaches are analyzed by application to simple test problems and by simulation of a supersonic hydrogen-air diffusion flame using a chemical scheme consisting of seven elementary reactions. © 2001 by The Combustion Institute

INTRODUCTION One of the severest problems when simulating turbulent reacting flows by means of ensembleaveraged species transport equations is closure of the mean chemical production term. Due to its strong non-linearity, neglect of turbulent fluctuations, i.e., setting the mean source term equal to the production rate calculated from mean values, may result in significant error. This approach is termed laminar chemistry. Norris and Hsu [1] showed that source terms determined by this simplification and those deduced from a probability density function (pdf) simulation can differ strongly in magnitude and even in sign. Because chemical reaction rates exhibit no spatial dependence, this problem is easily circumvented by introducing the single-point joint pdf. If all the variables on which chemical reaction depends are included in the pdf, integrating its product with the corresponding chemical source term over the entire composition space yields the mean reaction rate. There*Corresponding author. E-mail: [email protected] COMBUSTION AND FLAME 124:519 –534 (2001) © 2001 by The Combustion Institute Published by Elsevier Science Inc.

fore, pdf methods derive their importance from the fact that chemical reactions appear in closed form in this formulation. Two distinctive classes of pdf approaches are currently in use. In assumed pdf methods [2, 3], as the name indicates, the mathematical shape is presumed so that it is only necessary to solve transport equations for a few coefficients fully describing the function (usually some higher moments of the pdf). Although the virtue of this method lies in the fact that the numerical effort increases only moderately compared to laminar chemistry methods, the limitations are equally evident: the choice of mathematical functions recognizing boundedness of mass fractions is small, the forced shape may possibly not correspond to the physical one, and to take temperature-composition correlations into account proves to be non-trivial [4]. In the more elaborate (and more expensive) approach adopted here, a transport equation for the pdf is solved. No restrictions apply on the shape of the pdf, and all correlations are fully represented. The pdf of scalars (species mass fractions and enthalpy in the present case) is investigated here which makes possible clo0010-2180/01/$–see front matter PII S0010-2180(00)00207-8

520 sure of chemical sources. On the other hand, this approach must resort to the gradient diffusion assumption for modeling turbulent transport as velocity components are not included in the probability density function. Due to the usually high dimensionality of the pdf, its transport equation is commonly solved using a Monte Carlo scheme in which the pdf is represented by stochastic particles. Of the many aspects being relevant to Monte Carlo pdf methods, transport of these computational particles in physical space is the focus of this work. Transport can either be viewed from Eulerian or Lagrangian perspective. Both approaches have been widely used in the past for simulation of practical turbulent reacting and mainly subsonic flows. For applications of the former method see, e.g., Jones and Kakhi [5], for examples of the latter Tsai and Fox [6], and for a multitude of calculations employing both schemes consult Ref. [7]. In the Eulerian approach each particle represents a fraction of the pdf of its host cell. Particles do not have specific locations within the cell. The number of particles stays constant in time, and particles have equal significance. Due to transport processes, fractions of the pdf are replaced by the pdf’s of neighbouring cells. Interpreted in the sense of Monte Carlo schemes, this implies that a number of particles of a given cell is replaced by particles from surrounding cells that are selected at random. From the Lagrangian point of view, stochastic particles move through physical space independently of each other. They are assigned spatial coordinates and represent mass. Due to the stochastic nature of motion, the number of particles present in a certain cell changes in time. As can easily be derived for constant density flows, the average number of particles per cell is proportional to the cell volume if particle masses are equal. In order to prevent particle accumulations in large cells and to keep small cells from running empty, particles are ascribed a relative weight, which may vary throughout the domain and, as a consequence of movement, among the particles within a cell. The particles move in such a way that the evolution of their pdf equals that of the physical transport equation. The stochastic particles therefore need not necessarily behave as fluid particles (however useful this picture may be).

¨ BUS ET AL. H. MO In fact each of them represents a single realization of the turbulent flow. It should again be stressed that all simulations in this paper are based on a scalar pdf description. In even more sophisticated Lagrangian methods velocity is also treated as random variable and is incorporated into the pdf. This removes the need for modeling turbulent convection by the gradient diffusion assumption at the cost of higher dimension and further unclosed terms originating from the momentum equations. This approach is envisioned as further extension of the current work. In the next section, a brief overview of the fundamentals of the pdf transport equation is given. The following section dedicated to the Eulerian approach discusses the dependence of numerical diffusion on the CFL number for convection diagonal to the grid, which leads us to the development of an improved transport operator. Subsequently, the Lagrangian scheme is described shortly with emphasis on nonstandard features such as local time stepping for convergence acceleration. Finally, the paper presents test cases that demonstrate the differences between the applied methods. PDF TRANSPORT EQUATION In methods making use of the single-point joint scalar probability density function, the variables describing the thermochemical field are considered to be random, i.e., they may assume a different value in each repetition of a specific turbulent experiment. ⌽ ⫽ (Y1, . . . , YNk⫺1, h)T is a vector consisting of these random variables (species mass fractions and enthalpy in the present compressible case), and ⌿ ⫽ (⌿ 1 , . . . , ⌿ N k⫺1 , ⌿ h ) T denotes the corresponding vector of independent variables in sample space which comprises all possible states of ⌽. N k is the number of species. The pdf P is then defined such that P(⌿; x, t)d⌿ describes the probability of the event ⌿ ⬍ ⌽ ⬍ ⌿ ⫹ d⌿ at a particular location in physical space x and at time t. In Monte Carlo pdf methods the transport equation for P is solved which can be derived from conservation laws of species masses and energy by standard techniques [8] and takes the form

EULERIAN AND LAGRANGIAN PDF METHODS

␳៮

521

˜ ˜ ˜ ⭸P ⭸ ⭸ ⭸P ⭸ ⭸P ⭸ ˜⫺ ⫹ ␳៮ u ˜k ⫹ ␳៮ 具u ⬙k兩⌿典P ␳៮ D ⫽ ⫺␳៮ ⭸t ⭸ xk ⭸ xk ⭸ xk ⭸ xk ⭸⌿ ␣ ⭸⌿ ␤ ⫺ ␳៮

冓

冏

冔册

冏

⭸⌽ ␣ ⭸⌽ ␤ ˜ ⌿ P ⭸ xk ⭸ xk

D

⭸ S␣ ⭸ ⭸p ⭸p 1 ˜ ⫺ ␳៮ ˜ P ⫹ ⌽D ⌿ P ⫹ uk ⭸⌿ ␣ ␳ ⭸⌿ h ⭸t ⭸ xk ␳

˜ (⌿; x, t) ⫽ in terms of the Favre-weighted pdf P ␳ (⌿) P(⌿; x, t)/ ␳៮ . Unclosed terms in this formulation can be identified by angular brackets 具 . . . 典 signifying expectations conditional upon the event ⌽ ⫽ ⌿. In Eq. 1 ␳ denotes density, u k the velocity component in direction k, and p the pressure. An equal diffusion coefficient D is assumed for all species. S ␣ represents the chemical source term of species ␣, ⬃ denotes Favre averages, and ⬙ fluctuations about these, whereas – identifies conventional averages. The terms on the left-hand side of Eq. 1 repre˜ , mean convection, sent the change in time of P turbulent convection, and molecular diffusion, respectively. Those on the right-hand side stand for turbulent mixing, reaction, and the last one comprises compressibility effects and dissipation. By using an operator splitting technique, Eq. 1 can be factored into ˜ 共t ⫹ ⌬t兲 ⫽ 共I ⫹ ⌬tCD兲 䡠 共I ⫹ ⌬tM兲 P ˜ 共t兲 䡠 共I ⫹ ⌬tR兲 䡠 共I ⫹ ⌬tC兲 䡠 P

(2)

allowing consecutive treatment of the different processes which are spatial exchange (CD), turbulent mixing (M), chemical reaction (R), and compressibility (C). I represents the identity operator. The process of particular concern in this paper is spatial exchange due to convection and diffusion given by ⭸ ⭸ ⭸ ˜⫹ ˜⫹ ˜ ␳៮ P ␳៮ u ˜ P ␳៮ 具u ⬙k兩⌿典P ⭸t ⭸ xk k ⭸ xk ⫺

冔

冋冓

˜ ⭸ ⭸P ␳៮ D ⫽ 0. ⭸ xk ⭸ xk

(3)

Because velocity is not included in the pdf, turbulent convection remains unclosed in this formulation. It is commonly modeled by the familiar gradient diffusion assumption [9], leading to D being replaced by an effective diffusion coefficient that is the sum Deff ⫽ D ⫹ DT of molecular and turbulent (DT) diffusion coefficients. Because of the high dimensionality of the pdf,

(1)

the transport equation must be solved by a Monte Carlo scheme in which the pdf of each cell is represented by a large number N of stochastic particles [8]

冉 冘w N

P N共⌿兲 ⫽

n⫽1

写

Nk⫺1

冊

共n兲

⫺1

䡠

冘 冋w N

n⫽1

共n兲

册

␦ 共Y 共n兲 ␣ ⫺ ⌿ ␣兲 .

␣⫽1

䡠 ␦ 共h 共n兲 ⫺ ⌿ h兲

(4)

In the Eulerian approach, all particle weights w(n) are equal. The evolution of Eq. 1 in time is simulated by suitable manipulation of these ensembles. Again, our investigation is centered on the spatial exchange of particles resulting from convection and diffusion whereas the remaining processes take place in thermochemical space. Variables that cannot be extracted from the Monte Carlo particles, such as mean velocity, turbulent diffusion coefficient or turbulent frequency, are provided by a conventional finite volume flow solver with which the pdf simulation is coupled. For more details on our TASCOM solver see Refs. [10, 11], and for the q- ␻ -turbulence model consult Ref. [12]. For description of chemistry, the so-called abridged Jachimowski scheme of Baurle et al. [2] listed in Table 1 is employed. It comprises seven species and seven elementary reactions. The applicable third body effectivities for this mechanism are given by Baurle [13]. The change in composition space by chemical reaction is determined by directly integrating the chemical source term for each particle in time by an Euler semi-implicit scheme [14]. EULERIAN METHOD The essentials of the Eulerian scheme as it is employed here are given in a previous paper [14]. Therefore, we concentrate on the comparison of Eulerian and Lagrangian approaches. There are two principal deficiencies of the Eulerian pdf method: numerical diffusion and

¨ BUS ET AL. H. MO

522 TABLE 1 Reaction scheme [2] # 1 2 3 4 5 6 7

A

Reaction H2 ⫹ O2 º OH ⫹ OH H ⫹ O2 º OH ⫹ O OH ⫹ H2 º H2O ⫹ H O ⫹ H2 º OH ⫹ H OH ⫹ OH º H2O ⫹ O H ⫹ OH ⫹ M º H2O ⫹ M H ⫹ H ⫹ M º H2 ⫹ M

13

1.70 䡠 10 1.20 䡠 1017 2.20 䡠 1013 5.06 䡠 104 6.30 䡠 1012 2.21 䡠 1022 7.30 䡠 1017

b

Ea

0.0 ⫺0.91 0.0 2.67 0.0 ⫺2.0 ⫺1.0

48000 16500 5150 6290 1090 0 0

A in multiples of cm3/(mole s), E a in cal/mole.

spatial accuracy. The serious effects of numerical diffusion when dealing with non-linear chemistry have been demonstrated before along with proof that local time stepping is able to cure this problem in a simple one-dimensional case [14]. In the following discussion, the analysis is extended to two spatial dimensions with convection diagonal to the computational grid, and the findings stimulate the development of an improved convection operator. For the second problem, the first order accuracy of the simple spatial upwind discretization of convection, no workable cure is known. (Accuracy of physical diffusion processes is less critical because it is discretized by second order central differences.) Although conventional finite difference or finite volume methods may readily be extended to second order and higher accuracy, the positivity of the number of particles does not allow the application in Monte Carlo methods. Anand et al. [15] proposed a scheme in which a blend between second order central differences and first order upwinding is used by adding spatial exchange resulting from convection and diffusion. This method only results in second order accuracy if diffusive transport is at least of the same magnitude as convective transport, but in most practical problems diffusion is much less intense. Therefore, grid resolution would have to be enhanced in order to make both processes equally strong. At this point, however, spatial accuracy is no longer of concern anyway.

to three dimensions or axisymmetric problems (neglecting higher order effects) is straight forward. The number of particles that are transported across computational cell boundaries is usually determined from a finite difference discretization of Eq. 3 [16, 17]. Here, a first order upwind finite volume method is used instead, which principally leads to the same results but is somewhat more instructive in the development that follows. For an arbitrarily formed, quadrilateral cell i, j such as the one depicted in Fig. 1, the number of particles being convected into it over cell face cf (cf ⫽ AB, BC, CD, DA) ⌬N cf ⫽

⌬m cf 䡠 N, m i, j

(5)

with ⌬m cf max 关⫺␳ cfuជ cf 䡠 nជ cf, 0兴⌬t ⫽ m i, j ␳ i, jV i, j

(6)

Improved Treatment of Diagonal Convection In the following discussion, equations are given for the planar two-dimensional case. Extension

Fig. 1. Convection from diagonal cells in 2-D. Shaded areas represent mass transported into cell ABCD within time step ⌬t.

EULERIAN AND LAGRANGIAN PDF METHODS

523

and in cell i, j ⫹ 1 with probability

␤⫽

can easily be determined by help of the finite volume technique. Here, uជ ⫽ (u, v) T denotes the velocity vector, nជ cf ⫽ (⌬y cf , ⫺⌬x cf ) T the outward-pointing normal vector on cell face cf. The increments ⌬x and ⌬y are defined in counter-clockwise direction. ⌬m cf is the mass convected into cell i, j, whereas ␳ i, j V i, j represents the mass contained in cell i, j. In this procedure only adjoining cells ((i ⫹ 1, j), . . . ) are considered while no particles enter from diagonal cells ((i ⫹ 1, j ⫹ 1), . . . ). This results in serious deficiencies when convection diagonal to the computational grid occurs. In order to further appreciate this problem, a previous analysis [14] of numerical diffusion inherent in first order upwind schemes is repeated here for the convection problem depicted in Fig. 2. The Cartesian velocity components u and v are positive and constant as are grid spacings ⌬x and ⌬y and density. In Eulerian methods, a particle’s location is discrete and completely identified by cell indices i, j. For a first order method, a particle residing in cell i, j at time t will be found in cell i ⫹ 1, j after time step ⌬t with probability

T

⫺1

u ⌬t ⌬x

t VT ⫽ 2 ␰ ⫹ ␩2

(8)

The time step is limited by ␣ ⫹ ␤ ⫽ CFL ⱕ 1. Although we only consider mean convection in this problem, because of random selection the particle’s position is random with corresponding sample space variables ˆ, ı ˆ. j For a particle being located in cell i ⫽ 0, j ⫽ 0 at time t ⫽ 0 the pdf of its location after n time steps

Fig. 2. Model problem with diagonal convection.

␣⫽

v ⌬t. ⌬y

P共ıˆ, ˆ;n兲 j ⫽

n n ⫺ ˆı 冉 再 ˆı 冊冉 ˆj 冊共1 ⫺ ␣ ⫺ ␤兲

n⫺ıˆ⫺jˆ

ı ˆⱕ j n 䡠 ␣ˆı 䡠 ␤ˆj ˆ⫹

0

ˆ⫹ ı ˆ⬎ j n (9)

can be derived for constant ⌬t. Analytical solutions show that the expectations 具i典 and 具j典 assume the correct values of n␣ and n␤, respectively. More importantly, the matrix of variances is 具i⬘ 典 具i⬘j⬘典 冉 具i⬘j⬘典 具j⬘ 典 冊 具i典 䡠 共1 ⫺ ␣ 兲 ⫺n ␣␤ ⫽冉 . ⫺n ␣␤ 具j典 䡠 共1 ⫺ ␤ 兲 冊 2

V⫽

2

(10)

To gain deeper insight into the nature of this matrix, it is re-expressed in terms of ␰ ⫽ u/⌬x, ␩ ⫽ v/⌬y and the CFL number and transformed by T⫽

␰ 冉 ␩ 冑␰ ⫹ ␩ 1

2

2

⫺␩ ␰

冊

(11)

so that the first axis is parallel to uជ yielding

(7)

冉

␰ 4 ⫹ ␰ 3␩ ⫹ ␩ 3␰ ⫹ ␩ 4 ⫺ CFL共 ␰ 2 ⫹ ␩ 2兲 2 ␰␩ 共 ␩ ⫺ ␰ 兲 ␰⫹␩ ␰␩ 共 ␩ ⫺ ␰ 兲 ␰␩ 共 ␩ ⫹ ␰ 兲

The spatial discretization error manifests itself in the variance normal to convection and is independent of the CFL number. As in the one-dimensional case [14], a particle’s position

冊

.

(12)

along the direction of convection is subject to statistical error resulting from randomly selecting particles for exchange. It is possible to reduce this variance by increasing the time step

¨ BUS ET AL. H. MO

524 towards its limit CFL ⫽ 1. But unlike in the one-dimensional case, it can only become zero if the ratio of grid increments equals the ratio of velocity components (␰ ⫽ ␩). Obviously, this problem can be cured by including diagonal cells into the convection operator. The approach, borrowed from 2- or ⌬ND ⫽

冋

3-D unsplit methods [18, 19], can be explained as tracing back the mass being convected into cell i, j to its origin before the time step. By graphical consideration of Fig. 1, the number of particles originating, for example, in cell i ⫺ 1, j ⫺ 1 is found to be

册

冋

册

共sជDAi⫺1, j ⫻ uជ CD 兲 䡠 eជ3 共sជCDi, j⫺1 ⫻ uជ DA 兲 䡠 eជ3 1 ⌬mDA 1 ⌬mCD 䡠 max ⌬t, 0 ⫹ 䡠 max ⫺ ⌬t, 0 2 mi, j 共sជCDi, j⫺1 ⫻ ជsDAi, j 兲 䡠 eជ3 2 mi, j 共sជDAi⫺1, j ⫻ ជsCDi, j 兲 䡠 eជ3

(13) with a small error that vanishes for uniform velocity and density. Vectors along cell faces are denoted by ជs, e.g., ជsCDi, j⫺1 ⫽ (⌬xCD, ⌬yCD)T of cell i, j ⫺ 1, while eជ3 is the unit vector pointing out of the planar computational domain. Convection from diagonal cells is evidently a phenomenon proportional to ⌬t2, but, as was demonstrated above, the time step must be made as large as possible, thus giving rise to the need of considering second order effects. At this point, the analysis for the simple convection problem is repeated using the new convection operator. The probabilities of a par-

ticle being transported from cell i, j to cell i ⫹ 1, j or to i, j ⫹ 1 are now

␣⫽

u v ⌬t ⫺ ␥ and ␤ ⫽ ⌬t ⫺ ␥ , ⌬x ⌬y

(14)

respectively. Moreover, particles are convected to the diagonal cell i ⫹ 1, j ⫹ 1 with probability

␥⫽

u v ⌬t 2. ⌬x ⌬y

(15)

The time step restriction changes to ␣ ⫹ ␤ ⫹ ␥ ⫽ CFL ⱕ 1. After considerable manipulation, the pdf at time n⌬t is found to be

ˆ, ˆ兲 j ⫹k 冘 冉max共ıˆ,nˆ兲j ⫹ k冊冉max共ı min共ıˆ, ˆ兲 j 冊 l

P共ıˆ, ˆ; j n兲 ⫽

k⫽0

䡠

冉min共ıkˆ, ˆ兲j 冊 䡠 ␣

max共ıˆ⫺jˆ, 0兲⫹k

j ˆ, 0兲⫹k j j 䡠 ␤max共ˆ⫺ı 䡠 ␥min共ıˆ, ˆ兲⫺k 䡠 共1 ⫺ ␣ ⫺ ␤ ⫺ ␥兲n⫺max共ıˆ, ˆ兲⫺k

(16)

for ˆı ⱕ n ⵩ ˆj ⱕ n and zero elsewhere. The abbreviation l ⫽ min(ıˆ, ˆ, j n ⫺ ˆ, ı n ⫺ ˆ) j was used therein. For lack of an analytical solution for the moments of this pdf, Eq. 16 was integrated numerically to obtain expectations and variances. It was found that the expectations assume the same values as before (as must be) and that variances decrease as the time step increases, see Fig. 3. However, due to their sharp gradients, variances become small only for CFL numbers very close to unity. Moreover, it is only possible to reduce both variances to zero if the ratio of Cartesian velocities equals the ratio of cell dimensions. This result fully agrees with intuition and proves that numerical diffusion can

Fig. 3. Variance of particle’s position with new convection operator at t ⫽ 10 s. —— 具i⬘2典 ⫽ 具j⬘2典 for u/⌬x ⫽ v/⌬y ⫽ 1 s⫺1, – – – 具i⬘2典, – 䡠 – 䡠 – 具j⬘2典, the latter two for u/⌬x ⫽ 1 s⫺1, v/⌬y ⫽ 0.5 s⫺1.

EULERIAN AND LAGRANGIAN PDF METHODS

525

Fig. 4. Search and locate algorithm.

be completely suppressed by inclusion of diagonal cells into the convection operator. LAGRANGIAN METHOD In the Lagrangian method, the particles’ spatial locations evolve independently of each other by [8]

冉

x共 p兲共t ⫹ ⌬t兲 ⫽ x共 p兲共t兲 ⫹ ⌬t 䡠 u ⫹ ⫹ 共2⌬tD eff兲 1/ 2 䡠 ␰,

1 ⵜ共 ␳៮ D eff兲 ␳៮

冊

(17)

in which ␰ is a standardized joint normal random vector. Particularly in axisymmetric flows, where a particle’s weight w ( p) is initially proportional to its distance from the axis (if cells are of equal size and if density is constant), it is indispensable to adjust weights as the simulation proceeds. The method described by Howarth and El Tahry [20] was found to be very successful to that end. Search and Locate Algorithm In the course of the simulation it is necessary to continually identify the cells in which the particles are located. Here, a modified search and locate algorithm for arbitrary quadrilateral cells derived from the one presented by Allievi and Bermejo [21] is introduced. Instead of transforming the complete cell into a unit square, which results in an iterative process, the cell is split into two triangles as illustrated in Fig. 4. Subsequently, a particle’s location is mapped to unit triangles I and II as also described in Ref. [21], which can

be achieved directly without iteration. From the transformed coordinates, p I , q I , p II , and q II , the direction in which the search is to be continued can readily be deduced. For example, if ⫺1 ⱕ q I ⬍ 0 and all other values are positive, the search is continued in cell i, j ⫺ 1. It is terminated when all transformed coordinates are positive indicating that the cell hosting the particle has been identified. After locating all particles the whole array is re-sorted according to increasing cell indices. This saves a lot of indirect addressing in the program and allows recycling parts of the Eulerian code. Interpolation In the Lagrangian method it is necessary to calculate several flow quantities at the particles’ locations. The flow field is obtained from the finite volume solver with which the pdf module is coupled and is therefore given in discrete form at cell centres. From these, values at cell corners and at mid-points of cell faces are calculated by averaging four or two adjacent cell centre values, respectively. At solid walls or at symmetry lines physical boundary conditions (such as the no slip condition for velocity components) apply. Values at locations within the cell are obtained by biquadratic interpolation with basis functions given by Allievi and Bermejo [21]. In contrast to bilinear interpolation, it is then possible to prescribe the 9 values discussed above (at four corners, four mid-points of faces plus cell-centre and thereby to accurately fulfill physical boundary conditions imposed at cell

¨ BUS ET AL. H. MO

526 faces. From the computational point of view, biquadratic interpolation is not significantly more expensive than its bilinear counterpart. However, the difference observed between interpolation and simply assuming the flow variables to be constant within the computational cell (zeroth order interpolation) was very small and hardly visible for the grids used in our computations.

Fig. 5. Numerical diffusion caused by turbulent mixing.

Local Time Stepping In simulations of complicated geometries, fine grid resolution is often essential for resolving all features of the flow. If the time step is determined globally, i.e., the smallest cell limits temporal advancement, convergence of the simulation is slowed down drastically for the complete field, even in regions where characteristic residence time scales are long enough to allow larger time steps. This may finally lead to the point where simulations become impractical due to excessive computational requirements. When considering statistically stationary problems, local time stepping, which is a valuable remedy for this problem, has been applied successfully in Eulerian calculations before [14, 15, 22]. In this case, the time step for each cell is determined individually so that the CFL number is constant throughout the domain. From the Lagrangian point of view, implementation of local time stepping is not as straight forward as in the Eulerian method. In the latter case, after reaching stationarity there is an equilibrium between the number of particles being transported into the cell and the remaining processes that have no spatial dependence. In Lagrangian methods, however, convection and diffusion are realized principally without reference to the computational grid. Hence, it is impossible to control the number of particles brought into a certain cell by its own time step. Nevertheless, we found it sufficient to adjust the particles’ weights according to the ratio of time steps of their original and their destination cells w 共 p兲共t ⫹ ⌬t兲 ⫽ w 共 p兲共t兲 䡠

⌬t共x共 p兲共t ⫹ ⌬t兲兲 . ⌬t共x共 p兲共t兲兲

(18)

This of course does not eliminate the error of the particles’ locations due to the change of time steps between cells, which makes them travel too far or too short a distance into destination cells (depending on the time step ratio). This error, however, decreases with the CFL number. The computational effort, on the other hand, was reduced dramatically, and many calculations were rendered possible only by this measure. Mixing Models The adaption of the Modified Curl’s Model to unequal particle weights has been described by Nooren et al. [23]. In principle, the IEM mixing model [24], which is defined by 1 dY 共 p兲 ⫽⫺ 共Y 共 p兲 ⫺ 具Y典兲, dt 2␶␾

(19)

with ␶␾ denoting a time scale, can be applied to the Lagrangian scheme without alterations. However, in contrast to the Eulerian method, turbulent mixing constitutes a source of numerical diffusion in Lagrangian schemes. All mixing models must and do have the property of reducing scalar variance or, equivalently, moving particles towards the ensemble average ⬍Y⬎ in composition space. If there is a gradient in composition within a cell of finite size ⌬x, turbulent mixing will always tend to lower this gradient as illustrated in Fig. 5. Of course, this error vanishes with the grid size tending towards zero, similar to discretization errors. To circumvent grid refinement, which increases the computational cost, it proved to be very efficient in the course of the

EULERIAN AND LAGRANGIAN PDF METHODS present study to perform a bilinear least squares fit of all particles contained in a certain cell. The least squares fit must be weighed with the individual particle weights. Subsequently, the ensemble average in Eq. 19 is replaced by the local mean at the parti- cle’s location as estimated from the least squares fit. The least squares fit is nonsingular if more than two particles are present in a given cell and if not all of them are positioned on a straight line, both of which can be identified by the determinant of the least squares fit matrix being zero. In order to prevent bounded variables from exceeding their allowed range, the local mean is limited by the maximum and minimum mean values of the eight surrounding cells. This procedure does not introduce numerical diffusion as these local extrema are artifacts of the linear fitting. Finally, the changes necessary to apply the Continuous Mixing Model of Hsu and Chen [25] to particles of unequal weights are presented. As in the Eulerian case, all particles residing in a certain cell are paired randomly. In the absence of chemical reaction, the properties of a pair ( p), (q) evolve in time by

␾ 共 p兲共t ⫹ ⌬t兲 ⫺ ␾ 共 p兲共t兲 ⫽

527

Nw 共q兲 d ␾ 共 p兲 䡠 C 䡠 共 ␾ 共q兲 ⫺ ␾ 共 p兲兲, ⫽ N dt w 共i兲

(20)

Nw 共 p兲 d ␾ 共q兲 ⫽ 䡠 C 䡠 共 ␾ 共 p兲 ⫺ ␾ 共q兲兲, dt N 共i兲 w

(21)

冘

冘

i

i

with coefficient C ⫽ C⬘

␰ . ␶␾

(22)

Constant C⬘ must take a value of one to yield the correct variance decay. ␰ denotes a uniformly distributed random variable ranging from 0 to 1, and ␶␾ is a characteristic time scale of scalar dissipation that is usually linked to the time scale of turbulence ␶ through [8]

␶␾ ⫽

␶ C␾

(23)

with an assumed constant C ␾ ⫽ 2. Without chemical reaction, the analytical solution to Eqs. 20 and 21 yields

冋 冉

w 共q兲 N 䡠 共w 共 p兲 ⫹ w 共q兲兲 共q兲 共 p兲 共 ␾ 共t兲 ⫺ ␾ 共t兲兲 䡠 1 ⫺ exp ⫺ 䡠 C 䡠 ⌬t N w 共 p兲 ⫹ w 共q兲 w 共i兲

RESULTS Convection Test Case To demonstrate the effects of numerical diffusion associated with the Eulerian method, a problem of pure mean convection is considered first. At the inlet on the left boundary of the domain there are two separate streams having mixture fraction values of zero in the centre and one in the surrounding flow. These streams with horizontal velocity should not mix in the absence of diffusion. This problem is solved on the highly distorted grid shown in Fig. 6a. The distribution of mixture fraction along a cut through the exit plane depicted in Fig. 6b reveals a striking amount of numerical diffusion in the conventional Eulerian solution. The plot further shows that the inclusion of diagonal cells into the convection operator is not able to cure

冘

i

冊册

. (24)

the problem completely because convection is not exactly diagonal to the individual cells. Nevertheless, this method is able to visibly reduce numerical diffusion. For the sake of completeness, Fig. 6b also proves that the Lagrangian approach is free of numerical diffusion. It should be repeated that there indeed is a different source of numerical diffusion in Lagrangian methods stemming from processes that need to make use of the computational grid, namely turbulent mixing. Turbulent mixing has the effect of reducing variances and consequently smoothes gradients as was discussed above. Convection Diffusion Test Cases By comparison with finite difference reference solutions, the second test case should both serve as proof of spatial accuracy and as a demonstration of the efficiency of local time stepping. The

¨ BUS ET AL. H. MO

528

Fig. 8. Mixture fraction distribution for convection diffusion test case along bottom line of symmetry: a. Whole domain, b. Close-up view. Fig. 6. Mixture fraction distribution for 2-D convection test case: a. Computational Grid, b. Cut through exit plane.

relevance of spatial accuracy when dealing with reacting flows is then explored by a similar setup simulated with first and second order finite difference codes. We resorted to these conventional methods for this purpose to save computational resources. As sketched in Fig. 7, the inert planar test

Fig. 7. Convection diffusion model problem.

case consists of two co-flowing streams of constant velocity with equal and constant pressure and density initially having mixture fractions of zero and one. Turbulent mixing is not considered here. Axial velocity of both flows is 200 m/s and the uniform effective diffusion coefficient is taken to be D eff ⫽ 0.2 m2/s. The problem has been solved by both pdf methods as well as by use of upwind finite differences of first and second order accuracy. The grid is made up of 20 ⫻ 20 cells covering 4 m ⫻ 0.2 m. Figure 8a shows the mixture fraction distribution along the bottom line of symmetry, y ⫽ 0. As was to be expected, the Eulerian pdf approach yields exactly first order accuracy. In contrast, the Lagrangian scheme is even better than the second order method on the same grid; it can hardly be distinguished from the grid-independent finite difference solution (obtained on a 80 ⫻ 20 grid) which is

EULERIAN AND LAGRANGIAN PDF METHODS given as a close-up view in Fig. 8b. Also included in this figure is the solution obtained with the algorithm proposed by Anand et al. [15] that was designed to enhance spatial accuracy. The conclusion from the discussion above is underlined by this test case: in many practical problems the treatment proposed by Anand et al. will achieve little improvement since exchange by diffusion in the main flow direction is small compared to convection. To prove the concept of local time stepping in the context of Lagrangian schemes, a grid of 80 ⫻ 80 cells was constructed of which the corner cells posess the same dimensions as on the previous 20 ⫻ 20 grid. The grid is refined in both directions towards the centre of the domain resulting in a ratio of largest to smallest cell dimension of about 57 in each coordinate direction. Consequently, the time step in the computation varies by a factor of approximately 2,600. This might seem large but is a value quite typical of other simulations. Figures 9a and b display the centreline evolution of mixture fraction and a cut through the exit plane for CFL numbers of 0.5 and 0.9. Evidently, even for the larger CFL number the results comply very well with the finite difference solution although small discrepancies can be perceived. What is remarkable is that these solutions converged to statistical stationarity within 7,500 iteration cycles at the lower CFL number. We intended to present a comparison with global time stepping. This calculation, however, was so far from convergence after 20,000 time steps that continuation till stationarity was deemed too expensive. One further observation is to be made from these results. The statistical error inherent in Monte Carlo schemes is reduced by time-averaging the solution over n avg iteration steps (2,500 in the present case). In an extensive discussion Xu and Pope [26] explain that such time series contain a characteristic scale, ␶avg, over which events are correlated. Moreover, the stochastic error reduces by a factor of (2 ␶avg/ (n avg⌬t)) 1/ 2 . Because of the flow field’s homogeneity, ␶avg is uniform across the domain. Consequently, as ⌬t varies whereas n avg is constant, the statistical error is larger where the time step is smaller. This explains the deviation of the results obtained with local time stepping from

529

Fig. 9. Mixture fraction distribution for convection diffusion test case. Lagrangian scheme with local time stepping on 80 ⫻ 80 stretched grid: a. Along bottom line of symmetry, b. Cut through exit plane.

the reference solution near the centre of Fig. 9b (at y ⬇ 0.1 m). To determine the influence of spatial accuracy on ignition, a similar reactive laminar test problem was devised and studied by applying finite difference schemes of first and second order accuracy. As was seen above, these solutions correspond closely to the pdf results but can be obtained at a much smaller computational expense. The axial velocity is now set to 500 m/s, the diffusion coefficient to 0.0005 m2/s. Inlet temperatures of both streams amount to 1050 K. The lower stream is made up of equal mass fractions of oxygen and nitrogen, whereas the upper flow consists of hydrogen (5% of mass) diluted in nitrogen. Under constant pressure and adiabatic conditions, a considerable igni-

530

¨ BUS ET AL. H. MO with the number of discrete volumes. Moreover, the maximum time step reduces linearly with the cell dimension in main flow direction. Thus, the computational cost to reach a converged solution eventually increases quadratically with decreasing cell size. Although at a first glance the Lagrangian approach is much more expensive, taking these considerations for the Eulerian method into account it may even end up being the cheaper alternative. Supersonic Diffusion Flame Experiment

Fig. 10. Distribution of temperature and of mass fraction of OH along centreline for convection diffusion test case with reaction. Note the offset in x-direction.

tion delay is perceived in the computational domain of 0.08 m ⫻ 0.004 m. Figure 10 shows temperature and mass fraction of the OH radical along the centreline ( y ⫽ 0.002 m), which are good indicators of ignition. Obviously, the prediction of first order accuracy obtained on the initial 80 ⫻ 40 volume grid is far from satisfactory, and even a grid five times finer in main flow direction results in just as good a solution as the second order calculation on the coarsest grid. To get close to gridindependence, 160 ⫻ 40 cells are necessary in case of the second order scheme, whereas ten times as many cells in x-direction are needed for the first order method. This clearly demonstrates the need for high spatial accuracy to be able to reliably predict ignition. If the stochastic error of the Monte Carlo simulation is not to increase, the number of particles per cell has to be constant. In this case the total number of particles increases linearly

Finally, the supersonic hydrogen-air diffusion flame experiment of Cheng et al. [27] is investigated, which has been used for examination of the Eulerian code before [14]. In this axisymmetric burner experiment a stream of pure hydrogen at 545 K is injected at sonic speed into a vitiated air flow (Y H2O ⫽ 0.1749, Y N2 ⫽ 0.5802 and Y O2 ⫽ 0.2449) with a Mach number of 2 and 1250 K static temperature. A lifted flame develops under ambient conditions. Experimental data of means and rms fluctuations of major species (N2, O2, H2, H2O) and of the OH radical are available in seven downstream planes. Here, we concentrate on grid-dependence of the simulations. For further details, e.g., implementation of boundary conditions as well as modeling of compressibility effects, and for an extended comparison with the experiment, see Ref. [14]. The results shown below were obtained with the Continuous Mixing Model discussed above. Originally, an orthogonal grid of 68 ⫻ 124 cells was used. The Eulerian calculation on this grid was still found to be grid-dependent [14]. In the present paper, the Eulerian simulation is repeated on a grid with twice the number of volumes in the streamwise direction. Lagrangian calculations were performed on both the original and the refined grid. Figure 11 depicts the maximum molar fraction of the OH radical in planes of constant distance x downstream of the injector. As in the simple test case above, the main influence of spatial accuracy is on the prediction of ignition. The difference between the solutions obtained on the two grids is strongest for the finite volume calculation using laminar chemistry. With the Lagrangian scheme, only a slight shift of the ignition point is perceived so that this simulation may be close to

EULERIAN AND LAGRANGIAN PDF METHODS

Fig. 11. Maximum radial molar fraction of OH as function of axial location for supersonic hydrogen-air diffusion flame: a. Whole domain, b. Ignition point.

grid-independence. Even on the finer grid the Eulerian method deviates from the Lagrangian solution. Consequently, the Eulerian simulation requires further grid refinement for convergence. The relatively large difference in finite volume solutions on the two grids can be attributed to the effects of artificial dissipation used for shock-capturing. This measure reduces the scheme locally to first order accuracy, but artificial dissipation diminishes with enhanced grid resolution too. Unfortunately, results for both pdf methods exhibit much lower levels of OH than the experiment, whereas the laminar chemistry method seemingly does a better job. On the other hand, it must be mentioned that the experimental data were deteriorated by problems achieving true axisymmetry. Cheng et al. [27] cite difficulties with centering the hydrogen jet tube in the apparatus and further influences

531

of unequal thermal expansion as sources of obvious unsymmetry of their measurements. In supersonic flows, however, the consequences are hard to estimate as shock patterns may change significantly and in turn influence the temperature field and flame strongly. Therefore, the comparison of calculations and experimental results should be viewed with caution, and differences should not be overestimated. As this is the only known experiment with measurements of radical species and of fluctuations, it is nevertheless a valuable basis for discussion of topics relevant to the present investigations. The two pdf schemes are further investigated by plotting profiles at the downstream position x/D ⫽ 21.5, which is located shortly after ignition. Here, the flame exhibits the ring-like structure typical of diffusion flames. Concluding from Fig. 12, neither the profiles for mean temperature nor those for inert nitrogen, reaction product H2O or radical OH change dramatically from one pdf scheme to the other. Differences are most pronounced in the area of chemical activity (indicated by the presence of OH) and become negligible in the outer region of the flow. Compared to the difference between pdf calculations and finite volume results with laminar chemistry, the deviation of the two pdf schemes is indeed very small. This view is underlined by results for rms fluctuations plotted in Fig. 13. Again, the differences are pretty small. To weigh costs against benefits, the time consumed in each of the methods was measured. One set of calculations was performed on a NEC SX-4 vector supercomputer. 200 and 100 particles per cell were used in case of the Eulerian and Lagrangian methods, respectively. As to reduce stochastic errors and because the coupled finite volume solver proved to be very sensitive to even minor wiggles in the distribution, ensemble averages were further averaged over 2,500 iteration steps [26] and then spatially smoothed [14]. Because it is the same in both approaches, the cpu time for calculation of chemical reaction progress is referred to as 1 work unit. As given in Table 2, on this machine spatial exchange in the Lagrangian method takes nearly ten times as long as in the Eulerian version. This is attributed to a high degree of indirect addressing and relatively short vector

532

Fig. 12. Mean temperature and molar fraction profiles of N2, H2O, and OH at x/D ⫽ 21.5. —— Eulerian pdf calculation, – – – Lagrangian pdf calculation, – 䡠 – 䡠 – finite volume calculation, { experiment [27].

¨ BUS ET AL. H. MO

Fig. 13. RMS profiles of temperature and of molar fractions of N2, H2O, and OH at x/D ⫽ 21.5. —— Eulerian pdf calculation, – – – Lagrangian pdf calculation, { experiment [27].

EULERIAN AND LAGRANGIAN PDF METHODS TABLE 2 Computer Time on NEC SX-4 in Work Units Procedure Chemical reaction Spatial exchange

Eulerian

Lagrangian

1 0.121

1 1.159

533

major ingredients) on the scalar than on the vector machine. This is due to the outstanding performance of the entirely vectorizable reaction module on the vector computer. CONCLUSIONS

lengths in the Lagrangian scheme. Moreover, the procedure used to re-order particles (about 19% of the effort for spatial exchange) is not completely vectorizable, reducing performance considerably. Even though fully vectorized, 17% of the time for exchange is currently consumed by generating Gaussian random numbers. A future version of the program will feature an algorithm similar to the one presented by Hsu et al. [28] by providing one long sequence of random numbers from which subsets are randomly selected for individual cells. Although this will improve performance, spatial exchange in the Lagrangian approach will remain about as expensive as direct integration of chemical reaction progress of a seven species, seven elementary reaction scheme. As in another study on numerical efficiency [29], a slightly different picture emerges on a scalar SUN Sparc 60 workstation. On this platform, merely 50 particles per cell could be accommodated. It must be pointed out that the entire code has been optimized to run efficiently on vector computers. By relating times to the effort spent on calculating chemical reaction, we believe to obtain a realistic picture anyway. Table 3 demonstrates that the Lagrangian scheme takes about seven times as long for spatial exchange as the Eulerian approach on this scalar machine. Although the ratio is somewhat better than on the NEC SX-4, the difference is nonetheless considerable. Also noteworthy is that the reaction module constitutes a much larger portion of the total computational effort (of which reaction and exchange are TABLE 3 Computer Time on SUN Sparc 60 in Work Units Procedure Chemical reaction Spatial exchange

Eulerian

Lagrangian

1 0.032

1 0.228

Eulerian and Lagrangian approaches for the solution of the pdf transport equation are compared concerning numerical diffusion, spatial accuracy, and numerical efficiency. An analysis of numerical diffusion in Eulerian schemes for two-dimensional flows is presented which gives impetus for the development of an improved treatment of convection. A simple test case shows that this new convection operator is able to visibly reduce numerical diffusion. It is further pointed out that turbulent mixing constitutes a source of numerical diffusion in case of the Lagrangian approach and a cure for this problem is presented for the IEM model. With respect to numerical efficiency of the Lagrangian scheme, a simple though very efficient algorithm to realize local time stepping is devised. In many simulations of technical relevance this will reduce the immense computational costs considerably. A simple test case demonstrates that the discretization of the Eulerian approach is of first order spatial accuracy (as expected). The Lagrangian scheme, in contrast, is shown to coincide with the grid-independent solution of the mixture fraction transport equation for this problem. Moreover, it was deduced from a further artificial case including chemical reactions that spatial accuracy is extremely relevant to obtain a grid-independent prediction of ignition with a reasonable number of cells. Although the Lagrangian scheme increases the numerical expense, in this case it was deemed the cheaper choice for a grid-converged solution. Simulation of a supersonic hydrogen-air diffusion flame demonstrates the applicability of the Lagrangian scheme to practical problems and provides further comparison with both the Eulerian approach and the laminar chemistry method. The main influence of improved spatial accuracy is on the prediction of ignition. Elsewhere, the results of Eulerian and Lagrangian approaches are similar. However, the Lagrangian method

¨ BUS ET AL. H. MO

534 consumes between seven- and ten-fold the computer time for spatial exchange of particles compared to the Eulerian scheme. Therefore, the Lagrangian method is mainly seen as a stepping stone towards implementation of a full joint-pdf approach including velocity and turbulent frequency, which renders the gradient diffusion assumption superfluous at little extra cost. We would like to thank the Deutsche Forschungsgemeinschaft for funding of this project under contract No. BR 1713/1-2 and within the SFB 259 at the University of Stuttgart as well as for the provision of computational resources at the HLRS Stuttgart. We are also thankful to Prof. T. S. Cheng for making available the experimental data in electronic form. We further appreciate discussions with M. Rieber on numerical diffusion resulting from convection diagonal to the computational grid. REFERENCES 1.

2. 3.

4.

5. 6. 7.

8. 9.

Norris, A. T., and Hsu, A. T., 30th AIAA/ASME/SAE/ ASEE Joint Propulsion Conference, AIAA 94-3356, Indianapolis, IN, 1994. Baurle, R. A., Alexopoulos, G. A., and Hassan, H. A., J Propulsion Power 10:473– 484 (1994). Gerlinger, P., Mo ¨bus, H., and Bru ¨ggemann, D., 30th AIAA Fluid Dynamics Conference, AIAA 99-3775, Norfolk, VA, 1999. Baurle, R. A., and Girimaji, S. S., 37th AIAA Aerospace Sciences Meeting and Exhibit, AIAA 99-0928, Reno, NV, 1999. Jones, W. P., and Kakhi, M., Combust. Flame 115:210 – 229 (1998). Tsai, K., and Fox., R. O., AIChE J. 42:2926 –2940 (1996). Proceedings of the Third International Workshop on Measurement and Computation of Turbulent Nonpremixed Flames, Boulder, CO (1998). (Available at http://www.ca.sandia.gov/tdf/3rdWorkshop/Boulder.html) Pope, S. B., Prog. Energy Combust. Sci. 11:119 –192 (1985). Hsu, A. T., AIAA 22nd Fluid Dynamics, Plasma Dynamics and Lasers Conference, AIAA 91-1780, Honolulu, Hawaii, 1991.

10. 11. 12.

13.

14.

15.

16. 17. 18. 19. 20. 21. 22.

23.

24. 25.

26. 27. 28.

29.

Gerlinger, P., and Bru ¨ggemann, D., Int. J. Numerical Methods Fluids 24:1019 –1035 (1997). Gerlinger, P., Stoll, P., and Bru ¨ggemann, D., J. Computational Physics 146:322–345 (1998). Coakley, T. J., and Huang, P. G., 30th AIAA Aerospace Sciences Meeting and Exhibit, AIAA 92-0436, Reno, NV, 1992. Baurle, R. A., Modeling of Turbulent Reacting Flows with Probability Density Functions for Scramjet Applications, Ph.D. thesis, North Carolina State University, Raleigh, NC, 1995. Mo ¨bus, H., Gerlinger, P., and Bru ¨ggemann, D., 37th AIAA Aerospace Sciences Meeting and Exhibit, AIAA 99-0198, Reno, NV, 1999. Anand, M. S., James, S., and Razdan, M. K., 34th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, AIAA 98-3856, Cleveland, OH, 1998. Pope, S. B., Combust. Sci. Technol. 25:159 –174 (1981). Roekaerts, D., Comput. Fluids 21:97–108 (1992). Colella, P., J. Computational Physics 87:171–200 (1990). Saltzman, J., J. Computational Physics 115:153–168 (1994). Howarth, D. C., and El Tahry, S. H., AIAA J. 29:208 – 218 (1991). Allievi, A., and Bermejo, R., J. Computational Physics 132:157–166 (1997). Biagioli, F., in Advanced Computation and Analysis of Combustion (G. D. Roy, S. M. Frolov, and P. Givi, Eds.), ENAS Publishers, Moscow, 1997, pp. 450 – 465. Nooren, P. A., Wouters, H. A., Peeters, T. W. J., Roekaerts, D., Maas, U., and Schmidt, D., Combust. Theory Modelling 1:79 –96 (1997). Dopazo, C., Physics of Fluids 18:397– 404 (1975). Hsu, A. T., and Chen, J.-Y., Eighth Symposium on Turbulent Shear Flows, pp. 22-4-1–22-4-5, Munich, Germany, 1991. Xu, J., and Pope, S. B., J. Computational Physics 152:192–230 (1999). Cheng, T. S., Wehrmeyer, J. A., Pitz, R. W., Jarrett, O. Jr., and Northam, G. B., Combust. Flame 99:157–173 (1994). Hsu, A. T., Tsai, Y.-L. P., and Raju, M. S., 31st Aerospace Sciences Meeting and Exhibit, AIAA 930087, Reno, NV, 1993. Mo ¨bus, H., Gerlinger, P., and Bru ¨ggemann, D., in Computational Fluid Dynamics ’98 —Proceedings of the 4th ECCOMAS Conference (K. D. Papailiou, D. Tsahalis, J. Periaux, C. Hirsch, M. Pandolfi, Eds.) John Wiley & Sons, 1998, 1:162–168.

Received 18 October 1999; revised 20 August 2000