Direct numerical simulation of fluid flow and mass transfer in dense fluid-particle systems with surface reactions

a r t i c l e i n f o


a b s t r a c t
In this paper, an efficient ghost-cell based immersed boundary method is introduced to perform direct numerical simulation (DNS) of mass transfer problems in particulate flows. The fluid-solid coupling is achieved by implicit incorporation of the boundary conditions into the discretized momentum and species conservation equations of the fluid phase. Taking the advantage of a second order quadratic interpolation scheme utilized in the reconstruction procedures, the unique feature of this ghost-cell based immersed boundary method is its capability to handle mixed boundary conditions at the exact position of the particle surface as encountered in systems with interplay between surface reactions and diffusion.
A fixed Eulerian grid is used to solve the conservation equations for the entire computational domain. Following a detailed verification of the method in the limiting case of unsteady molecular diffusion without convection, we apply our method to study fluid-particle mass transfer for flow around a single sphere and a dense stationary array consisting of hundreds of spheres over a range of Damköhler numbers.

Introduction
Fluid-particle flows are frequently encountered in a wide range of processes in the chemical, petrochemical and energy 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. In the past decades, numerous empirical correlations for fluid-particle mass transfer have been proposed for both single particle and multi-particle systems (Gunn, 1978;Wakao and Funazkri, 1978;Hughmark, 1980;Wagner et al., 1997;Xu et al., 2000;Piché et al., 2001). Although these correlations are helpful for a quick and rough estimation of mass transfer rates for design purposes, they do not consider the influence of different reactor geometries and wider ranges of reaction conditions, especially variable reaction rates in complex particle configurations. In other words the interplay between reactivity and transport is not easily quantified on basis of these empirical mass transfer correlations. In this case, detailed information such as local values, spatial https://doi.org/10.1016/j.ces.2017.10.018 0009-2509/Ó 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). distribution and temporal evolution are hard to obtain. To overcome this problem, three-dimensional transient computer simulation has rapidly gained interest due to its capability to produce detailed quantitative information instead of an average mass transfer rate.
With the development of computational technology, DNS has emerged as a powerful tool to resolve all the details at the smallest relevant length scales and quantitatively derive microscale 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. Taking the advantages of efficient CPU/memory utilization and easy grid generation, IBM has been applied in various studies including complex geometries, moving particles and deformable immersed objects (Fadlun et al., 2000;Udaykumar et al., 2001;Tseng and Ferziger, 2003;Mittal et al., 2008;Wang et al., 2008;Takeuchi et al., 2010). Recently Das et al. (2017b) used IBM for flow through complex cellular porous structures. Following the fluid flow equations, additional equation for species transport can be added with relative ease using the same methodology.
The IBM was first introduced by Peskin (1977) for simulation of blood flow around the flexible leaflet of the human heart. The main idea of this method is to use a Cartesian grid for fluid flow simulation where the immersed boundary is represented by Lagrangian marker points. A forcing term is introduced in the transport equation to represent the interaction between the immersed boundary and the fluid, whose magnitude is taken such that the boundary conditions are fulfilled in an interpolated manner. A regularized Dirac delta function is used to distribute this singular force over the Eulerian cells surrounding each Lagrangian point. This method is categorized as a continuous forcing method (CFM), and many researchers have contributed to the further development of this method (Goldstein et al., 1993;Saiki and Biringen, 1996;Uhlmann, 2005;Wang et al., 2008;Vanella and Balaras, 2009;Di and Ge, 2015). The second category of IBM is referred to as discrete forcing method (DFM), which was first proposed by Mohd-Yusof (1997), and later extended by Fadlun et al. (2000), Tseng and Ferziger (2003), Marella et al. (2005), Ghias et al. (2007), Zhang and Zheng (2007), Haugen and Kragset (2010), Seo and Mittal (2011), Zeng and Farhat (2012) and Lee and You (2013). In this method the ghost cell approach is applied, where virtual variable values inside ghost cells (i.e. inside the immersed body) are calculated by applying boundary conditions to extrapolate fluid variables near the boundary to the ghost cell. DFM treats the immersed boundary as a sharp interface, and does not require the explicit addition of a discrete force in the governing equations, thus the stability conditions are the same as that without treatment of the immersed boundaries.
Although IBM has been widely used for studies of momentum transfer in fluid-solid systems, very few computational results are available in the field of mass transfer. Deen and Kuipers (2013) applied three-dimensional DFM with a directional quadratic interpolation scheme to mass transfer problems with infinitely fast reactions at the particle surface in dense fluid-particle systems. Feng and Musong (2014) studied heat and mass transfer of 225 spheres in a fluidized bed with the IBM method assuming Dirichlet conditions for the scalar variables at the surface of the spheres. The solid-fluid coupling is accounted for by forces located at the Lagrangian surface points that are distributed to the grids by regularized delta-functions. Gong et al. (2014) also used CFM for mass transfer across multiple deformable moving interfaces where the interfacial flux is calculated from the concentration jump and then interpolated to the Eulerian grids over the entire domain.  In this paper, an efficient ghost-cell based immersed boundary method is proposed for the simulation of reaction rate controlled mass transfer problems in fluid-solid system. The reconstruction procedures involve a second order quadratic interpolation scheme. As the unique feature, the general mixed boundary condition is realized and enforced exactly at the particle surface at the level of the discrete conservation equations. Taking advantage of the quadratic interpolation scheme, a surface reaction rate is incorporated into the general mixed boundary condition to study the interplay between chemical reaction and mass transfer processes in fluid-particle systems. It should be noted that no heat effect is considered in the current work.
The organization of this paper is as follows. First, the description of the model is given, including the governing equations, numerical solution method, fluid-solid coupling and incorporation of surface reactions. Following a detailed verification case of unsteady diffusion, both single particle and multi-particle dense systems with different reaction rates 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, the fluid-solid coupling as well as the representation of the surface reaction embedded in the general mixed boundary condition. 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 spherical particles, which are nonporous catalytic pellets so that only the external surface is reactive. Both fluid and solid phases have constant physical properties.

