Estimation of Numerical Errors Related to Some ... - ACS Publications

10588

Ind. Eng. Chem. Res. 2010, 49, 10588–10605

Estimation of Numerical Errors Related to Some Basic Assumptions in Discrete Particle Methods Sofiane Benyahia* and Janine E Galvin National Energy Technology Laboratory, Morgantown, West Virginia 26507

Discrete particle methods that track the motion of individual particles and their collisions are computationally very expensive. To accelerate these numerical simulations, some basic assumptions have been introduced and reported in the literature. This study investigates two of these common assumptions: (1) the use of computational parcels, or clouds, wherein many particles are lumped together so that only parcels and their collisions are tracked, and (2) the multiphase particle-in-cell, or MP-PIC, wherein the collision forces are replaced by a solids pressure term with the main purpose to avoid exceeding the maximum packing of the granular assembly. Using several cases relevant to the fluidization community, errors associated with these assumptions are computed. For these cases the magnitude of error in the time-averaged flow variables indicates that further research on the validity of these assumptions is warranted. Introduction Particulate flows occur commonly in nature, such as in dust storms and volcanic eruptions, as well as in many industries, such as chemical and energy conversion technologies. Traditionally, two approaches have been used to model particulate flow systems. The first is based on granular kinetic theory,1,2 which treats the granular media as a continuum (also called the Eulerian or two-fluid model). The second is based on discrete particle methods3,4 (Lagrangian), which describe the motion of each individual particle following Newton’s second law and resolves the collisions between particles and collisions between particles and wall boundaries using either a soft-sphere or a hard-sphere model. The discrete particle method is more accurate than the continuum description of particulate flows since it involves fewer assumptions in the model derivation. However, numerical simulations based on this method are extremely expensive for any practical system due to the large numbers of particles. To avoid tracking large numbers of particles in a discrete particle simulation, researchers have employed the concept of computational parcels (also referred to as clouds or numerical particles by others in the literature) which contain many particles moving at the same velocity.5-7 This concept is commonly used in discrete simulations of droplets to model sprays8,9 and will be referred to herein as the coarse particle model. The development of this technique in spray modeling was originally based on the assumption of noninteracting droplets. The possibility that droplets may collide was later incorporated using stochastic principals to determine the probability of collisions and their outcome9,10 (character of the resulting collision). It is worth noting that the probability of two parcels colliding takes into consideration the number of particles per parcel. Using computational parcels seems fairly reasonable for dilute flow conditions where collisions are scarce. The parcel approximation may also be reasonable in sprays for dense, or collisional, flow conditions if accompanied by modifications for predicting the incidence of collision. The accuracy/validity of the simulation for predicting the spray behavior is still subject to the precise model for describing the outcome of collision. More recently, the parcel approximation has been extended to model fluidized beds7,11-14 and other dense particle * To whom correspondence should be addressed. E-mail: [email protected]. 10.1021/ie100662z

flow,5,6,15-18 where both collisions and fluid-particle interactions play a major role in the transport mechanisms of particles. These simulations are most often found combined with a particle-in-cell technique, discussed shortly, where collisions between parcels are not treated explicitly, but are instead subsumed into a continuum pressure term. In either case (dense or dilute), the concept of computational parcels for gas-solids flows needs to be rigorously verified for as many practical applications as possible. The literature studies for particulate flows have focused on qualitative validation of this approach using experimental data and/or observations.7,11,16-18 Accordingly, the primary goal of the current study is to quantitatively verify the concept of computational parcels by comparing the numerical results obtained using different parcel sizes. The results using the smallest numbers of particles per parcels are assumed to be the most accurate (i.e., 1 particle per parcel) and are used to quantify the errors in simulations with larger numbers of particles per parcel. Two systems are examined. The first system is wall-bounded shear flow of particles in a vacuum. Various parcel sizes are used ranging from a single particle to 10 particles per parcel. The numerical error associated with the computational parcel assumption is assessed by comparing the time-averaged results (e.g., particle velocity, concentration, granular temperature). The parcel method is expected to perform poorly in this shear system since collisions are a major means of momentum and energy transport. Accordingly, the shear system may be considered a worst-case scenario for the parcel method. The second system, more relevant to fluidization, involves the flow of air and particles in a periodic vertical riser. As with the shear system investigated, simulations are conducted using several parcel sizes and the resulting time-averaged flow variables are compared. The second goal of this study is to investigate the impact of using a continuum solids pressure gradient in the particle equation of motion in place of the traditional collision forces associated with soft-sphere and hard-sphere discrete models (i.e., collisions are ignored). The continuum solids pressure is primarily employed to avoid exceeding the maximum packing of the granular assembly. This approach is commonly referred to as the multiphase particle-in-cell (MP-PIC) technique,5,7,11,14 and it also uses computational parcels as described earlier.

This article not subject to U.S. Copyright. Published 2010 by the American Chemical Society Published on Web 05/03/2010

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

The particle-in-cell (PIC) technique was originally designed for modeling fluid flows involving large distortions, such as, at colliding interfaces or shocks.19 The approach combines features from both the Lagrangian and Eulerian approaches resulting in a dual representation of the fluid. Essentially the PIC technique involves partitioning the fluid into discrete (Lagrangian) particles that are moved (advected) through an underlying Eulerian grid. Properties of the fluid (e.g., mass, velocity/momentum, energy) are assigned to the particles while interactions among the particles are solved on the computational grid. Generally, particles properties are interpolated to the grid, then the continuum equations for field variables are solved on the grid, and finally, the information is transferred back to the particles to update their velocity and position. Later variants of the PIC technique worked to improve the accuracy of the method by modifying the algorithm, for example, increasing and/or changing the information carried by the particles and/or improving the interpolation methods.20-22 In these later extensions, the grid may no longer carry permanent information but merely act as a computation convenience, or scratch-pad, for organizing particle data to calculate particle interactions.20 Besides fluid dynamic problems, the PIC method has been employed in simulating sea ice dynamics23 and in plasma physics applications24 where the method essentially consists of following the trajectories of charged particles in an electromagnetic field computed on a fixed mesh. Another technique related to the PIC approach is called the material point method (MPM).25 The MPM is primarily targeted at applications in soil dynamics to simulate large material deformations (or very dense/static type granular flows). Andrews and O’Rourke5 adapted the PIC method for simulating coupled gas-solids flow wherein the PIC model is used to describe situations with flowing solids. In their approach, called multiphase PIC (MP-PIC), interactions between solid particles are described by a continuum solids pressure term that depends on the total solids volume fraction. For evaluating the advanced time solids pressure, an implicit approximation to the advanced time solids volume fraction is calculated based on a Taylor series expansion of the expression for the volume fraction at a grid point. In 1997, the MP-PIC technique was extended to two dimensions11 and shortly after to three dimensions.16 The former effort also included improved interpolation methods (mapping particle properties to and from the Eulerian grid) while the latter offered some modifications to the application of the solids pressure term, specifically, a dependency on each solids’ velocity vector that eliminated the need for implicit calculation of the volume fraction. Other groups have also employed this MPPIC method for simulating gas-solids flow.7,14,15 To quantify the error associated with the MP-PIC assumption (i.e., the continuum solids pressure gradient), sets of discrete parcel simulations using either the solids pressure gradient term or the collision force based on a soft-sphere linear spring-dashpot model are conducted. The time-averaged flow profile results obtained using the MP-PIC model are compared to those obtained using the soft-sphere collision model. Differences between these results are considered to reflect the numerical error associated with the MP-PIC assumption. In this portion of the study, simulations of a bubbling fluidized bed and a periodic riser flow, similar to the riser system used earlier, are employed. The numerical errors associated with the computational parcel assumption and the MP-PIC approach are summarized for several basic particle-air flow systems of interest to the fluidization community (e.g., bubbling and transport fluidized

10589

beds). The computational speedup by invoking these assumptions in the simulations and the qualitative similarities in the numerical results, such as the core-annulus regime, are also reported. Description of Model The model equations used in the simulations of this work are outlined. The fluid phase is described via averaged equations of motion, specifically, the single-phase continuum Navier-Stokes equations for a Newtonian fluid, modified to include the fluid volume fraction as well as the fluid-particle interaction due to the drag force. Accordingly, the conservation of mass and momentum for the fluid flow are ∂(Fgεg) + ∇ · (Fgεgb Vg) ) 0 ∂t

(1)

and

[

]

∂(Fgεgb V g) f + ∇ · (Fgεgb V gb Vg) ) -εg∇Pg + ∇ · fτ ∂t NT φp b (x βp (V b)-b Vp) + εgFgb g φc g p p)1

∑

(2)

where εg is the fluid volume fraction, Fg is the fluid-phase f density, Pg is the fluid phase pressure, fτ is the fluid viscous stress tensor, b g is the specific gravity force, b Vg is the local bp) is the fluid velocity average velocity of the fluid phase, b Vg(x at the parcel location, b Vp is the velocity of the pth parcel (discussed below), and βp is the interphase momentum exchange coefficient. NT represents the total number of parcels in the system, and φp and φc are the volume of the pth parcel and the fluid cell, respectively. For a fully specified set of equations, closure models are needed for the fluid viscous stress tensor and the drag force. The stress tensor for a Newtonian fluid, in which the viscous stresses are proportional to the rate-of-strain, is given as f f 1 b 1 b ff ff b t τ ) 2µgfSg, where fSg ) (∇V g + (∇Vg) ) - ∇ · VgI 2 3

(3) ff where µg is the gas viscosity, S is the deviatoric part of the f rate-of-strain tensor, and If is the identity tensor. Throughout this study, the standard Wen-Yu26,27 interphase momentum exchange coefficient for the pth parcel in a fluid cell is used to close the drag force bg(x bp) - b Vp | -2.65 3 Fgεg |V βp ) CD εg 4 dpt

(4)

where dpt is taken as the particle diameter (discussed below) and CD is the drag coefficient expressed as

