Direct Numerical Simulation of Fluid Flow and Mass Transfer in Particle Clusters.

In this paper, an efficient ghost-cell based immersed boundary method is applied to perform direct numerical simulation (DNS) of mass transfer problems in particle clusters. To be specific, a nine-sphere cuboid cluster and a random-generated spherical cluster consisting of 100 spheres are studied. In both cases, the cluster is composed of active catalysts and inert particles, and the mutual influence of particles on their mass transfer performance is studied. To simulate active catalysts the Dirichlet boundary condition is imposed at the external surface of spheres, while the zero-flux Neumann boundary condition is applied for inert particles. Through our studies, clustering is found to have negative influence on the mass transfer performance, which can be then improved by dilution with inert particles and higher Reynolds numbers. The distribution of active/inert particles may lead to large variations of the cluster mass transfer performance, and individual particle deep inside the cluster may possess a high Sherwood number.

Introduction

Circulating fluidized bed reactors (CFB) are frequently applied in a wide range of processes in the chemical, energy and coating industries. Understanding mass transport processes as well as the fluid flow in such complex heterogeneous systems is of tremendous importance to improve performance and facilitate optimal design of process equipment. However, the occurrence of particle clusters is observed[1−4] and well accepted to have strong impact on the hydrodynamics[5−9] and consequently the mass transport behavior.[10−14] The occurrence of particle clusters was first reported by Yerushalmi et al.,[15] and later Tsuo and Gidaspow[16] studied the underlying principles for the formation and propagation of clusters. Although there is no firmly established definition for a cluster,[1,5,9,17−19] one feature is universally accepted that a cluster is composed of particles and possesses an internal solid volume fraction significantly larger than its surroundings (which is normally below 0.1 in CFB). Clusters can be characterized by duration of their existence, occurrence frequency and internal solids concentration, as studied by Sharma et al.[18] and other researchers.[19−22] In the past decades, numerous empirical correlations for fluid-particle mass transfer have been proposed for multisphere systems.[23−29] Although these correlations are helpful for a quick and rough estimation of overall mass transfer behavior for design purposes, it does not consider the influence of particle heterogeneities. Significant deviations of the Sherwood number in risers have been reported, and Breault[30] and Chalermsinsuwan et al.[31] pointed out these values even differ by several orders of magnitude. The reason behind this finding is explained by the presence of particle clusters. Such dense solid phases result from the particle–particle collisions as well as the particle–wall collisions, which are characterized by enhanced fluid bypass, low slip velocities and fluid back-mixing for the hydrodynamic phenomena. These reduced fluid performance in CFB will lead to poor fluid–solid contacting and negatively impact the interfacial mass transfer processes.

In clusters, the mass transfer to an individual “active” sphere is significantly influenced not only by the surrounding gas–solids suspensions but also by the nature of the surrounding particles, namely it is active or inert.[32] In this circumstance, dispersed active particles keep their intrinsic kinetics whereas two effects from the inert surrounding particles must be taken into consideration.[33,34] First, the inert particles decrease the volume available for mass transfer around the active particle. This makes the analog between heat and mass transfer invalid, as the particles provide an additional conductive pathway to the heat transfer process. Second, the flow field and the formation of the mass boundary layer around the active particle are altered due to the presence of the other particles. Mixtures of active particles diluted by inert ones are used to keep the system at a low conversion rate (for positive reaction orders)[35] and minimize complications caused by bed heterogeneities.[36] In catalytic reaction engineering, catalyst beds with dilution are often applied to study the kinetics isothermally in highly exothermic heterogeneous reactions[37] and various investigators have mentioned their considerable merits.[38−40] For mass transfer processes in diluted systems, researchers are interested in the influence of the dilution ratio (ratio of inert and total) and the influence of the distribution of inert and active particles.[32,34,37,41,42]

Because of the inherent difficulty to experimentally quantify the mass transfer behavior in a dense environment, several numerical studies have been performed for a cluster system to estimate its mass transfer coefficient. However, due to the high computational cost, simulations are restricted to large scale closure models. Dong et al.[43] applied the subgrid energy-minimization multiscale (EMMS/mass) model for mass transfer computation, in which the species concentration in each grid cell is divided into the mass fraction of the particle-rich dense phase and the mass fraction of the fluid-rich dilute phase. In this case, mass transfer performance varies at the subgrid level, and the authors reported that both gas velocity and solids flux have significant impact on the variation of the Sherwood number. Chalermsinsuwan et al.[31] also used the EMMS model to estimate the mass transfer coefficient in a PSRI riser. In their work, cluster sizes are computed first and then the Sherwood numbers are estimated. They found that the Sherwood number scales with the particle cluster diameter. Wang et al.[44] studied the mass and heat transfer of a spherical cluster by a modified k–ε turbulence model. Their cluster consisted of 65 stationary particles which are arranged regularly, and the corresponding porosity is 0.9. Based on their simulations, the influences of cluster porosity and Reynolds number are quantified. Kashyap and Gidaspow[45] analyzed the mass transfer coefficients in fluidized bed by using a kinetic theory based multiphase mode. The cluster diameter is suggested to replace the particle diameter in the conventional Sherwood number computation.

With the development of computational technology, DNS has become a powerful tool to resolve all the details at the smallest relevant length scales and quantitatively derive microscale transfer and transport coefficients to gain fundamental insight in fluid–solid interactions. In recent years, the immersed boundary method (IBM), as a branch of DNS, has received a lot of attention. IBM has distinct advantages such as efficient CPU/memory utilization and easy grid generation, and thus has been applied in various studies including complex geometries, moving particles and deformable immersed objects.[46−51] Next to the fluid flow equations, additional equations for species transport can be solved with relative ease using the same methodology.

