Development and assessment of speed-up algorithms for the reactive CFD–DEM simulation of fluidized bed reactors

Abstract

In this work, we propose a particle agglomeration (PA) speed-up algorithm to reduce the computational cost of the CFD–DEM framework based on the coupled solution (CP) of the gas–solid transport and heterogeneous reaction at the level of each particle in a fluidized bed reactor. The performances of PA are assessed and compared with the ones provided by the previously proposed operator-splitting (OS) and in situ adaptive tabulation (ISAT) speed-up algorithm. The analysis is carried out on methanation and steam reforming systems to investigate different reactions and transport characteristic times under different operating conditions. The results revealed that OS & ISAT is generally computationally more efficient by providing a higher chemical speed-up than CP & PA both in methanation (15 vs. 6) and steam reforming (13 vs. 4). However, OS requires a smaller simulation time step with respect to the CP to achieve convergent results when fast transport/reaction phenomena are considered. Hence, the overall computational cost might be even larger with respect to the CP despite the ISAT technique. Consequently, the selection of the best performing algorithm is not trivial and strongly depends on the characteristic times of transport and reaction which are a function of the local conditions experienced in the bed. Thus, we propose a strategy to select the most adequate algorithm (i.e., CP & PA or OS & ISAT). The strategy is based on the simulation of a test reactor characterized by a small computational cost but conceived to be representative of the chemical and fluid dynamic behavior of the target reactor. This enabled the simulation of a million particle methanation reactor experimentally investigated in the literature. A good agreement between simulation and experimental outlet concentrations (error up to 7% for the main species) is observed along with a significant 20-fold chemical speed-up.

---

Introduction

Fluidized bed technology is of great interest in the context of important catalytic processes, e.g. methanation^1–3 or methanol-to-olefin.^4,5 The description of the complex multiphase fluid dynamics coupled with reactivity is a crucial task for gaining detailed insights into the system behavior^6–8 and eventually for designing fluidized systems.⁹ In particular, reactive CFD–DEM is usually applied to track a few thousands of reactive particles in fluidized beds in the order of 10⁵ inert sand particles in the context of char combustion¹⁰ or biomass gasification.¹¹ In the case of catalytic processes, all the particles are catalytically active, and a computational cell basis for the evaluation of the catalytic chemistry is mandatory.¹² To overcome such a limitation, in our previous work, we proposed a reactive numerical framework,¹³ based on a computational fluid dynamics–discrete element method (CFD–DEM) algorithm,^14,15 able to couple the catalytic reactions at the active sites with the Lagrangian tracking of each particle in the fluidized system along with the solution of Navier–Stokes and transport equations at the reactor scale. However, the application of CFD–DEM to dense lab scale catalytic reactors (usually in the order of 10⁶ particles) is strongly hampered by the computational cost related to the chemistry. In this respect, we implemented an in situ tabulation speed-up technique^13,16,17 (ISAT) to reduce the number of required ODE integrations. In particular, we first implemented the solution of the species and energy balances in each particle by means of the operator-splitting^18,19 (OS) approach. By doing so, the chemistry and gas–particle transport phenomena were solved separately, within a time step for each particle in the system, enabling the efficient application of ISAT to the chemistry solution, and the analytical solution of the interphase transport phenomena.¹³ The combined OS & ISAT speed-up algorithm showed relevant speed-up factors up to 12 in a reactor of 10⁴ particles. However, the accuracy¹⁹ and the speed-up¹³ provided by the algorithm are dependent on the minimum time step used to split gas–particle transport and reactions. This was found to affect the maximum allowed simulation time step depending on the characteristic times of gas-particle transport and catalytic reactions. In the case that this maximum time step of splitting is very small, the computational cost of the simulation becomes the limiting factor.

In this work, we propose a numerical framework based on the coupled solution of transport and reaction terms (i.e. coupled approach, CP). In fact, by accounting the reaction and the transport at the same time, this method can be a valuable alternative when OS becomes too limiting in the selection of the maximum time-step. Moreover, to tackle the computational burden in the CP approach, we propose a particle agglomeration (PA) speed-up technique based on cell agglomeration (CA).^20–22 According to this speed-up algorithm, similar particles along the fluidized bed are binned at each time step in computational parcels. The parcels are solved by means of the integration of the ODE system describing transport and reactions coupled, and the results of the parcel are finally mapped back to the underlying particles. As a result, the number of ODE integrations are reduced with a consequent beneficial effect on the computational time.

To assess the performance of the CP & PA speed-up algorithm with respect to the previously proposed OS & ISAT, we applied the two algorithms to the reactive CFD–DEM framework¹³ under different operating conditions. The CP & PA algorithm has provided a very good accuracy and a relevant speed-up under the investigated conditions. However, when the selection of the maximum time step turned out not to be limiting, the OS & ISAT algorithm has provided higher speed-up factors.

The selection between the two procedures cannot be performed a priori, especially when complex mixtures are present. As such, we propose a methodology for selecting the best algorithm for different cases. First, a small test reactor is considered and employed to identify the best performing algorithm and to achieve an optimal selection of the parameters of the speed-up technique. Then, the resulting optimal algorithm is employed for the simulation of the target reactor, where the same numerical tests are very computationally demanding or even not sustainable.