{

0.687 ) Re < 1000 and Re ) CD ) 24/Re(1 + 0.15Re 0.44 Re g 1000 bg(x bp) - b Vp |dpt Fgεg |V

µg

(5)

The MFIX code (https://mfix.netl.doe.gov) is employed for all simulations in this effort. [The development version of the code, downloaded May 2009 from MFIX CVS https://mfix.netl.doe.

10590


gov/cgi-bin/cvsweb.cgi/mfix/, was used.] In MFIX, the above continuum conservation equations are discretized on a staggered Cartesian computational grid and are solved using a semiimplicit scheme based on the SIMPLE method with a variable time step. In a typical discrete particle model (DPM), the motion of each indiVidual particle is described using Newton’s law of motion. The first goal of the current effort, however, is to assess the use of computational parcels, a technique where a large number of particles are grouped together. Thus, rather than tracking individual particles, they are modeled in groups (or parcels). In this approach, a particle becomes a special case of a parcel made of one particle. The advantage of this approach is reduced computational time, since fewer particles and collisions need to be tracked and resolved. Accordingly, the position, linear velocity, and angular velocity of each indiVidual parcel (pth parcel) are calculated from Newton’s laws dx bp )b Vp dt

(6)

bp dV ∇Pg βp b (x )b g + (V b)-b Vp) + b Fc,p /mp dt Fp Fp g p

(7)

and (2) the errors inherent in the treatment of parcel collisions via the soft-sphere model. Typically, such particle collision models are replaced with indirect collision models, such as the stochastic collision model used in the droplet collision algorithm of O’Rourke and Amsden,28 or the empirical pressure term frequently employed within the MP-PIC approach, e.g., refs 5, 7, and 11. Note that Patankar and Joseph7 also ran simulations with a continuum solids pressure term for comparison. A brief overview of the traditional soft-sphere approach3 is presented here. For more specific details on the model and its implementation in MFIX, the reader is referred to the work of Garg et al.29 In the soft-sphere approach, a linear spring-dashpot bt) bn) and tangential (F model is used to express the normal (F b components of the elastic collisional force (Fc). The spring accounts for repulsion between colliding particles and the dashpot allows for dissipation of kinetic energy (due to inelastic collisions). Accordingly, the normal and tangential components of the contact force are decomposed into the spring (conservative) force and the dashpot (dissipative) force. The net contact bc,p) is the sum of these forces force on a given particle/parcel (F Fn + b Ft) over all collisions of a parcel with other (i.e., b Fc ) b parcels and walls. In this model, the normal component of the collision force is written as

(8)

b Fn kn ηn ) δb n b V mp mp n mp n

f dω p )b Tp Ip dt

where b xp is the position of the pth parcel, Fp is the density of Fc,p is the net contact the parcel, and mp is the mass of the parcel, b force due to collisions of parcel p with other parcels and walls, Tp is the sum of all ω b p is the angular velocity of the parcel, b torques acting on the parcel, and the other terms are as defined earlier. The forces acting on a parcel, presented on the righthand-side of eq 7, are the gravitational force, a generalized buoyancy force, the drag force, and the contact forces (due to collisions). In this approach, all of the particles in a computational parcel are assumed to have the same velocity as the parcel bp). Thus, the terms “particle velocity” and “parcel velocity” (V may be used interchangeably. While different methods may be used to describe the parcel density, the parcel density is typically set equal to the material density of an individual particle (Fpt), as is the case in this effort (i.e., Fpt ) Fp). Since the density of a parcel is taken as that of a particle (i.e., effectively no void within a computational parcel), an equivalent parcel diameter (deq) can be calculated as deq ) (6npφpt/π)1/3 ) np1/3dpt, where np is the number of particles per parcel and φpt is the volume of an individual particle. The parcel mass can be readily calculated from deq as mp ) Fpdeq3π/6. Similarly, the parcel volume, found in eq 2, can also be calculated, for example, φp ) npφpt or φp ) π/6deq3. In regard to the drag force, the interphase momentum exchange coefficient is described by eqs 4 and 5 above, where the particle diameter is used rather than the parcel diameter so that the fluid-parcel interaction is similar to the fluid-particle interaction. This procedure is equivalent to the approach used in the literature when computational parcels rather than particles are employed.7,11 When considering parcel-parcel collisions, the same ad-hoc approach of Patankar and Joseph7 is followed, that is, a softsphere collision model. The term ad-hoc is applied since a softsphere algorithm is still being used to treat parcel interactions even though collisions between parcels are not directly described using collision models designed for particles. Therefore, in this study two types of errors are being estimated: (1) the errors generated due to lumping particles in a computational parcel

(9)

where kn and ηn are the normal spring stiffness and dashpot damping coefficients, respectively; b n is the unit normal vector Vn are the normal overlap along the line of contact; and δn and b and normal relative parcel velocity during collision, respectively. Using the differential equations for a damped harmonic oscillator, the dashpot damping coefficient can be related to the normal restitution coefficient by the following:30,31 ln(e) ) -

π(ηn /mp)

√2(kn/mp) - (ηn/mp)2

(10)

In eqs 9 and 10, the collision parameters (kn, ηn) are divided by the mass of the colliding particles. In polydisperse systems where particles with different masses are considered, an equivalent mass (m1 + m2/m1m2) of the two colliding particles is typically used in place of the particle mass in eq 10.32 While the current study only deals with monodisperse systems, comparisons between cases using parcels of different masses are conducted. Therefore, the spring stiffness and viscous damping coefficients per unit mass of the parcel are kept the same for simulations using different parcel sizes. bt) is The tangential component of the contact force (F δt)/mp] expressed similarly to the normal force b Ft/mp ) [(ktb Vt)/mp] and is only limited by Coulomb frictional force (µFn, [(ηtb where µ is the friction coefficient), where kt and ηt are tangential spring stiffness and damping coefficients, respectively; b νt is the tangential component of the relative contact velocity; and b δt is the tangential displacement vector. In this effort, the tangential collisional parameters are related to the normal collisional parameters. Specifically, the ratio of normal to tangential spring constant (kn/kt) is taken as 2/7, which can be found for the case of uniform spheres by equating the normal and tangential periods of oscillations for a damped harmonic oscillator.31,33 The ratio of normal to tangential dampening coefficient (ηn/ηt) is taken as 1/2, which has been employed by Silbert et al.33 These ratios are also used for determining wall values of the tangential collisional parameters. Note that both normal and tangential


parameters can be directly estimated from material properties (see the work fo Di Maio and Di Renzo34 for example). Since the current study does not involve any comparisons with experimental data, such treatment is unnecessary. Unlike the continuum governing equations, solving the equations of motion for discrete parcels (or particles) does not require a computational grid (such an approach is considered “meshless”). In the soft-sphere approach, a constant time step is used to update parcel velocities through a first-order integration scheme for Newton’s law. The parcel position is then updated based on its velocity. The DPM time step is typically taken as equal to the minimum collision time divided by 50, where the minimum collision time is calculated at the beginning of the simulation based on the properties of all the parcels. In coupled simulations, the simulation essentially proceeds in two parts. First, the continuum fluid equations (i.e., eqs 1 and 2) are solved for a fluid time step, wherein the particles are stationary. Second, the discrete particle simulation ensues, wherein the particles are advanced in time (eqs 6-8) until the time elapsed equals that of a fluid time step. At this point, the fluid and particles are advanced to the same time and the process repeats. Recall that the second goal of this study is to assess the MPPIC approach, wherein particle-particle collisions are ignored, e.g., refs 5, 7, and 11. That is, the last term of eq 7 is generally replaced by a normal stress acting on each particle (or parcel) with the main purpose to avoid exceeding maximum packing of the granular assembly. Hence, for the MP-PIC approach, the equation of motion for each individual parcel becomes bp dV ∇Pg βp ∇Ps bg(x )b g + (V b p) - b Vp) dt Fp Fp Fpεs

(11)

where εs is the solids volume fraction, Ps is the solids pressure, and the other terms are as previously defined. The expression of solids pressure is taken from previous studies7,11 and is modified to account for the fact that solids volume fraction may exceed the maximum packing

Ps )

{

10591

Figure 1. Instantaneous (after 60 s) snapshot of parcel position and velocity in the pseudo-3D shear cell. Particles colored in red indicate velocity magnitude higher than 1 m/s, and blue indicates a near-zero velocity magnitude. Parcel diameters are drawn to scale.

outside the flow domain due to large values of solids pressure. Although the value of maximum particles packing was set as 0.7, the computed overpacking was typically 0.701-0.703. In two simulations, the value of maximum packing was reduced to that of random packing (i.e., εsmax ) 0.6) with no significant changes in the time-averaged computed results. Note that as εsmax is reduced, a larger pressure is felt throughout the system (i.e., pressure is felt at all solids packing fractions and not just at εsmax). Furthermore, pressure is the only force that stops εs from exceeding εsmax . Part 1. Verification of the Coarse-Particle (Parcel) Model

P*

εsβ εsmax - εs

Psold

εs < εsmax

(12)

εs g εsmax

where εsmax is the solids volume fraction at maximum packing, P* is a constant with units of pressure, and β is also a constant (not to be confused with the interphase momentum exchange coefficient βp). Like the algorithm of Patakar and Joseph,7 but unlike that of Snider et al.,11 the solids pressure gradient in eq 11 is implemented explicitly without an implicit approximation for the solids volume fraction on the Eulerian grid. However, the values assigned to P*, β, and εsmax are similar to those typically employed by Snider et al.11 and are 100 Pa, 2, and 0.7, respectively. Larger values of P* were found to increase dispersion even at low solids volume fraction, and larger values of β reduced dispersion, as in agreement with the report by Snider et al.11 The current simulations were also observed to be very sensitive to the values of these parameters. For example, an order of magnitude increase in the value of P* was observed to push parcels beyond the flow domain and cause the simulation to terminate prematurely. A major issue, also related to the solids volume fraction exceeding maximum packing, was resolved by setting the value of the current solids pressure to that of previous discrete particle time step (Pold s ). This adjustment allowed parcels to slightly overpack but stopped parcels from being forced