There are two categories of IBM: continuous forcing method (CFM)[52−57] and discrete forcing method (DFM).[58−63] The fluid–solid interaction is accounted for by a forcing term in the fluid governing equations and ghost variable values inside the immersed body, respectively. Mohaghegh and Udaykumar[64] performed a comparison study between CFM and DFM for simulations of particulate flows with both stationary and moving boundaries. IBM has been widely used for studies of momentum transfer in fluid–solid system; however, scarce studies have been reported for mass transfer processes to clusters composed of active and inert particles. In this paper, a previously developed ghost-cell based immersed boundary method[65] is applied for the simulation of mass transfer problems in particle clusters. The swarm effect together with the influence of the inert particles and the Reynolds number are investigated.

The organization of this paper is as follows. First, the description of the model is given, including the governing equations, numerical solution method and fluid–solid coupling. After that, our DNS model is verified by two test cases for which analytical or empirical solutions exist. As the main results two stationary cluster systems, a nine-sphere cuboid cluster and a random-generated spherical cluster, are considered and analyzed. Finally, the conclusions are presented.

Model description

In this part, we describe the governing equations that need to be solved in DNS, the numerical details involved in the finite difference scheme, as well as the fluid–solid coupling. For the model presented in this paper, the following main assumptions are applied:

The fluid phase is incompressible and Newtonian.

The solid phase consists of stationary spherical particles. In case of active catalyst particles, an external mass transfer limited chemical reaction proceeds at the external surface.

Both fluid and solid phases have constant physical properties.

The mass transfer is modeled as Fickian diffusion.

Governing Equations

The transport phenomena in the fluid phase are governed by the conservation equations for mass, momentum and species transfer, respectively given bywhere ρ is the fluid density and μ is the fluid viscosity, whereas D represents the species mass diffusivity in the fluid.

Numerical Solution Method

The governing equations are solved on a 3D staggered Cartesian grid with uniform grid spacing in all three directions. Following the work of Deen et al.,[66] the numerical solution of the governing equations is acquired utilizing second order discretization schemes. At the same time, small computational stencils are preferred for computational efficiency. First, the momentum equation is discretized in time as given by the following expression:where n is the time step index. The convective and diffusive momentum fluxes C and D are calculated by spatial discretization of:The convection term is spatially discretized by a second-order total variation diminishing scheme, whereas a standard second-order central differencing scheme is used for the evaluation of the diffusion term.

The solution of Equation is achieved by using a fractional step method where a tentative velocity field u̅** is first computed by using the pressure gradient at the old time step p. As the second step, the velocity field at the new time step n+1 is obtained based on the new pressure gradient calculated from Poisson equation at time step n+1. A robust and efficient parallel Block-Incomplete Cholesky Conjugate Gradient (B-ICCG) solver is used to obtain u̅** and p. For more information, we refer a more detailed description of this method to the work of Deen et al.[66] and Das et al.[67]

The species conservation equation is temporally discretized in the same way as for the momentum equation, namely the Adams–Bashforth scheme is applied for the convective transport whereas fully implicit Euler backward scheme is used for the diffusion term.with the convective and diffusive fluxes C and D respectively given byThe spatial discretization of species transport equation is the same as for the momentum equation.

The boundary condition is enforced at the exact position of the immersed object surface, for both velocity and concentration computation. It is handled at the level of the discretized momentum and species conservation equations and will be introduced in detail in the next section.

Fluid–Solid Coupling

After discretization, the differential conservation equations can be written as algebraic equations in the following generic form:where ϕ is the fluid variable for which we want to obtain a solution, namely the velocity and concentration field obtained from the momentum and species equation, respectively. This equation provides the relationship between the fluid quantity ϕ at any position of the simulation domain and its six neighboring cells indicated as ϕ. Due to the application of a nonbody conforming mesh, the boundary condition enforcement and consequently the solution of this equation are no longer straightforward at the fluid–solid interface. Ghost points are used to overcome this problem, which are inside the solid phase but possess at least one neighbor in the fluid phase. Every fluid cell is checked and the predefined boundary condition is applied if any of its six surrounding neighbors represents a ghost point.

The fluid–solid coupling is the key element of our DNS model, which involves a second order quadratic interpolation scheme. Taking the advantage of the Robin boundary condition imposed exactly at the object surface in a sharp interface manner, active catalysts are simulated by applying the zero-value Dirichlet boundary condition whereas inert particles are realized by enforcing the zero-flux Neumann boundary condition. In other words, mixed boundary conditions of individual spheres are handled consistently in our DNS model. The full methodology of the enforcement of the Robin boundary condition is introduced in our previously published paper, which we refer to for more details for the interested reader.[65] In this paper, we only outline the essential equations.

In the quadratic interpolation scheme, a generic variable ϕ in the vicinity of the immersed object surface can be approximated in terms of a second-order polynomial:where x, y and z are relative coordinates with respect to the origin located at the boundary point. This equation is in fact the approximation of ϕ using the Taylor expansion near the boundary point:For the enforcement of the Robin boundary condition:the first four coefficients are required:Nine neighboring fluid points and one image point are used to obtain the ten coefficients C in a second-order polynomial, where the resulting equation can be written as follows:where ϕ and C are the vectors for species concentration and coefficients respectively, and X is the Vandermonde matrix given byThe coefficients C are obtained by multiplication of the inversed matrix X–1 and the concentration vector ϕ, which can be written as a linear combination of ϕ values. With those required coefficients the value at the image point can be evaluated by satisfying the boundary condition, from which the value at the ghost point can be finally computed.

