Coupling CFD–DEM and microkinetic modeling of surface chemistry for the simulation of catalytic fluidized systems

A numerical framework is proposed to couple detailed microkinetic modeling and CFD–DEM for the simulation of gas–solid fluidized systems.

1 Closure models for the gas-particle, particle-particle and particle-wall interactions

Gas-particle heat and mass transfer
The Ranz-Marshall correlation has been selected to properly compute the Nusselt number Nu, and thus the coefficient h, to quantify the heat transfer rate between the catalyst and the gas phase in the energy balance of the catalytic particle (Eq. (1)): = 2 + 0.6 0.5 0.33 (S1) where = ‖ − ‖ is the particle Reynolds number dependent on the gas density , the gas dynamic viscosity , the norm of the gas-particle relative velocity ‖ − ‖, the particle diameter and the gas volumetric fraction in the computational cell hosting the particle p.
The gas-particle mas transfer rate has been modeled by means of the Ranz-Marshall correlation.
Therefore, the Sherwood number ℎ , and thus the coefficient , , has been computed according to Eq. (S2) to quantify the species mass transfer rate between the catalyst and the gas phase in the mass balance of the ℎ species in the catalytic particle (Eq.(2)) : ℎ = 2 + 0.6 1/2 1/3 (S2)

Gas-particle drag and buoyancy forces
The drag force contribution (Eq.(S3)) is computed proportionally to the gas-catalyst relative velocity ( − ) and to the particle volume : A combination of the empirical Ergun 1 and Wen-Yu 2 correlations has been selected for the drag coefficient , according to the Gidaspow 3  where is the gas dynamic viscosity, is the particle diameter, is the density of the gas phase at the particle position, − is the relative velocity between a generic particle p and the gas phase and is evaluated according to Schiller and Neumann 4 : The buoyancy force is characterized by means of the pressure gradient ∇ :

Particle-particle and particle-wall collisions
The normal and tangential component of the collisional force between two colliding particles a and b are computed according to Tsuji et al. 5 : where is the relative velocity between the two colliding particles a and b, is the unit vector connecting the centers of a and b, is the friction factor of the slider and and are the spring stiffness and the damping coefficient of the dashpot which can be derived from the Young modulus E, the Poisson ratio and the restitution coefficient of the particles as proposed in 5 . Subscripts n and t refer to the normal and tangential directions with respect to the plane normal to .
According to the selected approach, either two colliding particles a and b or particle a and wall b are allowed to slightly overlap of a quantity which is exploited to compute the normal (Eq.(S7)) and tangential (Eq.(S8)) components of the collisional force.
The normal and tangential components ( and ) of the overlapping vector are reported in Eqs.
Where and are the radii of the colliding particles a and b, ‖ − ‖ is the distance between the centers of a and b during the collisional event, and Δ is the time step for the solution of the governing equations of the particles, which must be selected lower than the characteristic time of particle collisional events, quantified as follows 6 : where = 1− 2 + 1− 2 is the effective Young Modulus and = 1/ + 1/ is the effective radius.

Computation of momentum, heat and mass transfer source terms for the gas phase balances
The momentum, heat and mass transfer rates between gas and a catalytic solid particle, are evaluated per unit of computational cell volume for each pellet p contained in a generic cell i, according to the drag force and Ranz-Marshall heat/mass transfer correlations, reported in the previous sections.
Finally, these contributions are summed up over all the particles contained in the i th cell of the computational domain, representing the reactor, to evaluate the source terms in the gas phase governing equations (Eqs.(7)-(9)): The resulting gas-particle exchange contributions , , ℎ, and , are added respectively to the source gas fields , ℎ and , corresponding to the momentum, heat and species mass transfer between the gas phase and the catalytic bed in the gas governing equations (Eqs. (6)-(9)).