Governing equations of fluid phase
The transport phenomena in the fluid phase are governed by the conservation equations for mass, momentum and species transfer, respectively given by: In the above equations, q f and l f represent the fluid density and the fluid viscosity respectively, and D f is the species mass diffusivity in the fluid.

Numerical solution method
The governing equations are solved by a finite difference scheme implemented for a staggered Cartesian grid. The grid is defined in three dimensions (3D) with a uniform grid spacing in all three directions. Building on the work of Deen et al. (2012), the numerical solution of the equations described in the previous section is acquired embedding high order discretization schemes as well as small computational stencils. We perform the time discretization of the momentum equation, which leads to the following expression: In this equation, n is the time step index. The convective and diffusive momentum fluxes C m and D m are calculated by spatial discretization of: The convection term is spatially discretized by a second-order total variation diminishing scheme, whereas the diffusion term is computed with a standard second-order central differencing scheme. The convection term is computed explicitly whereas the diffusion term is computed implicitly. The solution of Eq. (4) is achieved by using a two-step projection method where a tentative velocity field u ÃÃ is first computed by neglecting the pressure gradient contribution. As the second step, the velocity field at the new time step n þ 1 is obtained based on the new pressure gradient calculated from the Poisson equation at time step n þ 1. For the interested reader, we refer for a more detailed description of this method to the work of Deen et al. (2012).
The species convection-diffusion equation is temporally discretized in the same way as for the momentum equation, namely the Adams-Bashforth scheme is applied for the convective transport while a fully implicit Euler backward scheme is used for the diffusion term.
with the convective and diffusive fluxes C sp and D sp respectively given by: The spatial discretization methods used for the species transport equation is the same as for momentum equation. For both momentum and species equation, the boundary condition is enforced at exactly the immersed boundary surface, which is handled at the level of the discretized conservation equations and will be introduced in detail in the next section.