With the above-described procedure, the matrix coefficients in Equation can be updated. Altered coefficients within the original stencil are incorporated in the implicit scheme, whereas neighbors outside the original stencil are accounted for in an explicit way. This procedure is carried out for all ghost points to ensure that the desired local boundary condition applies everywhere at the immersed object surface. The pressure, velocity and concentration field are obtained for the entire computation domain in our DNS model.

Verification

In this section, two test cases are simulated by using our proposed DNS model. One is the Graetz–Nusselt problem, i.e., mass transfer from the tube wall to the fluid phase. The other one is the convective flow to a single stationary active/inert sphere. The results obtained from computer simulations are compared with the analytical/empirical solutions, through which our DNS model is verified.

Mass Transfer in Laminar Tube Flow

In the first test, the well-known Graetz–Nusselt problem is considered, which deals with forced mass convection combined with mass diffusion of a fully developed laminar flow in a tube. The boundary conditions at the cylindrical tube wall are treated with the methodology outlined in the previous section. The classical Graetz–Nusselt problem was investigated by Graetz[68] and later independently by Nusselt[69] for the case of a constant species concentration at the wall of the tube. Later, this problem was extended to the constant mass flux boundary condition, and Tao[70] and Tyagi[71] provided a mathematical solution to that. The amount of mass transfer at the wall is defined as the ratio of convective to diffusive radial mass transfer, which is commonly expressed by the Sherwood number Sh:where k is the mass transfer coefficient, D is the diameter of the tube and D is the mass diffusivity in the fluid. The analytical solutions are 3.668 and 4.364 under fully developed conditions, for the boundary condition of constant wall concentration and constant wall flux, respectively.[72]

By applying our DNS model, both boundary conditions can be realized. The fluid enters the tube with zero concentration, and a uniform velocity of 1 m/s is specified to maintain a laminar flow. The data used for the simulations are summarized in Table . The mass transfer coefficient is calculated as the mass flux divided by the driving force for the mass transfer process. For the case of constant wall concentration, the mass flux is averaged over the perimeter of the tube whereas the driving force is the difference between the prescribed wall concentration and the cup-average concentration of the fluid. For the case of constant wall flux, the mass flux is a prescribed value and the driving force is computed as the difference between the averaged wall concentration and the fluid cup-average concentration. The local concentration gradient (required for local mass flux calculation) and the local wall concentration are easily obtained from the values at the image and ghost point, whereas the cup-average concentration is defined byIn this equation, the integration is performed over a surface S perpendicular to the flow direction x, and u(x,y,z) is the x component of the fluid velocity at this certain point (x,y,z). It should be noted that the Sherwood number obtained from simulation is computed at locations far beyond both the hydrodynamic entrance region and the mass transfer entrance region. In other words, the Sherwood number is calculated under the circumstance of a fully developed flow. Simulations results are listed in Table , using different grid size. The computational domain, e.g., 100 × 10 × 10, indicates the grid numbers in the x, y and z directions, respectively. From the table, good agreements with the analytical results are observed for both boundary conditions.

Table 1

Data Used for the Simulations of Convective Mass Transfer for Laminar Flow in a Tube

Parameter	Value	Unit
Time step	1 × 10–4–5 × 10–5	s
Tube length	0.2	m
Tube diameter	0.02	m
Fluid density	1000	kg/m3
Fluid viscosity	1	kg/m/s
Species diffusivity	0.001	m2/s
Initial concentration	0	mol/m3

Table 2

Sherwood Numbers Obtained from Simulations for Graetz-Nusselt Problem, at Several Mesh Resolutions

	Computational domain
	100 × 10 × 10	200 × 20 × 20	400 × 40 × 40	800 × 80 × 80
Constant wall concentration	3.719	3.678	3.669	3.668
Constant wall flux	4.464	4.393	4.374	4.365

Single Stationary Sphere under Forced Convection

In the second test, we consider convective flow to a single stationary sphere in an enclosure. The sphere is located at the center of the domain laterally while it is positioned at a distance of two times the sphere diameter from the inlet in the flow direction. We consider two cases for the sphere: active and inert. For the former case, an external mass transfer limited reaction proceeds at the sphere surface, while the fluid solely flows over the sphere without any chemical conversion for the latter case. The data used for the numerical simulations are given in Table . The domain boundary condition in lateral direction is set to be free-slip and zero flux for velocity and concentration field computation, respectively. Fluid containing a single species with constant concentration of 10 mol/m3 flows into the system at a uniform velocity. At the outlet, pressure is prescribed as the standard atmospheric pressure and the species boundary condition is specified to be zero slope. The simulations are performed on a 160 × 160 × 160 grid with uniform grid spacing in all directions. The ratio of domain size to the particle size is 8 while the mesh resolution is 20. The mesh resolution N is defined as the ratio of the sphere diameter to the grid size. N = 20 is selected after a mesh convergence test, in which we used the case of Re = 200 for active catalyst particle. With N = 10, 20, 32 and 40, particle Sherwood numbers of 9.04, 10.46, 10.48 and 10.49 were obtained respectively, revealing a good mesh convergence. It should be noted that the Schmidt number is unity, indicating the same thickness of the momentum and mass boundary layers.

Table 3

Data Used for the Simulations of Single Stationary Sphere under Forced Convection