Thanks to this procedure, we were able to simulate the 1.2 million Geldart A particle methanation reactor, experimentally investigated in the literature by Li and Yang.²³ The results have shown good agreement between the predicted outlet compositions of the reactor and the experimental data (error up to 7% for the main species) with a 20-fold speed-up factor.

All in all, these methods make it possible to simulate fluidized bed with a particle resolved description of the chemistry and with the number of particles in the order of 10⁶. The possibility of accounting for the chemical phenomena at the particle level with such a number of particles paves the way towards the hierarchical refinement²⁴ of Euler–Euler models for the simulation of devices of technical relevance.

Coupled & particle agglomeration speed-up algorithm

The governing equations are reported in the ESI† – section 1. They consist of gas-phase and particle mass and energy conservation equations and the Navier–Stokes equations. These equations are sequentially solved for a given time step. Within the time step, the reaction and transport operators in the mass and energy equations at the particle can either be treated by means of operator-splitting (OS) or by means of a coupled approach (CP) – the two operators are simultaneously accounted for in the numerical solution. In our previous work, we presented and assessed the OS approach with the ISAT technique to speed-up the simulation. Here, instead, we propose and assess a particle agglomeration (PA) algorithm for the speed-up of the numerical solution of the ODE system within the CP approach.

PA is a runtime classification technique, based on the cell agglomeration (CA) algorithm,²² aimed at decreasing the computational cost through the reductions of ODE integrations by agglomerating similar particles along the fluidized bed. According to PA, an index (i_PA) is computed at each time step for each particle in the domain. This index is a measure of the similarity among the particles in terms of characteristic variables, i.e., composition, coverages, temperature and the Reynolds number of the particle. The gas thermo-chemical and fluid dynamic states are also considered in the index formula for PA, since they are key factors in the quantification of gas–solid transfer rates. Then, the particles characterized by the same index are grouped into a computational parcel which is integrated by the ODE solver. Finally, the result of the parcel is used to update all the particles that belong to that parcel.

The index i_PA is computed as follows:

where NV is the total number of the characteristic variables, and ϕ_i is the i-th normalized variable and tol_PA is the tolerance value. The ceiling of the ratio ϕ_i/tol_PA is considered during the computation of the index in eqn (1) to prevent negative values of the logarithmic function. Here, differently from the CA algorithm,^20,22 we consider the logarithmic operator in eqn (1) as a hash mapping function. This is due to the fact that when systems with high dimensionality (i.e. NV) are considered, the typical CA hash function²⁰ could result in numerical overflow in the evaluation of the index. In eqn (1), tol_PA affects the number of particles per computational parcel: the lower the tolerance value the stricter the particle similarity criterion, thus resulting in fewer particles per parcel, and, in turn, larger number of total parcels. It is a numerical parameter which is usually determined by means of a trial and error procedure for the different catalytic chemistries under investigation.

The normalized variables ϕ_i are computed as follows:

where y_i,max and y_i,min are the maximum and minimum values of the ith variable. In case the minimum value of the ith variable is zero, a value of y_i,min equal to 10⁻²⁰ is employed as the minimum value. In particular, a normalization of the variables based on a logarithmic scale ranging from the minimum to the maximum value of the variable calculated in the fluidized bed has been introduced with respect to the linear scale usually adopted in the case of CA.^20,22 We found this choice to be more robust especially in the cases where there is a strong uneven distribution of the species mass fractions in the composition space, for instance, during the simulation of the start-up of the fluidized bed reactor. In these situations, linear scaling is not able to properly discriminate the sharp transition, thus binning the particles characterized by different compositions, unless a very low tolerance is employed. By using a logarithmic scale, instead, the classification is refined close to low concentrations, thus improving the efficiency of the classification without reducing the tolerance.

Once the parcels are generated by means of eqn (1)–(2), the properties of each parcel are computed as the arithmetic average of the properties of the particles in the bin. Then, the parcel is integrated, and the results are mapped back to each particle by means of the following formula:²²

y^f_i,p = y⁰_i,p + (y^f_i,b − y⁰_i,b) (3)

where y⁰_{i, p} and y^f_{i, p} are the initial and final conditions of the particle, whereas y⁰_{i, b} and y^f_{i, b} are the initial and final conditions of the parcel where particles have been agglomerated.

Procedure for the selection of the test reactor

The selection of the two procedures (i.e., OS & ISAT and CP &PA) is not always obvious, especially when several species are present with significantly different diffusion properties. Here, we propose a methodology for selecting the best algorithm and for tuning its parameters. As shown in Fig. 1, a small test reactor – which retains the features of the target reactor – is considered and employed to identify the most adequate algorithm and to optimally tune the parameters of the speed-up technique. In particular, once the test reactor is selected (step 1), CP & PA and OS & ISAT are assessed by using this test reactor and for the process under investigation, thus avoiding the simulation of the massive number of particle reactor. By means of this step, the choice of the best performing speed-up algorithm and the tuning of its appropriate tolerance are carried out (step 2). Finally, the simulation of the lab scale reactor is performed adopting the selected and properly tuned speed-up algorithm (step 3).

