CFD Modeling and Simulation of the Axial Dispersion Characteristics of a Fixed-Bed Reactor.

In this study, the axial dispersion characteristics of a fixed-bed reactor with different packed structures were investigated via computational fluid dynamics (CFD) simulation. The discrete element method was employed to develop the physical model of a fixed bed. Then, CFD simulations were performed to investigate the flow resistance coefficient under different Reynolds numbers. The prediction values were in fair agreement with those calculated by the Carman equation, thereby validating the proposed CFD model. The tracer pulse method and the step method were employed to evaluate the residence time distribution characteristics in the fixed-bed reactors where the mean residence time and axial dispersion coefficient were calculated. The distribution characteristics of the tracer concentration and fluid velocity were also obtained and used to explain the mixing performance of the fixed bed. This simulation study can contribute to the optimization design and scaling up of reactors with porous packed structures.

Introduction

Fixed-bed reactors have been widely used in industrial production processes, including adsorption,[1,2] catalytic oxidation,[3] methane reforming,[4,5] and wastewater treatment.[6] Among these processes, the axial dispersion behaviors of species play an important role in determining the final reactor performance.[7,8] Owing to the complicated packed structure, however, it is challenging to obtain an in-depth understanding of these characteristics. Further investigations of the hydrodynamic characteristics of fixed-bed reactors based on a more effective research method are still required for their optimal design.

In earlier times, the two-dimensional pseudohomogeneous model was used to design fixed-bed reactors and a significant amount of attention was focused on global parameters such as bed porosity, pressure drop, and drag force.[9−11] However, the complicated flow structures in porous media regions were usually unclear. Subsequently, the flow patterns in fixed beds were visualized using advanced experimental technologies such as laser Doppler anemometer (LDA),[12,13] particle image velocimetry (PIV),[14,15] and magnetic resonance imaging (MRI).[16,17] The inertial flow structures were discussed in detail. Several valuable experimental data have been reported. These attempts and efforts provided a systematic research method for fixed-bed reactors and indicate the direction for their optimization design. However, it is expensive and almost impossible to examine all of the key influencing factors solely via experiments.

Computational fluid dynamics (CFD) simulations have recently been employed as a powerful tool for characterizing the flow structures in fixed-bed reactors. The discrete element method (DEM) has been widely used to reconstruct three-dimensional (3D) packed beds in recent years, providing detailed flow structures in the interstices between particles, such as velocity vectors and vortices, as well as global parameters.[18−20] In this respect, previous studies focused on investigating the pressure drop and drag coefficient,[21−23] evaluating the dependence between the accuracy of CFD results and the number of particles present,[24,25] and studying the effect of particle properties (i.e., particle shape and size distribution) on the flow patterns.[26,27] Knowledge of such detailed flow structures in the interstices can have significant implications for the estimation of the local heat and mass transfer rate.

There were other studies focusing on the mixing performance of fixed-bed reactors. Maier et al. simulated the velocity distributions in a column of glass beads using the lattice Boltzmann method, and they found that the longitudinal diffusion had a significant effect on the axial velocity.[28] Later, they simulated the tracer dispersion in the pore scale of regular and random spheres based on a random-walk particle-tracking method.[29] They also investigated the hydrodynamic dispersion in open-cell polymer foam and confined packed beds; their simulation results agreed well with the nuclear magnetic resonance spectroscopy experimental data.[30,31] Gutsche and Bunke developed a hydrodynamic model to predict the axial dispersion and external mass transfer in a fixed bed.[32] The proposed hydrodynamic model exhibited universality, simplicity, and good prediction capabilities. However, the model can potentially be only suitable for the creeping flow regime. Lima and Zaiat observed that there was an optimal value for the degree of back-mixing for hydrogen production in anaerobic fixed-bed reactors.[33] Dixon and Medeiros investigated the wall effects on the gas-phase radial dispersion in fixed beds via CFD simulations.[34] They reported that a three-parameter model based on a mixing length concept, where the transverse dispersion coefficient varies with the radial position, can significantly improve the ability to reproduce the near-wall concentration profiles. Mondal et al. investigated the mixing characteristics of a serpentine millichannel-based packed-bed device.[35] They examined the effect of the flow rate on the residence time distribution and proposed new correlations for the frictional resistance and axial dispersion. More recently, Petrazzuoli et al. predicted axial Péclet numbers using direct numerical simulations under single-phase laminar flow conditions. The proposed workflow produced faster data than an experimental approach.[36] Generally, the packed structure is dependent on the particle shape and size distribution, which in turn determines the flow characteristics and mass transfer rate. Moreover, it is still a challenge to investigate dispersion effects in fixed-bed reactors at high Re by means of experiments.