Parameter	Value	Unit
Time step	2 × 10–5–5 × 10–5	s
Grid size	2.5 × 10–4	m
Sphere diameter	0.005	m
Fluid density	1	kg/m3
Fluid viscosity	2 × 10–5	kg/m/s
Species diffusivity	2 × 10–5	m2/s
Initial concentration	10	mol/m3
Inlet concentration	10	mol/m3

For the case of an active catalyst particle, the verification can be done by comparing the particle Sherwood number Sh obtained from the simulation with the empirical value given by the well-known Frössling correlation:where Re is the particle Reynolds number and Sc is the Schmidt number, respectively defined as

The Sherwood number obtained from simulation work is computed by the following expression:where Φ is the mass transfer rate, with the normal pointing outward of the solid, calculated by the integration of the concentration gradient at sphere surface over the whole sphere (S indicating the external surface area of the sphere):With the flow velocity u0 varying from 0.04 to 1.6 m/s, the particle Reynolds number increases from 10 to 400. The comparison between the simulation results and the empirical values are listed in Table , where a good agreement is obtained.

Table 4

Comparison of Particle Sherwood Numbers Obtained from Simulation and Empirical Correlation for Various Particle Reynolds Numbers

Res	Empirical	Simulated	Relative error
10	3.90	3.68	–5.58%
20	4.68	4.51	–3.70%
40	5.79	5.67	–2.15%
60	6.65	6.56	–1.32%
100	8.00	7.94	–0.73%
200	10.49	10.46	–0.29%
300	12.39	12.26	–1.07%
400	14.00	13.75	–1.76%

For the case of an inert nonreactive sphere, the verification can be done by checking the mass balance of the system. In other words, the total amount of the species entering the system should equal the total amount of the species leaving the system if zero-flux Neumann boundary condition is imposed correctly at the sphere surface. The mass flux at inlet J and outlet J are computed by the following two equations respectively:where 1 and n indicate the first and the last grid point in the streamwise direction, respectively. It should be noted that the diffusive term is not included in the calculation of the outlet mass flux due to the enforcement of the zero-gradient Neumann boundary condition at the outlet for species computation. The mass balance was checked for all particle Reynolds numbers listed in Table , the outlet mass flux gave the same value as the inlet mass flux for all cases. For flows with Re ≤ 200, the deviation for each Re case is a constant with the magnitude of E-15 which is the computational accuracy we applied in the simulations. For higher Reynolds number flows, the deviation value is oscillating within the range of E-13 and E-12, which can be reasonably understood as the result of vortex dynamics behind the sphere. This verifies the excellent enforcement of the zero-flux Neumann boundary condition at the sphere surface.

Results

In the previous section, our proposed DNS model has been verified by two benchmark cases, where the simulation results reach a good agreement with well-established solutions. In this section, we will demonstrate results for two cluster systems simulated by applying our DNS model. These systems are closely related to engineering applications and reveal the strong points of our DNS model in the framework of fluid–particle systems. In Section , the mass transfer between the fluid phase and the solid phase are analyzed for different cluster proximities and configurations, by studying a gas flowing through a stationary cuboid cluster consisting of nine spheres. Subsequently, in Section , a more complex system is considered. A spherical cluster consisting of hundred spheres is generated by random packing, and the influence of inert particle dilution on the cluster mass transfer performance is determined.

Nine-Sphere Cuboid Cluster

In this system, we consider convective mass transfer to a cluster which is composed of nine spheres. The basic configuration of this cluster is a cube with eight spheres at its vertices, where its front face is perpendicular to the flow direction x. One single sphere is added into the cube and positioned in the center. All nine spheres are of the same size d = 0.005 m. The cluster proximity decreases gradually to assess the impact of clustering of particles on mass transfer behavior. The proximity of a cluster S is measured as the center-to-center spacing between the spheres in the front face. Five proximities are considered for this cluster system, with the detailed information listed in Table . The packing density is calculated asThe central sphere is fixed at the location centrally in y and z directions. In the flow direction, a minimal distance of 2d and 6d is maintained for the front sphere to the inlet and the back sphere to the outlet, respectively. In the lateral direction, a minimum distance of 3.5d is reserved for spheres to the domain boundaries. The computational domain and particle arrangement for the nine-sphere cuboid cluster is shown in Figure . It should be noted that the distance between the front face and the back face is also the cluster proximity S in this specific cube case. The data used for the simulations and the boundary conditions applied at the domain boundaries are the same as those in Section . Based on the mesh convergence test in the last section, the particle diameter is 20 times larger than the grid size.

Table 5

Sphere–Sphere Distance and Corresponding Packing Density for Five Cluster Proximities

Cluster proximity S	3ds	2ds	1.8ds	1.6ds	1.4ds
Packing density η	0.074	0.175	0.215	0.268	0.341

Figure 1

Simulation domain and particle configuration for the nine-sphere cuboid cluster.

In the simulations, three particle Reynolds numbers 10, 50 and 200 are used to assess the influence of different flow patterns on the mass transfer behavior. For each simulation (5S × 3Re), we consider two cases. In one case, all nine spheres are active, whereas in the other case, the only active sphere is positioned in the center with eight inert spheres positioned at vertices. Figure shows the concentration distribution in the central plane of the nine-sphere cuboid cluster with different proximities for both all-sphere-active case and central-sphere-active case at Re = 200. From the figure, steady and symmetric concentration fields are observed, and we can clearly identify the clustering effect. The wake behind the central sphere becomes wider with smaller cluster proximity, and finally separates into two pieces at S = 1.4d due to the strong blockage effect of the back spheres. The influence of eight spheres at vertices is also demonstrated in the figure. For the same cluster proximity, the wake is always thinner if the eight surrounding spheres are inert, as the wake will include part of the mass transfer effects of these eight spheres if they are active. This difference is especially distinct at the smallest cluster proximity, the back spheres in the inert case solely split the wake without further reducing the species concentration.