Fluid-solid coupling
The fluid-solid coupling constitutes the key element of our model. In order to impose variable boundary conditions in a sharp interface method, a ghost-cell based immersed boundary method is used. The discretization of the momentum and species equation leads to algebraic equations of the following generic form: where / is the fluid variable for which we want to obtain a solution, namely velocity and concentration field for momentum and species equation, respectively. This equation provides the relationship between any fluid quantity / c and its six neighboring points indicated as / nb . Since the surfaces of the immersed objects do not coincide with the mesh boundaries, we need a special treatment for the nodes near the fluid-solid interface. In our method, the first step is to identify all ghost points, which are defined as points inside the solid phase but possessing at least one neighbor in the fluid phase.
These points are used to obtain a solution of the governing equations exterior to the solid phase. The second step is to check every fluid node whether any of its six surrounding neighbors represents a ghost point. If this is the case, a boundary condition has to be applied. Fig. 1 demonstrates two different interpolation schemes used in our model (examples given in a 2D domain for illustration purpose. The shaded cells represent the initial computational stencil containing four neighbors). On the left, Fig. 1(a) the directional quadratic interpolation scheme is shown, which has the advantages of second order accuracy and a compact computational stencil. However, this interpolation scheme is limited to the application of the Dirichlet boundary condition, thus it is only used for enforcement of the no-slip condition in our model. The numerical velocity field is tested to be divergence free (Das et al., 2017b). After imposing the boundary condition, / 0 is eliminated from the stencil. For the interested reader, we refer for more details of this method to the work of Deen and Kuipers (2013).
For the general mixed boundary condition encountered in species transport calculation, we need to resort to quadratic interpolation, which is presented in Fig. 1(b), to describe the gradient in the normal direction to the surface of the immersed object. As illustrated in Fig. 1(b), / 0 is removed from the computational stencil after imposing the boundary condition, while / 1 and / 2 are added to it. This method is inspired by Luo et al. (2008), who applied a high-order polynomial to describe the boundary condition in acoustic problems. In their work and the following work of Xia et al. (2014), a third order polynomial is used and consequently 20 coefficients need to be determined. A spherical region of adaptively chosen radius R around the boundary point is searched to provide adequate fluid data points, and subsequently a weightedleast square error minimization is applied to determine these coefficients. In our model, a second order polynomial is used, meanwhile, a quadratic reconstruction procedure is chosen to obtain the corresponding coefficients.
In this quadratic interpolation scheme, a generic variable / in the vicinity of the immersed object surface can be approximated in terms of a second-order polynomial as follows: where x, y and z are relative coordinates with respect to the origin located at the boundary point.
Eq. (11) is in fact the approximation of / using the Taylor expansion near the boundary point: In the 3D case, the number of coefficients for a second-order polynomial is ten. In order to determine these coefficients c ijk , we need / values from nine neighboring fluid points and one image point. The image point is defined as the mirror point of the ghost point through the boundary in the normal direction, which has the same distance to the boundary point as the ghost point. For the 2D case, only five fluid points plus one image point are required for computation of c ij coefficients, as indicated in Fig. 1(b). It should be noted here that the resulting linear problem may be ill-conditioned in case the image point is close to the sphere surface. In this circumstance, the ghost value may be correct algebraically but problematic numerically. To remedy this problem, a minimal distance of half of the mesh size is maintained between the image point and the sphere surface in case the original distance falls below this value.
With adequate data points (ten values at fluid/image point for ten coefficients in a second-order polynomial), the resulting equation for coefficients c ijk can be written as follows: where / and c are the vectors for species concentration and coefficients respectively, and X is the Vandermonde matrix given by: To solve Eq. (13), the Vandermonde matrix is inverted by applying LU decomposition using the Crout algorithm. The coefficients c ijk are obtained by multiplication of the inversed matrix X À1 and the concentration vector /, which can be written as a linear combination of / values. For the treatment of the general mixed boundary condition the following coefficients are required: Therefore, for the general mixed boundary condition at the immersed object surface: the image point value can be evaluated by satisfying the boundary condition at the boundary point: where M m is defined according to the following equation, with the components of the normal unit vector indicated as n x , n y and n z respectively: Considering the general correlation between image point and ghost point, the following two equations can be obtained: where L I is the distance from the image point to the sphere surface, L G is the distance from the ghost point to the sphere surface, and DL is the mutual distance between the image point and the ghost point which equals L I þ L G . The value at the ghost point can be finally computed as: With Eq. (24), the matrix coefficients in Eq. (10) can be updated. Altered coefficients within the original stencil are incorporated in the implicit scheme, while neighbors outside the original stencil are accounted for in an explicit way. The procedure described above needs to be carried out for all ghost points to ensure that the desired local boundary condition applies everywhere at the immersed boundary surface. Note that the pressure, velocity and concentration field are obtained for the entire computation domain, i.e. also for the cells inside the particles, with our DNS model proposed above.

Surface reaction incorporation
In this work, we consider the situation where a first order irreversible surface reaction occurs on the external surface of spherical particles. The corresponding rate law is written as follows: where k is the surface reaction rate coefficient with the unit meter per second. At steady state there is no accumulation on the surface, so that the diffusion flux to the sphere surface (Eq. (26)) equals the rate of consumption due to reaction. In this circumstance, an equilibrium equation is obtained (Eq. (27)).
where D f is the species mass diffusivity and n is the normal direction pointing outward of the sphere. Eq. (27) is in fact of the functional form of the general mixed boundary condition (Eq. (19)), by setting the following parameters: We will now focus on a system with a spherical particle for which two dimensionless variables are defined based on the initial species concentration in the fluid c f ;0 and sphere radius R s : through which Eq. (27) is non-dimensionalized leading to the Damköhler number (Da) for surface reactions in fluid-particle systems: In case of a small Da, diffusion occurs much faster than the reaction thus the process is reaction limited. For a large Da, the reaction rate is much greater than the diffusion rate thus the overall conversion rate is limited by mass transfer. As a dimensionless number, Da provides a quick estimation of the degree of conversion that can be achieved in fluid-particle reactions (Fogler, 1999). It is important to realize that, at extremely high reaction rate and extremely low reaction rate the boundary condition actually becomes zeroboundary-value Dirichlet boundary condition and zero-flux Neumann boundary condition respectively. It should also be noted here, although we assume a first order irreversible surface reaction in current work, it is flexible to apply our model to reactions in any order. Eq. (27) is a straightforward modification of the general mixed boundary condition, where c f on both left and right hand sides are computed and incorporated into the discretized species conservation equation implicitly through fluid-solid coupling. In case of an nth order reaction, the reaction kinetics could be linearized to keep the format of a first order reaction, with c nÀ1 f integrated into the reaction rate coefficient k which is then calculated explicitly.