The first part of this study focuses on assessing the accuracy of the parcel method, wherein a number of individual particles are grouped together. Two systems are employed to evaluate the errors associated with the parcel assumption: (1) a walldriven (Couette) granular shear flow under zero gravity and (2) a pressure-driven (Poiseuille) flow of gas and particles in a periodic vertical riser. Description of Granular Shear Flow Simulation (Part 1 Case 1). In this effort, parcels are sheared between two parallel walls moving with opposing velocities (see Figure 1 for reference). The frictional contact between the parcels and the moving walls drives the flow with momentum transferred away from the walls by subsequent parcel-parcel collisions. Thus, the walls act as an energy source and the momentum and energy of the parcels is lower away from the walls due to the inelastic nature of the parcel-parcel collisions. In this system, external forces are neglected and no fluid is present (i.e., in a vacuum under zero gravity). The input parameters for the simulation, including particle/ parcel properties and the system geometry, are summarized in Table 1. In the first portion of this work, four cases are examined, each with a different number of particles per computational parcel (np): 1, 2, 5, and 10 particles/parcel. Recall that the density of a parcel is set to that of a particle (no void fraction within a parcel), so an equivalent parcel diameter (deq)

10592


Table 1. Simulation Conditions for a Granular Shear Flow System quantity particle ) parcel density particle diameter particles per parcel parcel diameter

symbol

units

value

Fpt ) Fp

kg/m

2900

dpt np deq

m

0.01 1, 2, 5, 10 0.01, 0.0125, 0.0171, 0.0215 5730, 2865, 1146, 573 0.9 0.5 5.3 × 106

3

m

total number of parcels nt restitution coefficient friction coefficient spring constant/ parcel mass wall velocity system geometry

e µ k/mp

1/s2

Uw m/s (12.5 Lx × Ly × Lz m × m × m 0.2 × 1.0 × 0.1

Table 2. Relative Simulation Speeds for Granular Shear Flow System particles per parcel

1

2

5

10

simulation speedup

1

3.6

26

102

is calculated as deq ) (6npφpt/π)1/3 ) np1/3dpt. Note that, although the parcel size changes, the total solids volume fraction of the system (εs ) 0.15) is maintained. The results of the discrete simulation are sensitive to the value of the spring stiffness (kn, kt). To remove the effect of the spring constant value on the simulation results, the value of the spring constant per unit mass of parcel is kept constant between all four cases. Another way to deal with this issue, although computationally very expensive, is to use relatively large values of the spring stiffness such as those used by Ye et al.35 who found good agreement between the soft- and hard-sphere collision models. To accommodate the larger volume of a parcel (compared to a single particle), pseudo-3D geometry with a depth of several particles diameters in the z direction is necessary. Such an approximation of the flow geometry has been used elsewhere; for example, see the work of Lohse et al.36 The quantities Lx, Ly, and Lz represent the simulation domain length in the x, y, and z directions, respectively. The simulation is bounded on the left (x ) 0) and right (x ) Lx ) 0.2 m) by solid walls moving in opposing y-directions with speed equal (Uw. The value of Uw is of comparable magnitude to that used by Karion and Hunt32 for a similar problem setup. Initially the parcels are randomly placed on a nearly cubic lattice and are assigned velocities from a Gaussian distribution with zero mean and a standard deviation of 0.2 m/s. The initialization ensures that the parcels will begin to move into contact with the moving walls so that a shear flow develops as a consequence. The parcel-wall collision properties are set equal to those used for the parcel-parcel collision properties (i.e., the same k, e, and µ). The four remaining sides, namely the upper, lower, forward, and backward boundaries, are standard periodic conditions (e.g., any particle that leaves the top side will re-enter from the bottom side). Simulation Results for Granular Shear Flow System (Part 1 Case 1). The main motivation for employing the parcel method is the significant reduction in simulation run time compared to a full particle simulation. Table 2 gives the relative speed of the three simulations using the coarse parcels compared to that using 1 particle per parcel (i.e., using real particles). As is evident, the parcel technique yields faster simulations with a speedup of two orders of magnitude for the simulation with 10 particles per parcel. Before any new technique is considered reliable, however, it should be carefully verified for accuracy, which is one goal of the current effort.

Figure 2. Time-averaged (30-90 s) flow variables across the width of the 3D shear cell. The DPM results (line) using real particles are compared with those using computational parcels containing 2 or more particles (symbols). These graphs are scaled using dp ) 0.01 m, Uw ) 25 m/s, and Lx ) h ) 0.2 m.

A snapshot of the shear system is presented in Figure 1. As is evident, the simulation results for all four cases show some qualitatively similar behavior for momentum and density. In particular, the parcels near the moving walls have a relatively high velocity (magnitude) due to the frictional contact, and the parcels in the core of the shear cell are lower due to dissipative parcel-parcel collisions (both normal and tangential collisions are dissipative). Accordingly, higher parcel concentration is observed at the core of the shear cell where the parcel energy is lower. At first glance, the behavior between the four cases appears to be similar; however, Figure 1 also shows that momentum of the parcels increases as the number of particles/ parcel increases (going from left to right). To better assess the quantitative differences between the four cases, the time and spatially averaged lateral profiles in solids


volume fraction (εs), y component of velocity (U), and granular temperature (T) are reported in Figure 2. To generate these figures, the data was collected into cells (or bins/strips) for a period of time and then averaged. Specifically, the domain was discretized using 10 computational cells along the x-direction between the shear walls (Lx) and 1 cell along the length (Ly) and depth (Lz) of the system. Since these simulations do not involve the solution of any fluid flow continuum equations, the computational cells serve no purpose except for data collection and subsequent post-processing (averaging the discrete data). All data was collected every solids time step in each strip starting at 30 s and ending at 90 s after which the temporal averages were calculated for each strip. The average solids volume fraction within each strip was determined by including the volume of parcels whose center resides in the given strip at the instant of measurement. Similarly, the y-component of velocity and granular temperature in the computational strip were found by including only the value of y-velocity component and temperature of parcels whose center resides in the strip at the instant of measurement. For example, the strip averaged N Vp,y/N, where y-velocity component was calculated as U ) ∑p)1 Vp,y is the y-component of the parcel velocity, and N is the number of parcels whose center resides within the strip. In this effort, the granular temperature is defined as T ) 1/3Cp2, where Cp is the magnitude of the fluctuating particle velocity relative to the local mass average velocity. The local mass average velocity in the computational cell was evaluated similarly, that is, by including only those parcels’ velocities whose centers reside in the cell at the instant of measurement. In all four cases, the solids volume fraction is lowest near the moving walls, which as previously mentioned, is due to energy acquired through frictional parcel collisions with the walls. The high-energy parcels then collide with surrounding parcels, effectively pushing them to the center of the system where they subsequently “cool” due to dissipative collisions. This behavior is typical for granular Couette flow where particle motion is driven by a frictional wall (see the work of Benyahia,37 for example). Thus, parcel velocity and granular temperature are highest near the walls and lowest at the center of the shear cell. While the solids volume fraction profiles between the four cases are quantitatively similar, as the parcel size increases, the solids volume fraction becomes slightly larger at the walls and lower at the center. Unlike the solids volume fraction profiles, significant differences are evident in the granular temperature and y-velocity profiles between the four cases. Note that a log-scale is used to more clearly represent the temperature profiles. The granular temperature and the magnitude of velocity are highest across the domain for the coarsest parcel case (10 particles/parcel). These differences may be explained by considering what happens as particles are grouped together into parcels. Namely, as the number of particles per parcel increases, the total number of parcels decreases (while the total solids volume fraction is maintained), and this decrease will affect the collision frequency. With fewer collisions the overall amount of dissipation due to parcel-parcel collisions will be lower, and conversely, the granular energy will be higher. Since collisions between parcels are a major mechanism of momentum transfer for this walldrive shear flow, the difference in the computed granular temperature and velocity profiles between the 10-particle-perparcel case and the 1-particle-per-parcel case is not unexpected. Table 3 gives a quantitative measure of the differences in the flow field variables between the four cases, specifically, the fractional error associated with the parcel method. The relative

10593

Table 3. Average and Maximum Errors (Fractional) in Three Parcel Size Cases Relative to the 1-Particle-Per-Parcel Case for Granular Shear Flow particles per parcel error

2 avg

5 max

avg

10 max

avg

max

volume fraction 0.018 0.037 0.036 0.077 0.065 0.12 velocity 0.36 0.95 3.60 4.90 11.0 16.0 granular temperature 1.80 2.0 67.0 75.0 615.0 740.0

error is computed using the simulation results for the 1-particleper-parcel case as a baseline since it is considered the most accurate simulation. The smallest errors are in parcel concentration, which is not surprising since the averaged solids volume fraction is the same for all simulations. The largest relative errors are in granular temperature followed by the y-component of the velocity. Even for the 2-particles-per-parcel simulation, large relative errors are computed. The results not only demonstrate the importance of collisions for granular shear flow but also indicate that, for such cases, the technique of grouping particles together into parcels induces large errors that may be unacceptable depending on the end use. It is worth noting that both of the smaller parcel cases (1- and 2-particle-per-parcel cases) were still losing energy during the time averaging interval of up to 90 s and so were not at a steady state (the larger parcel cases had reached a steady state). If the two cases had been continued until steady state, however, the subsequent comparisons with the larger parcel cases would reflect even greater errors. Finally, these simulations were repeated with an increase in total solids volume fraction (εs ) 0.30), wherein all four cases achieved a steady state within the time interval and the errors followed similar trends (results not shown) as those with the lower initial solids volume fraction (εs ) 0.15). Description of Periodic Riser Flow Simulation (Part 1 Case 2). In this section, the parcel assumption is evaluated in the context of a more complex flow system, namely, the flow of Geldart group B particles in air through a periodic riser (see Figure 3 for reference). The flow of particles is driven against the gravitational force by a pressure drop in the gas phase (Poiseuille flow). The input parameters for this simulation, including the physical properties of the gas phase and parcels as well as