In this study, we aim to evaluate the axial dispersion characteristics of a fixed-bed reactor with different packed structures. Physical models of the fixed-bed reactors were developed using the DEM method. Then, CFD simulations were performed to investigate the flow resistance coefficients, which were compared with those calculated by the Carman and Ergun equations. The tracer pulse method and step method were employed to evaluate the residence time distribution characteristics in a fixed-bed reactor. The mean residence time and axial dispersion coefficient were calculated based on the residence time distribution density curves. Much attention was focused on the axial dispersion characteristics at high inlet velocity (2 m/s) with the local Re of up to 4100. Furthermore, the distribution characteristics of the tracer concentration and fluid velocity were obtained and used to explain the mixing performance of the fixed bed.

Fixed-Bed Reactor Mesh Model

The diameter and height of the fixed-bed reactor corresponded to 60 and 120 mm, respectively. The fixed-bed reactor was randomly filled with spherical particles of different size distributions. The mean particle size was set to 5 mm, and the standard deviations were set to 0.25, 0.75, 1.5, and 2.5. Then, particles with a normal size distribution were generated. Additionally, a fixed bed packed with a constant particle size was developed for the control study. In this study, these packing structures were modeled using the commercial EDEM software. Spherical particles were generated at a specific height (H = 120 mm) and allowed to naturally fall under the action of gravity to simulate the real particle filling process. The physical parameters of granular material include Poisson’s Ratio (0. 25), solids density (2400 kg/m3), and shear modulus (1 × 107 Pa). The contact parameters include the coefficient of restitution (0.5), coefficient of static friction (0.5), and coefficient of rolling friction (0.01). Generally, we only cared about the final filling structure (i.e., bed voidage) but ignored the filling process, which only took 0.2 s. About 6000 particles were inserted every second. The particles filling height was 60 mm, and the number of filled particles varied for different size distributions. Figure shows the particle size distributions for the different fixed beds generated by DEM simulations.

Figure 1

Fixed-bed physical model and size distribution of generated particles via EDEM.

Then, the commercial software ANSYS ICEM was employed to mesh the developed 3D-packed models. An O type segmentation method was adopted to better solve the mesh distortion at the vertex of the block of arc or other complex shapes while generating an ideal boundary layer mesh near the wall surface. Tetrahedral elements were used to generate an unstructured mesh in the porous media region. To minimize backflow effects, an empty bed of 30 mm was added at the entrance and exit of the fixed bed. Namely, the fixed-bed reactor had a height of 120 mm, but the height of the packed region was only 60 mm. Moreover, the packed region was located in the middle of the fixed bed. The two empty bed regions were discretized via a hexahedral mesh with a maximum mesh element of 0.8 mm. Finally, the three parts of the mesh model were assembled and used for the CFD simulations.

CFD Modeling and the Numerical Simulation Method

CFD Modeling

Based on the developed physical mesh model, CFD simulations were conducted to solve the flow fields of velocity, pressure, and tracer concentration. The detailed governing equations included the continuity, momentum, and species balance equations. The laminar, transition, and turbulent flow regimes were often encountered in fixed-bed reactors. Therefore, the laminar model (Reα < 10), the k–ω turbulent model with low Reynolds number corrections (10 < Reα < 300), and the RNG k-ε model (Reα > 300) were used according to the specific Reynolds number. The governing equations are expressed as followswhere p is the pressure, τ is the stress tensor, and ρ and are the gravitational body force and external body forces, respectively. In eq , μ is the molecular viscosity and μt is the turbulent viscosity, which is computed by the turbulent model.