Results
In this part, we will show results obtained for three systems simulated by using our proposed 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 2.1 simulation results are compared with analytical solutions for the limiting case of unsteady molecular diffusion, which serves as a verification case of our DNS model. After that, our DNS model is applied to both a single particle system (Section 2.2) and a dense multi-particle array system (Section 2.3), to study the mass transfer behavior and chemical conversion for several reaction rates.

Unsteady diffusion around a single sphere
In this test case, we consider the unsteady diffusion of a chemical species to a sphere, followed by a first order irreversible chemical surface reaction. The sphere is positioned in the center of a large pool of quiescent fluid. The governing equation for unsteady mass diffusion in the fluid phase is described by Eq. (3), with u set to zero. The initial condition is given by: where R s is the sphere radius and c f ;0 is the initial value for the species concentration. The boundary condition at the confining walls is a prescribed concentration, which has the same value as the initial species concentration.
The boundary condition at the sphere surface is a special case of the general mixed boundary condition where the diffusion rate equals the chemical reaction rate.
D f @c f @r This simulation results are compared to the analytical solution of the spherical symmetric problem (so only r dependence), which is defined by the following expression: In the analytical problem there are no bounding walls. However, for short enough times no significant difference is expected as long as the diffusion fronts have not reached the domain walls in simulation. The solution of the instantaneous Sherwood number and time dependent concentration profile around the sphere are expressed by the following two equations. For a detailed derivation, we refer to Appendix A.
with s and x defined as: where Fo, the Fourier number, is the dimensionless time: For the DNS the particle is located in the center of a cubic box with a size of 0.04 m. Simulations are run for short enough times such that the diffusion fronts are kept far away from the domain boundary. The data used for the numerical simulations are listed in Table 1. It should be noted here that a sufficiently small time step is used to reduce temporal discretization errors. Mesh convergence tests were performed by using grids of 80 Â 80 Â 80, 160 Â 160 Â 160, 256 Â 256 Â 256 and 320 Â 320 Â 320 for infinitely large Damköhler number. These grids correspond to the mesh resolution of 10, 20, 32 and 40 respectively, which is defined as the ratio of the sphere diameter to the grid size. The simulation results of mesh convergence tests are given in Table 2, in which the instantaneous Sherwood number is compared with the analytical solution. In Fig. 2 the relative errors at Fo ¼ 0:032 are plotted as a function of the mesh resolution and a second order convergence is clearly observed. The Sherwood number is computed from the following expression: where U f !s is the mass transfer rate, with the normal pointing outward of the solid. This quantity is calculated by the integration of the concentration gradient at sphere surface over the whole sphere (with S s the external surface area of the sphere): The driving force in Sherwood number calculation is defined as the difference between the wall concentration c f ;0 (namely the initial concentration) and the surface concentration c B . For general mixed boundary condition, c B is the average species concentration at sphere surface. The calculation of c B will be described in the next paragraph. From Table 2, excellent spatial convergence as well as a decreasing relative error along with the time development are observed. The larger errors at initial time steps can be explained by the extremely steep concentration profile near sphere surface, which cannot be fully captured by the quadratic interpolation scheme. Taking both accuracy and computational cost into consideration, mesh resolution N = 20 is considered to offer an acceptable trade-off. Consequently all subsequent simulations were computed on a 160 Â 160 Â 160 grid with a uniform grid size of 2.5 Â 10 À4 m in all directions. Six reaction rates varying from slow reaction to fast reaction are imposed at the sphere surface, which correspond to the Damköhler numbers of 0.01, 0.1, 1, 10, 100 and infinity. From the simulation, the species concentration in the fluid phase around the sphere is obtained directly. The concentration value at the exact position of the sphere surface could also be obtained through the definition of our ghost-cell based immersed boundary method by post-processing. The local surface concentration and the local concentration gradient are computed based on Eqs. (22) and (23), respectively. The value of the overall average sphere surface concentration c B is the surface average of local ones. It should be noted that, the local concentration gradient here is in fact the one used in the calculation of the mass transfer rate in Eq. (44). The concentration values at the first ten grid points near the sphere surface as well as the one at the exact sphere surface at Fo ¼ 0:032 are plotted in Fig. 3 and compared with the analytical solutions. As indicated before, the diffusion fronts are kept away from the confining walls to maintain the boundary condition for simulations. Even for the infinitely fast reaction, the species concentration is still very close to the initial value at the distance of one sphere radius from the sphere surface. As clearly demonstrated in the figure, very accurate results are obtained for the concentration profiles for all Damköhler numbers, including the steep ones at high reaction rates. The evolution of the concentration profiles around the sphere are shown in Fig. 4, for a Damköhler number Da ¼ 100. In this figure, the surface value and the values of the first five fluid points are presented and reveal a good agreement with the analytical solution. As expected the biggest deviations are obtained for small values of Fo because of the existence of steep concentration profiles during the early stage of the process. Fig. 5 shows the development of the instantaneous Sherwood number (calculated by Eq. (43)) for six different Damköhler numbers, including the limiting case of infinite fast surface reaction, versus the Fourier number. The profiles are bounded by the profile for an infinitely fast reaction as the lower limit and the slowest reaction as the upper limit. The main difference between simula-tion results and analytical solutions exist in the initial time steps. Large deviations are found there due to the incomplete capture of the concentration profile by our interpolation scheme, but later the deviation decreases and the simulation results converge to the analytical solution well.
The increased Sherwood number for Da ! 0 might be unexpected, because this limit is associated with infinitely slow reactions. However, the Sherwood number does not indicate the magnitude of the mass flux but the ratio of the flux to the driving force (i.e. difference between bulk and surface concentrations). For Da ! 0 both tend to zero. At small Damköhler number the (very small) flux to the particle surface is fully dominated by the surface reaction rate as it occurs at the fixed bulk concentration. Therefore the limit Da ! 0 gives the Sherwood number for the case of imposed flux. Both the simulations and the analytical solution show that the Sherwood number at constant flux is larger than the Sherwood number at constant concentration.

