DFT calculation is conducted for H₂ reducing fully oxidized CMTF8341. The reaction path is illustrated in Figure 2, consisting of (three elementary steps) as Equations 4–6. First, one H₂ molecule undergoes a dissociative adsorption onto the surface; one H atom is combined with an active O atom on a Mn site, forming a hydroxyl (OH) radical; the other H atom is more likely to be placed on a Ca atom, rather than another Mn site. With Mn sites represented by * signs, and Ca sites by # signs, the elementary reaction is expressed as H 2 g + O * + # ⇌ OH * + H # (4)

Subsequently, the single H atom is associated with the OH radical, forming one H₂O molecule on the Mn site. OH * + H # ⇌ H 2 O * + # (5)

The last step is the H₂O molecule desorbing into the gas phase, leaving a vacant Mn site on the surface. H 2 O * ⇌ H 2 O g + * (6)

The energy diagram, having undergone zero-point correction, is plotted in Figure 3. The vibration frequencies of the species and transition states, along with other basic parameters, are listed in Table 1.

The rate equation of a grain involves not only the solid conversion, but also the gas partial pressure at the surface, which is not equal for all grains due to their spatial distribution. Grains at the particle surface obtain the most concentrated gas from the atmosphere, while less gas is received by inner grains as it diffuses through the pores, simultaneously consumed by outer grains. Assuming that the grains are evenly distributed throughout the particle’s volume, a radial gas distribution is governed by an equation considering both reaction and diffusion, as Equation 28 ∂ p ∂ t = 1 r 2 ∂ ∂ r r 2 D eff ∂ p ∂ r − ω p ∂ p ∂ r r = 0 = 0 p r = r p = p s (28) where r p is the particle radius, and p s is the surface gas partial pressure. The reaction term is derived from the gas–solid stoichiometry, as Equation 29. ω = 1 − α α ρ t R OC R T M O · L i · k + 1 1 − X p (29)

The diffusion term is described by the Fick’s Law; the diffusivity shall combine molecular and Knudsen diffusion because the pore radius is smaller than the molecular mean free path length. The particle conversion is the average of all grains, given by Equation 30. d X p d t = ∫ 0 r p L i · k + 1 1 − X g p · 4 π r 2 d r / 4 3 π r p 3 (30)

The governing equation, Equation 28, requires numerical solution due to both p and X g varying along r , which is not applicable to numerous particles in a reactor. Alternatively, a further simplification based on Thiele modulus (Sedghkerdar and Mahinpey, 2015; Yang et al., 2016), proposed in previous studies (Wang et al., 2017; Wang et al., 2021), offers an approximate analytic solution as Equation 31 d X p d t = η · L i · k + 1 1 − X p p s (31) where η is the effectiveness factor of intraparticle gas diffusion, given by Equation 32. η = 3 ϕ 2 ϕ ⁡ coth ⁡ ϕ − 1 ϕ = r p ω D eff (32) the Thiele modulus ϕ representing the influence of reaction versus diffusion.

The conversion rate of a particle, as Equation 31, is based on its surface gas partial pressure. Each particle in a fluidized bed moves along a unique trajectory under the fluidization by the inlet gas flow, and is therefore exposed to a gas partial pressure different from other particles. The fluidization behavior is described by the Eulerian–Lagrangian CFD–DEM model in this study, coupling the single-particle reaction kinetics with the particle–fluid motion. The thermal effects are omitted due to the low reduction heat generally observed from various oxygen carriers (García-Labiano et al., 2005; Hallberg et al., 2011), as well as the stable temperature measured during the experiment as Section 3.1.

Validation experiment is conducted on a micro-fluidized-bed thermogravimetric analyzer (MFB–TGA) (Li et al., 2019). A fluidized bed reactor, 3 cm in diameter and 10 cm in height, is placed on an electronic balance, measuring the real-time solid mass change under fluidization throughout the reaction process. More information is provided in the Supplementary Material.

