Direct Numerical Simulation of Reactive Fluid-Particle Systems Using an Immersed Boundary Method.

In this paper, direct numerical simulation (DNS) is performed to study coupled heat and mass-transfer problems in fluid-particle systems. On the particles, an exothermic surface reaction takes place. The heat and mass transport is coupled through the particle temperature, which offers a dynamic boundary condition for the thermal energy equation of the fluid phase. Following the case of the unsteady mass and heat diffusion in a large pool of static fluid, we consider a stationary spherical particle under forced convection. In both cases, the particle temperatures obtained from DNS show excellent agreement with established solutions. After that, we investigate the three-bead reactor, and finally a dense particle array composed of hundreds of particles distributed in a random fashion is studied. The concentration and temperature profiles are compared with a one-dimensional heterogeneous reactor model, and the heterogeneity inside the array is discussed.

Introduction

Fluid–particle systems are commonly encountered in a wide range of industrial applications, such as synthesis of chemicals by heterogeneous catalysis, combustion of pulverized coal, and coating on particle surfaces. Simultaneously with the mass transfer, these processes are often accompanied with significant heat effects, which introduces additional complexity to the system. Accurate predictions of the gas–solid interactions are of great help to improve process performance and facilitate equipment design. Therefore, it is important to understand the heat- and mass-transport processes in a coupled manner in such complex heterogeneous systems.

In the part decades, extensive experimental investigations have been conducted for fluid–particle systems, from which various correlations have been proposed for heat and mass transfer.[1−6] These correlations are very helpful to provide a quick estimation of the average heat and mass-transfer rates for engineering purposes; however, detailed information such as local variation and temporal development cannot be easily quantified. Three-dimensional transient simulation of complex multiphase flows has attracted considerable interest because of the fast development of computational capabilities, through which detailed quantitative information can be produced instead of an average value in heat- and mass-transfer processes. As the most detailed level of the multiscale modeling approach,[7] DNS is a powerful tool to resolve all the details at the smallest relevant length scales to gain fundamental insight in fluid–particle interactions and quantitatively derive closures for applications in more coarse-grained models. In recent years, the immersed boundary method (IBM), as a means to enforce the proper boundary conditions on immersed objects, has attracted a lot of attention. Taking the advantages of simple grid generation and efficient memory utilization, IBM has been successfully applied in different studies including moving particles, complex geometries, and deformable immersed objects.[8−13] After solving the momentum field, additional equations for species and thermal energy transport can be computed using the same methodology as the fluid flow equations.

Fluid–solid coupling schemes can be roughly categorized into two types for IBM: (1) continuous forcing method (CFM); (2) discrete forcing method (DFM). In CFM a Cartesian grid is used for fluid field simulation, whereas the immersed object is represented by Lagrangian marker points.[14−18] The interaction between the solid boundary and the fluid is accounted for by explicitly introducing a forcing term in the governing equations, which is the result of the distributed singular force over the Eulerian cells surrounding each Lagrangian point by a regularized Dirac delta function. In DFM the immersed boundary is treated as a sharp interface.[19−23] Rather than a forcing term determined by the feedback mechanism, the predefined boundary condition on the immersed object surface is enforced by extrapolating the surrounding fluid variables to the ghost cells (cells inside the immersed object but possessing at least one fluid cell neighbor). The virtual ghost value is then directly incorporated into the discretized governing equations of the fluid phase. DFM is also known as the ghost cell approach, and it is the one utilized in the current paper.