Single stationary sphere under forced convection
In this system, we consider convective mass transfer to a single stationary sphere in an enclosure, where a first order irreversible chemical reaction proceeds at the sphere surface. Fluid flows into the system containing a single species with constant inlet concentration of 10 mol/m 3 . The sphere is located at the center of the domain laterally while it is positioned at a distance of two times of the sphere diameter from the inlet in the flow direction. The data used for the numerical simulation are summarized in Table 3. 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 eight to remove the confining effect from the simulation boundaries. Free slip boundary conditions are applied  at the domain boundaries for the calculation of velocity field, while a Neumann boundary condition (non-penetrating walls) is used for the species conservation equation. For the simulation a uniform fluid velocity is imposed at the inlet, whereas the pressure at the outlet is prescribed as the standard atmospheric pressure. The spe-cies outlet boundary condition is specified as a zero slope condition. For this single sphere system, simulation results using either a Dirichlet boundary condition or a Neumann boundary condition at the sphere surface were presented in our previous work (Lu et al., 2016). In that work, the Sherwood number for both boundary conditions were calculated, and we found an improved mass transfer performance in case of the Neumann boundary condition. In this work, we consider two groups of simulations, with only a reactive surface boundary condition. In one group the effect of the particle Reynolds number is assessed at fixed reaction rate whereas in the other group at fixed Reynolds number the effect of the rate of the surface reaction is assessed.
For the first group, an intermediate reaction rate (Da ¼ 1) is enforced at the sphere surface, and five representative particle Reynolds numbers are used to assess the influence of different flow patterns on the interplay between convective mass transport and surface reaction. The Sherwood number, total fluid-solid mass transfer rate and the species concentration at sphere surface   Reynolds number, which indicates a better interfacial mass transfer behavior between the fluid phase and the solid phase. It is interesting to find the species concentration at sphere surface is also increasing with higher Reynolds number, which is an expected behavior because the species supply rate increases whereas the kinetics of the chemical reaction remains constant. In this circumstance, the driving force for mass transfer is actually decreasing, thus the increased mass transfer rate from the fluid to the particle is fully due to the reason of larger fluid flow velocity. The concentration distribution around the sphere for each Reynolds number is shown in Fig. 6. From left to right, the Reynolds number is 20, 50, 100, 200 and 400 respectively. The boundary layer becomes thinner at higher Reynolds number, whereas the wake region behind the sphere narrows. In the wake region, a circulation appears at Re ¼ 200 and it further develops into vortices at Re ¼ 400. In contrast to the distinct concentration difference developed under Dirichlet boundary condition (Lu et al., 2016), a less pronounced concentration change is observed which is in fact due to the smaller reaction rate. At Da ¼ 1 and increasing Sh the mass flux to the particle surface becomes more and more reaction dominated. The Sherwood numbers listed in Table 4 are computed from Eq. (43), by replacing the initial concentration c f ;0 by the inlet concentration c f ;in . However recall the definition of the first order irreversible surface reaction, in which the external mass flux to the sphere surface is equal to the consumption rate on the surface: From this equation, the average surface concentration c B can be calculated and consequently the overall effective mass transport coefficient could be defined as: In our case (Da ¼ 1), the surface reaction rate coefficient k is specified as 0.008, whereas the external mass transfer coefficient k m is given by the empirical Frössling correlation: where Re s is the particle Reynolds number and Sc is the Schmidt number, respectively defined as: Sc By comparing these two values, the external mass transfer coefficient is only a few times the surface reaction rate coefficient, thus this process is neither in a mass transfer limited regime nor in a reaction limited regime. In Table 5, the overall effective mass transport coefficients obtained from the simulations are compared with the values calculated by Eq. (46), and they are in a good agreement. For the second group of simulations, the particle Reynolds number is specified to be 100, whereas variable reaction rates are imposed at the sphere surface which lead to Damköhler numbers increasing from 0.001 to 1000. The particle Sherwood number, concentration gradient and concentration value at the sphere surface computed from our simulations are given in Table 6. The gradient and the surface concentration are the overall averaged values over the whole external sphere surface. The Sherwood number is observed to increase with higher reaction rate up to Da ¼ 1 and then decrease with further increase of the Damköhler number. However, it is interesting to see that the Sherwood numbers for different Damköhler numbers are very close, i.e. within a deviation Table 4 Computed Sherwood number, total fluid-solid mass transfer rate and sphere surface concentration for a single stationary sphere with first order irreversible surface reaction (Da ¼ 1) for five particle Reynolds number.  of 5%. As expected, species build up along the boundary and the surface concentration approaches the inlet concentration in case of slow reactions, whereas large concentration gradient and close-to-zero concentrations are observed at the sphere surface for fast reactions. The overall effective mass transport coefficient for each reaction rate obtained from the simulations is also listed in Table 6 and compared with the empirical values computed from Eq. (46). Good agreement is found for all Damköhler numbers, except for extremely slow reaction (Da ¼ 0:001). However considering the negligible overall mass transport behavior or in other words the nearly zero-flux Neumann boundary condition in this case, the calculation of the overall effective mass transport coefficient does not make much sense any more. The overall mass transport behavior goes from the reaction limited regime to the mass transfer limited regime in this group of simulations.