The measurement accuracy is 1 mg, compared with the maximum mass change of 16 mg. Gas inlet and outlet are connected to the reactor with soft tubes to reduce mass fluctuation introduced by the gas circuit. The temperature is maintained constant by an electric furnace and monitored by a K-type thermocouple, which do not contact the reactor. Consequently, the signal from the electronic balance shall accurately present the total solid mass change in the reactor.

The bed inventory consists of silica sand and oxygen carrier particles, whose fluidization properties are listed in Table 2. A bubbling fluidization regime is observed under a gas flux of 1.2 NL/min. Continuous redox cycles are performed by switching the gas among inert (N₂), oxidizing (O₂) and reducing (H₂) components, under different H₂ concentrations (5 vol%, 10 vol% and 20 vol%) and temperatures (750 °C, 800 °C, 850 °C and 900 °C).

The simulating conditions in CFD–DEM are identical to those in the MFB–TGA experiment, computational settings and reaction parameters listed in Table 3. Initialization is done by injecting given numbers of sand and fully oxidized oxygen carrier particles, which are immediately packed under gravity. Subsequently, an inert fluidizing gas stream is introduced from the bottom at a given temperature. The gas is switched to a reducing composition after a stable fluidization is achieved, and the reaction simultaneously begins. A mesh-refinement test showing grid-independence is presented in the Supplementary Material. Every case is parallelized into 24 CPU cores and run on a supercomputing center, taking 8 h of real time to progress 1 s of simulation time.

The open source software CFDEM®coupling (Kloss et al., 2012) is selected for CFD–DEM computation, which serves as an interface alternately progressing CFD and DEM time steps. CFD is solved by OpenFOAM^® with the PISO algorithm; discretization is accomplished through the finite volume method. Particle motion and collision are solved by LIGGGHTS^®, employing a first-order Euler time integration scheme.

The reaction kinetics is not implemented in the original software; instead, it is programmed by the authors into an extended package, whose architecture is shown in Figure 4. The overall reaction is described by a class on the top layer, which outputs the particle conversion rate as Equation 31 to CFDEM®coupling; a list of elementary reactions is recorded as its component. The elementary reaction class calculates the rate constants as Section 2.2.2, relying on the reactant and product species along with the transition state. The species class aggregates submodels calculating the equation of state and partition functions; respective model types are selected for every species. DFT calculation is separately performed on the VASP software (Kresse and Furthmüller, 1996) before CFD–DEM computation.

All inputs required by the model, including measured properties, DFT results, reaction equations and submodel types, are provided in a text file during runtime; various overall reactions can be constructed by modifying the inputs, which results in different implementations of the abstract classes.

The conversion of all oxygen carrier particles in the reactor is derived from both computation and experiment as Figure 5. The simulated result represents the average conversion of all particles, corresponding to the experimentally measured conversion of the whole bed inventory. An approximately linear conversion growth is observed in the initial fast stage under every condition, which is transformed to a slow stage as the conversion approaches 0.8. Such phenomenon is affected by the mass transfer resistance of the reactor; as particles near the inlet are converted, more gas is supplied to farther particles, thus compensating the reaction deceleration caused by increased conversion. Therefore, restricting the residence time of particles, to maintain the solid conversion below 0.8, would be beneficial to chemical looping which requires fast redox cycles rather than complete conversion.

The conversion curves under different inlet H₂ concentrations, at a temperature of 900 °C, are plotted in Figure 5A, which shows that the reaction is accelerated as the inlet concentration increases. An average conversion of 0.8 is reached at ∼8 s for 20 vol% H₂, while the time is prolonged to ∼15 s for 10 vol% H₂, and ∼30 s for 5 vol% H₂. Above results indicate that the fast-stage rate is approximately proportional to the inlet gas concentration, proved by the particle rate equation, Equation 31, where the gas partial pressure appears as an independent factor.