In this study, CFD simulations were performed to investigate the flow fields in a fixed-bed reactor at different inlet velocities, which ranged from 0.01 to 5 m/s. According to the calculated Re, three flow regimes, such as the laminar flow, transitional flow, and turbulent flow regimes, were involved in this study. During the CFD simulations, when the flow rate was very slow (i.e., v = 0.01 m/s), the laminar flow model was used. With the increase in flow rate, we found that the k–ω model with low Re corrections could help the CFD simulation to converge quickly. With the further increase in the flow rate, the k–ω model failed and the RNG k–ε model was successfully used.

The k–ω turbulent modelwhere the effective diffusivities are given byThen, the turbulent viscosity was computedWhen a low Reynolds number correction was consideredwhereThe RNG k–ε modelFor the high Reynolds number, the turbulent viscosity is computed asWhen the low Reynolds number effects are considered, the differential relation is employed to calculate the turbulent viscositywhere

In this study, the dispersion of tracer in a fixed-bed reactor was considered. The species transport equation is as followswhere Y is the mass fraction of tracer, D is the diffusion coefficient of tracer, Sc is the source term, which equals to zero in this study, and DT is the turbulence diffusivity, which was evaluated from the turbulent Schmidt numberIn the Ansys Fluent, the default ScT is 0.7. Note that turbulent diffusion generally overwhelms laminar diffusion, and the specification of detailed laminar diffusion properties in turbulent flows is generally not necessary.

Numerical Simulation Details

In this study, a single-phase two-component system (air and tracer, ρ = 1.225 kg/m3, μ = 1.7894 × 10–5 Pa/s–1, D = 7.84 × 10–5 m2/s) was selected to investigate the effects of the inlet velocity and particle size distribution on the mixing performance of a fixed-bed reactor. With the exception of the velocity inlet and pressure-outlet boundary conditions, all of the other boundary conditions were set to “wall.” Namely,

for the inlet boundary (z = 0 mm)for the outlet boundary (z = 120 mm)for all of the wall boundaries

First, the steady-state CFD simulation was performed and convergence was reached when the scaled residuals for each transport equation were less than 1 × 10–4. Then, the tracer was introduced, and the steady-state CFD simulation was switched to transient simulation. In the tracer pulse method, the tracer (50 wt %) was patched at the entrance of the first layer of particles because we mainly focused on the mixing characteristics in the packed region. Generally, the patch height of the tracer should be small enough. The selected patch height (5 mm) was equal to the particle size. The residence time density function was obtained by recording the mass fraction of the tracer at the exit of the packed bed. In the step method, the fluid region below the packed bed was patched with the tracer (50 wt %), and the inlet boundary was also switched to the mixture feed (50% tracer, wt). In this manner, the residence time distribution function was obtained by recording the mass fraction of the tracer at the exit. When the residence time density function and the residence time distribution function were obtained, the transient CFD simulations manually stopped. All CFD simulations were executed on a 2.2 GHz Intel 2 Central Processing Unit (16 cores) with 128 GB RAM, and about two days were required for each simulation.

Residence Time Distribution Analysis

Based on the time-dependent curves of the tracer concentration, the density distribution function is defined asThen, the mean residence time (tm) and the variance (σ2) can be calculated in turnA lower variance value indicates a narrow distribution. Based on the calculated mean residence time and variance, the Péclet number can be available. The relation between the Péclet number and the variance for the open–open boundary conditions is expressed as follows[37]Finally, the axial dispersion coefficient can be calculated based on the Péclet numberwhere uave is the average flow rate in the porous media region and dp is the particle size.

Results and Discussion

Study of Grid Independence and Model Validation

For the CFD simulation, grid accuracy plays a significant role in determining the predicted results. In this study, grid independence analysis is performed by comparing the velocity profiles. In Figure , the velocity is the area average velocity. We select 20 cross sections and calculate their average velocity. As shown, when the maximum grid size is less than 0.6 mm, the grid accuracy seems to have little effect on fluid velocity. There is reasonable consistency between the results with a grid size of 0.4 and 0.6 mm. Although the agreement may not be good enough to confirm the resolution is sufficient, it is likely to be close. Generally, the fluid velocity is significantly affected by the bed porosity. Lower bed porosity is usually responsible for the higher velocity. Therefore, a verification study about the voidage is also necessary when the maximum grid size is 0.6 mm. There are two ways to determine the bed voidage. Based on the physical modeling process, the volume of all of the packed particles can be counted. Thus, the true value of the bed voidage is available. Additionally, the bed voidage can be obtained by calculating the volume of the fluid in the porosity media region via ANSYS Fluent. Table displays the voidage for each fixed-bed reactor when the maximum grid size is 0.6 mm. As shown, the simulated voidages agree well with the true values. From Table , it is observed that the bed voidages for each fixed bed with different particle size distributions are basically the same when the average particle size remain unchanged.