Reynolds number
Finally, it is verified that in the steady state the calculated species consumption rate ðkc B Þ equals the calculated mass flux to the sphere surface D f @c f @r , for all simulation cases above. It should be mentioned that both c B and @c f @r here are the overall averaged values over the whole sphere surface. This illustrates the correct enforcement of the rate-controlled reactive boundary condition (general mixed boundary condition) at the particle surface.

Dense stationary array
In this section, our proposed DNS model is applied to a dense stationary array, which is a more physically complex system containing a relatively large number of particles. In this simulation, fluid flows through a stationary random array of particles, and a first order irreversible chemical reaction proceeds at the surface of these particles. The particle array is created by the hardsphere Monte-Carlo method and distributed in a random configuration over the computational domain with a predefined solid phase packing density of 0.3. The data used for the simulations are listed in Table 7. The simulations are computed on a 3D domain with a length of 0.3 m in the flow direction and a length of 0.025 m in the lateral direction. In total, 716 spheres are positioned inside the simulation domain in a random fashion. It should be noted here that part of the 716 spheres seem cut by the boundaries, but actually these are not cut as periodic boundary conditions are applied in the cross-sectional directions. The first 100 grids and the last 100 grids in the flow direction are utilized as inlet and outlet regions respectively, hence no sphere is positioned in these two sections. The incorporation of these two empty sections is essential to avoid problems for inlet flow development and outflow recirculation respectively, at high particle Reynolds numbers. For the simulation a prescribed uniform fluid velocity is imposed at the inlet, with the value of 0.32 m/s, 0.64 m/s and 0.96 m/s corresponding to the particle Reynolds number of 80, 160 and 240 respectively. At the inlet, the fluid enters the system with a uniform species con-centration of 10 mol/m 3 . For lateral domain boundaries, the periodic boundary condition is applied for both velocity field and species concentration calculation. The pressure at the outlet is prescribed as the standard atmosphere pressure, whereas a zero concentration gradient (Neumann boundary condition) is set there for the species computation.
In Fig. 7, the particle configuration of 716 spheres is shown, together with the computed velocity distributions at the particle Reynolds number of 160 for the central plane of the array of particles. The right sub-figure is a small array used for mesh convergence test, which is randomly cut from the big particle array. The lateral size is maintained the same while the streamwise length is one tenth of the original one, which gives a cubic array with the dimension of 0.025 Â 0.025 Â 0.025 m 3 . In this small particle array, there are in total 84 spheres. In Fig. 7 the periodically distributed particle configuration across the domain boundaries is demonstrated clearly in both 2D and 3D views, whereas the preferred flow pathways in the array is observed in the velocity map. Inside the particles the computed velocity field is zero due to the particular enforcement of the no-slip boundary condition at the sphere surface. For the reaction proceeding at the sphere surface, four reaction rates are specified, corresponding to the Damköhler number of 0.1, 1, 10 and infinity. Although mesh convergence has been tested in previous two sections, here we use the case of Da ¼ 1 with an even higher particle Reynolds number of 250 to check the potential particle-particle mutual influence in the array. Two finer meshes are used, corresponding to the mesh resolution of 32 and 40. The total mass transfer rate (in the unit mol/s) of the system is considered to be a good quantitative indicator for the mesh convergence test. With mesh resolution of 20, 32 and 40, simulations give the results of À6.366EÀ04, À6.475EÀ04 and À6.504EÀ04 respectively. If we extrapolate the results of the performed simulations to the limiting case of infinite mesh resolution (grid size approaching zero) by a second order polynomial, the value is estimated to be À6.583EÀ04. In this circumstance, the  deviation is 3.30%, 1.64% and 1.20% for aforementioned increasing mesh resolutions respectively. Previous work from our group (Das et al., 2017a) applied the same way for the mesh convergence test in large domain simulations. This illustrates a good mesh convergence behavior, and indicates the proper selection of the mesh resolution of 20 for all simulations. For industrial mass transfer processes, the evolution of the cup-average concentration is of high interest, which is defined by: In this equation, the integration is performed over a surface S f 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, S f only accounts for the part occupied by the fluid inside this cross-sectional surface, whereas the points inside the solid phase are not considered in calculation. In Fig. 8, the cup-average concentration is shown as a function of the axial coordinate of the array, for the range of used particle Reynolds numbers and enforced Damköhler numbers. In this figure, the inlet and outlet regions can be clearly recognized, which are characterized by a constant cup-average concentration (for x < 0:025 and x > 0:275). For the same reaction rate (Damköhler number) proceeding at the sphere surface, higher Reynolds numbers will slow down the concentration decay rate and hence less reactant is converted in the array. This influence is especially distinguished in the cases of small and intermediate Damköhler numbers (Da ¼ 0:1 and 1, respectively), whereas the difference is mainly observed in the developing region for large Damköhler numbers (Da ¼ 10 and Da ¼ 1) as most reactant is finally consumed in these circumstances. For small Damköhler number a nearly linear decay of the cup-average concentration profile is observed, while the decay curves become sharper for the whole array domain at intermediate Damköhler number. For large Damköhler numbers, the concentration decays quickly in the front part of the array, and finally a full conversion of the reactant is achieved. The development of the cup-average concentration profile results from the interplay between reactivity and transport. The concentration distributions of these four Damköhler number cases at the particle Reynolds number of 240 for the central plane of the particle array are shown in Fig. 9, where the concentration difference resulted from the variable reaction rate of the first order irreversible chemical reaction proceeding at the particle surface could be clearly visualized. The computed concentration field inside the particles is zero due to the assumption of nonporous catalytic pellets with and only with reactive external surface. For the range of Reynolds numbers considered here (with Sc ¼ 1) axial dispersion is negligible compared with transport by convection. Therefore the differential equation describing the evolution of cup-average concentration is expected to the following expression (i.e.: 1D heterogeneous plug flow model): where k eff is the average overall effective mass transport coefficient defined by Eq. (46) and a s is the specific fluid-particle mass transfer surface given by: It should be noted that k eff is used and should be used here to consider not only the physical phenomenon of external mass transfer from the bulk fluid to the sphere surface but also the chemical process of species conversion due to the surface reaction. In this case, the corresponding driving force is ðc f À 0Þ. By integrating with the boundary condition at the inlet, c f ¼ c f ;0 at x ¼ 0, the following equation is obtained: In this equation, the Damköhler number is a predefined value, whereas the Sherwood number is unknown. The average Sherwood number of the full particle array can be computed from the slicebased Sherwood number, which describes the spatial distribution of the local Sherwood number in the surface perpendicular to the flow direction. This local Sherwood number in each slice is computed by using the following equation: where rc f Á n is the averaged concentration gradient at the sphere surface and c B is the averaged surface concentration for all particle sections located in the selected surface perpendicular to the flow direction. The cup-average concentration of this cross-sectional sur- face is used for calculating the local driving force, which is discussed earlier in this section (Eq. (50)). In Fig. 10, the distributions of the slice-based Sherwood number in the streamwise direction are presented, for four Damköhler numbers at Re p ¼ 160 and three Reynolds numbers at Da ¼ 1, respectively. For all cases, the slicebased Sherwood number is oscillating in the whole packed region and the average values are listed in Table 8. For the case of the same Reynolds number, larger oscillating amplitude is observed as reaction rate increases, and the profile rises slightly. If the Damköhler number is maintained, higher Reynolds number will shift the profile upwards and meanwhile increase the amplitude slightly. In Table 8, average Sherwood numbers estimated by the empirical Gunn correlation (Gunn, 1978) are also listed and compared with our DNS results.
Sh s ¼ ð7 À 10e þ 5e 2 Þ 1 þ 0:7Re where e is the void fraction of the particle array which is calculated as ð1 À gÞ, and Re s and Sc are the particle Reynolds number (based on the inlet fluid superficial velocity u 0 ) and the Schmidt number respectively. Both the average Sherwood number obtained from the slice-based computation and the Gunn correlation are used in the 1D plug flow model (Eq. (53)), which is compared with DNS results for the steady-state axial profiles of the cup-average concentration. As presented in Fig. 11, the comparison is performed for the complete packing section of the array, excluding the inlet and outlet   regions. The agreement between the simulation results and the 1D model profiles is quite reasonable, and the trend of the development of the profiles reaches a good agreement. The difference between DNS and the empirical model increases with higher Reynolds number, which could reasonably result from the inhomogeneous flow pattern inside the random packing at high Reynolds number.
Taking the advantages of all detailed information obtained from DNS, we can even go one level deeper, which is the Sherwood number of individual particles. The local particle Sherwood number is computed by using Eq. (43), with the driving force in the denominator redefined as the difference between the local cup-average concentration of the cross-sectional surface at the sphere center and the overall average surface concentration. The local cupaverage concentration is the same one we used in the calculation of the slice-based Sherwood number. It should be noted that the sphere center has high possibility to locate in a non-grid-point position due to the random Monte-Carlo distribution. In this case, a weighted averaging based on the neighboring two surfaces is used to give a value for the driving force calculation. In Fig. 12, the distribution of the individual particle Sherwood numbers is given for Re p ¼ 240 with increasing Damköhler numbers. The shape of the curve is a sharp log-normal distribution that becomes wider with increasing Damköhler numbers. The expectation value of each distribution is computed as 17. 08, 17.15, 17.67 and 18.40 for increasing Damköhler numbers, and the standard deviation is 5. 16, 5.25, 5.96 and 7.17 respectively. In addition, the distribution of the individual particle Sherwood numbers is also demonstrated for Da ¼ 1 with increasing Reynolds numbers in Fig. 13. Almost the same distribution pattern is observed, which shifts towards larger Sherwood number with the maximal probability maintained around 0.09. The expectation values of 12.45, 15.38 and 17.15, with corresponding standard deviation of 6.72, 5.87 and 5.25, are obtained for the case of Re p ¼ 80, Re p ¼ 160 and Re p ¼ 240 respectively. The decrease of the standard deviation is due to the reduction of the tails at higher Reynolds number. From both figures, we observe that a significant variation in the local particle Sherwood number exists in the particle array. This behavior can be reasonably Fig. 12. Distribution of the local particle Sherwood numbers for Rep ¼ 240 with Damköhler number of 0.1, 1, 10 and infinity. Fig. 13. Distribution of the local particle Sherwood numbers for Da ¼ 1 with particle Reynolds number of 80, 160 and 240. explained by flow features within the dense particle array where the flow distribution is strongly irregular. Stagnant zones may occur behind particles, whereas flow channeling may take place between particles. The local flow velocity is strongly influenced by the local porosity and higher velocities are observed at places of low porosity. As a result of the locally varying flow velocity, gas-solid mass transfer fluctuates strongly around individual particles.