Figure 3. Instantaneous parcels position and velocity in the pseudo-3D periodic riser. Particles colored in red indicate velocity magnitude higher than 3 m/s, and particles colored in blue indicate a near-zero velocity magnitude. Parcels diameters are drawn to scale. These images represent the complete pseudo-3D domain (i.e., not a sample or a slice).

10594


Table 4. Simulation Conditions for Periodic Riser Flow System quantity

symbol

gas density and viscosity particle ) parcel density particle diameter particles per parcel parcel diameter total number of parcels restitution coefficient friction coefficient spring constant/parcel mass superficial gas velocity system geometry Eulerian grid Eulerian & DPM time steps

Fg and µg Fpt ) Fp dpt np deq nt e µ k/mp Vsf,g Lx × Ly × Lz i×j×k dt and dts

numerical parameters, are summarized in Table 4. Four cases are examined, each with a different number of particles per computational parcel (np) ranging from 10 to 100 particles/ parcel. The smallest parcel considered in this study contains 10 particles due to computational time limitations for running a more refined simulation. Although a case with only one particle/parcel was not used, this study still clearly demonstrates the effect of increasing the parcel size on the predicted flow behavior. The same total solids volume fraction (εs ) 0.05) is maintained for all five cases. As done earlier, the parcels start by being randomly placed on a nearly cubic lattice and their velocities are assigned from a Gaussian distribution with a zero mean and a standard deviation of 0.2 m/s. The value of the spring constant per unit mass of parcel is kept constant between all cases. As in the granular shear flow case, pseudo-3D geometry is employed; however, the term pseudo-3D flow system takes a slightly different meaning here than in the previous case due to the gas-solids coupling. Namely, for this case the solution for the fluid flow is obtained in a 2D system because only one computational point lies in the z-direction. The simulation is bounded on the left (x ) 0) and right (x ) Lx ) 0.1 m) by walls. Again the parcel-wall collision properties (k, e, and µ) are set equal to those used for the parcel-parcel collision properties. A no-slip boundary condition is used at the walls for the gas phase. The remaining four sides are treated as standard periodic boundaries. The main flow of gas and particles occurs in the positive y-direction (upward) against the gravitational force that acts in the negative y-direction (downward). The gas flow rate is held constant with a superficial velocity of 6.0 m/s. To maintain a constant gas flow rate, the pressure drop is automatically adjusted every iteration. Simulation Results for Periodic Riser Flow (Part 1 Case 2). Figure 3 shows a snapshot of the periodic riser flow system. All four cases show some qualitatively similar behavior. That is, large clusters of parcels form mostly near the walls of the riser, although occasionally clusters are observed in the core of the riser, and larger concentrations of parcels are present in the annulus (near the walls). Additionally, parcels with the highest velocity magnitude (colored yellow-red) tend to flow in the core region of the riser where the gas momentum is the highest. Animations of the simulations illustrate how a cluster may occasionally form near one wall, flow downward due to gravity, and after a few seconds begin to disappear and form on the opposite side of the riser. As evident in the figure, two clusters may also exist simultaneously at both walls of the riser. Figure 3 also reveals some differences between the simulations using different parcel sizes (e.g., compare the panels in Figure 3 for the 10 particles/parcel case to the 100 particles/parcel case). As the parcel size becomes smaller (or the total number density

units 3

kg/m and Pa s kg/m3 µm µm

1/s2 m/s m×m×m s

value 1.093 and 1.95 × 10-5 2500 500 10, 20, 50, 100 1077, 1357, 1842, 2320 30576, 15286, 6112, 3059 0.95 0.1 1.8 × 107 6.0 0.1 × 0.4 × 0.01 15 × 60 × 1 1.0 × 10-4 and 1.0 × 10-5

of particles becomes larger), the clusters become more dense. Moreover, simulations with smaller parcel size generally show parcels in all flow regions of the domain, while the simulations with coarser parcel size show regions with no parcels (more open/void space). To assess the quantitative differences between the different cases, the time-averaged lateral profiles in solids volume fraction (εs), y-component of gas and solids velocity (Vg and Vs), and granular temperature (T) profiles are reported in Figure 4. To obtain these figures, the data was collected into cells for a period of time and then averaged. The Eulerian grid indicated in Table 4 was employed for data collection (both discrete and continuum). The discrete (parcel) data, y-component of solids velocity, and granular temperature are spatially averaged in each Eulerian computational cell, and the calculations of the ycomponent of solids velocity and granular temperature follow the same procedure used for in Figure 2 for the granular shear flow system, whereas, the calculation of solids volume fraction differs. In particular, the solids volume fraction in a computational cell includes only the parcels volume that resides within the cell (rather than the entire volume based on the parcel center). While minor variations in the flow variables occur in the flow-wise (y) direction, the data is spatially averaged over that direction so only changes across the lateral direction are presented. In some instances, parcels may not be present in a given computational cell. In this event, the spatial average (over y cells) for y-velocity and granular temperature may be handled differently. Two possible options are to set the solids y-velocity and granular temperature to zero in those cells, or to ignore these cells, as they occur, in the averaging process. The latter was considered in Figure 4. If the cells with zero values were used in the spatial averaged data, then the magnitude of solids y-velocity and granular temperature would be slightly less. As observed in the granular shear case, as the number density of particles increases in the riser by considering smaller parcels (20 particles/parcel or less), parcels will exist in any given Eulerian cell. In such cases, the two approaches will yield the exact same averaged results. After an initial transient of 5-10 s, all data was collected every 0.01 s in each cell for a sufficient time (30-50 s). At this point, the spatial, then temporal, averages were calculated. The averaging time period is usually considered sufficient when symmetric profiles along the width of the riser are obtained, as those shown in Figure 4. The solids volume fraction profiles, presented in Figure 4, show similar trends for all four cases using different parcel sizes and reflect what was observed qualitatively in Figure 3. Namely, core-annulus flow is predicted, wherein high solids concentration occurs near the walls and relatively dilute conditions occur at the core of the riser. This flow regime occurs under a wide range of flow conditions,38 and similar flow behavior (core-


10595

Table 5. Average and Maximum Errors (Fractional) for Three Parcel Size Cases Relative to 10-Particles-Per-Parcel Case for Periodic Riser Flow particles per parcel

Figure 4. Time and spatially averaged (along the vertical y-direction) flow variables computed across the width of the 3D periodic riser. The width of the riser is made dimensionless using the particle diameter dpt ) 500 µm.

annulus flow) has been predicted using a continuum modeling approach for the particles based on granular kinetic theory.39

20

50

100

error

avg

max

avg

max

avg

max

volume fraction velocity granular temperature

0.056 0.49 0.057

0.14 0.62 0.11

0.088 0.93 0.12

0.18 7.1 0.27

0.113 1.75 0.41

0.25 11.0 0.56

While the results are generally similar, the concentration of parcels near the walls is lower in the coarser parcel simulations than the smaller parcel simulations (i.e., the clusters are more dense in the simulations with smaller parcel sizes). The difference in solids concentration near the wall may be partly due to a possible sensitivity in the measurement to the computational cell volume, which remained constant as the parcel size increased. This sensitivity could be reduced by collecting particles-related data on a grid independent of the mesh used for solving the Eulerian equations. Unlike the solids volume fraction profiles, significant differences are evident in the profiles of granular temperature and the y-component of velocity between the simulations with differing particle sizes. The profiles corresponding to the coarsest parcel case (100 particles/parcel) display the greatest deviation from the smallest parcel case (10 particles/parcel). All of the granular temperature profiles exhibit a maximum away from the walls except for the simulation with the coarsest parcel (100 particles/parcel), which shows a maximum near the walls. As the parcel size is reduced, an overall increase in granular temperature is obtained (e.g., compare the 20 particles/parcel case to the 10 particles/parcel case). This result may be nonintuitive since granular temperature is generally expected to decrease as clusters become denser due to increased dissipation. This behavior can be explained by considering what is happening across the riser as the cluster density increases. Recall that as the parcel size decreases the total number density increases, and as we have seen in Figure 3, the predicted clusters are more dense. The denser clusters induce higher downward velocities and steeper vertical velocity profiles (as seen in Figure 4), which in turn, produce higher values of granular temperature due to shear. Similar results were previously computed using the continuum approach based on granular kinetic theory.39 In addition, the large downflow of clusters near the walls entrain gas; therefore, the gas velocity is lower in this flow region, with the lowest gas velocity observed in the simulation with the smallest parcel size. Conversely, since a constant average gas flow rate is maintained in these simulations, the highest gas velocities are computed at the core of the riser, and are higher as the parcel size is reduced. Table 5 gives the relative errors obtained using the coarse parcel simulations with respect to the simulation with the smallest parcel (10 particles/parcel). The largest errors are in parcel velocity, and the smallest, in parcel concentration. Although the errors using coarse parcels can be quite large, these errors are significantly lower than those computed for the granular shear flow (see Table 3) where collisions are an important mechanism of momentum transport. In this periodic riser flow system, gas-particle (parcel) drag force plays a significant role on parcel transport and may reduce the relative importance of collisions, which are not fully accounted for in coarse parcel simulations. Table 6 gives the relative speed of the coarser parcel simulations compared to the smallest parcel size simulation

10596


Table 6. Relative Simulation Speeds for Periodic Riser Flow System particles per parcel

10

20

50

100

simulation speedup

1

3.3

11

15