Figure 2

Grid independency analysis: the velocity distribution in the axial direction of the fixed bed packed with spherical particles of different size distributions. (A) Constant size; (B) σ = 1.5.

Table 1

Packed Structure Parameters and Voidages for Each Fixed-Bed Reactor when the Maximum Grid Size Was 0.6 mm

standard deviation	particle numbers	voidage (true value)	voidage (CFD)
0	1500	0.4575	0.4612
0.25	1570	0.4294	0.4330
1.5	1580	0.4234	0.4250
2.5	1590	0.4206	0.4232

Subsequently, the selected meshes were employed to simulate the pressure drop at different inlet velocities. In fixed-bed reactors, the well-known pressure drop correlation was defined aswhere f is the flow resistance coefficient. The classical Ergun[11] and Carman[38] resistance coefficients are defined as eqs and 28, respectively.where the particle Reynolds number (Reα) is defined as follows

The simulated pressure drops were then employed to calculate the flow resistance coefficient. As shown in Figure , cases with constant size and size distribution were considered. When the packed particles exhibit a size distribution, the small particles could be filled between large particles, and this led to a more compact bed. We found that when the average particle size remained unchanged, a wider particle size distribution was responsible for a smaller voidage, which decreased from 0.4575 to 0.4206. Hence, the simulated flow resistance coefficients met higher values. In this study, the investigated Reynolds numbers ranged from 6 to 4000. According to the comparison results, it was found that the simulated flow resistance coefficients were in fair agreement with those calculated by the Carman correlation over the entire range of Reynolds numbers, including laminar and turbulent flow regimes. However, the Ergun correlation overestimated the flow resistance coefficients at high Reynolds numbers. Similar results have been reported by many other researchers.[39,40] In this study, a simple single-phase flow model was employed to solve the flow fields. The rationality of the CFD model had already been verified elsewhere.[41,42] The fair agreement between simulation results and calculated data validated the proposed CFD model.

Figure 3

Comparison results of CFD results and classical correlations in terms of flow resistance coefficient under different Reynolds numbers.

The axial pressure distribution and contours of the static pressure in fixed-bed reactors are shown in Figure . Although the local pressure can fluctuate owing to the change in bed porosity, the facet average pressure is almost linear along the bed height. The linear scaling rule for the pressure drop is commonly used for an actual fixed-bed reactor. Furthermore, there are no entrance and exit effects because constant pressures are observed.

Figure 4

Axial pressure distribution and contours of static pressure in fixed-bed reactors (inlet velocity = 2 m/s).

In summary, the effectiveness of the physical mesh models was confirmed. The accuracy of the model was further verified based on the comparison results for the flow resistance coefficient. Further CFD simulations were conducted to evaluate the mixing performance of fixed-bed reactors.

Residence Time Distribution Characteristics

The mixing performance of fixed-bed reactors was investigated in this section. The tracer pulse method and step method were used to determine the residence time distribution via CFD simulations. The tracer with a mass fraction of 0.5 was initially patched at the entrance of the packed structure. Then, transient simulation was performed, and the mass fraction of the tracer at the exit of the packed structure was monitored and recorded. Figure shows the residence time distribution density function (5A and 5B) and the residence time distribution function (5C and 5D) for two different packed structures. It was observed that the mean residence time decreased when the particle size distribution was considered. This phenomenon was determined by the bed structure. Typically, the increased flow velocity was responsible for the decreased mean residence time. As discussed above, more small particles were filled between large particles, which in turn led to a decrease in bed porosity and an increase in the flow rate. In this study, the height of the packed bed was 60 mm. When the inlet velocity was 0.01 m/s, the mean residence time was approximately 3.7 s, which decreased to 0.017 s when the inlet velocity increased to 2 m/s. When the bed porosity was 0.425 and inlet velocity was 0.01 m/s, the residence time of ideal plug flow behavior was 2.55 s, which was lower than that calculated via CFD simulation (3.7 s). This indicated that the actual fluid flow in the packed region significantly deviated from the ideal plug flow. The stagnant zones might be responsible for the increased residence time. There were many tracer concentration enrichment zones in the packed region, delaying the outflow of tracer (Figures and 8).