Figure 2

Concentration distribution in the central plane of the nine-sphere cuboid cluster with different proximities at Re = 200. From top to bottom, the cluster proximity decreases from 3d to 1.4d. The all-sphere-active case and the central-sphere-active case are shown in the left and right column, respectively.

In Figure , the development of the particle Sherwood numbers, calculated by Equation , of all spheres for the all-sphere-active case is shown. Due to the symmetric geometry, as expected, the Sherwood numbers of the four spheres in the front/back face are the same and indicated by “front sphere” and “back sphere” in the figure. For the front sphere, higher Reynolds number results in a more pronounced improvement of its mass transfer behavior at smaller cluster proximities. The central sphere has similar behavior, but the minimum Sherwood number always occurs at the smallest cluster proximity. Particle clustering has a large influence on the mass transfer performance of the back spheres. With increasing Reynolds numbers, the Sherwood number increases less compared to the ones of the front and the central sphere.

Figure 3

Development of particle Sherwood numbers of all spheres along with the cluster proximity. Solid lines are of the front sphere, dashed lines are of the central sphere and dotted lines are of the back sphere. Simulations with the Reynolds number of 10, 50 and 200 are indicated by blue-plus, red-circle and green-star lines, respectively.

The comparison of the particle Sherwood numbers of the central sphere between all-sphere-active case and central-sphere-active case is shown in Figure . It is clearly observed that for all cluster proximities the central sphere has a better mass transfer performance in case the eight surrounding spheres are inert. However, this difference can be dramatically reduced by using higher fluid velocities.

Figure 4

Comparison of the particle Sherwood numbers of the central sphere at varying cluster proximities. Solid lines are of the central-sphere-active case whereas dashed lines are of the all-sphere-active case. Simulations with the Reynolds number of 10, 50 and 200 are indicated by blue-plus, red-circle and green-star lines, respectively.

To study the influence of cluster geometry on the mass transfer behavior, four more configurations are considered for the nine-sphere cuboid cluster: 45° rotation in z plan, 45° rotation in z plan following 45° rotation in y plan, elongation of the front-back face distance to 2 times of the cluster proximity and elongation of the front-back face distance to 4 times of the cluster proximity. These four extended cluster configurations, shown in Figure , are noted as Rotated 1, Rotated 2, Half and Quarter, respectively. The basic configuration is indicated by Normal in the following text.

Figure 5

Four extended configurations for the nine-sphere cuboid cluster. From left to right, the configuration is Rotated 1, Rotated 2, Half and Quarter, respectively.

In Figure , we show the concentration distribution in the central plane of the nine-sphere cuboid cluster for all five configurations. The examples are given for both all-sphere-active case and central-sphere-active case at the Reynolds number of 50 and the smallest cluster proximity S = 1.4d. The Normal and the Rotated 1 configurations have quite similar behavior, with the wake significantly split into two pieces. The Rotated 2 configuration behaves more like an isolated big particle, where the inert particle enforcement can be clearly visualized at the front, top, bottom and back spheres. With larger distance between the front and back faces, spheres have more independent behavior. In the Half configuration, no clear “split” behavior is observed in the wake region, whereas each sphere behaves almost as an isolated one in the Quarter configuration. For all configurations, a much lower concentration field is detected for the cluster in case all nine spheres are active.

Figure 6

Concentration distribution in the central plane of the nine-sphere cuboid cluster with different configurations at Re = 50 and S = 1.4d. From top to bottom, the cluster configuration is Normal, Rotated 1, Rotated 2, Half and Quarter, respectively. The all-sphere-active case and the central-sphere-active case are shown in the left and right column, respectively.

The Sherwood numbers of the central sphere in all five cluster configurations are plotted as a function of the cluster proximity in Figure for both all-sphere-active and central-sphere-active cases. Similar behavior of the Sherwood number profiles are observed, but with a much denser distribution in the central-sphere-active case. At the same fluid velocity, the central sphere mass transfer performance is pronouncedly improved at smaller proximities in case the eight vertex spheres are inert. Among five cluster configurations, Rotated 2 configuration has the worst mass transfer performance, which is more distinct at higher Reynolds numbers. A unique profile is noticed for the Rotated 2 configuration at Re = 200 that the Sherwood number increases at small cluster proximities.

Figure 7

Development of the particle Sherwood number of the central sphere along with the cluster proximity for all-sphere-active case (left) and central-sphere-active case (right). Simulations at Reynolds number 10, 50 and 200 are represented by solid lines, dashed lines and dotted lines, respectively. Blue-plus, red-circle, green-star, yellow-cross and black-triangle lines indicate the cluster configuration of Normal, Rotated 1, Rotated 2, Half and Quarter, respectively.

For this nine-sphere cuboid cluster system, in total 150 simulations were performed. Results are introduced and discussed qualitatively in the current section; nevertheless, we refer for detailed simulation results in the Supporting Information.

Random-Generated Spherical Cluster

