Article pubs.acs.org/IECR
Direct Numerical Simulation of Fluid Flow and Mass Transfer in Dense Fluid−Particle Systems Niels G. Deen* and J. A. M. Kuipers Multiphase Reactors Group, Department of Chemical Engineering and Chemistry, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands ABSTRACT: In this paper, a novel simulation technique is presented to perform direct numerical simulation (DNS) of fluid flow and mass transfer in dense fluid−particle systems. The fluid−solid coupling is achieved via direct (i.e., implicit) incorporation of the boundary condition (with a second-order method) at the surface of the particles at the level of the discrete momentum and species conservation equations of the fluid. A fixed (Eulerian) grid is utilized to solve the Navier−Stokes equations for the entire computational domain. Dissipative particle−particle and/or particle−wall collisions are taken into account via a hard-sphere discrete particle (DP) approach, using a three-parameter particle−particle interaction model that accounts for normal and tangential restitution, as well as tangential friction. Following the verification of our method using wellknown empirical expressions for the Sherwood number, we apply our method to study fluid−particle mass transfer in dense multiparticle systems involving random arrays of stationary particles.
■
INTRODUCTION Dense gas−particle flows are frequently encountered in a variety of industrial processes. More specifically, these flows are encountered in gas-fluidized beds that are frequently applied in the chemical, petrochemical, metallurgical, environmental, and energy industries in large-scale operations involving physical (coating, drying, and granulation) and chemical (synthesis of fuels and base chemicals) transformations. Very often, (simultaneous) mass and heat transport between the fluid and the particles takes place. In the literature, numerous empirical correlations for fluid−particle mass transfer have been presented both for isolated particles and particles in (dense) swarms, as encountered in fluidized suspensions. Although these correlations are very useful for design purposes and allow engineers to make quick estimates of the average rates of mass and heat transfer, it is not possible to obtain information on the (spatial) distribution of these quantities in dense arrays. Moreover, for systems involving polydisperse and/or nonspherical particles, these correlations are not well-established. For the accurate prediction of dense gas−particle flows in engineering scale equipment, which, in practice, can only be achieved with continuum models, it is essential that both the fluid−particle interactions as well as the particle−particle interactions are accurately taken into account. This requirement has led to the adoption of a multiscale modeling approach (see Figure 1) for these flows.1−3 At the most fundamental level, direct numerical simulation (DNS) can be used, which offers the possibility to directly compute the fluid−particle interaction and the associated drag closure laws, which are required at the more-coarse-grained levels of modeling. Because of computer processing unit (CPU) and memory constraints, typically O(103) particles can be treated simultaneously in this type of simulation. At the intermediate level, the discrete element model (DEM) model is used,4−8 in which individual particles are tracked in the computational domain, taking into account two-way coupling and dissipative collisions between the particles and/or confining walls. In the DEM model, the flow © 2013 American Chemical Society
Figure 1. Multiscale approaches for dense fluid−particle flows: direct numerical simulation (DNS), the discrete element model (DEM), and the two-fluid model (TFM) based on the Kinetic Theory of Granular Flow (KTGF).
field computation is based on the volume-averaged Navier− Stokes equations, where the control volumes typically contain a large number of particles and, as a consequence, the details of the flow in the vicinity of the particles is lost, necessitating the incorporation of a fluid−particle drag closure. The DEM model offers the advantage that rather-complicated particle−particle interaction models, including normal and tangential restitution and tangential friction, can be incorporated to assess their impact on the flow structure. The disadvantage of this type of Special Issue: Multiscale Structures and Systems in Process Engineering Received: Revised: Accepted: Published: 11266
December 11, 2012 January 14, 2013 January 17, 2013 January 17, 2013 dx.doi.org/10.1021/ie303411k | Ind. Eng. Chem. Res. 2013, 52, 11266−11274
Industrial & Engineering Chemistry Research
Article
model is given by the fact that typically only O(106) particles can be treated simultaneously, which is mainly due to CPU constraints. Finally, for the simulation of engineering-scale gasfluidized beds, a two-fluid model (TFM) or continuum model,9−11 based on the Kinetic Theory of Granular Flow (KTGF), is used. In this type of model, both the continuous phase and the particulate phase are described as continuous media with mutual interaction. In the continuum models for dense particulate flows, usually only the phase coupling due to drag is considered, whereas for the dissipative particle−particle interactions, only normal restitution is considered. With the advent of increasing computational resources, DNS has become a powerful tool for quantitatively deriving microscale transport coefficients, which are required to gain fundamental insight in gas−particle interactions, as illustrated with the above examples. Many authors applied immersed boundary (IB) methods to model the transport around static or moving particles. The IB method12−15 uses a fixed Eulerian grid to solve the conservation equations of mass, momentum, and thermal energy of the continuous phase. Lagrangian markers are placed on the surface of the immersed body and are advected along with the particle. The boundary conditions on the surface of the particle are enforced, by introducing appropriate source terms in the transport equations. In this way, no-slip conditions, as well as fixed temperature or flux conditions, can be applied. After the restoring force is calculated at the positions of each of the markers, the Eulerian force density is computed via distribution to the Eulerian grid. The IB method has been widely used to study fluid−structure interactions; it was pioneered by Peskin12 and applied to cardiac flow problems. Subsequently, the range of applications of this powerful method has expanded considerably. For excellent reviews, the interested reader is referred to the works by Peskin14 and Mittal and Iaccarino.15 One advantage of the IB method is its flexibility, with respect to the degree of rigidity (from elastic to rigid) of the bodies. Moreover, this method is relatively easy to implement. Disadvantages include the explicit treatment of the fluid−solid interaction, which leads to stiffness problems for rigid particles. In addition, appropriate values for the fluid−solid interaction parameters (such as the spring stiffness) must be determined for each particular class of problems. A different class of IB technique is the direct forcing method (DFM), first introduced by Mohd Yusof.16 Contrary to the vast literature on immersed boundary techniques to enforce the boundary conditions for the momentum equations, considerably fewer studies have focused on the application of these techniques to enforce the boundary conditions for scalar transport equations. In the DFM, the force densities are computed directly from the discretized Navier−Stokes equations, taking into account the no-slip conditions at the particle surface, at Lagrangian points (or surface elements) that are distributed over the particle surface. Many different versions of the DFM have been presented since its introduction. For example, Kim et al.17 introduced a mass source/sink in the continuity equation that leads to a decreased magnitude of the error in the solution. However, their spatial discretization error remained second-order. An elegant method that combines the strong points of the method of Peskin12 with the advantages of the DFM for moving rigid particles has been reported by Uhlmann.18 The current work builds on an alternative IB method, which we reported recently.19 This involves a fully implicit method that is second-order accurate. We have
demonstrated the merits of the model to predict gas−solid particle heat-transfer coefficients for static and moving (i.e., fluidized) particle arrays. Most DNS studies of particulate systems have been limited to momentum transfer in static particle arrays.2,20,21 Only recently, these studies have been extended to also include heat transfer.19,22−27 Dixon et al.,28 in their study of fixed-bed heat transport, included particle−particle heat transport through contact points. Little work has been reported so far on combined momentum, heat, and mass transfer, also including chemical reaction.29,30 Yang et al.31,32 studied convective heat transfer in structured beds of nonspherical particles. The organization of this paper is as follows. First, the description of the model is given. Subsequently, the numerical solution method is discussed and the results are presented. Finally, the conclusions are presented.
■
MODEL DESCRIPTION Our model consists of two main parts: one part involves the solution of the fluid phase equations, accounting for the presence of the solid particles, whereas the other part addresses the solution of the solid-phase equations, accounting for all forces acting on them and the possible nonideal collisions between the particles themselves and/or confining walls. First, the main assumptions are given: (a) The fluid phase is incompressible and Newtonian. (b) All physical properties of both phases are constant. (c) The solid phase consists of spherical catalyst particles. (d) At the surface of the particles, an infinitely fast chemical reaction occurs. Fluid-Phase Equations. The transport phenomena in the fluid phase are governed by the conservation equations for mass, momentum, and chemical species, which are respectively given by (∇· u ̅ ) = 0
∂ρf u ̅
+ (∇·ρf uu) = −∇p + μf ∇2 u ̅ + ρf g ̅
(2)
+ (∇·cA,f u ̅ ) = DA,f ∇2 cA,f
(3)
∂t
∂cA,f ∂t
(1)
In these equations ρf, μf, and DA respectively represent the density, viscosity, and diffusivity of the fluid. Solid-Phase Equations. The translational and rotational motion of the suspended solid particles is governed by the Newtonian equations of motion, which are respectively given by
mp
Ip
dwp̅ dt
dωp̅ dt
= mpg ̅ + Ff̅ →s
= Tf̅ →s
(4)
(5)
where mp and Ip respectively represent the mass and the moment of inertia of the particle. The final terms on the righthand sides in eqs 4 and 5 account for the fluid−particle interaction (respectively, drag and torque) and are given by Ff̅ →s = −
∬S (τf ·n̅ + pn̅ ) dS p
11267
(6)
dx.doi.org/10.1021/ie303411k | Ind. Eng. Chem. Res. 2013, 52, 11266−11274
Industrial & Engineering Chemistry Research Tf̅ →s = −
Article
Δt u ̅ n + 1 = u ̅ ** − ∇p n + 1 αρf
∬S ( r ̅ − rp̅ ) × (τf ·n̅ + pn̅ ) dS p
=−
∬S ( r ̅ − rp̅ ) × (τf ·n̅ ) dS p
Since u̅n+1 must satisfy the incompressibility constraint, upon taking the divergence of eq 14, the pressure Poisson equation is obtained:
(7)
where τf represents the viscous stress tensor, which is given by τf = −μf [(∇u ̅ ) + (∇u ̅ )T ]
α 1 2 n+1 ∇p = (∇· u ̅ **) ρf Δt
(8)
It should be noted that the buoyancy term is not included in an explicit manner in the equation of motion for the particles, since it is included in the calculation of the pressure force acting on the particle surface. The mass transfer to the surface of the particles is calculated using the following equation: cA,f |particle−surface = 0
∬S (Df ∇cf ·n̅ ) dS p
(9)
αcA,f n + 1 = βcA,f n + γcA,f n − 1 − ΔtCs̅ n + Δt[DA,f ∇2 cA,f n + 1]
(10)
(16)
To evaluate the total fluid−particle force (F̅f→s), the torque (T̅ f→s) and mass-transfer rate (Φf→s), as well as the velocity gradients and the concentration gradients, must be known. The evaluation of these quantities will be explained in detail later.
where the net convective species molar flux (Cs) is given by Cs = (∇· uc̅ A,f )
■
α[ρf u ̅ n + 1] = β[ρf u ̅ n] + γ[ρf u ̅ n − 1] − Δt ∇pn + 1 (11)
where n indicates the time index and C̅ m is the net convective momentum flux, which is given by Cm̅ = ρf (∇·uu)
Δt wp̅ n + 1 = wp̅ n + g ̅ Δt + [Ff̅ →s n + Ff̅ →s n + 1] 2mp
(12)
For the discretization of the time derivative, we use a secondorder three-point backward discretization scheme (with α = 1.5, β = −2, and γ = 0.5), requiring three levels of storage for the velocity components. The solution of eq 11 is achieved via a two-step projection method, where, as a first step, the tentative velocity field u** ̅ is computed from
ωp̅ n + 1 = ωp̅ n +
Δt [Tf̅ →s n + Tf̅ →s n + 1] 2Ip
(18)
(19)
This implies that the total fluid−particle interaction force and the torque must be evaluated for each particle at the beginning and the end of each time step. Once the new particle translational and rotational velocities are obtained, an eventdriven hard-sphere collision model is invoked. In this model, it is assumed that the interaction forces are impulsive, and, therefore, all other finite forces are negligible during collision. The closure of this collision model involves three micromechanical parameters: the coefficients for normal and tangential restitution and the tangential friction coefficient, which in principle can be obtained from impact experiments. Fluid−Solid Coupling. The fluid−solid coupling constitutes the key element of our method. The discretization of the species and momentum conservation equations lead to algebraic equations of the following generic form:
α[ρf u ̅ **] = β[ρf u ̅ n] + γ[ρf u ̅ n − 1] − ΔtCm̅ n + Δt[μf ∇2 u ̅ **] + ρf g ̅ Δt
(17)
For the spatial discretization of the convection and diffusion terms in eq 16, the schemes used are exactly the same as those used for the solution of the momentum equations. Solid-Phase Equations. The source terms appearing in the Newtonian equations of motion are treated as known (explicit) terms and, therefore, the integration of these equations can be conducted in principle with any integration technique for ordinary differential equations (ODEs). For the simulations reported in this paper, we have used a second-order trapezoidal rule, producing translational and rotational velocities and at the new time level computed, respectively, as follows:
NUMERICAL SOLUTION METHOD Fluid-Phase Equations. The fluid-phase equations are solved in three dimensions (3D) on a Cartesian grid with uniform grid spacings in all directions. First, we perform time discretization of the momentum equation to obtain
− ΔtCm̅ n + Δt[μf ∇2 u ̅ n + 1] + ρf g ̅ Δt
(15)
which again is solved with a robust and efficient ICCG algorithm to obtain the pressure at the new time level. Finally, the velocity field at the new time level is obtained from the second (correction) step (eq 14). This correction step has a negligible effect on the correct enforcement of the no-slip condition at the fluid−particle boundaries. It should be stressed here that the pressure and velocity are obtained for the entire domain (i.e., also for the grid points interior to the particles). Subsequently, the species conservation equation is solved using the same temporal discretization as that used for the solution of the momentum equations:
implying a completely mass-transfer-limited chemical reaction at the surface of the particles. The mass-transfer rate from the fluid to the particle is given by Φf→s = −
(14)
(13)
The Laplace operator is approximated with standard central second-order finite difference representations, whereas the convection terms are computed with a second-order fluxdelimited Barton scheme.33 The enforcement of the no-slip condition at the surface of the immersed boundary is handled at the level of the discretized momentum equations and will be explained later. We use a robust and very efficient incomplete Cholesky conjugate gradient (ICCG) algorithm to solve the resulting sparse matrix equation for each velocity component. The velocity field at the new time level n + 1 is related to the tentative velocity field as follows:
acϕc =
∑ anbϕnb + bc nb
11268
(20)
dx.doi.org/10.1021/ie303411k | Ind. Eng. Chem. Res. 2013, 52, 11266−11274
Industrial & Engineering Chemistry Research
Article
where ϕ corresponds to the species concentration or one of the velocity components in the fluid phase. This equation is obtained for all “fluid” nodes denoted as “c” exterior to the particles (see Figure 2) and relates quantity ϕc to six (3D) or
Figure 3. Incorporation of the boundary condition for a general fluid quantity ϕ. The five indicated nodes are used for the discrete representation (in 2D) of the PDE governing fluid quantity ϕ. For the nodes residing inside the (moving) particles a and b, the connectivity of these nodes to particle a or particle b must be established during each time step. Because of the utilization of a staggered grid, this connectivity must be determined for each velocity component and species concentration separately. For the situation shown, the boundary condition (see Figure 2) must be applied for three neighbor nodes that reside inside the particle a (1 node) or particle b (2 nodes).
Figure 2. Incorporation of the boundary condition for a general fluid quantity ϕ. The five indicated cells are used for the discrete representation (in 2D) of the PDE governing fluid quantity ϕ (fluid velocity component or fluid temperature). The “solid” node corresponding with ϕ0 resides inside the particle, whereas the “fluid” nodes corresponding to ϕ1 and ϕ2 reside in the fluid. The value of ϕ0 is expressed in terms of a second-order polynomial, according to ϕ = aξ2 + bξ + c, where the values of the coefficients a, b, and c are obtained from the known values of ϕ at the three nodes (ξ = 0: ϕ = ϕ0; ξ = 1: ϕ = ϕ1; and ξ = 2: ϕ = ϕ2), leading to eq 21. Note that ξ = 0 corresponds to the location of the “solids” node.
the local force exerted by the fluid on the particle and the local fluid-to-particle mass flux over the external surface of each particle and require the evaluation of, respectively, the velocity gradient and concentration gradient at the particle surface. For a general quantity ϕ, the (local) dimensionless gradient at the particle surface is given by the following expression: dϕ dξ
four (2D) neighboring nodes, indicated with “nb”. For moving particles, the detection of these nodes must be carried out during each time step, taking into account the location of the node in the staggered computational grid (velocity nodes are located at the faces of the cells, the node for scalar quantities, such as the species concentration is located at the center of the cells). In that case, a boundary condition must be applied where the value of ϕnb = ϕ0 at that particular “solid” node is expressed as a one-dimensional (1D) linear combination of the ϕ values of the relevant “fluid” nodes ϕ1 and ϕ2: ⎛ 2ξs ⎞ ⎛ ξ ⎞ ϕ0 = −⎜ ⎟ϕ1 + ⎜ s ⎟ϕ2 ⎝ 1 − ξs ⎠ ⎝ 2 − ξs ⎠ ⎤ ⎡ 2 +⎢ ⎥ϕ ⎣ (1 − ξs)(2 − ξs) ⎦ p
ξ = ξs
⎛ 2 − ξs ⎞ ⎛ 1 − ξs ⎞ =⎜ ⎟ϕ1 − ⎜ ⎟ϕ ⎝ 1 − ξs ⎠ ⎝ 2 − ξs ⎠ 2 ⎡ ⎤ 3 − 2ξs −⎢ ⎥ϕ ⎣ (1 − ξs)(2 − ξs) ⎦ p
(22)
For the evaluation of F̅f→s, the contribution due to the pressure also must be evaluated. Note that the viscous and pressure contribution of the drag are obtained separately, unlike the Immersed Boundary Method (IBM) proposed by Uhlmann.18 We have implemented two different methods that work equally well. In the first method, the pressure at the surface of the particle is obtained by (linear) extrapolation of the pressure at the fluid nodes in the vicinity of the particle surface. In the second method, the surface integral of the pressure is converted to a volume integral (see eq 23):
(21)
Ff̅ →s = −
where ϕp represents the desired value of ϕ at the boundary of the particle and ξs is a dimensionless distance that is known from the intersection of the surface of the particle and the grid line. Expression 21 is used to eliminate ϕ0 from eq 20, which leads to an algebraic equation similar to eq 20 but with modified constants. This procedure must be carried out for all “solid” nodes to properly account for the local boundary condition to be met at the surface of the particle. It is possible that, for fluid node “c”, boundary conditions must be applied for multiple neighboring nodes (temporarily) coinciding with different (moving) particles; therefore, it is essential to keep track of the “connectivity” of nodes to particles (see Figure 3). For the representation of the fluid−particle interaction, the force F̅f→s, the torque T̅ f→s, and the mass-transfer rate Φf→s must be evaluated. These quantities are obtained by integrating
∬S pn̅ dS = −∭V ∇p dV p
p
(23)
The volume integral is evaluated by summing the values of the pressure gradient calculated at all nodes interior to the particle under consideration. The components of the pressure gradient are calculated at the velocity nodes in the staggered grid using standard (second-order) central finite-difference approximations. In the limiting case of no motion (a fluid− particle system at rest), the pressure gradient in eq 23 follows, from hydrostatics ∇p = ρfg ̅ (see eq 2) and the total fluid− particle ineraction force is given by the well-known buoyancy force: Ff̅ →s = −
∭V ∇p dV = −∭V ρf g ̅ dV = −ρf Vpg ̅ p
11269
p
(24)
dx.doi.org/10.1021/ie303411k | Ind. Eng. Chem. Res. 2013, 52, 11266−11274
Industrial & Engineering Chemistry Research
■
Article
RESULTS Convective Mass Transfer to a Stationary Sphere. In our first test, we consider convective mass transfer to a single stationary sphere in an enclosure. The data used for the numerical simulations are given in Table 1. Laterally, the sphere
Nevertheless, the implementation of AMR is a major task in itself and clearly goes beyond the scope of this paper. Convective Mass Transfer to Stationary Arrays of Spheres. In this section, we will consider the flow through a stationary random array of particles (see Table 3 for the data
Table 1. Data Used for the Simulations of Convective Mass Transfer to a Stationary Sphere
Table 3. Data Used for the Simulations of the Stationary Random Particle Configuration
parameter
parameter
value 0.16 m × 0.16 m × 0.16 m 1.0 × 10−4−2.5 × 10−4 s 0.020 m 1.00 kg/m3 2 × 10−5 kg/(m s) 2 × 10−5 m2/s
domain size time step sphere diameter fluid density fluid viscosity fluid diffusivity
was placed at the center of the domain, whereas in the flow direction, the center of the sphere was positioned at a distance 2dp from the inlet. For these simulations, a 1603 grid was used with uniform grid spacing in all directions. For a single sphere, the expression for the particle Sherwood number (Shp) is given by the well-known empirical Froessling equation: Shp =
k md p DA,f
= 2.0 + 0.6(Rep)1/2 (Sc)1/3
a
ρf u0d p μf
and
Sc =
μf ρf DA,f
(26)
The comparison between the computed Sherwood numbers and the results obtained from the Froesling equation is quite good (see Table 2), despite the fact that the ratio of the particle Table 2. Computed Mass-Transfer Coefficients for a Stationary Sphere for Several Particle Reynolds Numbers and Comparison with the Froessling Equation Mass-Transfer Coefficient, km (m/s) Rep
simulated
Froessling equation
relative error (%)
10 20 40 60 100 200 300 400
0.01520 0.01852 0.02319 0.02676 0.03245 0.04338 0.05180 0.05850
0.01559 0.01873 0.02318 0.02659 0.03200 0.04190 0.04960 0.05600
−2.50 −1.12 +0.04 +0.64 +1.41 +3.53 +4.44 +4.46
150 5× 1× 3×
× 150 × 600 = 13.5 × 106 10−4 m 10−4 s 10−3 m
0.0 < x-position < 0.075 m 0.0 < y-position < 0.075 m 0.0 < z-position < 0.200 m 1.0 kg/m3 2.0 × 10−5 kg/(m s) 2.0 × 10−5 m2/s 60−400 m/s
Random configuration (3000 spheres).
used in the simulation). The spheres were distributed in a random fashion over the computational domain (excluding the top section with z > 0.20 m), using a standard Monte Carlo procedure to produce a predefined voidage of 0.70. The incorporation of the empty top section was found to be essential to avoid recirculation over the outflow boundary at higher particle Reynolds numbers (Rep > 180). For the simulation, a prescribed uniform velocity and a prescribed pressure were imposed at the bottom and top boundaries, respectively. For the other boundaries, the no-slip condition was imposed. At the bottom wall, the fluid entered with a uniform species concentration of 10.0 mol/m3. At all remaining walls, a zero-concentration gradient (Neumann boundary condition) was imposed. The simulations were performed on a single-node standard desktop PC and required ∼4 weeks to be completed. The particle Reynolds number changed from 60 to 420, corresponding to a superficial velocity ranging from 0.2 m/s to 1.4 m/s. In Figure 4, the particle configuration is shown (3000 spheres in a domain with a vertical coordinate of