Figure 5

Residence time distribution density function and residence time distribution function under two different packed structures: (A) and (C) inlet velocity is 0.01 m/s and (B) and (D) inlet velocity is 2 m/s.

Figure 7

Contour of tracer concentration at different times when the inlet velocity is 2 m/s.

Figure 8

Radial distribution of tracer concentration at different bed heights when the inlet velocity is 2 m/s and flow time is 0.005 s.

Additionally, the bed structure appeared to affect the half-peak width of the residence time distribution density function in a limited manner. This effect increased slightly as the inlet velocity increased from 0.01 to 2 m/s. Mondal et al. reported that the tracer residence time increased at a lower flow rate, which also broadened the residence time distribution curve.[35] This was due to the fact that the diffusive process dominated at a lower flow rate. Furthermore, the residence time distribution density function curve did not exhibit a trailing tail, which implied that there was almost no dead zone in the packed region.

Based on the measured density function curves, the mean residence time, variance, and axial dispersion coefficient for different inlet velocities were calculated using the user-defined MATLAB program. It should be noted that small differences in the residence time distribution parameters were observed for two different bed structures. Therefore, Table only displays the residence time distribution parameters for the particles with size distribution (σ = 1.5).

Table 2

Mean Residence Time and the Axial Dispersion Coefficient in the Fixed-Bed for Different Inlet Velocities

inlet velocity (m/s)	average velocity in packed region (m/s)	particle Reynolds number, Re	mean residence time, tm (s)	variance σ2 (s2)	axial dispersion coefficient, Dax, (m2/s)
0.01	0.027	9.3	3.66	1.04	4.67 × 10–6
0.05	0.14	46.6	0.71	5.11 × 10–2	2.97 × 10–5
0.1	0.27	93.4	0.41	3.12 × 10–2	9.69 × 10–5
0.5	1.39	477.4	0.068	2.51 × 10–4	1.71 × 10–4
1	2.79	954.1	0.034	5.93 × 10–5	3.32 × 10–4
2	5.55	1901.0	0.017	4.40 × 10–5	1.77 × 10–3

The effect of the inlet velocity on the mean residence time is shown in Figure A. At the lower inlet velocity, the tracer residence time and variance increased, thereby indicating a broadening of the residence time distribution curve. The correlation was determined via least-squares fitting as followsThere was an almost inverse proportional relationship between the mean residence time and inlet velocity. Similarly, Mondal et al.[35] investigated the mean residence time in the millichannel-based fixed-bed device for different flow rates. Based on their experimental data, the relationship between mean residence time and flow rate was

Figure 6

Dependency of the residence time on the inlet velocity (A) and the axial dispersion coefficient on the particle Reynolds number (B) for the two different bed structures.

As shown, eq is similar to eq obtained by us. The difference in parameters may be caused by different operating conditions. In Mondal et al.’s experiments, the velocity ranges from 5.8 × 10–4 to 2.67 × 10–3 m/s, which are much smaller than those used in our study.

The axial dispersion coefficient was calculated and is shown in Figure B, where it increased as the particle Reynolds number increased. Typically, dispersion dominated diffusion at a high Reynolds number. In this study, the highest axial dispersion coefficient was obtained as 1.77 × 10–3 m2/s at a higher inlet velocity (2 m/s). However, for the low Reynolds number (Re < 100), the axial dispersion coefficient was less than 9.69 × 10–5 m2/s, thereby resulting in a dominating diffusive process. Furthermore, the same order of magnitude of the axial dispersion coefficient was observed by Mondal et al., who investigated the hydrodynamics and mixing characteristics of a millichannel-based serpentine fixed-bed device.[35] The axial dispersion coefficient also varied for the two-bed structures. A higher axial dispersion coefficient was observed when the particle size distribution was considered. The deviation in the axial dispersion coefficient continued to increase as the particle Reynolds number increased.