considered in this study (10 particles/parcel). As in the granular shear flow system, simulations are faster with coarser particles and larger speedups are obtained as coarse parcels are used. Part 2. Verification of the Multiphase Particle-in-Cell (MP-PIC) Approach The second part of this study focuses on verification of a basic assumption in the MP-PIC approachsthe use of a continuum solids pressure term in place of resolving individual particle-particle collisions (referred to in this effort as the MPPIC assumption). As a result, the particle (or parcel) force balance includes a term proportional to the gradient of solids pressure whose main function is to prevent the simulation from exceeding maximum particle packing. For the simulations of this study, the concept of parcels (instead of particles) is employed to reduce the computational times. Coarse-grained simulations (using parcels) are commonly used together with the MP-PIC assumption. To evaluate the impact of the MPPIC assumption, results from simulations based on the MPPIC approach are compared to those using a DPM technique where collisions between parcels are considered. In conducting these comparisons, the DPM method is assumed to be correct, meaning that the DPM model has the right physical models and predicts the right quantitative flow behavior. While nothing is absolutely correct, the MP-PIC approach is, in itself, an approximation of the DPM approach, and its main assumption is that parcel-parcel (or particle-particle) collisions are unimportant and can therefore be ignored. The purpose of this effort is to estimate the errors incurred in neglecting collisions. Two systems are employed in this study: (1) a 2D bubbling fluidized bed using a uniform air distributor and (2) gas-solids flow in a 2D periodic vertical channel. Description of Bubbling Fluidized Simulation (Part 2 Case 1). The 2D IIT bubbling fluidized bed used in this study was first simulated by Ding and Gidaspow1 and has also been simulated by Snider et al.11 While the results of the current MPPIC simulations are evaluated against those of the DPM simulations, they may also be compared against those presented by Snider et al.11 The input parameters for this simulation, including the physical properties of the gas phase and parcels, the numerical parameters, and system geometry, are summarized in Table 7. As indicated, a single test case characterized by 8000 total parcels with 668 particles per parcel is examined using both

simulation methods (DPM and MP-PIC). The initial parcel distribution is automatically generated using MFIX-DPM, wherein the particles are seeded on a nearly cubic lattice starting near the origin then filling in layers first along the x-direction followed by the y-direction until the specified solids volume fraction is attained. The solids volume fraction (εs) is set to 0.4 corresponding to an initial bed height of 0.5 m. The parcel velocities are initialized to zero. As mentioned in the first portion of this study, simulations employing parcels with different sizes should be conducted in 3D to maintain an equivalent volume fraction of particles in the system and also to allow for realistic particle-particle interactions (i.e., spherical particles imply a 3D system). In this part of the study, however, the results from only two simulations with the same numerical settings (i.e., same number of parcels and particles/parcel) are being compared, so the specific geometry is inconsequential. Thus a 2D geometry, which was also employed by the earlier studies,1,11 is acceptable. The simulation is bounded on the left (x ) 0) and right (x ) Lx ) 0.4 m) by walls. Note that while the MP-PIC assumption involves neglecting parcel-parcel collisions, a soft-sphere collision algorithm is still used for parcel-wall contacts, as is done in DPM, to avoid parcels moving across the walls. The parcel-wall collision properties for e and µ are equal to those used for the parcel-parcel collision properties, and the spring constant is increased by a factor of 100 for the walls. A no-slip boundary condition is used at the walls for the gas phase. At the inlet, the bed is uniformly fluidized with air at 1 m/s, which is equivalent to about four times the minimum fluidization velocity (Umf).1 At the outlet, atmospheric pressure is specified. Simulation Results for 2D IIT Bubbling Fluidized Bed (Part 2 Case 1). Figures 5 and 6 show a series of successive snapshots, spaced 0.5 s apart, of parcel position in the 2D fluidized bed using the soft-sphere DPM approach and the MPPIC approach, respectively. Two differences are clearly noticeable when comparing the two figures. First, the MP-PIC predicts greater bed expansion than the DPM method with the number of particles ejected being higher in the freeboard. Second, bubbles are not as clearly defined (more diffuse) in the MPPIC simulation. Similar bubble behavior was predicted in the work by Snider et al.11 (see Figure 36 of Snider’s work, page 63), which is not surprising since the same solids pressure model is used here as in their study. Qualitatively, the MP-PIC model, with its solids pressure term, does not seem to predict the same bubbling behavior observed in the DPM collision model. The diffuse behavior of the MP-PIC technique may be attributed to the occurrence of excessive parcel overlap in the simulation, which in turn stems from the neglect of parcel-parcel

Table 7. Simulation Conditions for 2D IIT Bubbling Fluidized Bed quantity

symbol

gas density and viscosity particle ) parcel density particle diameter particles per parcel parcel diameter total number of parcels restitution coefficient friction coefficient parcel-parcel spring constant/parcel mass empirical parameters maximum packing superficial gas velocity system geometry Eulerian grid Eulerian and DPM time steps

Fg and µg Fpt ) Fp dpt np deq nt e µ k/mp P* and β εsmax Vsf,g L x × Ly i×j dt and dts

units 3

kg/m and Pa s kg/m3 µm µm

1/s2 Pa and m/s m×m s

value 1.093 and 1.95 × 10-5 2500 500 668 4370 8000 0.99 0.1 1.8 × 107 100 and 2 0.7 1.0 0.4 × 1.0 15 × 30 ∼2.0 × 10-4 and 1.64 × 10-5


Figure 5. Instantaneous parcel position in the uniformly fluidized bubbling bed using the soft-sphere discrete particle method (DPM) approach. Large bubbles form as is expected when group B powder is fluidized with a high gas velocity.

Figure 6. Instantaneous parcel position in the uniformly fluidized bubbling bed using the multiphase particle-in-cell (MP-PIC) approach. Bubbles are diffuse and their shape not clearly defined. These results are qualitatively similar to the predictions of Snider et al.11 (see their Figure 36).

(or particle-particle) collisions. The extent of parcel-parcel overlap is illustrated in Figure 7, wherein a small portion of the MP-PIC simulation domain is magnified. As a result of this overlap, the locally averaged volume fraction is over predicted, and, in turn, the solids pressure becomes artificially large. Note that the MP-PIC solids pressure cannot prevent parcels from overlapping, since it is derived from the average solids concentration in each Eulerian grid. The only forces that prevent parcels (or particles) from interpenetrating are collisional in nature, which are neglected in the MP-PIC model. Figure 8 corroborates the observations made from Figures 5 and 6, that is, the overall fluidization behavior between the DPM and MP-PIC methods differ, as indicated by the differing amplitudes in the transient void fraction profile taken at a specific location in the bed. As is evident, the minimum specified void fraction (εg ) 0.3), or conversely the maximum solids packing (εsmax ) 0.7), is occasionally reached in the MP-PIC simulation. The mean void fraction corresponding to the DPM and MP-PIC data in Figure 8 is 0.69 and 0.76, respectively.

10597

Figure 7. Instantaneous snapshot (t ) 5.0 s) of parcel position in the uniformly fluidized bubbling bed using the MP-PIC technique. The figure to the right is a magnification of the parcels in the red square showing parcel clustering and significant overlapping.

Accordingly, the MP-PIC approach results in a more dilute (homogeneous) flow, which is attributed to the high solids pressure that tends to push particles apart and eject them high in the freeboard (as illustrated earlier in Figure 6). While Figure 8 is for a single point, this behavior is also reflected later in Figure 9 where the MP-PIC exhibits higher voidage in the bottom region and an overexpansion of the bed. The corresponding standard deviations in the void fraction data of Figure 8 are 0.136 and 0.143 for the DPM and MP-PIC methods, respectively. Thus, at this location in the bed, the errors in the mean and standard deviation of void fraction associated with the MP-PIC method relative to the DPM method are 10.1% and 5.1%, respectively. Another interesting result that can be obtained from Figure 8 is the bubble frequency which is extracted from the void fraction fluctuations by counting the number of bubbles that occur during the time interval 5-15 s. In this effort, a bubble is defined as flow regions of 80% void fraction and higher. Thus, the DPM and MP-PIC techniques predict bubble frequencies of 3 and 4.7 Hz, respectively. The time-averaged (from 5-15 s) gas pressure and void fraction profiles along the height of the fluidized bed are presented in Figure 9. The pressure profiles for the DPM and MP-PIC approaches are similar; however, the MP-PIC simulation predicts a slightly higher pressure drop, which may be related to the ejection of particles into the freeboard. That is, particles will lose energy by hitting the top boundary of the bed (while the top boundary is an outlet for the fluid phase, it acts as a wall for the particles as they are not allowed to leave the system). This overexpansion, previously noted in comparisons between Figures 5 and 6, is also evident in the void fraction profiles. Namely, the DPM approach predicts denser flow behavior near the bottom regions of the bed, while the MP-PIC simulation predicts some solids present all the way up to the outlet (y ) 1 m). As mentioned earlier, the primary goal of this second effort is to assess the MP-PIC approach using the DPM approach as a baseline. Table 8 gives the average and maximum errors in the MP-PIC approach relative to the DPM approach for the data presented in Figure 9. The errors in void fraction are calculated using the time-averaged data along the entire height of the bed (from y ) 0 to 1 m). In contrast, the errors in gas pressure are

10598


Figure 8. Transient void fraction from inside the 2D bubbling fluidized bed at x ∼ 0.2 and y ∼ 0.3 m using the DPM (with collisions) and MP-PIC methods (with continuum solids pressure and no collisions).

Figure 9. Time (5-15 s) and spatially (across the width of the bed) averaged gas pressure and gas volume fraction along the height of the 2D fluidized bed using the DPM and MP-PIC techniques. Table 8. Average and Maximum Errors (Fractional) for MP-PIC Approach Relative to DPM Approach for the 2D Bubbling Fluidized Bed quantity

average error

maximum error

volume fraction gas pressure bed expansion

0.04 0.24 0.67

0.09 0.54