In this section, our proposed DNS model is applied to a dense stationary cluster, which is a more physically complex system containing 100 random-generated spherical particles. All spheres are of the same size d = 0.005 m. The cluster is created by the hard-sphere Monte Carlo method and distributed in a spherical configuration with a predefined solid phase packing density η of 0.3. It should be noted this value is selected on one hand according to the significant clustering effect revealed at such a packing density in the previous section, on the other hand on the basis of a typical value of realistic clusters observed in riser flows. This gives that the diameter of the cluster d is about seven times the particle diameter d. The simulations are computed on a 3D domain with a length of 0.21 m (6 × d) in the flow direction and a length of 0.07 m (2 × d) in the lateral direction. In the flow direction, the center of the cluster is located at 1.5d from the entrance. The empty entrance and exit sections are essential to allow for flow development and avoid outflow recirculation, especially at high particle Reynolds numbers. It should be noted that the cluster Reynolds number will reach a high value even with intermediate particle Reynolds number due to the fact that d = 7d. The computational domain and particle arrangement are shown in Figure . The same configuration is used for all following simulations. It should be noted that the number of particles (i.e., 100) contained in the current cluster is not meant to be a large enough representative sample for proper statistical analysis.

Figure 8

Simulation domain and particle configuration for the random-generated spherical cluster.

In the simulations, fluid flows through the cluster with a prescribed uniform inlet velocity, with the value of 0.08, 0.2 and 0.4 m/s corresponding to the particle Reynolds number Re of 20, 50 and 100 (cluster Reynolds number Re of 140, 350 and 700), respectively. At the inlet, the fluid enters the system with a uniform species concentration of 10 mol/m3. For lateral domain boundaries, free-slip boundary condition is enforced for the velocity field computation, whereas zero-flux Neumann boundary condition is applied for species concentration computation which implies an isolated system. The pressure at the outlet is prescribed as the standard atmospheric pressure, and a zero concentration gradient is set there for the species computation. The same parameters are used for the simulations as in Section and Section . For the current work, there are around 70 million grid cells in total and the parallel computation were performed for around 3 months on 24 Intel Xeon E5-2690 processors. Studies for higher fluid velocities are possible, however might significantly increase the computational time and make the simulation extremely expensive.

The concentration distribution together with the computed velocity field are shown in Figure for the central plane of the spherical cluster at three different Reynolds numbers. In these simulations, all particles are active catalysts with an external mass transfer limited chemical reaction proceeding at the external surface. The mesh resolution, which is the ratio of the particle diameter to the grid size, applied in the simulations is 20. This value is selected according to the mesh convergence test done in the previous verification case and our published work simulating a dense particle array with the same solid phase packing density.[65] Inside the particles both computed concentration and velocity fields are zero due to the assumption of nonporous catalysts with reactive external surface and the enforcement of the no-slip boundary condition at the sphere surface, respectively. In Figure , it is observed that the cluster behaves like a large isolated sphere with the diameter of the cluster size d. A wake is observed at the rear of the cluster. At particle Reynolds number of 20, corresponding to cluster Reynolds number of 140, the streamlines converge much more slowly behind the cluster than they diverge before the cluster. The flow and concentration fields are steady and approximately axisymmetric. At particle Reynolds number of 50, the cluster-based Reynolds number is 350, and two circulating eddy rings are found behind the cluster. With further increased Reynolds number Re = 100 (Re = 700), alternating vortex shedding appears in the wake of the cluster. Inside the cluster, both concentration and velocity values are much lower than those outside. However, as fluid flow has stronger penetrating capability at higher fluid superficial velocities, the velocity and concentration are less depleted with increasing Reynolds numbers. Inside the cluster, preferred flow pathways are clearly observed at locations with low packing density.

Figure 9

Concentration distribution (left) and computed velocity magnitude (right) at the central plane of the spherical cluster for the case of all active spheres. Panels a, b and c are of the particle Reynolds number 20, 50 and 100, with corresponding cluster Reynolds number 140, 350 and 700 respectively.

For industrial mass transfer processes, the influence of inert particle dilution is of high interest. The spherical cluster is therefore composed of a mixture of active catalysts and inert particles. The inert particles are randomly selected and distributed inside the cluster. In the case of an inert particle, fluid flows over it without any chemical reaction. This is simulated by enforcing the zero-flux Neumann boundary condition at the sphere surface. As a unique feature of our DNS model, variable boundary conditions of individual particles are considered consistently. In our work, five active ratios AR = 0.1, 0.3, 0.5, 0.7 and 0.9 are used to assess the dilution effect, together with the limiting case of all active spheres. The active ratio is defined as the ratio of the number of active particles to the total particle number:For all cases with AR < 1, three different distributions of inert and active particles are applied in the simulations to reduce the influence of relative locations among active catalysts (thus in total 45 + 3 = 48 simulations were performed). In Figure , concentration distributions for the central plane of the spherical cluster are shown for six active ratios at the particle Reynolds number of 50. The mixture of active catalysts and inert particles, namely the enforcement of two different boundary conditions, are clearly visualized in the figure. For example in panel a, only four spheres are active in the current plane identified by the surrounding concentration gradient whereas the other spheres are inert. With increasing active ratio, corresponding number of particles switch from inert sphere to active catalyst and a lower species concentration is observed inside the cluster and also in the wake region. It should be noted that the concentration distribution pattern in the wake region remains unchanged during this AR change process, namely two symmetrical standing vortices are always attached to the cluster.

Figure 10

Concentration distribution for the central plane of the spherical cluster with variable active ratios at Re = 50 (Re = 350). AR = 0.1, 0.3, 0.5, 0.7, 0.9 and 1.0 are presented by figures from panels a to f, respectively.

