Coupling Euler–Euler and Microkinetic Modeling for the Simulation of Fluidized Bed Reactors: an Application to the Oxidative Coupling of Methane

We propose a numerical methodology to combine detailed microkinetic modeling and Eulerian–Eulerian methods for the simulation of industrial fluidized bed reactors. An operator splitting-based approach has been applied to solve the detailed kinetics coupled with the solution of multiphase gas–solid flows. Lab and industrial reactor configurations are simulated to assess the capability and the accuracy of the method by using the oxidative coupling of methane as a showcase. A good agreement with lab-scale experimental data (deviations below 10%) is obtained. Moreover, in this specific case, the proposed framework provides a 4-fold reduction of the computational cost required to reach the steady-state when compared to the approach of linearizing the chemical source term. As a whole, the work paves the way to the incorporation of detailed kinetics in the simulation of industrial fluidized reactors.


Euler-Euler Closure Models
Table S1-Closure models adopted in this work to describe the gas-solid drag force, the heat and mass transfer coefficients, the convective and diffusive fluxes and, the KTGF properties in the multiscale Euler-Euler framework Quantity Closure Model Ref.
Gas-particle drag force Dissipation of granular energy related to the inelastic particleparticle collisions,

Chemical and Mass Transfer Regime
The Multiphase Operator Splitting approach has been tested in the chemical and mass transport regime in order to assess the validity of the method with different regimes. This assessment has been performed by considering a fixed bed reactor filled with particles having a diameter of 500 µm. A simple kinetic mechanism composed of a heterogeneous first-order irreversible reaction (Eq. (S1)-(S2)) has been employed in the analysis: The operating conditions adopted for the numerical tests are a reactor temperature of 623 Damkohler numbers (Da). In particular, a Da equal to 10 3 has been adopted to simulate the system operated in external mass transfer limited regime, while a Da equal to 10 -2 has been selected as representative of the chemical regime.  An excellent agreement has been obtained assessing the capability of the MOS to accurately predict the evolution of the system in both the conditions. Hence, the hereby proposed splitting method is able to deal with different regimes.

OCM Microkinetic Mechanism
The elementary steps composing the gas phase and heterogenous microkinetic schemes are reported in Table S2 and Table S3, respectively. The reaction rate is calculated by means of the modified Arrhenius equation: where is the temperature-independent pre-exponential coefficient of the i-th reaction, T is the temperature, β i is the temperature exponent of the i-th reaction, is the activation energy of the ith reaction and, R is the ideal gas constant.

Configuration
To show the effect of the mesh on the lab-scale configuration, a mesh convergency analysis has been performed. In particular, the analysis consists of simulating the lab-scale system by using three different grids characterized by a cell to particle ratio in the range 10-20. A ratio lower than 10 has not been simulated due to the not accurate KTGF predictions when a fine computational grid is adopted 1 , while a ratio higher than 20 has not been tested due to the small number of cells in the cross section area that in our opinion is not suitable for accurate predictions. The size of the computational cells adopted to divide the lab-scale computational domain is smaller than the fluid dynamic structure present inside the units 10 . Therefore, the correlation proposed by Gidaspow 1 can be accurately used to describe the gas-solid interactions. The analysis has been carried out in terms of average methane and ethane content in the catalytic bed. Figure S2 shows the temporal trend of the volume averaged mass fraction of methane and ethane in the catalytic bed (Eq (21) of the manuscript) obtained with the different computational grids.

Figure S2 -Temporal trend of the volume average mass fraction of methane (a) and ethane (b) in the catalytic bed obtained with different computational domains in the lab-scale configuration.
The profiles of the two species show that the discretization has a negligible influence on the chemistry. In particular, a maximum deviation lower than 2% is obtained by using the three different computational grids both during the transient and at steady state. Consequently, the computational domain with a cell size of 10 dp can be properly adopted to validate the predictions of the MOS-based Euler-Euler framework against the experimental data.

Configuration
To show the effect of the mesh on the industrial scale configuration, a grid convergency analysis has been performed. In particular, the computational domain has been divided with three different computational grids characterized by a maximum cell to particle ratio in the range 250-750. Figure   S3 shows the maps of the steady-state void fraction obtained with the three different computational grids by considering the Gidaspow 1 model for the drag force. Figure S3 -2D slices of the steady-state solid fraction obtained by adopting a computational grid having maximum cell size equal to 250 dp (a), 500 dp (b) and 750 dp (c).   The decrement of the specific surface area of the catalyst leads to an increment of the volume average methane mass fraction due to the reduction of the catalyst activity. Coherently, the increment of the specific surface area of the catalyst leads to a decrement of the volume average methane mass fraction, as expected.  Table S4 reports the maximum absolute deviation with respect to the base case (i.e. 0 ) for each change of the catalytic loading. As shown in the table, a maximum deviation lower than 2.5% has been obtained even by modifying the specific surface area of the catalyst by 10%. Therefore, the S11 overall effect of the catalytic loading is sub-linear and thus, the adopted value is robust and representative of the system. S12

Source Term
The comparison between the MOS approach and the linearization of the chemistry source term performed in the lab-scale reactor revealed the presence of a delay in the reactor dynamics in the latter case. We also observed that the delay is a function of the time step and it is the result of numerical artifacts introduced by the linearization of the chemistry source terms. We ascribe the numerical delay to the implicit treatment of the linearized source term required to obtain a numerically stable reactive simulation. This approach employs a Taylor expansion as reported in eq (S8): where is the i-th species mass fraction, � ( )� is the i-th species production rate related to the i-th species mass fraction at the beginning of the time step and � ( +1 )� is the i-th species production rate related to the i-th species mass fraction at the end of the time step. In doing so, the derivative term during a time step is considered constant and equal to the derivative term at the beginning of the time step when the maximum concentration of reactant is experienced. Since the derivative of the reaction rate with respect to the reactant is negative, the reaction rate experiences during a long time step an excessive decrement since the derivative is based on the maximum reactant condition. Consequently, a lower average reaction rate is experienced during the time step introducing the delay in the dynamics of the species. Coherently, the longer is the time step, the more relevant is the error and delay introduced since the less adequate is the linear reaction rate approximation.