Conclusions
In this paper a ghost-cell based immersed boundary method is presented for direct numerical simulation of fluid-solid mass transfer processes. In this method, a boundary condition is enforced at the exact position of the particle surface and incorporated into the fluid phase governing equations at the discrete level implicitly. Taking the advantage of a second order quadratic interpolation scheme utilized in the reconstruction procedures, the unique feature of our ghost-cell based immersed boundary method is its flexible application of the general mixed boundary condition. Assuming a first order irreversible surface reaction at the external surface of a sphere, the general mixed boundary condition could be modified into such a form that the Damköhler number can be controlled by specifying variable reaction rate coefficients. With our proposed DNS model, the complex heterogeneous systems involving the interplay between reactivity and transport encountered in industrial mass transfer processes can be handled with high efficiency.
The first simulation case is a detailed verification of the method. A single sphere with unsteady diffusion is considered, where a first order irreversible chemical reaction proceeds at the sphere surface. The species concentration profile and the development of the instantaneous Sherwood number obtained from DNS have a good agreement with the analytical solutions, for different Damköhler numbers. After that, forced convection to a stationary sphere is considered, and simulations are performed for two parameter sets. In one set the same reaction rate is considered while varying the Reynolds number. Both Sherwood number and species concentration at the sphere surface increase with higher Reynolds number. The other set applies variable reaction rates for the surface reaction with a constant Reynolds number. In this case, the Sherwood number is quite independent of the Damköhler number. The overall mass transport process goes from the reaction limited regime to the mass transfer limited regime. For both sets, the effective mass transport coefficients obtained from DNS agree well with the empirically calculated values. In the last simulation, a stationary array of particles is considered, for which the concentration distribution is computed for several particle Reynolds numbers and Damköhler numbers. The axial profiles of the cup-average concentration are compared with the 1D plug flow model and reach a good agreement. The slice-based Sherwood number and the particle-based Sherwood number are computed, and both demonstrate a considerable variation in the particle array. It is interesting to notice that the Sherwood number is influenced by the Damköhler number in this particle array case. With higher reaction rate, the mass transfer performance is improved to some extent. This is guessed as the consequence of a more heterogeneous concentration field at larger Damköhler numbers. @g @s À @ 2 g @x 2 ¼ 0