As an isolated structure consisting of multiply particles encountered in, for example, CFB riser flows, the mass transfer behavior of the cluster is worth being investigated in its entirety. The Sherwood number of the spherical cluster is calculated by using the following equation:where Φ is the total mass transfer rate from the fluid phase to the spherical cluster which results from all active particles, d is the diameter of the spherical cluster and A = 4πr2 is the external surface area of the equivalent solid particle of the cluster size. The driving force for the cluster Sherwood number is defined as the difference between the inlet concentration c and the averaged species concentration inside the cluster c. By applying this Sherwood number definition, the mass transfer of the reactant from the bulk fluid to the cluster surface, the diffusion inside the cluster together with the dilution influence of inert particles are accounted for integratively. The direct influence of more active catalysts contained in the cluster on the increased mass transfer rate is canceled out by taking the lower reactant concentration inside the cluster into consideration. In Figure , the Sherwood number of the spherical cluster is plotted as a function of active ratio at three Reynolds numbers. It should be noted that the average value of three active/inert particle distributions is used in this figure and the figures discussed hereafter, with the error bar representing the 90% confidence intervals calculated by the Student’s t-distribution. For the same active ratio, the cluster Sherwood number increases with higher Reynolds numbers. This is a similar behavior to the case of a single sphere, and the improved gas-cluster mass transfer is due to the higher species supply rate at increased fluid superficial velocities. For all three particle Reynolds numbers, the cluster Sherwood number increases with larger active ratios. In other words, a better interfacial mass transfer performance between the fluid and the cluster is obtained when the cluster contains less inert particles. From the figure, the influence of the distribution of inert and active particles can also be observed from the size of the error bars. In other words, different relative locations of active catalysts may result in considerable variations in the cluster Sherwood number. It has less influence with increasing active ratio as a more homogeneous system is obtained, and its influence is enlarged at higher Reynolds numbers. We also calculated the Sherwood numbers with the Frössling correlation applied to the cluster as a whole, which are indicated by dashed lines in the figure. The computed cluster Sherwood numbers are much higher than the empirical values, especially at increased Reynolds numbers. This can be reasonably explained by the flow through the cluster.

Figure 11

Influence of active ratio and Reynolds number on the cluster Sherwood number.

Because the spherical cluster is a porous-like structure, it is interesting to introduce the concept of “effectiveness factor” to account for the effect of particle clustering and also the effect of inert particles on the overall mass transfer efficiency of the cluster. In catalytic reaction engineering, the effectiveness factor Ω of a porous pellet is defined as the ratio of the actual overall reaction rate to the reaction rate that would result if the entire surface were exposed to the bulk concentration.[73] Due to the assumption of an infinite fast surface reaction proceeding at the surface of active catalysts, the mass transfer rate of the reactant from the bulk fluid to the cluster is equal to the rate of reaction consumption of the reactant of the cluster, and thus the effectiveness factor of the spherical cluster can be computed asIn this equation, Φ is the same quantity used in Equation and Φ is the “ideal” total mass transfer rate if the same number of particles are active and exposed to the bulk fluid phase in an extreme dilute configuration.In Equation , k is the mass transfer coefficient obtained from the empirical Frössling correlation, as given in Equation , A is the external surface area of individual particle and Δc is the concentration driving force which is defined as (c – 0) for current computation. As all particles are of the same size in the current simulation work, Equation can be rewritten aswith defined asand Sh defined in Equation . The ratio of the overall average particle Sherwood number to the “ideal” Sherwood number computed by the Frössling correlation is actually defined as contact efficiency γ by Venderbosch et al.[74] in their pioneering studies of mass transfer in riser reactors. For riser flows, especially in the coarser scale models (DPM and TFM), gas–solid contact efficiency is widely employed to quantify the deviation of a riser flow from an idealized plug flow,[75−77] which does not only consider the exposed area of a single particle to the surrounding gas phase but also includes the effect of the depleted reactant inside the gas pocket resulted from the particle cluster. In Figure , the influence of active ratio and particle Reynolds number on the effectiveness factor (contact efficiency) of the spherical cluster is shown. These results demonstrate that low efficiency (inefficient gas–solid contacting) is due to the agglomeration of particles into a cluster. Even with few active catalysts contained in the cluster, its efficiency is much below 1.0. However, this situation can be considerably improved by increasing the particle Reynolds number which gives a better gas permeance through the cluster. At a fixed gas superficial velocity, the cluster has a lower efficiency at increasing active ratios. In other words, the gas inside the cluster becomes more depleted of reactant as the cluster contains more active catalysts, which leads to a poorer gas–solid contacting. All these finding are consistent with both experimental and DPM simulation results published by Venderbosch et al.[74] and Carlos Varas et al.,[77] respectively.

Figure 12

Influence of active ratio and particle Reynolds number on the effectiveness factor (contact efficiency) of the spherical cluster.

For a highly reactive catalyst inside a cluster, which is our case, the whole mass transfer process is controlled by the external resistance from the bulk fluid to the catalyst surface. By using Equation for the cluster Sherwood number calculation, the additional mass transfer resistance attributed to the formation of a cluster is accounted for. We are now interested in the essential mass transfer behavior over the gas film around the active catalysts in the cluster, which we name as the effective particle Sherwood number . Due to the consistent mass transfer rate between the fluid phase and the solid phase, a relationship among the overall particle Sherwood number , the cluster Sherwood number Sh and the effective particle Sherwood number is obtained, which is written asIn this equation, N is the number of active particles, and A, d and A, d are the external surface area and the diameter of the particle and the spherical cluster, respectively. In Figure , the effective particle Sherwood number is plotted at six active ratios and three Reynolds numbers. At each Reynolds number, the data can be perfectly fitted by a linear function. This behavior indicates that the essential mass transfer performance of the active catalyst contained in the cluster can be linearly improved by increasing the number of active particles. It is understood as the result of a more heterogeneous fluid concentration field inside the cluster. For all active ratios, higher effective particle Sherwood number is obtained at increased Reynolds numbers.