In contrast to numerous researches of interfacial momentum exchange in fluid–particle systems, much less IBM computational results are reported in the field of mass and heat transfer.[24−29] Only a few studies on reactive systems involving coupled mass and heat transfer processes have been reported. Due to the focus of the current paper on coupled heat and mass transport, only related literature are reviewed; however, for literature review of heat and mass-transport studies without such a coupling, we refer interested readers to our earlier papers.[30,31] Dierich et al.[32] applied 2D modeling to study the behavior of ice particles melting in the water, where the interfacial heat transfer as well as the phase-change phenomena play a significant role. The Dirichlet boundary condition is applied on the ice surface. Kedia et al.[33] introduced a second-order method for simulations of low-Mach number chemically reacting flow around heat-conducting immersed solid objects in combustors. A buffer zone is used to account for the fluid–solid interactions, where the conjugate heat transfer and the zero species penetration are applied. Deen and Kuipers[34] applied the directional quadratic interpolation scheme to study the coupled heat and mass transfer in a dense particle array, where they assume an infinitely fast chemical reaction proceeding at the particle surface. Also in the combustion field, Abdelsamie et al.[35] studied low-Mach number turbulent reacting two-phase flows in arbitrary geometries. In their work, CFM and DFM are implemented for moving particles and static boundaries respectively, and the source terms for species and thermal energy are directly incorporated in the governing equations. Dierich et al.[36] applied a variant of CFM, namely, a fixed-grid method, to study the simultaneous change of the solid interface and the solid porosity of a moving char particle during heterogeneous chemical reactions. The Boudouard reaction takes place at the outer surface as well as the inside of the porous particle, and the whole setup is two-dimensional and axi-symmetric. Luo et al.[37] developed a DFM-based IBM to study the combustion process of a single char particle. Different reconstruction schemes are applied to enforce the boundary conditions of different variables, and the effect of surface reactions on the boundary conditions is accounted for by balancing the mass and energy at the gas–solid interface. It should be mentioned that, although not using IBM, Dixon and co-workers have developed a particle-resolved CFD methodology to study methane stream reforming and ethylene partial oxidation processes using realistic microkinetics.[38−41] In their work, flow, heat, and species transport in fixed-bed reactors are modeled, as well as the intraparticle variations accounted for by the solid-particle method.[42]

Building on our previous study,[30,31,43] a DNS methodology which is based on an efficient ghost-cell based IBM is extended to the simulation of reactive fluid–particle systems. The coupling of heat and mass transfer arises as a consequence of an exothermic chemical reaction proceeding at the exterior surface of the particles that the particle temperature is increased by the liberated reaction heat and subsequently transfers the thermal energy to the fluid phase. As the Robin boundary condition is realized at the exact fluid–solid interface through a second-order quadratic interpolation scheme in our IBM, a surface reaction rate is incorporated to describe the interplay between chemical transformations and external mass-transport processes. For heat-transfer processes, the Robin boundary condition simply switches into the Dirichlet boundary condition with the solid temperature determined by the particle thermal energy equation, which serves as a dynamic boundary condition for the fluid thermal energy equation.

This paper is organized as follows. First, the DNS methodology is described, including governing equations, numerical solutions and fluid–solid coupling. Second, four reactive fluid–particle systems with increasing complexity are studied: unsteady diffusion around a single sphere, forced convection to a single sphere, the three-bead reactor and a dense particle array. Finally, the conclusions are provided.

DNS Methodology

In this part, we present the governing equations which need to be solved in DNS, the numerical details of the solution methods and the fluid–solid coupling. For the methodology described in the current paper, we assume the following main assumptions:

Both fluid and solid phase have constant physical properties.

The fluid phase is incompressible and Newtonian.

Diffusion is Fickian.

The solid phase is composed of spheres, on which external surface a chemical reaction proceeds.

The temperature gradients inside particles are negligible, and no heat effect on the reaction rate is considered.

Governing Equations

The following conservation equations for mass, momentum, species, and thermal energy are solved to describe the transport phenomena in the fluid phase:In above equations, ρf is the fluid density, μf is the fluid viscosity, Df is the species mass diffusivity in the fluid, and Cp, f and λf are the heat capacity and thermal conductivity of the fluid phase, respectively.

The particle temperature is governed by the following equation assuming of a uniform particle temperature:In this equation, Vs is the particle volume and Cp, s is the volumetric heat capacity of the solid phase. The first term on the right-hand side is the fluid–solid heat-transfer rate while the second term represents the rate of reaction heat liberated from a chemical reaction. The heat- and mass-transfer rates, with the normal unit vector pointing outward of the particle, are computed by the following two equations, respectively:It should be noted that in the present work, we consider a single heterogeneous reaction, and hence, eq is valid. In the case of multiple reactions, the second term on the right-hand side should be the summary of the heat effect from all reactions. Considering an exothermic chemical reaction proceeding at the exterior surface of the particles, the heat liberation is assumed to be rapidly transported to the interior of the particle with a negligible intraparticle temperature gradient. The coupling between the fluid thermal energy equation and the fluid species equation is fulfilled through the solid-phase thermal energy equation. In other words, the particle temperature of individual particle offers a dynamic boundary condition for the thermal energy equation of the fluid phase.

Numerical Solution Method