Fig. 1 Schematic representation of the three step strategy for the selection of the best performing algorithm along with the tuning of its parameters. First, a test reactor is designed to provide similar fluid dynamic and chemical behaviour of the target reactor. Then, the most adequate algorithm is found by comparing the accuracy and speed-up evaluated for the test reactor. Finally, the selected approach is employed to carry out the simulation of the target reactor.

Crucial for the reliability of the procedure is the selection of the test reactor. The geometry of both the test bed and particles has to be carefully set to reproduce the target reactor behavior. Moreover, the inlet composition and temperature along with the inlet superficial velocity must be imposed to be representative of the target reactor. In doing so, it is crucial to ensure that the test reactor shows the same chemical and interphase transport phenomena occurring in the target reactor. In particular, the geometry of the bed can be the same shape of either the lab scale reactor or a pseudo-2D parallelepipedal reactor, typically used in the literature for the numerical investigation of small fluidized beds.^25–27 The ratio between the height and width of the test reactor must be chosen to avoid the occurrence of slugging behavior.²⁸ The mechanical characteristics of the particles must be the same as the ones used in the target reactor since they affect the fluidization behavior of the bed.²⁹ The temperature and composition of the target reactor are adopted, since they lead to the similar characteristic times for both the interphase species mass transfer and catalytic reactions. The inlet superficial velocity and the number of particles in the bed (and, in turn, the height of the bed) are finally selected to achieve fluid dynamic and chemical similarities with respect to the target reactor. The fluid dynamic similarity is obtained by comparing the fluidization ratio of the two reactors, defined as the ratio between the adopted inlet superficial velocity and the minimum fluidization. The chemical similarity is achieved by assuring a comparable gas-hourly-space-velocity (GHSV) in the test and target reactors. Both fluid-dynamics and GHSV are affected by the number of particles. Thus, a trade-off between the two constraints must be achieved, and it can be reached by using a maximum ratio of 2 between the GHSV of the two reactors and between the fluidization ratio of the two reactors to assure an adequate similarity between the test and lab scale reactors.

Computational domains and numerical settings

We focus our analysis on the methanation million-particle fluidized bed reported by Li and Yang,²³ that we use as a showcase in this work (target reactor). In this section, we present the numerical details of the test reactor and of the target reactor used for the analysis.

Million-particle lab scale reactor (target reactor)

The target reactor is characterized by a cylindrical geometry with a diameter of 1.5 cm and a height of 9 cm, as shown in Fig. 2. The bed consists of 1.2 million particles, which are initially packed at the bottom of the reactor. The geometrical and mechanical characteristics of the particles are reported in Table 1.

Fig. 2 Cylindrical reactor employed for the simulation of the lab scale reactor and packed bed configuration adopted as a starting point for the reactive simulations of the methanation lab scale reactor reported by Li and Yang.²³

Table 1 Geometrical features of the methanation lab scale reactor reported by Li and Yang²³ along with the particles' mechanical properties

Geometrical properties
Particle diameter DP [μm]	125
Particle density ρP [kg m^−3]	1684
Particle number NP [−]	1 244 130

---

Mechanical properties
Young modulus E [Pa]	10^6
Poisson ratio ν [−]	0.3
Restitution coefficient e [−]	0.8
Friction factor μ [−]	0.3

---

The following numerical settings have been adopted for the simulation of the catalytic process. An atmospheric pressure has been imposed at the top of the reactor, while a zero-gradient condition has been assumed for the lateral walls and for the bottom (i.e., inlet). At the bottom of the reactor, three superficial velocities equal to 0.162, 0.270 and 0.323 m s⁻¹ have been investigated. A no-slip condition has been imposed for the lateral walls of the cylinder. The temperature and composition of the gas phase have been imposed at the inlet of the reactor equal to the operating conditions of the inlet feed stream. A zero-gradient condition has been imposed on the remaining boundaries. The simulations are carried out under isothermal conditions. The computational grid has been generated ensuring a minimum cell-to-particle-diameter ratio equal to 3.³⁰ Additional information on the mesh and on the numerical schemes is reported in the ESI† – section 2.

Test reactor

The pseudo-2D parallelepiped geometry employed in our previous work¹³ has been chosen. Its dimensions have been adapted to obtain the reactor configuration (i.e. 6 mm width, 3 cm height and 500 μm depth) containing 10⁴ particles with the same mechanical properties reported in Table 1. The number of particles is equal to the one simulated in our previous work, which provided representative fluidization behaviors of the bed and an affordable computational time. Thus, the inlet superficial velocity has been selected to be equal to 0.152 m s⁻¹ in order to keep a maximum ratio of 2 between the fluidization ratio and the GHSV of the two reactors, as reported in Table 2, for the case of the superficial velocity equal to 0.162 m s⁻¹ assumed to be representative of all the investigated conditions. Finally, the same composition and temperature adopted in the million-particle reactor have been used. The simulations are carried out under isothermal conditions. Additional information on the mesh and on the numerical schemes is reported in the ESI† – section 2.