The effect of temperature under an equal H₂ concentration (10 vol%) is shown in Figure 5B, a significant rate increase observed at a higher temperature. Temperature has a decisive impact on the rate constant k + 1 ; two temperature-dependent factors compose the rate constant as Equation 12, including an exponential function of the energy barrier, and a quotient of partition functions. Changes of the rate constant and the two factors against temperature are plotted in Figure 6. The partition function quotient is approximately a power function of the temperature, which grows slower than the exponential factor in the given temperature range; consequently, the rate constant changes in a similar trend to the exponential factor. The impact level of temperature is determined by the energy barrier, Δ ε 0 , + 1 ; the reduction rate by H₂ in this study is more sensitive to the temperature compared with CO reduction in a prior study (Wang et al., 2025), which is because a higher energy barrier demands more energy for the reactants to become transition states, thus magnifying the effect of temperature.

The error margin of solid conversion is determined by the measurement accuracy of MFB–TGA (1 mg), resulting in a solid conversion error of 0.0625. The comparison between computational and experimental results are displayed in Figure 7, showing that the multiscale computation is accurate at most data points. However, exceeded errors are observed under high H₂ concentration (20 vol%) or low temperature (750 °C) in the middle of the reaction processes. The concentration-related error is mainly introduced by mass transfer, which has a greater impact when a higher inlet concentration leads to a non-uniform distribution; unresolved CFD–DEM omits the fine structure of the particle surface boundary layer, tending to underestimate the mass transfer resistance (Derksen, 2014; Srinivasakannan et al., 2012). Integrating a surface film model into the particle-scale diffusion is a possible way to show the effect of boundary-layer mass transfer; interparticle gas diffusivity in the dense phase should be more accurately modeled considering the sub-scale velocity distribution, which may need validation by particle-resolved DNS–DEM before application to CFD–DEM (Wang et al., 2023). The temperature-related error depends on the energy barrier Δ ε 0 , + 1 , which is generated during DFT calculation; an increased energy barrier results in a more temperature-sensitive reaction rate. The energy error mainly depends on the functional selection; as different functionals are suitable for different systems, accuracy can be improved by choosing an appropriate functional for calculation (Kim et al., 2013). Besides, various energy correction schemes have been developed against systematic errors, based on error analysis versus experimental data, or functional dependence on energies (Christensen et al., 2015).

The areal density of Mn sites ( Λ ) and the number of neighboring sites ( d * ) are hypothetic parameters based on Langmuir’s adsorption model, which are not yet measurable; therefore, they are regarded as adjustable parameters to be assigned from the experimental results. The Mn site areal density (4.60 × 10⁻⁸ mol/m²) has been determined in a previous study (Wang et al., 2025), where the single-site mechanism of CMTF8341 reduced by CO, occurring on the same Mn site, was modeled; Λ appeared as the only unknown parameter, and was adjusted according to the CO-reduction results. The number of neighbors is valued as 1 so that the computational results fit those from the experiment. This result indicates that every Mn site is paired to exactly one nearest Ca site, and every dual-site reaction occurs on a predetermined pair of sites.

The microkinetics as Equations 18–20 is originally solved as an ordinary differential equation set; an example condition of 20 vol% H₂, 1 vol% H₂O and 900 °C is adopted in this section. The coverages of all species on the primary (Mn) site vary as Figures 8A,B. The major components are the reactant ( O * ) and product ( * ), which undergo substantial conversion throughout the reaction process; coverages of the minor intermediates ( OH * and H 2 O * ) grow in the same trend, but are restricted in an order of 10^–9 and 10^–5, respectively. The forward and reversed rates of every elementary reaction change as Figure 8C. Every elementary reaction starts with a positive net rate, converting O * to OH * , H 2 O * , and * in order. An equilibrium is established for reaction Equation 6 immediately after the reaction begins, owing to a high turnover frequency of the H₂O adsorption–desorption process. As the reaction further progresses, Equation 5 also approaches equilibrium at ∼1 s as its reversed rate increases. Consequently, Equations 5, 6 can be accurately described by the partial equilibrium assumption. Above results have proved that Equation 4 is the rate-determining step, which is far from equilibrium throughout the reaction. The approximated coverages of surface species in Equation 24 are thus verified, leading to an analytic grain conversion rate as Equation 26.