A finite difference scheme is used to solve the aforementioned governing equations on a 3D staggered Cartesian grid with a uniform grid spacing in all directions. Following our previous work,[30,43] the numerical solution of these equations is acquired by embedding second-order discretization schemes as well as compact computational stencils. The momentum, species, and thermal energy conservation equations are discretized temporally by applying the Adams–Bashforth scheme for the convective transport and the Euler backward scheme for the diffusive transport:In these equations, n is the time step index. The convection terms for momentum m, species Csp, and thermal energy Cth are, respectively, given by:and the diffusive momentum fluxes m, molar flux Dsp, and heat flux Dth are computed as:

Spatially, a second-order total variation diminishing scheme and a standard second-order central differencing scheme are applied for the discretization of the convection and diffusion terms, respectively. The solution of eq is obtained by applying a two-step projection method where an intermediate velocity field u̅** is first computed by using the pressure gradient at the old time step.[12,24] Subsequently, the velocity field is updated using the new pressure gradient computed from the Poisson equation at the new time step n + 1.

The governing equation of the solid phase is solved after the governing equations of the fluid phase, and the trapezoidal rule is applied for the time integration:

Fluid–Solid Coupling

After discretization, the differential equations are transformed into algebraic equations providing the relationship between the any fluid variable inside the simulation domain and its six neighbors. It should be noted that the six neighbors are the implicitly computed cells (i.e., diffusion term), while the convection term, which involves more neighbors, is explicitly accounted for. At the fluid–solid interface (i.e., the central point is the fluid phase but any neighbor point is the solid phase), IBM is invoked to incorporate the predefined boundary condition. The fluid–solid coupling is the most important element in our DNS methodology, which enforces the Robin boundary condition exactly at the immersed object surface and implicitly at the level of the discretized equations. The IBM method, which possesses a second-order accuracy, is published in our earlier paper,[30] and we refer the interested reader to that paper for a full description.

As indicated in the Introduction, the Robin boundary condition switches to the Dirichlet boundary condition straightforwardly for the calculation of the temperature field in the fluid phase, whereas the rate constant of the chemical reaction proceeding at the particle surface is incorporated into the Robin boundary condition for the mass-transfer calculation. The Damköhler number is used to quantify the ratio of the reaction rate to the diffusion rate:where rs is the particle radius. The original mass balance equation as well as the subsequent nondimensionalization procedure were also reported in the aforementioned literature.

Results and Discussion

In this part, four cases with increasing complexity are presented for coupled heat and mass transfer in fluid–particle systems, with an exothermic first-order irreversible chemical reaction proceeding at the particle surface. Following the comparison between DNS results and analytical (empirical) solutions for two single-particle systems, simulations are performed for systems with multiple particles.

Unsteady Heat and Mass Transport

Here the unsteady mass diffusion to a sphere, followed by the unsteady conduction of thermal energy to the surrounding fluid, is considered. We assume the spherical particle is located in an infinitely large pool of static fluid. The governing equations for unsteady mass and heat diffusion in the fluid phase are described by eq and 4, respectively, with convection terms set to zero. The boundary conditions are:at the boundaries of the simulation domain, andat the sphere surface. cf,0 and Tf,0 are the initial conditions for species and thermal energy equation, respectively, and eq ) and 20 are valid providing the simulation time is short enough to keep the diffusion fronts away from the confining walls. The particle temperature Ts in eq serves as a dynamic boundary condition for the fluid phase thermal energy equation, and governed by the particle thermal energy equation (Equation ).

In DNS, the sphere is positioned in the center of a cubic box with a length of 0.04 m. Five reaction rates varying from the case of reaction limiting to the case of mass-transfer limiting are imposed at the sphere surface, corresponding to the Damköhler number of 0.01, 0.1, 1, 10, and 100. The data used for the numerical simulation are given in Table .

Table 1

Parameter Settings for the Simulation of Unsteady Heat and Mass Transport

parameter	value
time step [s]	5 × 10–5
grid size [m]	2.5 × 10–4
sphere diameter [m]	0.005
fluid density [kg/m3]	1.0
fluid thermal conductivity [W/m/K]	0.025
fluid heat capacity [J/kg/K]	1000
species diffusivity [m2/s]	2 × 10–5
solid volumetric heat capacity [J/m3/K]	1000
reaction enthalpy [J/mol]	–10–5
fluid initial concentration [mol/m3]	1.0
fluid initial temperature [K]	293
particle initial temperature [K]	293

An “exact” solution for the particle temperature is obtained by solving the spherically symmetric model using a standard finite differencing technique with second-order accuracy. In this case, the governing equations for unsteady mass and heat diffusion in the fluid phase are only r dependence and, respectively, described as:It should be noted that a relatively high mesh resolution, defined as the number of grid points distributed in the radial direction, was used to obtain this “exact” solution.