Table 2 Fluidization ratio and space velocity of the lab scale reactor (ϕ = 1.5 cm, H = 9 cm) and the pseudo-2D test reactor (0.6 cm × 0.05 cm × 3 cm)

	Test reactor	Lab scale reactor
Fluidization ratio [−]	20	40
GHSV [1/h]	70 728	40 416

---

Average quantities, error estimation and speed-up calculation

The average composition of the bed 〈ω〉 as a function of the simulation time has been selected as the reference variable representative of the reactor behavior:^13,25,26where k is the k-th simulation time step and NP is the total number of particles.

The accuracy has been assessed by means of the evaluation of the mean global error:¹⁷

where Y_p is the vector containing the composition and temperature of the p-th particle in the fluidized bed, and k is the k-th simulation time step. Moreover, the temporal average of the mean global error has been used:where NT is the total number of time steps composing the simulation.

The speed-up of the solution of the catalytic particles, i.e., the chemical speed-up factor (SP_C), has been computed to measure the computational gain provided by each speed-up algorithm:

where τ^speed-up_k and τ^{no speed-up}_k are the computational cost, i.e., the wall clock time, required for the solution of the interphase transport and reaction in all the particles of the bed at the k-th time step with and without speed-up algorithms. In the case of parallel simulations, they represent the maximum computational cost among the processors at the k-th time step.

Results and discussion

Here, we first present the assessment of the performance of CP & PA and OS & ISAT in the test reactor configuration for two different processes, i.e., methanation and steam reforming. The assessment is performed in terms of both accuracy and speed-up. The operating conditions are reported in Tables 3 and 5. The test reactor for the methanation process for the simulation of the target reactor reported by Li and Yang²³ has been described in the previous section. For the case of steam reforming, we have considered the same reactor geometry, particle dimensions and number of particles as well as isothermal conditions of the test reactor for the methanation case. Moreover, the inlet velocity has been set equal to 0.128 m s⁻¹ to obtain the same fluidization ratio of the methanation test reactor under the operating conditions selected for the methane steam reforming process (Table 5).

Table 3 Operating conditions of the methanation process²³

Methanation
Temperature [K]	673.15
Pressure [bar]	1.01325
Inlet mole fractions [−]
Carbon monoxide	0.2
Hydrogen	0.6
Argon	0.2

---

Then, the simulation of the target reactor (lab-scale reactor with 10⁶ particles) in the case of the methanation process is performed by employing the algorithm and the parameter obtained and optimized at the test reactor level. The accuracy of the CFD–DEM and the computational gain in this system characterized by such a massive number of particles are assessed and quantified with respect to the experimental data from the literature.²³

Analysis of the speed-up algorithms in the test reactor

Methanation reaction. First, the methanation system is considered. The rate equation kinetics proposed by Li and Yang²³ has been adopted. The operating conditions are reported in Table 3.

CP & PA algorithm. First, the CP & PA algorithm is analyzed. A time step of 5 μs is employed, since the convergence of the results predicted with the coupled approach (CP) has been obtained for this time step. Three PA tolerances have been investigated, i.e. 0.2, 0.1 and 0.05, to investigate their effect on the accuracy and speed-up. Fig. 3 shows the results of the simulations in the test reactor for CP & PA in terms of speed-up and average composition in the reactor. Fig. 3a shows the comparison of the chemical speed-up factors obtained with all the tolerances investigated. The CP & PA algorithm is capable of providing a satisfactory chemical speed-up factor between 2.5 and 10.5 according to the tolerance. In particular, a stricter (i.e., lower) tolerance decreases the computational speed-up due to the larger number of parcels that have to be solved. The performance of the CP & PA algorithm is then analyzed in terms of accuracy by means of the evaluation of the temporal average error as compared with the solution without the PA speed-up technique. A very good agreement between the simulations carried out with and without speed-up techniques is obtained with all the tolerances (maximum error 1.2% – Fig. 3a). The optimal (i.e., trade-off between speed-up and accuracy) value of the tolerance is found to be equal to 0.1 at which a speed-up of around 6 along with an average error of 0.9% is obtained. Below such a tolerance, indeed, a decrement of the speed-up without no substantial improvement in the accuracy of the algorithm is observed.

Fig. 3 Speed-up (a) and accuracy (b) obtained with the CP & PA algorithm in the case of the methanation process. The speed-up is evaluated by comparing the speed-up factor related to the solution of the particle ODE systems obtained for the three PA tolerances investigated. The average speed-up value over the whole simulation is also reported for each tolerance (dashed line). The accuracy is assessed by comparing at each time step the average composition of the bed obtained without speed-up (solid line) and the one obtained with CP & PA (open symbols) for the optimal PA tolerance.

Fig. 3b shows the average composition of the bed as a function of the simulation time for the main species involved in the process obtained with the coupled algorithm with and without the PA speed-up technique. The comparison of the results makes it evident that the excellent temporal averaged error corresponds to a good agreement between the simulations with and without the speed-up algorithm at each time step (maximum mean global error ε equal to 2.49%), thus proving the accuracy provided by the algorithm both during the dynamics and the pseudo-steady state regime of the process.