The overall reaction Equation 17 has an equilibrium point when elementary reactions Equations 4, 5, 6, 14 are all equilibrated, as required by the principle of detailed balance. Let Equations 4, 5, 6, have zero net rates in Equations 5, 8 and 9, combined with the ion-diffusion equilibrium as Equation 15, resulting in Equation 41. p H 2 O C V O p H 2 C O O = q H 2 O 0 V q * 0 q H 2 0 V q O * 0 exp − Δ ε 0 , overall k B T (41)

Note that the solid concentrations are involved because crystal defects are treated as non-localized systems, which is not the case for perfect crystal species in other common reactions. The equilibrium solid conversion changes against gas concentrations as Figure 9, suggesting that the solid is almost fully converted under normal conditions when both p H 2 and p H 2 O have an order of kPa–MPa. Only under extreme conditions when p H 2 / p H 2 O < 10 − 13 , typically when the gas reactant is depleted due to insufficiency, can the solid achieve incomplete conversion.

The number of neighboring sites ( d * ), in every dual-site elementary reaction rate, is a positive integer depending on the surface site distribution. With (4) being the rate-determining step, d * appears as an independent factor in the grain’s rate as Equation 26. Figure 10 shows the grain’s conversion curve at 20 vol% H₂ and 900 °C, as d * increases from 1 to 4. Rates proportional to d * are observed in the early stage when X g < 0.5 . Consequently, the proportional effect of d * propagates to the particle and reactor scale, given that all grains and particles share the same d * . A proper value of d * should fit every concentration and temperature, as it remains constant under all operating conditions.

The evolution of particle distribution at 900 °C, 20 vol% H₂, is illustrated in Figure 11. The moment of 0 s is the starting point of reaction, when a stable bubbling fluidization has been established. Particles in the center are elevated by the gas stream, and subsequently fall along the walls under gravity; the bed height fluctuates around 3 cm as the void fraction of the bed inventory zone changes. The conversions of oxygen carrier particles are marked by color in Figure 11, where only a minor difference exist among all particles, suggesting that an even mixture is accomplished by fluidization.

Both the computational and experimental results in Section 4.1 display the overall solid conversion, which equals the average of all particles. However, the conversion of a single particle may deviate from the average. The average, maximum and minimum conversions among all particles are plotted in Figure 12A, along with three typical particles, at 20 vol% H₂, 900 °C. Particle 3 undergoes an average-rate reaction, while Particle 1 is faster and Particle 2 slower. These particles are exposed to particular H₂ concentrations at every moment, determined by their unique trajectories during fluidization. The positions of the three particles, sampled during the initial fast stage of the reaction (0–5 s), are presented in Figures 12B,C. Particle 1 has the lowest average position in height, and Particle 2 the highest, which indicates that Particle 1 is exposed to a higher concentration of H₂ than Particles 3 and 2.

A deviation among particle conversions is accumulated during the initial fast stage, with the fastest rate being 2–3 times the slowest rate, which subsequently decreases as the particles are converted. Histograms in Figure 13A display the particle conversion distributions at different average conversions, which all show a normal distribution scheme. Besides, the particle conversion distribution is biased to the lower side in the beginning, which is because most H₂ is consumed by the few particles near the inlet, while other particles only obtain a low H₂ concentration.