Numerical implementation of the operator-splitting
The governing equations solved for each sub-step of the operator-splitting algorithm discussed in section 3.1.2 are reported in the following. In particular, in the first sub-step of the algorithm, i.e. the first transport step, only the gas-particle transport of heat and mass are accounted for in the ODE system composed of the species, site species and energy balances on the particle (Eq. ((1)-(3))), allowing for their analytical solution, reported in Eqs. (S15)-(S17). In the second sub-step, i.e. the reaction step, the catalytic reactions are accounted for, neglecting the gas-particle transport contribution, thus requiring the solution of the ODE system reported in Eq. (S18), adopting the results of the first sub-step as initial conditions. Finally, in the third sub-step, i.e. the second transport step, 3 Additional details about the fluid dynamic predictions of the proposed framework

Assessment of the expansion dynamics of the fluidized bed
The temporal evolution of the bed height has been considered to investigate the entire bed expansion dynamics. The arithmetic mean of the height of all the particles (Eq. (S22)) has been adopted to enable the comparison with experimental data 7 . Figure S1 reports the results of this analysis. Even though the instantaneous experimental data of the bed height are underestimated by both the proposed framework and the numerical results of Goldschmidt, the average behavior of the bed is well reproduced. In particular, a temporal averaged bed height of 10 cm has been obtained, which is in good agreement with both the experimental and numerical value reported by Goldschmidt, i.e. 11.4 cm and 9.7 cm, respectively.

Selection of the simulation time step
The mean global error 〈 〉 introduced by operator-splitting has been assumed as the indicator for the selection of the simulation time step to ensure the accuracy of the proposed techniques: where , is the vector of the state variables derived with the operator-splitting approach of the p th catalyst particle at the i th time step, , is the vector of the state variables derived with the coupled approach, NT is the total number of simulation time steps and NS is the number of species in the catalytic particles. In particular, , vector represents the composition of the catalytic particle in case of rate equation kinetics, and both the composition and coverage of the particle in case of microkinetic modeling.
The description of the Catalytic Partial Oxidation (CPO) of methane by means of rate equations 8 has been adopted to select the proper simulation time step Δ . Table S1 shows the deviations introduced by the operator-splitting technique in function of the parameter under investigation. A Δ equal to 5 • 10 −6 s has been selected since further decrements of the simulation time step does not relevantly affects the accuracy of the operator-splitting. In addition to the selection of the time step, we performed also a scalability test by means of the adopted time step to properly choose the number of processors, since the proposed methodology is parallelizable by means of the Message Passing Interface (MPI) protocol. We computed the parallelization speed-up for an increasing number of CPUs and the operating conditions for the methane CPO reported in Table 1. The microkinetic model proposed by Maestri et al. 9 has been selected to account for the heterogeneous chemistry. The computational domain has been distributed over the CPUs by means of vertical slices parallel to the flow direction. In doing so, the relevant computational cost related to the particles has been properly subdivided over the cores and no particlefree CPUs has been found during the simulations, thus obtaining the best parallelization efficiency.  In this study, we selected 16 CPUs for the reactive tests reported. In fact, it is the highest number of cores characterized by an acceptable overhead (i.e. parallelization efficiency above 75%), introduced both in the solution of gas equations and in the DEM algorithm due to the inter-processor particle transfers.