OS & ISAT algorithm. The performance of the OS & ISAT algorithm has then been analyzed in the same reactor and for the same operating conditions. First, an analysis of the effect of the time step on the predictions is carried out.

Fig. 4a shows the comparison of the average composition of the bed as a function of the simulation time for the main species involved in the process. The time step has a strong effect on the predictions of the species. We found that a time step equal to 1 μs is capable of providing a solution which is independent of the time discretization. In contrast, a time step equal to 5 μs, i.e. the same employed in the case of CP & PA, provides inaccurate results since an overestimation of the hydrogen content in the order of 15% is observed, as shown in Fig. 4a. In this case, a violation of the atomic balances is observed due to the misprediction of the hydrogen concentration. The reason of the deviation is ascribed to the OS algorithm which accounts for the competitive and simultaneous transport and reaction phenomena in separate terms. To guarantee the convergence to the solution, it is necessary that the simulation time step is lower than the characteristic time of the transport of the species.

Fig. 4 Comparison of the average composition of the bed as a function of the simulation time obtained by solving the transport and reaction according to operator-splitting (OS) and coupled (CP) approaches. The effect of the time step is assessed by comparing the results obtained with 5 μs (triangle) and 1 μs (square) (a) along with a comparison of the performance of the OS (square symbol) with respect to the CP (continuous line) in the case of a time step equal to 1 μs (b).

To better clarify this issue, we have computed the characteristic times of the gas–particle transport of all the reactive species involved in the process according to eqn (8):³¹

where K_c,j is the mass transfer coefficient of the j-th species and S_v is the specific geometric external area of the catalytic particle. The characteristic time of species transport has been compared with the time step used in operator splitting (i.e. the simulation time step). The minimum, maximum and average characteristic times of the species gas–particle transfer for all the particles in the fluidized bed have been evaluated (Table 4). The average time has been obtained as the arithmetic average of coefficients computed for all the particles in all the investigated simulation times. The characteristic time of the hydrogen gas–particle transfer is lower than the splitting time, i.e., 5 μs. Thus, the loss of accuracy of OS is related to the decoupled solution of the hydrogen interphase transfer from the catalytic reactions when the time step is larger than the maximum allowed. A smaller time step, i.e. lower than the hydrogen transport characteristic time, avoids the inaccuracies. Fig. 4b shows the average composition of the bed for a simulation carried out with a time step equal to 1 μs, i.e., lower than the characteristic time of the hydrogen transport. An excellent agreement is observed with respect to the coupled approach along the entire simulation time confirming the capabilities of the OS to predict the species profiles when the appropriate time step is considered.

Table 4 Minimum, maximum and average mass transfer characteristic time for all the species present in the methanation system. Minimum and maximum characteristic times are computed from minimum and maximum mass transfer coefficients calculated in the bed (eqn (8)), whereas the average time is computed from the arithmetic mean of the mass transfer coefficient over all the particles and all the investigated simulation times

	τ min [μs]	τ average [μs]	τ max [μs]	τ max [μs]
CH4	7.18	9.77	11.8	11.8
H2	3.33	3.84	4.10	4.10
CO	7.59	10.5	12.8	12.8
CO2	9.20	12.9	16.1	16.1
H2O	6.08	8.61	10.5	10.5

---

However, the smaller time step required by the OS dramatically increases the number of time steps necessary to simulate the same real time resulting in an additional computational burden. To improve the computational efficiency, we assessed the effect of the ISAT technique in reducing the simulation computational cost.¹³ The performance of the ISAT speed-up technique has been tested by considering different tolerances. In particular, three values equal to 1 × 10⁻⁵, 5 × 10⁻⁶, and 1 × 10⁻⁶ have been employed. ISAT boosts the simulation by reducing the computational burden of the chemical sub-step. An average chemical speed-up equal to 15.5, 14.9 and 9.0, respectively, has been found resulting in an overall simulation speed-up¹⁷ in the order of 5.5, 5.4 and 4.5, respectively. The simulation speed-up provided by ISAT is able to only partially recover the additional cost due to the reduced time step. However, the overall computational cost of the OS & ISAT method is similar to the one of the sole CP algorithm. Hence, the overall performance of the OS & ISAT is worse with respect to the CP & PA method for these operating conditions.

Steam reforming. An analogous analysis has been simulated in the same reactor in the case of steam reforming. Steam reforming has been selected as the second case study since it involves the same species of the methanation reactor, obtaining the same complexity of the mixture. The operating conditions are reported in Table 5 and they are selected to span different ranges of both characteristic times for the transport and reaction due to the higher pressure and temperature with respect to the methanation case. The rate equation kinetics proposed by Donazzi et al.³² have been employed.

Table 5 Operating conditions of the methane steam reforming process³³

Methane steam reforming
Temperature [K]	1048.15
Pressure [bar]	29
Inlet mass fractions [−]
Methane	0.1911
Steam	0.7218
Carbon monoxide	0.0001
Carbon dioxide	0.0200
Hydrogen	0.0029
Nitrogen	0.0641