Figure 13

Influence of active ratio and Reynolds number on the effective particle Sherwood number, with corresponding linear fitting profiles.

From the DNS results, further detailed information can be obtained such as the Sherwood number of individual particles. The local particle Sherwood number is computed by using Equation , with the inlet species concentration c substituted by the local average concentration c. In this case, the driving force is defined locally as the difference between the surface concentration of the active catalyst which is zero by our assumption and the local average species concentration c which is obtained from the following expression:The integration is performed over the local fluid volume surrounding each particle, and L is the distance between the local fluid point and the center of the particle. A cubic box of the size 5d with its center coinciding with the center of the particle is selected to be the local fluid volume. The method for computing the local average concentration c together with the selection of the local fluid box were published in the previous work of our group,[66,78] and we refer to these papers for a detailed description. For the computation of the local particle Sherwood number, the mass transfer rate Φ, calculated by the integration of the concentration gradient at sphere surface, automatically goes to zero for inert spheres. In Figure , the Sherwood numbers of individual particles for the case of AR = 1 are plotted as a function of the dimensionless radial distance to the cluster center at three Reynolds numbers. At Re = 20 the particles at the boundary of the cluster (r/R > 0.8) have higher Sherwood numbers than most of the inner particles. However, this behavior does not hold at higher Reynolds numbers. The increased Reynolds number leads to a larger improvement of the Sherwood number of the inner particles than the particles at the boundary shell, especially those close to the center of the cluster. In other words, high values of the local particle Sherwood number do not predominately occur at the cluster boundary, while they could also occur inside the cluster at high Reynolds numbers. This phenomenon can be reasonably explained by the nature of clusters that the gas bypassing is enhanced at higher fluid velocities and simultaneously the fluid has stronger convection through the cluster. The existence of preferred fluid pathways inside the cluster, which is a pronounced property of the porous medium, will further amplify the latter effect. Although this figure reveals a better mass transfer performance for almost all particles at higher Reynolds numbers, few exceptions are observed which can be explained by the same reason above. The distributions of the individual particle Sherwood numbers in radial direction are shown in Figure for increasing active ratios at Re = 50. This figure confirms the previous finding that particles deep inside the cluster might also have a good fluid–solid mass transfer behavior regardless of the active ratio. The particle Sherwood numbers are observed to have a wider spreading in values as active ratio increases. Although in a previous discussion we concluded that the overall average particle Sherwood number decreases with increasing active ratios, it is interesting to observe in this figure that not all individual particles have smaller Sherwood numbers at larger active ratios. Some particles have a completely opposite behavior. This is due to the varying heterogeneous concentration field resulting from the switch of particles from inert to active in the neighborhood.

Figure 14

Distribution of the local particle Sherwood numbers of the spherical cluster for the case of AR = 1, at three Reynolds numbers.

Figure 15

Distribution of the local particle Sherwood numbers of the spherical cluster for five active ratios, at Re = 50.

Conclusions

In this paper a previously introduced ghost-cell based immersed boundary method is applied to study fluid-particle mass transfer for flow passing through particle clusters consisting of active catalysts and inert particles. Taking the advantage of a second order quadratic interpolation scheme utilized in the reconstruction procedures, mixed boundary conditions, Dirichlet and Neumann boundary condition for active and inert particle respectively, can be realized consistently at the exact position of the particle surface.

For the nine-sphere cuboid cluster, it is found that for almost all cases the mass transfer performance of the central sphere decreases due to the formation of the cluster. The only exception is the Rotated 2 configuration at Re = 200 for which the Sherwood number increases with smaller cluster proximities. The mass transfer performance of the central particle is improved if it is surrounded by inert particles, especially at small cluster proximities. Higher Reynolds number will increase the Sherwood number under any circumstance. For the random-generated spherical cluster, the cluster effectiveness factor, also known as contact efficiency in large scale DPM and TFM models, is much below 1.0 and further decreases as the cluster contains more active catalysts. The cluster Sherwood number and the effective particle Sherwood number, which describe the mass transfer from the bulk fluid to the inner cluster and the mass transfer over the gas film around the active catalysts respectively, increase with more active catalysts contained in the cluster. Regarding the local particle Sherwood numbers, it is found that high values can also occur deep inside the cluster at increased Reynolds numbers and the value may increase with higher active ratios. The distribution of active/inert particles may lead to large variations of the mass transfer behavior, which decreases with higher active ratios and increases with higher Reynolds numbers.

Based on the current work, it is evident that our DNS model is a powerful tool for a systematic study of the mass transfer behavior of particle clusters. Besides the active ratio and Reynolds number, more variables of a cluster such as solid volume fraction, shape, size and orientation can be investigated. The particle level Sherwood number reveals the inaccurate description of the local heterogeneity by using mass transfer correlations developed for homogeneous systems (which is widely used in DPM and TFM simulations). The Sherwood number and the effectiveness factor at the cluster level can be parametrized and incorporated into mass transfer closures for aforementioned coarser scale models. This will significantly improve the prediction of the global mass transfer phenomenon.