Speed-up trend related to the operator-splitting for the rate equation kinetics
To furtherly investigate the trend of the speed-up provided by the operator-splitting reported in Figure   6, the maps of the chemical speed-up in the fluidized bed have been studied, for the simulation times 0.1 and 1.2 s, i.e. before and after the beginning of the syngas production. First, 10 time-steps have been simulated with and without the operator-splitting starting from coupled approach results obtained for the selected times. Then, the computational costs evaluated for each particle at the last  The highlighted dependence of the speed-up factor from the chemical composition in the reactor can be explained considering the characteristic times of the involved chemical phenomena. In fact, a slowdown of the simulation caused by the operator-splitting approach has been found whenever the characteristic time of consumption or production of at least one of the involved species resulted smaller than the adopted simulation time step. In particular, in the specific CPO case study, the aforementioned situation is generated by fast combustion reactions on Rhodium catalyst and in particular by the CO and H2 combustions which are two and three orders of magnitude faster than the methane oxidation 8 , whose characteristic time is of the same order of magnitude of the simulation time step, thus explaining the key role of oxygen and syngas.
In order to better describe the complex behavior of the chemical speed-up provided by the operator-splitting technique, the stiffness of the ODE system composed of the species mass balances in the catalytic pellet has been analyzed. The rate of change of the transport and reaction phenomena and of the sole catalytic reactions has been investigated for coupled approach ODE system and the operator-splitting one related to the reaction step. Since these phenomena are proportional to the first time derivative of the composition for a generic particle p, the vector of the second order time derivatives of the species mass fraction , (Eq. (S24)) has been selected as the indicator of the stiffness of the system 10 : where is the Jacobian matrix of the ODE system to be solved. In fact, the higher are the second derivatives, the faster are the changes of the reaction rates during the simulation time and the lower is the time step of the ODE solver required to correctly describe the phenomena, resulting in a higher computational cost. c Figure S4 reports the norm of the second derivatives vector in function of the integration time of the particles (i.e. until 5 • 10 −6 s), for four catalytic pellets composing the fluidized bed, each one representative of one of the four efficiency regions reported in Figure S3b for the pseudo steady state.
In particular, the analyses of the ODE system related to both the coupled approach and the reaction step of the operator-splitting has been reported to define the influence of species transport on the stiffness of the CPO kinetic scheme.
In the fresh feed region, the higher methane oxidation reactivity is not smoothed by the gas-particle mass transfer when the operator-splitting approach is adopted, since gas-solid transport phenomena are not accounted for in the reaction step. Consequently, an increment of the ODE system stiffness with respect to the coupled approach is observed in Figure S4a, thus explaining the slow-down of the solution of the particles reported in Figure S3b. In the total oxidation region, the reactivity is smoothed by the lower concentration of the oxygen. Thus, the ODE system in the reaction step can be integrated more easily with respect to the system in the coupled approach. In fact, the interplay between transport and reaction phenomena, characterized by relevantly different time scales, must not be solved, thus the coupled approach is order of magnitudes stiffer than the operator-splitting, as reported in Figure   S4b. Therefore, computational gain is achieved for the methane total oxidation reactivity ( Figure   S3b). In the intermediate region, the simultaneous presence of non-negligible amounts of oxygen and syngas makes not negligible the stiffness introduced by the side reactions of CO and H2 combustions ( Figure S4c), causing a relevant slow-down of the simulation of the particles ( Figure S3b). In fact, the reactivity of the side syngas combustion is no more controlled by the lack of syngas and the mitigating effect of gas-solid transport, which lowers the concentration of syngas produced in the particle, becomes pivotal to control the characteristic times of the chemical phenomena. Therefore, a relevant decrement of the solution step adopted by the ODE solver is required in the reaction step of operator-splitting to properly describe the fast combustion kinetics which causes the complete depletion of oxygen at 10 −7 s. As expected, after this time, the absence of oxygen stops the syngas oxidation and the reaction step ODE system becomes less stiff than the coupled approach one ( Figure  S4c), coherently with the syngas production slice of the bed ( Figure S4d) where no oxygen is available ( Figure S3b).

Numerical explanation of the speed-up trend related to ISAT for the rate equation kinetics
The distribution of the Chemical speed-up factor has been studied for the simulation time 1. As expected, Chemical speed-up factors between 20 and 30 are achieved in almost the whole catalytic bed, since the solution of the heterogeneous chemistry of these catalytic particles is retrieved from the ISAT table, as reported in Figure S5b. Only few particles still require the solution of the ODE system for the growth and addition operations of the ISAT algorithm and, thus, still experience the