Flow Field Characteristics

As discussed above, the residence time distribution characteristics of fixed-bed reactors were determined by their hydrodynamic performance. Therefore, it was necessary to provide a more detailed discussion of the flow field in fixed-bed reactors for an in-depth understanding of the mixing performance. Figure shows the contour of the tracer concentration at different times when the inlet velocity was 2 m/s. As shown, a nonideal flow phenomenon that deviated from slug flow was observed in the fixed-bed reactor. At the initial time, the maximum concentration of tracer was 0.5, and it decreased to 0.26 after 0.005 s. Additionally, there was almost no dead zone in the fixed-bed reactor. This result had already been discussed previously. This could be reflected by the residence time distribution density function curve, which did not exhibit a trailing tail.

Additionally, the radial distribution of the concentration is extremely important for characterizing the mixing performance of fixed-bed reactors. Figure displays the radial distribution of the tracer concentration at different bed heights when the inlet velocity was 2 m/s and the flow time was 0.005 s. In a fixed-bed reactor, the particles are filled to form a porous medium region where the flow channel is curved. When the fluid flows through the porous medium area, it is continuously dispersed, just like through multiple layers of distributors, which leads to the radial distribution of the fluid in the reactor. This implies that radial concentration distribution is inevitable, and this phenomenon can be explained by the uneven velocity distribution in the fixed-bed reactor. As shown in Figure , the radial distribution of tracer concentration in the fixed-bed reactor is uneven and a large concentration gradient in the local position is observed. This concentration distribution is typically determined by the velocity distribution. The low local flow rate led to tracer enrichment.

Figure displays the axial and radial velocity distributions in the fixed-bed reactor when the inlet velocity is 2 m/s. It is evident that a large concentration gradient is usually located in the low-velocity regions. When the flow velocity decreases to zero, a stagnant zone is formed. In Figure , a higher flow velocity is observed at the narrow clearance between particles. However, wall effects only appear at the upper left of the fixed-bed reactor. In this study, the tube-to-particle diameter ratio (D/d) is higher than 10. This leads to a uniform fixed bed, especially when the packed particles exhibit a size distribution. As shown in Figure , many small particles are packed close to the wall of the reactor, thus weakening the wall effect. Many researchers also applied the criteria of D/d > 10 for neglecting wall effects.[21,43−45]

Figure 9

Axial and radial velocity distributions in the fixed-bed reactor when the inlet velocity is 2 m/s.

Figure shows the contours of different velocities in the two different bed structures. As shown, although the inlet velocity was 2 m/s, the flow velocity in the clearance between the particles was more than 12 m/s. When the particle size distribution was considered, the fixed bed was more compact, and this yielded a higher flow velocity. Hence, the mean residence time decreased. Generally, fluid velocity played a significant role in determining the local mass and heat transfer rate. Therefore, further studies are required to investigate the effect of the packed structure on the hydrodynamic performance of fixed-bed reactors.

Figure 10

Contour of different velocities in two different bed structures.

Conclusions

In this study, the DEM was employed to develop different packed structures, which were then used for CFD simulations for investigating the axial dispersion characteristics of a fixed-bed reactor. The accuracies of the proposed model and the numerical simulation method were validated by the flow resistance coefficients, which were calculated by the Carman equation. The tracer pulse method and step method were employed to evaluate the residence time distribution characteristics in a fixed-bed reactor. It was determined that the mean residence time decreased when the particle size distribution was considered. The tracer residence time and variance increased as the inlet velocity decreased. An almost inverse proportional relationship between the mean residence time and inlet velocity was observed. Higher axial dispersion coefficients were observed when particle size distribution was considered.

The distribution characteristics of the tracer concentration and fluid velocity were obtained and used to explain the mixing performance of the fixed bed. The nonideal flow phenomenon that deviates from slug flow was mainly determined by the uneven velocity distribution, which was also responsible for the concentration distribution characteristics in the clearance between particles. When the fixed bed was packed with particles with a size distribution, the wall effects were further weakened. A more compact fixed bed can yield a higher flow velocity, thereby leading to a decrease in the mean residence time.