The comparison of the temperature difference between the particle and the bulk fluid for the “exact” solutions and the DNS results at 3 s are listed in Table . From the table, a good agreement is observed for all reaction rates. In Figure , the particle temperature evolution profiles are plotted against time, to compare the DNS results and the “exact” solutions. As clearly demonstrated in Figure , these two solutions are in good agreement. For all reaction rates, the species flux is comparatively high initially, so that the heat liberated from the exothermic reaction rapidly heats the particle up from the initial temperature. After that, a temperature difference between the particle and the surrounding fluid is established, so that the particle transfers the thermal energy to the fluid phase through unsteady heat conduction. At the final stage, the heat removal rate is approaching the heat liberation rate, and therefore, a low reaction rate results in a low particle temperature. The final stage is achieved faster with higher reaction rates because of the dominating role of the Damköhler number. In other words, the solid temperature evolution is fully controlled by the unsteady diffusion, namely, the Sherwood number, for the Da = 0.01 case.

Table 2

Comparison between “Exact” Solutions and DNS Results for the Temperature Difference (in K) between the Particle and the Bulk Fluid Far from the Particle at 3 s

Da	0.01	0.1	1	10	100
“exact”	0.66	6.19	36.89	72.55	80.21
DNS	0.67	6.25	37.75	73.78	81.16

Figure 1

Comparison of particle temperature evolution profiles between the “exact” solutions and the DNS results, which are indicated by the solid lines and the dashed lines respectively.

Single Sphere under Forced Convection

Building on the last case, we now consider a single stationary sphere under forced convection. In lateral direction, the sphere is positioned at the center of the domain, whereas in the flow direction, it is located at a distance of two times of the sphere diameter from the inlet. Fluid flows into the system with constant inflow concentration of 1 mol/m3 and constant inflow temperature of 293 K. For this system, we consider 3 reaction rates corresponding to the Damköhler number of 0.1, 1, and 10. Besides the same data used in the last case, fluid viscosity is also required for the current numerical simulation which is specified to be 2 × 10–5 kg/m/s. The simulations are computed in a cubic box with 160 grids in each direction. At the inlet, uniform fluid velocity is imposed, which is varied to give the particle Reynolds number of 20, 40, 60, 100, 200, and 300. Free slip boundary conditions are applied at the transversal domain boundaries for velocity calculation, whereas the standard atmospheric pressure is set at the outlet. Homogeneous Neumann boundary condition is used for both concentration and thermal energy equations at all domain boundaries except the inlet.

According to Equation , the particle temperature at steady state is described by:The particle temperature can be calculated from the mass-transfer coefficient km and heat-transfer coefficient α.where Ss is the particle surface area. Equation can be rearranged to obtain the following expression:with the external mass- and heat-transport coefficient predicted by the well-known empirical Frössling and Ranz–Marshall correlations:where Res is the particle Reynolds number, Sc is the Schmidt number, and Pr is the Prandtl number, respectively, defined as:In eq , Le is the Lewis number defined as the ratio of the thermal diffusivity to the mass diffusivity:In the current work, Sc = 1 (and Pr = 0.8) is used, which means that the momentum and mass (temperature) boundary layers have a similar thickness.

In Table , the comparisons between the simulation results obtained from DNS and the empirical values calculated from eq are shown. As observed from the table, the particle temperature increases with higher reaction rates and decreases with larger Reynolds numbers. The former behavior is understood as the result of more heat liberated from the chemical reaction, whereas the latter behavior is well explained by the stronger convective heat transfer from the particle to the fluid. All results are in good agreement. At the final steady state, the calculated ratio of the heat and mass-transfer rates is also verified to equal the reaction enthalpy (−100 kJ/mol).

Table 3

Comparison of the Particle Temperature at Steady State between DNS Results and Empirical Values (Indicated as EMP)

	Da = 0.1		Da = 1		Da = 10
	DNS	EMP	DNS	EMP	DNS	EMP
Res = 20	296.7	296.4	319.7	318.0	363.6	360.6
Res = 40	296.0	295.8	315.7	314.5	360.6	358.1
Res = 60	295.6	295.5	313.4	312.5	358.3	356.2
Res = 100	295.1	295.1	310.5	309.9	354.7	353.4
Res = 200	294.6	294.6	306.9	306.6	348.8	348.7
Res = 300	294.4	294.4	305.0	304.8	344.9	345.6