calculated using only the time-averaged data up to y ) 0.6 m. Above that height, the errors in the time-averaged gas pressure become significantly larger as the gas pressure approaches atmospheric (zero gauge) for the DPM simulation. As illustrated in Figure 9, the DPM approach predicts a lower bed expansion (about 0.6 m), while the MP-PIC approach predicts a bed expansion as high as the maximum height (1 m). In lieu of resolving the forces involved in individual parcel-parcel (or particle-particle) collisions, the MP-PIC approach employs a continuum solids pressure. The expression for the solids pressure model used in the MP-PIC model of this effort (see eq 12) and that of Snider et al.11 requires setting a value for the maximum solids packing as an input. As indicated by Table 7, this value was set to 0.7. Figure 8 demonstrates that the MP-PIC simulation will, at least at times, predict regions characterized by this maximum packing. In contrast, the transient profile in void fraction corresponding to the DPM simulation shows lower amplitudes (less deviation) in the solids packing. ) 0.7 was based on the previous study of Snider The value εmax s et al.;11 however, a lower value would be a more reasonable representation of random close packing (and more congruent with the behavior of the DPM simulation). For this reason, the MP-PIC bubbling bed simulation was repeated using εsmax )

0.6 and the resulting transient void fraction profile is shown in Figure 10 along with the profile when εsmax ) 0.7. As expected, the amplitudes in void fraction are generally smaller in the simulation with the lower maximum packing (εmax s ) 0.6), which is attributed to the correspondingly larger pressure force that acts on the granular assembly at all solids volume fractions. A bubble frequency of 3.3 Hz was obtained in this case, which is closer to the value predicted using the DPM technique (Figure 8). The gas pressure profile obtained with the lower packing value was similar to that shown in Figure 9 and will not be reproduced. While the empirical constants P* and β may also be adjusted, as mentioned earlier in this paper, they were maintained constant as their impact on bed behavior is known and their values were selected based on previous research in this area.11 Description of 2D Periodic Channel Flow Simulation (Part 2 Case 2). In this section, the MP-PIC assumption is evaluated in the context of a 2D periodic channel. The premise is similar to the periodic riser flow employed in part 1 of this effort. In particular, a flow of particles is driven against a gravitational force by a pressure drop in the gas phase (Poiseuille flow). For reasons explained in the description of the bubbling fluidized bed simulation, the geometry here is 2D rather than pseudo-3D used in part 1 (assessing the MP-PIC assumption rather than effect of parcels). The current system employs many of the same physical and numerical parameters as used in the 2D bubbling fluidized bed system, and Table 9 summarizes the input parameters that differ from Table 7 (similarities are not reproduced). Accordingly, the system geometry is doubled in size in the vertical direction in order to decrease the average solids concentration in the system by half (i.e., εs ) 0.1). The restitution coefficient of inelastic collisions is decreased, which will induce higher dissipation than present in the 2D bubbling bed case. The superficial gas velocity is set to 6 m/s, and, as before, the pressure drop is automatically adjusted each iteration to maintain constant gas flow rate. Simulation Results for 2D Periodic Channel Flow (Part 2 Case 2). Qualitative differences between the MP-PIC and DPM approaches are clearly represented in Figures 11 and 12, which show a series of successive snapshots of parcel position spaced 1.0 s apart. The DPM simulation shows a large cluster that forms on one side of the channel and moves to the other side. The intermittent formation of large clusters in particle laden airflow through a 2D periodic channel has been reported previously in the literature.39 This behavior is not clearly captured in the MP-PIC simulation. In particular, the shape of the cluster is different in the MP-PIC approach, which is due to the diffusive nature of the solids pressure applied to the


10599

Figure 10. Transient behavior of gas void fraction inside the 2D bubbling fluidized bed at x ∼0.2 m and y ∼ 0.3 m obtained using the MP-PIC technique for two different values of maximum solids packing: 0.6 and 0.7. Table 9. Simulation Conditions for 2D Periodic Channel Flow quantity restitution coefficient superficial gas velocity system geometry Eulerian grid

symbol e Vsf,g Lx × Ly i×j

units m/s m×m

value 0.95 6.0 0.4 × 2.0 15 × 50

particles force balance as seen previously in the bubbling bed case. The appearance of a more diffuse cluster near the walls in the MP-PIC results may be misleading due to the overlapping of particles also observed in this simulation and indicated previously in Figure 7. To better assess the differences between the two approaches, the time-averaged lateral profiles in solids volume fraction (εs), y-component of gas and solids velocity (Vg and Vs), and granular temperature (T) profiles are reported in Figure 13. To obtain these figures, the data was collected into cells for a period of

Figure 11. Instantaneous parcel position in the 2D periodic channel flow using the DPM soft-sphere approach. A large cluster of particles forms at one side of the channel and then moves to the other side.

Figure 12. Instantaneous parcel position in the 2D periodic channel flow using the MP-PIC approach. The formation of a core-annulus flow is qualitatively captured, but the shape of clusters captured by the DPM approach (Figure 11) is not qualitatively reproduced with the MP-PIC approach.

time starting at 10 s and ending at about 60 s and then averaged. (See the periodic riser flow system presented in part 1 for details on calculating these average results.) The 50-s time interval is long enough to obtain fairly symmetric results along the channel width, although some profiles (e.g., solids volume fraction) may show better symmetry with additional time-averaging. It is worth noting that the major differences between the MP-PIC and DPM results cannot be explained by the duration of time-averaging. Another important fact to mention is that similar to previously obtained continuum predictions based on granular kinetic theory,39 the discrete approach also demonstrates that core-annulus flow arises from time-averaging the solids concentration in the riser. That is, cluster formation along the side walls is a dynamic process, but when time-averaged, the clusters that occasionally form on either side of the channel shape a profile with a larger concentration of solids near the walls of the channel and a low concentration in the core region. Figure 13 reveals that the MP-PIC approach predicts a more dilute core region, which produces larger gas and solids velocities in this flow region and is slightly denser in regions away from the narrow core region. The MP-PIC approach also predicts negative solids velocity near the walls. However, this behavior cannot be explained by the solids concentration profiles since similar solids concentrations are also predicted by the DPM approach in this flow region. To help explain the phenomenon of negative solids velocity near the wall in the MP-PIC simulation, the transient profile of solids volume fraction in the region next to the wall is examined in Figure 14. In particular, the results show that in the DPM simulation large clusters occasionally form (see the larger peaks around 15 and 30 s and also Figure 11), but the solids concentration near the wall is relatively lower for most of the simulation. On the other hand, the MP-PIC approach exhibits relatively smaller fluctuations (i.e., steadier concentration of solids), which in turn lead to steadier down flow of solids near the wall. The resulting time-averaged solids velocity near the wall reflects this behavior. In addition to different axial velocity profiles, Figure 13 reveals unlike granular temperature profiles for the DPM and MP-PIC techniques. The MP-PIC approach is characterized by significantly larger values of granular temperature than those computed by the DPM technique, which is attributed to the correspondingly larger solids velocity gradients present in that system. The velocity and granular temperature profiles predicted by the simulation (DPM or MP-PIC) will depend on the specific particle and gas properties employed, as well as the operating conditions. Regardless, the DPM approach will yield more accurate results, and for a given system, any differences between the two can be considered as errors resulting from the MP-PIC

10600


assumption. Table 10 gives the average and maximum errors in the MP-PIC approach relative to the DPM approach for the data presented in Figure 13. Although both techniques predict the core-annulus flow regime reasonably well, large deviations are evident in the MP-PIC simulation, wherein parcel-parcel collisions are neglected. While Tables 8 and 10 reflect large errors in the MP-PIC approach relative to the DPM technique for the two systems examined, the value of these errors may depend on the specific settings for several parameters. As discussed earlier, the MPPIC predictions depend on the continuum solids pressure, which requires setting values for P*, β, and εsmax in the empirical expression for the continuum solids pressure. Unfortunately, little measured data is available for the empirical parameters P* and β, and in this effort, their values were selected to match those used by Snider et al.11 who have already optimized them as explained earlier. Another quantity that will impact the relative error is the value of the restitution coefficient. The restitution coefficient is a measure of the particle hardness (a material property), and its value is defined by the type of particles being simulated. The current simulations, however, are not being employed to model a real system with a prescribed restitution coefficient. Instead, the focus of this effort was to assess the MP-PIC assumption using the DPM technique as a baseline. Thus, the value of the restitution coefficient may be considered rather arbitrary. One trouble with this view, though, is that the DPM approach does resolve individual collisions while the MP-PIC approach does not. Thus, the DPM results are likely to be sensitive to the specific value of restitution coefficient. In contrast, the MPPIC approach does not depend on e except for describing particle-wall interactions. From previous experience,39 lower values of restitution coefficient are known to result in larger clusters and more solids down flow. The behavior is affirmed in the time-averaged DPM results presented in Figure 15, in which the DPM simulation was repeated with a lower restitution coefficient. As illustrated, the simulation with the lower restitution coefficient is characterized by slightly denser solids concentration in the wall region and lower solids velocities in the annulus. The latter is attributed to the formation of denser/ larger clusters of particles in the wall region. Thus, in this case the errors in the MP-PIC simulation relative to the DPM simulation will be less with the lower restitution coefficient. While a better match between the two simulations may be obtained, it should be kept in mind that the restitution coefficient is not a fitting parameter but a physical property of the particle. One final point worth noting is that parcels, rather than particles, are being employed in these simulations, which is also an approximation. Thus, the restitution coefficient based on particle properties does not necessarily translate to a restitution coefficient for parcels of the same particle type. This observation is further explored in the following section. Part 3. Simple Improvement to the Coarse-Particle Model

Figure 13. Time and spatially averaged (along the vertical y-direction) flow variables across the width of the 2D periodic channel using the MP-PIC and DPM technique.

In the last part of this study, a simple improvement to the coarse-particle (parcel) model is proposed. The modification is based on the observation that simulations using the parcel assumption tend to predict larger values of granular temperature compared to those using particles in systems with zero gravity (Figure 2), while the opposite behavior (i.e., lower values of granular temperature) is observed in systems under a gravitational force (Figure 4). A reduction in inelasticity would result in the opposite behavior in either system, which will be


10601