---

CP & PA algorithm. We start with an analysis of the performance of the CP & PA algorithm. The same PA user tolerances used in the methanation process have been herein investigated. Fig. 5a shows the comparison of the chemical speed-up factors obtained with all the tolerances investigated. CP & PA provides a speed-up between 2 and 7 according to the tolerance. A very good accuracy has been observed for all the tolerance with average global error below 1%. Fig. 5b shows the comparison of the temporal trend of the average composition of the bed obtained with and without the CP & PA speed-up algorithm for a tolerance of 0.1.

Fig. 5 Speed-up (a) and accuracy (b) obtained with the CP & PA algorithm in the case of the steam reforming process. The speed-up is evaluated by comparing the speed-up factor related to the solution of the particle ODE systems obtained for the three PA tolerances investigated. The average speed-up value over the whole simulation is also reported for each tolerance (dashed line). The accuracy is assessed by comparing at each time step the average composition of the bed obtained without speed-up (solid line) and the one obtained with CP & PA (open symbols) for the optimal PA tolerance.

OS & ISAT algorithm. Then, the performance of the OS & ISAT algorithm has been investigated by employing the same three ISAT tolerances of the methanation case. Moreover, the time step (i.e., the splitting time) selected for the CP & PA analysis has been employed. The convergence of the OS solution with this time step was verified. Fig. 6 shows the results of the OS & ISAT algorithm in comparison with the results obtained by means of the reactive CFD–DEM framework without any speed-up algorithm. Fig. 6a shows the comparison of the chemical speed-up factors obtained with all the tolerances investigated and the corresponding temporal averaged error. The OS & ISAT algorithm is capable of providing a relevant chemical speed-up factor up to 13. Moreover, it has been found that the tolerance has a minor influence on the computational gain since all the speed-up profiles are quasi-superimposed. Nevertheless, it can still be noticed that a reduction of the tolerance decreases the speed-up as expected. In this specific case, the optimal tolerance (i.e., trade-off between speed-up and accuracy) has been found to be equal to 10⁻⁵, below which no substantial improvements in the accuracy provided by the algorithm are observed (temporal averaged errors between 0.66% and 0.65%) at the expense of a detrimental effect of the speed-up. As a further insight, the comparison of the average composition of the bed as a function of time for the main species involved in the process is reported in Fig. 6b in the case of the optimal tolerance. An excellent agreement has been obtained during the entire simulation for all the species. This is ascribed to the proper accounting of the gas–particle species mass transfer phenomena. Steam reforming is operated at both high temperature and pressure. The effect of temperature and pressure on the transport characteristic time is opposite. An increment of temperature reduces the characteristic times due to the higher diffusivity, whereas an increment of pressure increases the characteristic times due to a decrement of the diffusivity. In this specific case, the dominant effect is due to the pressure, resulting in a slower gas–solid transport. Table 6 lists the characteristic times of transport of all the species. In particular, it is evident that the characteristic time of the hydrogen transport is around 30 μs. This characteristic time is higher than the splitting time step (5 μs) in line with the fact that the solution is converged. The analysis of steam reforming in the test reactor revealed that the OS & ISAT algorithm surpasses the performance of the CP & PA since a three times higher speed-up factor is obtained.

Fig. 6 Speed-up (a) and accuracy (b) obtained with the OS & ISAT algorithm in the case of the methane steam reforming process. The speed-up is evaluated by comparing the speed-up factor related to the solution of the particle ODE systems obtained for the three ISAT tolerances investigated. The average speed-up value over the whole simulation is also reported for each tolerance (dashed line). The accuracy is assessed by comparing at each time step the average composition of the bed obtained without speed-up (solid line) and the one obtained with OS & ISAT (open symbols) for the optimal ISAT tolerance.

Table 6 Minimum, maximum and average mass transfer characteristic times for all the reactive species present in the methane steam reforming system. Minimum and maximum characteristic times are computed from minimum and maximum mass transfer coefficients calculated in the bed (eqn (8)), whereas the average time is computed from the arithmetic mean of the mass transfer coefficient over all the particles and all the investigated simulation times

	τ min [μs]	τ average [μs]	τ max [μs]
CH4	79	119	160
H2	31	43	54
CO	83	140	160
CO2	99	176	189
H2O	77	125	127

---

As a whole, the analysis revealed that different speed-up strategies may have different performances even in similar reacting environments elucidating the importance of the selection of the most adequate algorithm for the efficient simulation of reactors characterized by a large number of particles.

Simulation of the million particle reactor: methanation

The results of the analysis of the test reactor enabled us to properly define both the numerical algorithm and the parameters to carry out a numerical simulation of a lab scale reactor. The selection of the best performing speed-up algorithm and its tuning is key to enable the simulations of million particle reactors whose computational cost would be unaffordable without any speed-up technique. Thus, we employed the algorithm and the optimal parameters obtained in the analysis of the test reactor. Thus, the methanation lab scale reactor (“the target reactor”) experimentally investigated by Li and Yang²³ has been simulated by means of the CP & PA algorithm, which has been the most effective for this specific process.