Three-Bead Reactor

In this section, we consider three spheres positioned in a line, the so-called three-bead reactor. The spheres are positioned in a cuboidal domain with the length of 0.10 m in the flow direction and 0.01 m in the lateral direction. At the domain boundaries, velocity obeys the free slip condition and the system is isolated and adiabatic for species and thermal energy calculation. At the outlet, the standard atmospheric pressure is set, and zero slope boundary condition is set for both concentration and temperature. The first sphere is located at a distance of 2 times that of the sphere size from the inlet in the flow direction, and the other two spheres are located in such a way that the mutual distance between all sphere centers is 1.5 times that of the sphere size. In this simulation, the reaction rate at the sphere surface is maintained at Da = 1, indicating the equivalent time scale for reaction and diffusion, whereas 3 particle based Reynolds numbers Res = 60, 100, and 200 are applied to the system. All other simulation parameters are the same as those in the previous case.

The particle temperature evolution profiles are presented in Figure . There are two contributors to the rise of the particle temperature: liberated reaction heat and convective heat transfer. From the figure, one can observe that a thermal energy wave is propagating through the in-line array of three spheres. For all Reynolds numbers, the particle temperature increases from the first sphere to the third sphere. All spheres rapidly heat up due to the exothermic chemical reaction proceeding at the surface. Because of a temperature difference between the particle and the surrounding fluid, the thermal energy is transferred from the solid phase to the fluid phase and further transported downstream by the fluid flow. For the first sphere, the particle temperature soon reaches a constant value as the removed heat equals the generated reaction heat. For the second and third spheres, the unconverted reactant is partly converted on their surface and the thermal energy in the fluid phase gives additional heating to these spheres. Because of this reason, the last sphere takes the longest time to reach the steady state. From the figure, it is also clear that with higher Reynolds numbers, the final solid temperature is lower and achieved faster for all particles.

Figure 2

Evolution profiles of the particle temperature for the three-bead reactor. Particle Reynolds number of 60, 100, and 200 are presented by the solid lines, dashed lines, and dotted lines, respectively.

The evolution profiles of the fluid phase cup-average concentration and temperature in the streamwise direction are shown in Figure , which demonstrates the relative contribution of the individual spheres to the overall reactant conversion and temperature rise. The cup-average value ϕ (namely concentration or temperature in the current case) is defined by:In this equation, the integration is performed over the area occupied by the fluid in a cross section Sf perpendicular to the flow direction, and u is the axial component of the fluid velocity. From the figure, as expected, it is observed that three spheres have the same contribution to the overall species conversion and temperature rise. This is due to the low reaction rate specified at the sphere surface, and hence, all three spheres possess almost equivalent reaction consumption. Higher Reynolds number will decrease the species conversion at the sphere surface, and consequently lead to lower fluid temperature in the system. This corresponds well with the lower particle temperature in Figure .

Figure 3

Fluid phase cup-average concentration and temperature profiles along the flow direction. Three simulation cases with the particle Reynolds number 60, 100, and 200 are presented by the solid lines, dashed lines, and dotted lines, respectively, whereas the blue and red color indicate the concentration and temperature, respectively.

From the simulation, the overall conversion of the reactant in the three-bead reactor is obtained, which is 0.055, 0.035, and 0.019 for Res = 60, Res = 100, and Res = 200, respectively. The theoretical fluid outlet temperature can be computed from the adiabatic temperature rise:The fluid outlet temperatures given by the simulations are 298.69, 296.59, and 294.95 K for increasing Reynolds numbers, which are in good agreement with the theoretically calculated values: 298.51, 296.48, and 294.91 K.

Dense Particle Array

For the last case, DNS is performed to study a dense array composed of a relatively large number of stationary particles with the same size. The particle array is created by the hard-sphere Monte-Carlo method, with periodic boundary conditions in all three directions. In other words, the spheres that cross the boundaries are duplicated at the other side. In total, 573 spheres are distributed in a random configuration over a 3D domain with the packing height of 0.15 m in the flow direction and a length of 0.025 m in the lateral direction, with a predefined void fraction of 0.6. For the simulation three prescribed uniform fluid velocities are imposed at the inlet, leading to the particle Reynolds number of 60, 100, and 200, respectively. A reaction rate equaling the diffusion rate, namely, Da = 1, is specified at the particle surface. Other parameters used in the simulation are the same as those in the second case. The two parameters Schmidt number and Prandtl number, describing the relative thickness of the momentum boundary layer to the mass and thermal boundary layer, are 1.0 and 0.8, respectively. The periodic boundary condition is applied at the lateral domain boundaries for all velocity, concentration, and temperature calculations. The pressure at the outlet is set as the standard atmospheric pressure, and a zero slope boundary condition is set there for species and thermal energy equations.