Figure 14. Transient behavior of spatially averaged (along the vertical y-direction) solids volume fraction in the region next to the wall (x ) 1.3 cm) using the MP-PIC and DPM techniques for the 2D periodic channel flow. Table 10. Average and Maximum Errors (Fractional) for MP-PIC Approach Relative to DPM Approach for the 2D Periodic Channel Flow quantity

average error

maximum error

volume fraction parcel velocity gas velocity granular energy

0.07 1.25 0.58 3.90

0.20 3.25 0.99 6.70

demonstrated in this section. Recall that granular temperature is a byproduct of interparticle (interparcel) collisions, as well

as from the random velocity of the particles (parcels) themselves.40 In the parcel approach, collisions between parcels are described with the same equations as collisions between particles, while collisions between particles within the parcel are ignored. Thus in this section, the equations to describe the collision forces acting on a computational parcel are revisited and modified to account for the particles within a parcel. Simulations of the two systems originally examined in part 1 (granular shear flow and gas-particle flow in a periodic riser) are repeated using the modified model. Finally, the results and numerical errors associated with the parcel concept using the modified model are presented. Derivation of Modified Model. The normal collision force between two parcels in contact, eq 9, is rewritten as follows b Fn ) knδn - ηnVn

(13)

Equation 13 clearly does not consider any collisions that may occur between the particles within a parcel (np) during a parcel-parcel collision. One simple remedy to this issue is to multiply both sides of eq 13 by np and then define a new collision force b Fn′ as b Fn′ ) npb Fn ) npknδn - npηnVn ) kn′δn - ηn′Vn

(14)

Essentially, when two parcels collide, all particles within one parcel are considered to collide with all particles of the other parcel. If the number of collisions considered increases, then the dissipation of granular energy should also increase. One of the primary mechanisms for dissipation of granular energy is through the restitution coefficient. To account for the extra dissipation due to the increased number of collisions, a modification to the restitution coefficient, evaluated using eq 10, is proposed based on the new collision constants (see eq 14): ln(e') ) -

Figure 15. Effect of reducing the restitution coefficient (e) on solids concentration and axial velocity profiles showing denser concentration and lower velocity due to formation of denser clusters of particles.

π(ηn′ /mp)

√2(kn′ /mp) - (ηn′ /mp)

)2

√npπ(ηn /mp)

√2(kn/mp) - np(ηn/mp)2

(15)

The new restitution coefficient (e′) is related to the original restitution coefficient (e) as follows (divide eq 15 by eq 10)

10602


ln(e') ) √np ln(e)

√2(kn/mp) - (ηn/mp)2 √2(kn/mp) - np(ηn/mp)2

(16)

After some manipulation, a closed form solution to eq 16 for e′ in terms of e and np can be obtained, which is expressed here as ln(e') ) √np ln(e)

√1 - (ln(e))2 /[(ln(e))2 + π2]

√1 - np(ln(e))2/[(ln(e))2 + π2]

(17)

For a limiting case where the spring constant is much larger than the damping coefficient, eq 16 can be reduced as np ln(e') ≈ n , for cases where kn /mp . (ηn /mp)2 ln(e) √ p 2

(18)

As is evident, eq 17 has a singularity for np g [(ln(e))2 + π2]/ (ln(e))2. One example to avoid this issue is where np should be less than 3752 for e ) 0.95. Equation 18 also shows that essentially e′ f 0 as np f ∞. This means that for finite values of e, the modification is limited to parcels with relatively small number of particles. For the values of e and np used in this study (see Tables 1 and 4), eqs 17 and 18 yield essentially the same modified restitution coefficients. Simulation Results with the Modified Restitution Coefficient. The restitution coefficient used in the simulations of the granular shear flow system described in Table 1 is e ) 0.9. For a case with 10 particles per parcel, the corresponding modified restitution coefficient is e′ ) 0.72, which can be obtained using either eq 17 or 18. With the modified value of restitution coefficient, the simulation with the 10 particles per parcel case is repeated. The averaged lateral profiles in solids volume fraction (εs), y-component of velocity (U), and granular temperature (T) profiles are shown in Figure 16 along with results from the original 10-particle-per-parcel and 1-particle-per-parcel cases. See part 1 for details on the averaging technique. All profiles demonstrate better agreement with the single particle case using the modified (lower) restitution coefficient compared to the original restitution coefficient. The solids volume fraction profile displays a small improvement, whereas fairly good agreement with the single particle case is found for the particles velocity profile. The latter is greatly influenced by the amount of dissipation due to inelastic collisions in the shear cell. As expected, the magnitude of granular temperature is lower in the modified (lower) restitution coefficient case, and the profile is in much better agreement with the results obtained in the single particle per parcel. Table 11 gives the relative average and maximum fractional error for the 10-particles-per-parcel case with the two different values of the restitution coefficient relative to the 1-particleper-parcel case. With the modified restitution coefficient, the results are remarkably closer to the single particle case than before. The parcel velocity and granular temperature are now more accurate than those obtained with a parcel of only two particles but with the original restitution coefficient of e ) 0.9 (see Table 2 for comparison). These improvements suggest that the proposed modification, in which all particles effectively collide during a parcel-parcel impact, may be a useful tool for enhancing the accuracy of the parcel approach. For the pseudo-3D periodic riser flow system described in Table 4, with e ) 0.95, the modified restitution coefficient for a case with 100 particles per parcel is e′ ) 0.85, which may be

Figure 16. Effect of lowering the restitution coefficient on the time-averaged flow variables computed across the width of the pseudo-3D granular shear cell. The scaling factors can be found in Figure 2. Table 11. Effect of Restitution Coefficient on the Average and Maximum Fractional Errors in the Granular Shear Flow System for the 10-Particles-Per-Parcel Case Relative to 1-Particle-Per-Parcel Case restitution coefficient

0.72

0.9

error

average

maximum

average

maximum


0.05 0.50 0.13

0.15 0.63 0.39

0.065 11.0 615

0.12 16.0 740

obtained using either eq 17 or 18. Note in that effort, the smallest parcel size was 10 particles per parcel, which is what the results of the 100 particles per parcel case are compared against. Accordingly, the value of np used in eq 17 or 18 is 10 (np ) 100/10 ) 10). Similar to the granular shear results discussed above, the numerical results obtained in this periodic riser flow using the modified restitution coefficient and coarse parcel are


10603

Table 12. Effect of Restitution Coefficient on Average and Maximum Errors (Fractional) for Periodic Riser Flow for 100-Particles-per-Parcel Case Relative to 10-Particles-per-Parcel Case restitution coefficient

0.85

0.95

error

average

maximum

average

maximum


0.078 0.83 0.25

0.18 5.59 0.39

0.113 1.75 0.41

0.25 11.0 0.56

restitution coefficient was lowered, the magnitude of granular temperature actually increased in this system. This behavior is attributed to the formation of denser clusters, which are then pulled downward due to the gravitational force; in turn this creates larger production due to shear, as observed in the resulting parcel velocity profiles. This behavior is not observed in the granular shear flow study due to the absence of a gravitational force in that system. Table 12 reflects the significant improvements in the simulation results due to a lower restitution coefficient with parcel velocity showing the greatest benefit. Although better improvements were obtained in the granular shear flow system, this modification is still worth considering for gas-particle riser flows since only a very simple adjustment to the restitution coefficient is required. Summary and Conclusions

Figure 17. Effect of lowering the restitution coefficient on time-averaged flow variables computed along the width of the 3D periodic riser.

remarkably closer to those using the smallest parcel size (10 particles per parcel), as shown in Figure 17. Although the

In this study, two basic assumptions used commonly in discrete particle methods are investigated. The first is the coarsegrained assumption in which many identical particles are grouped together into a parcel, or cloud. The second assumption is referred to herein as the MP-PIC technique in which the discrete collision forces are replaced with a solids pressure gradient. Both techniques have been employed in the literature with the main purpose of reducing the computational time associated with tracking a large number of individual particles and evaluating their corresponding forces. Without such techniques, discrete particle model simulations of large-scale reactors would be impractical to run. The current study also observes a significant computer time saving by making these basic assumptions. However, this saving is accompanied by large numerical errors that may not be acceptable to many engineers and scientists. Two different particulate systems are used to test the coarsegrained assumption, including wall-bounded shear flow of particles in a vacuum (granular flow) and the flow of air and particles in a periodic vertical riser. The numerical results obtained using different parcel sizes are compared with the results using the smallest numbers of particles per parcels. The various cases using different parcel sizes show some qualitatively similar behavior for momentum and solids density (comparison of snapshots) but are quantitatively fairly different in regard to the reported time-averaged volume fraction, velocity, and granular temperature profiles. The MP-PIC technique is evaluated using simulations of a bubbling fluidized bed and a periodic riser flow, similar to the riser system used for evaluating the coarse-grained assumption. In these simulations the parcel assumption is also employed. The time-averaged flow profile results obtained using the MPPIC model are assessed by comparing them to those obtained using the soft-sphere collision model. In the first case, the MPPIC model did not predict the same qualitative bubbling behavior observed in the DPM collision model; however, agreement between the two models could be improved by

10604