The kinetics proposed by Li and Yang²³ (already used for the assessment of the algorithms in the test reactor) have been adopted, since it has been validated against the experimental data²³ for the operating conditions investigated in this study.

The operating conditions employed in the experiments have been adopted (see Table 3). Three different inlet velocities equal to 0.162, 0.270 and 0.323 m s⁻¹ have been selected to investigate the effect of different fluid dynamic conditions.

The following procedure has been followed to perform the simulation of the reactor. First, the packed bed configuration (Fig. 2) has been generated by injecting the particles from the top of the reactor and allowing them to settle. A particle flow rate of 8 × 10⁶ particles per second was employed. Then, the system has been fluidized by using an inert stream of Argon at 673.15 K injected from the bottom with a velocity of 0.162 m s⁻¹. Once a steady fluidization is achieved, the syngas has been fed from the bottom, and the reactor has been simulated until a pseudo-steady state of the process has been achieved, i.e. reproducing a real time higher than 8 times the residence time in the bed. Then, the value of the relevant quantities, i.e. pressure drops and concentrations, has been evaluated by temporal averaging the last two residence times. The other simulations have been carried out starting from the pseudo-steady state results at 0.162 m s⁻¹ by changing the inlet velocity value.

First, we assessed the fluid dynamic predictions of the reactive CFD–DEM framework by observing the temporal trend of the pressure drops inside the reactor. Fig. 7 shows the pressure profile which is characterized by an oscillating behavior caused by the expansion and contraction behavior of the fluidized bed as usually observed in fluidized bed reactors. The temporal average of the pressure drops is in excellent agreement (deviation below 5%) with the theoretical value obtained as the ratio between the weight of the bed divided by the cross-sectional area of the reactor.

Fig. 7 Pressure drops in the methanation fluidized bed reactor as a function of the simulation time (solid line). The dashed line represents the theoretical value of pressure drop in the fluidized regime.

Then, the accuracy has been assessed by a comparison of the outlet composition measured in the experiments with the one evaluated in the simulation. The reactor outlet composition has been evaluated at each time step by means of the cup-mix average over a section at a height of 4 cm from the bottom, i.e. slightly above the end of the expanded bed.

Fig. 8 shows a parity plot between the concentrations measured in the experiments²³ and the results of numerical simulations for all the investigated flow velocities. The results show a good agreement with the experiments, since all the species are within the experimental confidence interval (±5% for the main species and ±10% for CO₂). These results confirm the accuracy of the CFD–DEM framework and of the CP & PA algorithm for the simulation of catalytic processes in fluidized bed reactors.

Fig. 8 Parity plot comparing the numerical results with the experimental data by Li and Yang²³ for the outlet reactor compositions predicted with the reactive CFD–DEM framework and the CP & PA speed-up algorithm proposed in this work at 0.162 m s⁻¹ (square), 0.270 (triangle) m s⁻¹ and 0.323 (circle) m s⁻¹ inlet velocities.

Finally, the computational gain offered by the CP & PA speed-up algorithm in the case of the million-particle reactor has been computed. The entire simulation without any speed-up technique is unfeasible due to the computational cost. Hence, the evaluation of the computational gain is not carried out along the entire simulation but locally every 0.2 s according to the following procedure. Two time steps have been simulated with and without the speed-up algorithm, starting from the results obtained by means of CP & PA. Then, the ratio between the computational times required by the reactive CFD–DEM framework with and without the speed-up technique has been computed to evaluate the speed-up factor according to eqn (7). Fig. 9 shows the evolution of the chemical speed-up factor during the simulation. A chemical speed-up-factor up to 27 and an average value around 20 has been quantified, corresponding to a significant computational gain for the entire simulation.

Fig. 9 Trend of the chemical speed-up factor provided by the CP & PA algorithm as a function of the simulation time in the case of the methanation million particle reactor.

Conclusions

A speed-up algorithm based on the coupled ODE integration of the gas–particle transfer and catalytic chemistry (CP) and the PA speed-up technique (i.e. CP & PA) has been implemented in the existing CFD–DEM framework and compared to the operator-splitting & ISAT algorithm developed in our previous work.¹³ Then, a comparison between the performance of CP & PA and OS & ISAT has been carried out in a test reactor for both a methanation and a steam reforming system.

The analysis revealed that the OS algorithm requires a careful definition of the simulation time step which has to be lower than the minimum characteristic time of transport phenomena to achieve a converged solution. In contrast, the CP algorithm partially overcomes this limitation due to the simultaneous accounting for both the reaction and the transport. At the same time the performances of the speed-up techniques are different. The ISAT generally outperforms the PA providing usually higher computational gain. In this view, the selection of the approach, i.e., CP & PA or OS & ISAT, turned out to be strongly dependent from the system and the operating conditions investigated. Hence, the analysis carried out on the test reactor enables the assessment of the different approaches along with the tuning of the algorithm parameters.