In Figure , the computational domain as well as the particle configuration are demonstrated. It should be noted that two empty regions are reserved for the inlet and outlet, with a length of 0.015 and 0.05 m, respectively. These two regions are incorporated in order to avoid problems for the development of the inflow and the recirculation of the outflow. For our simulations, the mesh resolution N, defined as the ratio of the particle diameter to the grid size, is 20. This value is selected by performing a mesh convergence test in a small subarray, following the methodology published in our previous paper.[30] The deviation of the total heat-transfer rate in the current work using N = 20 is 4.14%, which is slightly higher than the one in our earlier work and can be explained by the increased solid phase packing density. The computed velocity field for all three Reynolds numbers are plotted in the longitudinal and transversal cross sections. The longitudinal cross section is the central plane of the particle array, whereas the three transversal cross sections are located at the packing height of 0.005, 0.075, and 0.145 m, respectively. In velocity maps, the periodic flow field in the lateral directions and the preferred flow pathways in the array are clearly observed. Note that in our simulations, the velocity field is solved not only for the fluid phase but also inside the particles; however, it is zero there due to the accurate enforcement of the no-slip boundary condition at the particle surface. The complexity of the flow structure is depicted by the transversal planes along the packing that the local increase of the fluid velocity significantly varies with the local porosity (caused by the variation of the local particle structure). With increased Reynolds numbers, besides a more heterogeneous flow field inside the array, a more unsteady wake is observed for the flow leaving the array which develops from straight streamlines into vortex. The periodically distributed particles, which are “cut” by the domain boundaries, are clearly visualized in Figure .

Figure 4

Particle configuration for the dense array simulation (left) and the computed velocity field (right) for the central plane and three lateral cross sections with the particle Reynolds number of 60, 100, and 200 (from left to right).

The evolution of the cup-average concentration and temperature along the particle array is of high interest for industrial applications, as it may indicate an approximate reactor length for the design purpose or an appropriate feedstock amount for the production process. The calculations of the cup-average concentration and temperature are given in the preceding section (eq ). In Figure , the cup-average concentration and temperature profiles of the fluid phase are shown as a function of the domain coordinate in the flow direction. The empty inlet and outlet regions are clearly visible as here the profiles are constant. Within the packing range, the species concentration decreases due to the chemical reaction proceeding at the external surface of the particles and the fluid temperature increases due to the particles that are heated up by the heat liberated from the exothermic reaction. As expected, higher Reynolds number will decrease the decay rate of the concentration and consequently lead to lower species conversion, which will accordingly reduce the temperature rise of the system. Similarly to the three-bead reactor, the computed fluid temperature at the outlet, 381.4, 370.2, and 352.8 K for the particle Reynolds number of 60, 100, and 200 respectively, can be compared with the value predicted by eq . For increasing Reynolds numbers, the overall conversion obtained from the simulations are 0.908, 0.796, and 0.620, and the sequential calculation gives the values of 383.8, 372.6, and 355.0 K respectively. The concentration and temperature distribution are shown in Figure for the central plane of the particle array. It should be noted that not only the information of the fluid phase but also the information of the solid phase is visualized in this figure. The concentration inside the particles is zero due to the assumption of the reactive external surface, whereas the temperature inside individual particles is its real solid temperature computed using eq . This is a very important point, as it may help us to find out where local high temperatures are located. Although the cup-averaged values in Figure show that the fluid temperature increases along the flow direction, from Figure it can be clearly observed that both fluid and particle temperature vary in the transversal cross section. High solid temperature occurs at regions with local dense particle configuration, especially for particles suffering from a blockage effect. The increased solid temperature consequently heats up the surrounding fluid and results in the nonuniform temperature distribution in the fluid phase. This heterogeneity is intensified at higher Reynolds numbers. In Figure , it is clearly observed that the system’s thermal energy continuously accumulates along the flow direction. In other words, the heat produced by the exothermic chemical reaction is stored in the particles and propagates to downstream particles by the convective heat transfer of the fluid.

Figure 5