adjusting parameters in the MP-PIC solids pressure model (e.g., εsmax). In regard to the second case, both techniques predict the core-annulus flow regime reasonably well; however, large deviations are evident in the MP-PIC simulation. In particular, the MP-PIC simulation did not clearly capture the behavior predicted by the DPM simulation, namely, the intermittent formation of large clusters on one side of the channel that then move to the other side. Worth noting is that soft-sphere collision model depends on the restitution coefficient while the MP-PIC model is independent of this parameter. A better match between the two simulations may be obtained by adjusting the restitution coefficient in the former. While the restitution coefficient is a physically measurable property of the particle, it does not truly represent the physics involved in a parcel collision. Finally, a simple modification to the restitution coefficient is presented that is based on allowing all particles contained in a parcel to collide with other particles during a parcel-parcel contact. This simple approach results in significant reductions in the estimated errors for the parcel technique in both the granular shear flow system and periodic gas-particle flow in a riser. The authors hope this study will motivate more research in this area to further verify and validate the basic assumptions commonly used in discrete particle models. Nomenclature CD ) drag coefficient dpt, deq ) diameter of a particle and a parcel, respectively e ) particle-particle restitution coefficient b Fc,p ) net contact force acting on a parcel during collision b F n, b Ft ) normal and tangential component of the net contact force b g ) acceleration of gravity f If ) identity tensor kn, kt ) normal and tangential spring stiffness constants, respectively b n ) unit normal vector along the line of contact between two colliding parcels Nt, np ) total number of parcels and number of particles/parcel, respectively mp ) mass of a parcel Pg ) gas pressure P* ) empirical constant with units of pressure in the MP-PIC model Re ) Reynolds number ff S ) deviatoric part of the rate-of-strain tensor g

T ) granular temperature b Tp ) sum of all torques acting on a parcel b Vg, b Vp ) velocity vector of the gas-phase and a parcel, respectively b V n, b Vt ) contact velocity vector in the normal and tangential directions, respectively b xp ) position vector of a parcel Greek Letters β ) dimensionless empirical constant in MP-PIC model βp ) interphase momentum exchange coefficient between the gas and a single parcel δn ) overlap or displacement in the normal direction b δt ) overlap or displacement vector in the tangential direction εg, εs ) volume fraction of the gas and solids-phase, respectively φpt, φp, φc ) volume of a particle, a parcel, and a computational cell, respectively ηn, ηt ) normal and tangential dashpot damping coefficients, respectively µ ) Coulomb friction coefficient µg ) gas-phase shear viscosity Fg, Fp ) constant density of the gas-phase and a parcel (or particle), respectively

ff τ ) gas-phase stress tensor ω b p ) angular velocity of a parcel Indices g ) gas phase max ) maximum packing n ) normal to a surface old ) value at previous time-step used in the explicit MP-PIC formulation p ) particle or parcel s ) solids phase t ) tangential to a surface w ) wall

Literature Cited (1) Ding, J.; Gidaspow, D. A bubbling fluidization model using kinetic theory of granular flow. AIChE J. 1990, 36, 523. (2) Lun, C. K. K.; Savage, S. B.; Jeffrey, D. J.; Chepurniy, N. Kinetic theories for granular flow - inelastic particles in Couette-flow and slightly inelastic particles in a general flowfield. J. Fluid Mech. 1984, 140, 223– 256. (3) Cundall, P. A.; Strack, O. D. L. A discrete numerical model for granular assemblies. Geotechnique 1979, 29, 47–65. (4) Tsuji, Y.; Tanaka, T.; Ishida, T. Discrete particle simulation of twodimensional fluidized bed. Powder Technol. 1992, 71, 239. (5) Andrews, M. J.; O’Rourke, P. J. The multiphase particle-in-cell(MPPIC) method for dense particulate flows. Int. J. Multiphase Flow 1996, 22, 379–402, no. 2. (6) Snider, D. M.; O ’Rourke, P. J.; Andrews, M. J. Sediment flow in inclined vessels calculated using a multiphase particle-in-cell model for dense particle flows. Int. J. Multiphase Flow 1998, 24, 1359–1382. (7) Patankar, N. A.; Joseph, D. D. Modeling and numerical simulation of particulate flows by the Eulerian-Lagrangian approach. Int. J. Multiphase Flow 2001, 27, 1659–1684. (8) Dukowicz, J.K. A particle-fluid numerical model for liquid sprays. J. Comput. Phys. 1980, 35, 229–253. (9) Schmidt, D. P.; Rutland, C. J. A new droplet collision algorithm. J. Comput. Phys. 2000, 164, 62–80. (10) O’Rourke, P. J. Collective drop effects on vaporizing liquid sprays. Ph.D. Thesis, Princeton University, Princeton, NJ, 1981. (11) Snider, D. M.; O ’Rourke, P. J. ; Andrews, M. J. An incompressible two-dimensional multiphase particle-in-cell model for dense particle flows; LA-13280-MS, UC-700, Los Alamos Nat. Lab., June 1997. (12) O’Rourke, P. J.; Zho, P.; Snider, D. A model for collisional exchange in gas/liquid/solid fluidized beds. Chem. Eng. Sci. 2009, 64, 1784– 1797. (13) Snider, D.; Banerjee, S. Heterogeneous gas chemistry in the CPFD Eulerian-Lagrangian numerical scheme(ozone decomposition). Powder Technol. 2010, 199, 100–106. (14) Leboreiro, J.; Joseph, G. G.; Hrenya, C. M.; Snider, D. M.; Banerjee, S. S.; Galvin, J.E. The influence of binary drag laws on simulations of species segregation in gas-fluidized beds. Powder Technol. 2008, 184, 275– 290. (15) Patankar, N. A.; Joseph, D. D. Modeling and numerical simulation of particulate flows by the Eulerian-Lagrangian approach. Int. J. Multiphase Flow 2001, 27, 1685–1706. (16) Snider, D. M. An incompressible three-dimensional multiphase particle-in-cell model for dense particle flows. J. Comput. Phys. 2001, 170, 523–549. (17) Snider, D. Three fundamental granular flow experiments and CPFD predictions. Powder Technol. 2007, 176, 34–46. (18) Sakai, M; Koshizuka, S. Large-scale discrete element modeling in pneumatic conveying. Chem. Eng. Sci. 2009, 64, 533–539. (19) Harlow, F. H. The particle-in-cell computing method for fluid dynamics. In Fundamental Methods in Hydrodynamics, 3; Alder, B., Fernbach, S., Rotenberg, M., Eds.; Academic Press: New York, 1964. (20) Brackbill, J. U.; Ruppel, H. M. FLIP: a method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimension. J. Comput. Phys. 1986, 65, 314–343. (21) Burgess, D.; Sulsky, D.; Brackbill, J. U. Mass matrix formulation of the FLIP particle-in-cell method. J. Comput. Phys. 1992, 103, 1–15. (22) O’Rourke, P. J.; Brackbill, J. U.; Larrouturou, B. On particle-grid interpolation and calculating chemistry in particle-in-cell methods. J. Comput. Phys. 1993, 109, 37–52.

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010 (23) Flato, G. M. A particle-in-cell sea-ice model. Atmosphere-Ocean 1993, 31 (3), 339–358. (24) Dawson, J. M. Particle simulation of plasmas. ReV. Mod. Phys. 1983, 55 (2), 403–447. (25) Bardenhagen, S. G.; Brackbill, J. U.; Sulsky, D. The material-point method for granular materials. Comput. Meth. Appl. Mech. Eng. 2000, 187, 529–541. (26) Gidaspow, D. Multiphase Flow and Fluidization: Continuum and Kinetic Theory Description; Academic Press: San Diego, CA, 1994. (27) Agrawal, K.; Loezos, P. N.; Syamlal, M.; Sundaresan, S. The role of meso-scale structures in rapid gas-solid flows. J. Fluid Mech. 2001, 445, 151. (28) O’Rourke, P. J.; Amsden, A. A. The tab method for numerical calculation of spray droplet breakup; SAE paper 872089, 1987; pp 1-10. (29) Garg, R.; Galvin, J.; Li, T.; Pannala S. National Energy Technology Laboratory. Documentation of open-source MFIX-DEM software for gassolids flows. https://mfix.netl.doe.gov/documentation/dem doc 2010.pdf. (30) Tsuji, Y.; Kawaguchi, T.; Tanka, T. Discrete particle simulation of two-dimensional fluidized bed. Powder Technol. 1993, 77, 79–87. (31) Schafer, J.; Dippel, S.; Wolf, D. E. Force schemes in simulations of granular materials. J. Phys. 1 France 1996, 6, 5–20, no. 1. (32) Karion, A.; Hunt, M. L. Wall stresses in granular Couette flows of mono-sized particles and binary mixtures. Powder Technol. 2000, 109, 145– 163. (33) Silbert, L. E.; Ertas, D.; Grest, S. G.; Halsey, T. C.; Levine, D.; Plimpton, S. J. Granular flow down an inclined plane: Bagnold scaling and rheology. Phys. ReV. E 2001, 64, 051302-1051302-14.

10605

(34) Di Maio, F. P.; Di Renzo, A. Analytical solution for the problem of frictional-elastic collisions of spherical particles using the linear model. Chem. Eng. Sci. 2004, 59, 3461–3475. (35) Ye, M.; van der Hoef, M. A.; Kuipers, J. A. M. From discrete particle model to a continuous model of Geldart A particles. Chem. Eng. Res. Des. 2005, 83, 833–843. (36) Lohse, D.; Bergmann, R.; Mikkelsen, R.; Zeilstra, C.; van der Meer, D.; Versluis, M.; van der Weele, K.; van der Hoef, M.; Kuipers, H. Impact on soft sand: void collapse of jet formation. Phys. ReV. Lett. 2004, 93, 19803-119803-4. (37) Benyahia, S. Validation study of two continuum granular frictional flow theories. Ind. Eng. Chem. Res. 2008, 47, 8926–8932. (38) Karri, S. B. R. ; Knowlton, T. M. Flow direction and size segregation of annulus solids in a riser. In Fluidization IX; Fan, L. S., Knowlton, T. M., Eds.; Engineering Foundation: New York, 1998; pp 189196. (39) Benyahia, S.; Syamlal, M.; O’Brien, T. J. Study of the ability of multiphase continuum models to predict core-annulus flow. AIChE J. 2007, 53, 2549. (40) Campbell, C. S. Rapid granular flows. Annu. ReV. Fluid Mech. 1990, 22, 57–90.

ReceiVed for reView March 17, 2010 ReVised manuscript receiVed April 13, 2010 Accepted April 14, 2010 IE100662Z

Estimation of Numerical Errors Related to Some ... - ACS Publications

Recommend Documents