The standard deviation reaches its maximum at an average conversion of 0.3–0.5, under all tested conditions as Figure 13B; a greater standard deviation is observed for a condition with a faster reaction rate (given by Figure 5). The maximum standard deviation can be selected as a quantitative indicator of particle diversity, which is plotted against the theoretical reaction rate in Figure 13C. The theoretical reaction rate is calculated from Equation 31, when the particle has a conversion of zero and obtains an inlet gas concentration, acting as a superior limit which leads to a rate between 0% and 100% of this limit for every single particle. Results show that the maximum standard deviation is approximately linear against the logarithm of theoretical rate. As the fluidization process provides the spatial H₂ concentration distribution with a similar pattern under all conditions, consequently, a greater deviation is formed at a faster theoretical rate. Such diversity of particles decelerates the overall reaction, because more gas is involved in a slower reaction stage with the highly-converted particles; the utilization efficiency of the solid material is thus affected, which shall be compensated in the system design.

The influence of intraparticle gas diffusion is represented by the effectiveness factor η , defined as Equation 32, which changes against the particle conversion as Figure 14. The effectiveness factor at every temperature remains above 0.9 throughout the reaction process, indicating that the conversion rate of the overall particle is approximately equal to that of a grain at the particle surface in the studied cases. A larger Thiele modulus results in a smaller effectiveness factor; thus, η increases with solid conversion, given that the rate decreases during the reaction; meanwhile, η decreases with temperature, owing to the rate constant ( k + 1 ) positively correlated to the temperature. Generally, the reduction of CMTF8341 by H₂ has a relatively low rate compared with gas diffusion, so intraparticle gas diffusion has a minor impact among all processes in the multiscale framework.

The overall rate equation derived from quasi-steady approximation, Equation 26, is first-order for the solid species, as the solid conversion appears as a factor of 1 − X g . This feature is mainly determined by the rate-determining step, Equation 4, where only one solid reactant ion ( O * ) is involved. Contrarily, previous studies on H₂ reduction kinetics (Li and Li, 2024; Li, 2022; Li Z. et al., 2021) have generally adopted the symmetrical dissociation mechanism as Equations 42–44 H 2 g + 2 O * ⇌ 2 OH * (42) 2 OH * ⇌ H 2 O * + O * (43) H 2 O * ⇌ H 2 O g + * (44) where the H₂ molecule is dissociated on two identical sites ( * ) and combined with two oxygen atoms. Repeating the surface- and grain-scale derivations, with the rate-determining step being Equation 42, results in a second-order overall rate equation as Equation 45 d X g d t = L i · d * k + 1 1 − X g 2 p H 2 (45) which significantly varies from the first-order rate under dual-site dissociation Equation 4.

Experimental validation in Section 4.1 has agreed with the first-order rate equation, along with the dual-site dissociation mechanism. Dual-site microkinetics lead to asymmetrical impacts on bulk diffusion, which only occurs on primary sites (Mn) rather than secondary sites (Ca). Therefore, the atom-scale reaction mechanism may have crucial impacts on the mesoscale mass transfer and macroscopic kinetics.

The heavy computational cost has been an obstacle to scaling-up CFD–DEM simulation to industrial-scale reactors. Parallel computing accelerates computation by independently handling local interactions. A scale-up test is conducted based on the setup as Section 3.2, by varying the number of processors between 2 and 48, thus changing the number of particles per processor. The real time cost per second of simulation time is plotted in Figure 15. An approximately proportional acceleration is obtained in lowly-parallelized cases. However, extra cost is introduced by interprocess communication; the computational speed thus remains stable and subsequently slows down as more processors are employed.

Therefore, further acceleration approaches shall provide essential improvements, for either computation or the DEM model itself. Graphical processing units (GPU), which are especially applicable to batches of simple computation, can be used to handle the massive particle collision (Lu, 2022). Existing investigation on modeling improvements include: (1) the coarse-grain model, which significantly reduces particles in the system, by combining multiple particles into a parcel (Sakai et al., 2014); (2) machine learning algorithms, which replaces direct force calculation with a prediction–correction updating (Lu et al., 2021). These modeling treatments have partially sacrificed the accuracy of fluidization; however, given the minor effect of particle diversity (as Section 4.3), such simplification would be acceptable, provided that the fluidization regime is maintained.