Cup-average concentration and temperature profiles obtained from the dense particle array simulations. The case with the particle Reynolds number 60, 100, and 200 are, respectively, indicated by the solid, dashed, and dotted lines, whereas the concentration and temperature profiles are indicated by the blue and red color, respectively.

Figure 6

Concentration and temperature distribution for the longitudinal plane in the array center, with the particle Reynolds number of 60, 100, and 200 (from left to right). Note that the particle temperature is also indicated in the figure.

In more empirical approaches, the transport behavior of packed beds is commonly described by the classic one-dimensional heterogeneous model.[44] By incorporating the chemical reaction rate, the following equations are obtained for the fluid phase at steady state:where ε is the void fraction of the particle array and as is the specific fluid–particle contact surface area given by:Dax and λax are the axial dispersion coefficient and the axial conductivity respectively. The Damköhler number is already specified well in our simulations, whereas the Sherwood number can be computed from the empirical Gunn correlation:[6]In this equation, Res is the particle Reynolds number calculated with the inlet fluid superficial velocity uin. For more details of this one-dimensional heterogeneous model, we refer interested readers to the Supporting Information. In Figure , the comparisons of the cup-average concentration and temperature profiles between the DNS results and the phenomenological model are shown. As presented by the figure, good agreements are reached for the overall features of the profiles. The important discrepancy is that the DNS results give a faster and higher species consumption and consequently a similar trend for thermal energy than the ones predicted by the 1D model. This discrepancy further increases at higher Reynolds numbers, and the reason is thought to be the inhomogeneous flow pattern inside the particle array. In the figure, we also plot the 1D profiles without the axial dispersion effect (for both species and thermal energy transport). As expected, the axial dispersion plays a more pronounced role at low Reynolds numbers, and it is negligible for the range of the Reynolds number considered in the current work due to the dominating role of the convective transport.

Figure 7

Comparisons of the fluid cup average concentration (upper) and temperature (lower) profiles among the DNS results, the 1D model using the Sherwood number computed from the Gunn correlation and the 1D model without the axial dispersion effect using the same Sherwood number.

The average Sherwood number obtained from the slice-based computation is as well used in the 1D model. It is calculated by the following expression in planes perpendicular to the flow direction:where the numerator and denominator account for the average concentration gradient and the average driving force respectively, for all parts of the particle surface within the current plane. As shown in Figure , the slice-based Sherwood number depicts the local variation of the mass-transfer performance in the streamwise direction, which oscillates in the full packing region due to the varying fluid–particle interface area in each transversal plane. From the figure, it can be concluded that the Sherwood number is statistically homogeneous because it is independent of the streamwise coordinate except for a very small slab at the end of the array. The Sherwood numbers within this slab significantly deviate from the homogeneous value and hence is excluded in the average value calculation. In Table , the average Sherwood number is compared with the value obtained from the Gunn correlation. As predicted by the correlation, the Sherwood number increases with higher Reynolds numbers, however, accompanied by larger deviations. A similar behavior has been reported by other studies,[27,29,45], and the values obtained from their refitted correlations are also listed in the table. The utilization of the average value does not help to improve the prediction of the 1D model, as depicted by Figure . Because of the lower values applied in the 1D models, the concentration profiles are somewhat higher than the ones computed using the Gunn correlation and the temperatures lower.

Figure 8

Streamwise profile of the slice-based Sherwood number for all three Reynolds numbers.

Table 4

Comparison of the Average Sherwood Number between the Gunn Correlation and the DNS Results

Res	Gunn	Tavassoli[27]	Sun[29]	slice-based	particle-based
60	12.90	10.99	10.30	11.65	14.47
100	15.81	14.30	13.03	14.13	15.73
200	21.59	-	-	17.69	18.60

Figure 9

Comparisons of the fluid cup-average concentration (upper) and temperature (lower) profiles among the DNS results, the 1D model using the Sherwood number based on the slice-based average and the 1D model using the Sherwood number computed from the Gunn correlation.