The results of this analysis in the test reactor for the methanation system enabled the selection of the approach (i.e., CP & PA) and optimized the algorithm tolerance. The results have proven the capability of the framework with the speed-up strategies of simulating millions of bed particles. In particular, a good agreement with the experimental data (error up to 7% for the main species) is obtained along with a 20-fold chemical speed-up factor, which translates into significant reduction of the computational cost related to the solution of the heterogeneous chemistry which now turns to be comparable with the cost of the particle tracking.

As a whole, the application of the reactive CFD–DEM framework by means of speed-up algorithms has been effective both in small test fluidized beds and in massive number of particle lab scale reactors, thus, enabling both the fundamental investigation of lab scale reactors and their efficient simulation in the context of a hierarchical analysis of fluidized systems.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The project leading to this work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Project SHAPE - grant agreement No. 677423). Computational time at CINECA, Bologna (Italy) is gratefully acknowledged.

References

1. J. Li, H. Wang, Q. Zhu and H. Li, Chem. Eng. J., 2019, 357, 298–308.
2. L. Sun, K. Luo and J. Fan, Catal. Commun., 2018, 105, 37–42.
3. J. Liu, D. Cui, C. Yao, J. Yu, F. Su and G. Xu, Fuel Process. Technol., 2016, 141, 130–137.
4. S. Xu, Y. Zhi, J. Han, W. Zhang, X. Wu, T. Sun, Y. Wei and Z. Liu, Advances in Catalysis for Methanol-to-Olefins Conversion, Elsevier Inc., 1st edn., 2017, vol. 61.
5. P. Tian, Y. Wei, M. Ye and Z. Liu, ACS Catal., 2015, 5, 1922–1938.
6. N. Jurtz, M. Kraume and G. D. Wehinger, Rev. Chem. Eng., 2019, 35, 139–190.
7. B. Partopour and A. G. Dixon, Ind. Eng. Chem. Res., 2019, 58, 5733–5736.
8. M. Maestri, Chem. Commun., 2017, 53, 10244–10254.
9. M. P. Dudukovic, Science, 2009, 325, 698–701.
10. Y. Zhuang, X. Chen, D. Liu and C. Bu, in Clean Coal Technology and Sustainable Development. ISCC 2015, Springer, 2016.
11. S. Yang, H. Wang, Y. Wei, J. Hu and J. W. Chew, Energy Convers. Manage., 2019, 196, 1–17.
12. D. M. Snider, S. M. Clark and P. J. O'Rourke, Chem. Eng. Sci., 2011, 66, 1285–1295.
13. R. Uglietti, M. Bracconi and M. Maestri, React. Chem. Eng., 2018, 3, 527–539.
14. Y. Tsuji, T. Kawaguchi and T. Tanaka, Powder Technol., 1993, 77, 79–87.
15. B. H. Xu and A. B. Yu, Chem. Eng. Sci., 1997, 52, 2785–2809.
16. S. B. Pope, Combust. Theory Modell., 1997, 1, 41–63.
17. M. Bracconi, M. Maestri and A. Cuoci, AIChE J., 2017, 63, 95–104.
18. G. Strang, SIAM J. Numer. Anal., 1968, 5, 506–517.
19. Z. Ren and S. B. Pope, J. Comput. Phys., 2008, 227, 8165–8176.
20. G. M. Goldin, Z. Ren and S. Zahirovic, Combust. Theory Modell., 2009, 13, 721–739.
21. L. Zhou, L. Zhong, W. Qin, W. Zhao and H. Wei, Fuel, 2018, 234, 1313–1321.
22. S. Rebughini, A. Cuoci, A. G. Dixon and M. Maestri, Comput. Chem. Eng., 2017, 97, 175–182.
23. J. Li and B. Yang, AIChE J., 2019, 65, 1–15.
24. S. Rebughini, M. Bracconi, A. G. Dixon and M. Maestri, React. Chem. Eng., 2018, 3, 25–33.
25. Q.-F. Hou, Z.-Y. Zhou and A.-B. Yu, Ind. Eng. Chem. Res., 2012, 51, 11572–11586.
26. M. J. V. Goldschmidt, R. Beetstra and J. A. M. Kuipers, Powder Technol., 2004, 142, 23–47.
27. J. Li, R. K. Agarwal, L. Zhou and B. Yang, Chem. Eng. Sci., 2019, 207, 100031.
28. D. Kunii and O. Levenspiel, Fluidization engineering, Butterworth-Heinemann, 1991.
29. D. Geldart, Powder Technol., 1973, 7, 285–292.
30. Z. Peng, E. Doroodchi, C. Luo and B. Moghtaderi, AIChE J., 2014, 60, 2000–2018.
31. X. Gao, B. Kong and R. D. Vigil, Bioresour. Technol., 2015, 198, 283–291.
32. A. Donazzi, A. Beretta, G. Groppi and P. Forzatti, J. Catal., 2008, 255, 259–268.
33. K. R. Rout and H. A. Jakobsen, Can. J. Chem. Eng., 2015, 93, 1222–1238.

---

Footnote

† Electronic supplementary information (ESI) available: Reactive CFD–DEM framework; computational domain and numerical methods; parametric analysis of the effect of y_i,min. See DOI: 10.1039/c9re00440h

---