The previous results indicate that the variable flow field inside the array, which results from the absence of a homogeneous packing structure, has considerable influence on the local transport performance of individual particles. Monitoring the variation from particle to particle may help us to gain deeper understanding of the heterogeneity. For this purpose, the Sherwood number of individual particle is computed, which is defined as:The numerator is the average concentration gradient, whereas the denominator is the local driving force for individual particle which is calculated as the difference between the cup-average concentration at the center position of the particle and the average surface concentration of the particle. There are two advantages by applying this definition. On one hand, it is closely related to the slice-based calculation, as the average over particles around a certain position should be close to the Shslice there. On the other hand, rather than the overall value obtained from the slice-based calculation, this particle-based quantity is commonly used in the coarser scale models to characterize the interfacial transfer. In Figure , the values of the computed Sherwood number of individual particles are plotted as a function of the streamwise (x) and spanwise (z is shown here) directions. As wall effects are eliminated by applying the periodic boundary condition at the spanwise boundaries, no entrance region is observed for Shparticle, i in the flow direction (left panel), and no predominant high values are observed at the lateral boundaries (right panel). The average value over all particles are listed in Table as well, and the standard deviation decreases with increasing Reynolds number which is 12.14, 9.21, and 8.13, respectively. From the figure, it seems that the Sherwood number of individual particles has a quite uniform and homogeneous distribution in all directions. However, some details are worth being discussed. First, there is a decrease of the Sherwood number at high Reynolds number. Although most particles behave in the same way as the overall Sherwood number, namely, that the value increases with larger Reynolds numbers, some particles have a completely opposite behavior, namely, that the maximal value occurs at Res = 60. This indicates the preferred flow path inside the array that the fluid channels somewhere and leads to stagnation in other places, especially at high Reynolds number. Second, variation of the Sherwood number in the same transversal cross section. The local heterogeneity can be more pronouncedly revealed by looking at the particles at the same position along the flow direction, which can even differ up to 100%. Since Shparticle, i is defined as the ratio of the concentration gradient to the mass-transfer driving force, it can be concluded that at least one of them is not constant at the same cross section perpendicular to the flow. This is in contrast to the assumptions employed in the 1D model that both gradient and driving force mainly change in the flow direction, and this may explain the deviation in our previous results. In Figure , the concentration gradient and the concentration driving force are plotted at the particle Reynolds number of 200 in the streamwise direction. As expected, both gradient and driving force decrease along the particle array. In the figure, one can clearly find that the particles in the same transversal plane experience some differences for the concentration gradient; however, differences for the concentration driving force are more significant. To confirm our finding, the gradient and the driving force for thermal energy transfer is also plotted in Figure , which show the same behavior. This well explains the variations in the previous results and may indicate the improper usage of the cup-average concentration for the driving force calculation. In other words, particles in the same spanwise plane may not feel the same driving force; instead, they are more likely affected by a local driving force which is based on the local microstructure around each particle. However, this may trigger further questions on how to derive an overall value from individual values, which can be then applied in phenomenological models. This raises a challenge for future researches to analyze really local and detailed information, which is realizable by using DNS.

Figure 10

Distribution of the Sherwood number of individual particles in the flow (upper) and lateral (lower) directions.

Figure 11

Distribution of the mass (upper) and thermal energy (lower) gradient and the driving force of individual particles along the flow direction, at Res = 200.

Conclusions and Outlook

Direct numerical simulations based on a ghost-cell based immersed boundary method are performed for reactive fluid–particle systems in this paper, namely, coupled heat and mass-transfer processes. An exothermic first-order irreversible chemical reaction proceeds at the exterior surface of the particle, whose reaction rate is incorporated into the Robin boundary condition which is further incorporated into the fluid species equation at the discrete level implicitly. The liberated reaction heat heats up the particle, whose temperature offers a dynamic boundary condition for the fluid thermal energy equation.

Four fluid–solid systems are studied. The particle temperature obtained from DNS agrees with the value obtained from the well-established solutions for a single sphere under both unsteady and convective situation. In the three-bead reactor, the particle temperature increases from the first sphere to the third sphere and the species conversion decreases at higher Reynolds numbers. The adiabatic temperature rise obtained from DNS is in a good agreement with the theoretical value. For the dense particle array, the fluid concentration and temperature profiles are obtained and compared with the 1D heterogeneous model. The qualitative agreement is good; however, the quantitative discrepancy is thought to be caused by the heterogeneity inside the array which is further amplified at higher Reynolds numbers. This point is visualized in both velocity and temperature distributions. Further analysis of the Sherwood number of individual particles reveals that the particles in the same transversal plane experience large variations, which mainly come from the widely scattered driving force.

In this paper, we revealed the strong power of DNS in modeling of reactive fluid–particle systems for engineering applications. Future work can be expanded in two aspects. First, the current methodology needs to be extended toward more realistic systems, which considers the temperature dependent reaction rate and incorporates the realistic reaction kinetics. Second, the local transport behavior of individual particles needs to be studied, from which a better estimate of the transfer performance is expected by quantitatively understanding the effect of local parameters.